projects / unix-history / blame
? search:
summary | tags | clone url | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
This commit was generated by cvs2svn to track changes on a CVS vendor
[unix-history] / gnu / usr.bin / cc / lib / real.c
CommitLineData
b2254bec
PR		 1/* real.c - implementation of REAL_ARITHMETIC, REAL_VALUE_ATOF,
2and support for XFmode IEEE extended real floating point arithmetic.
3Contributed by Stephen L. Moshier (moshier@world.std.com).
4
5 Copyright (C) 1993 Free Software Foundation, Inc.
6
7This file is part of GNU CC.
8
9GNU CC is free software; you can redistribute it and/or modify
10it under the terms of the GNU General Public License as published by
11the Free Software Foundation; either version 2, or (at your option)
12any later version.
13
14GNU CC is distributed in the hope that it will be useful,
15but WITHOUT ANY WARRANTY; without even the implied warranty of
16MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17GNU General Public License for more details.
18
19You should have received a copy of the GNU General Public License
20along with GNU CC; see the file COPYING. If not, write to
21the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
22
23#include <stdio.h>
24#include <errno.h>
25#include "config.h"
26#include "tree.h"
27
28#ifndef errno
29extern int errno;
30#endif
31
32/* To enable support of XFmode extended real floating point, define
33LONG_DOUBLE_TYPE_SIZE 96 in the tm.h file (m68k.h or i386.h).
34
35To support cross compilation between IEEE and VAX floating
36point formats, define REAL_ARITHMETIC in the tm.h file.
37
38In either case the machine files (tm.h) must not contain any code
39that tries to use host floating point arithmetic to convert
40REAL_VALUE_TYPEs from `double' to `float', pass them to fprintf,
41etc. In cross-compile situations a REAL_VALUE_TYPE may not
42be intelligible to the host computer's native arithmetic.
43
44The emulator defaults to the host's floating point format so that
45its decimal conversion functions can be used if desired (see
46real.h).
47
48The first part of this file interfaces gcc to ieee.c, which is a
49floating point arithmetic suite that was not written with gcc in
50mind. The interface is followed by ieee.c itself and related
51items. Avoid changing ieee.c unless you have suitable test
52programs available. A special version of the PARANOIA floating
53point arithmetic tester, modified for this purpose, can be found
54on usc.edu : /pub/C-numanal/ieeetest.zoo. Some tutorial
55information on ieee.c is given in my book: S. L. Moshier,
56_Methods and Programs for Mathematical Functions_, Prentice-Hall
57or Simon & Schuster Int'l, 1989. A library of XFmode elementary
58transcendental functions can be obtained by ftp from
59research.att.com: netlib/cephes/ldouble.shar.Z */
60
61/* Type of computer arithmetic.
62 * Only one of DEC, MIEEE, IBMPC, or UNK should get defined.
63 */
64
65/* `MIEEE' refers generically to big-endian IEEE floating-point data
66 structure. This definition should work in SFmode `float' type and
67 DFmode `double' type on virtually all big-endian IEEE machines.
68 If LONG_DOUBLE_TYPE_SIZE has been defined to be 96, then MIEEE
69 also invokes the particular XFmode (`long double' type) data
70 structure used by the Motorola 680x0 series processors.
71
72 `IBMPC' refers generally to little-endian IEEE machines. In this
73 case, if LONG_DOUBLE_TYPE_SIZE has been defined to be 96, then
74 IBMPC also invokes the particular XFmode `long double' data
75 structure used by the Intel 80x86 series processors.
76
77 `DEC' refers specifically to the Digital Equipment Corp PDP-11
78 and VAX floating point data structure. This model currently
79 supports no type wider than DFmode.
80
81 If LONG_DOUBLE_TYPE_SIZE = 64 (the default, unless tm.h defines it)
82 then `long double' and `double' are both implemented, but they
83 both mean DFmode. In this case, the software floating-point
84 support available here is activated by writing
85 #define REAL_ARITHMETIC
86 in tm.h.
87
88 The case LONG_DOUBLE_TYPE_SIZE = 128 activates TFmode support
89 (Not Yet Implemented) and may deactivate XFmode since
90 `long double' is used to refer to both modes. */
91
92/* The following converts gcc macros into the ones used by this file. */
93
94/* REAL_ARITHMETIC defined means that macros in real.h are
95 defined to call emulator functions. */
96#ifdef REAL_ARITHMETIC
97
98#if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
99/* PDP-11, Pro350, VAX: */
100#define DEC 1
101#else /* it's not VAX */
102#if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
103#if WORDS_BIG_ENDIAN
104/* Motorola IEEE, high order words come first (Sun workstation): */
105#define MIEEE 1
106#else /* not big-endian */
107/* Intel IEEE, low order words come first:
108 */
109#define IBMPC 1
110#endif /* big-endian */
111#else /* it's not IEEE either */
112/* UNKnown arithmetic. We don't support this and can't go on. */
113unknown arithmetic type
114#define UNK 1
115#endif /* not IEEE */
116#endif /* not VAX */
117
118#else
119/* REAL_ARITHMETIC not defined means that the *host's* data
120 structure will be used. It may differ by endian-ness from the
121 target machine's structure and will get its ends swapped
122 accordingly (but not here). Probably only the decimal <-> binary
123 functions in this file will actually be used in this case. */
124#if HOST_FLOAT_FORMAT == VAX_FLOAT_FORMAT
125#define DEC 1
126#else /* it's not VAX */
127#if HOST_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
128#ifdef HOST_WORDS_BIG_ENDIAN
129#define MIEEE 1
130#else /* not big-endian */
131#define IBMPC 1
132#endif /* big-endian */
133#else /* it's not IEEE either */
134unknown arithmetic type
135#define UNK 1
136#endif /* not IEEE */
137#endif /* not VAX */
138
139#endif /* REAL_ARITHMETIC not defined */
140
141/* Define INFINITY for support of infinity.
142 Define NANS for support of Not-a-Number's (NaN's). */
143#ifndef DEC
144#define INFINITY
145#define NANS
146#endif
147
148/* Support of NaNs requires support of infinity. */
149#ifdef NANS
150#ifndef INFINITY
151#define INFINITY
152#endif
153#endif
154
155/* ehead.h
156 *
157 * Include file for extended precision arithmetic programs.
158 */
159
160/* Number of 16 bit words in external e type format */
161#define NE 6
162
163/* Number of 16 bit words in internal format */
164#define NI (NE+3)
165
166/* Array offset to exponent */
167#define E 1
168
169/* Array offset to high guard word */
170#define M 2
171
172/* Number of bits of precision */
173#define NBITS ((NI-4)*16)
174
175/* Maximum number of decimal digits in ASCII conversion
176 * = NBITS*log10(2)
177 */
178#define NDEC (NBITS*8/27)
179
180/* The exponent of 1.0 */
181#define EXONE (0x3fff)
182
183/* Find a host integer type that is at least 16 bits wide,
184 and another type at least twice whatever that size is. */
185
186#if HOST_BITS_PER_CHAR >= 16
187#define EMUSHORT char
188#define EMUSHORT_SIZE HOST_BITS_PER_CHAR
189#define EMULONG_SIZE (2 * HOST_BITS_PER_CHAR)
190#else
191#if HOST_BITS_PER_SHORT >= 16
192#define EMUSHORT short
193#define EMUSHORT_SIZE HOST_BITS_PER_SHORT
194#define EMULONG_SIZE (2 * HOST_BITS_PER_SHORT)
195#else
196#if HOST_BITS_PER_INT >= 16
197#define EMUSHORT int
198#define EMUSHORT_SIZE HOST_BITS_PER_INT
199#define EMULONG_SIZE (2 * HOST_BITS_PER_INT)
200#else
201#if HOST_BITS_PER_LONG >= 16
202#define EMUSHORT long
203#define EMUSHORT_SIZE HOST_BITS_PER_LONG
204#define EMULONG_SIZE (2 * HOST_BITS_PER_LONG)
205#else
206/* You will have to modify this program to have a smaller unit size. */
207#define EMU_NON_COMPILE
208#endif
209#endif
210#endif
211#endif
212
213#if HOST_BITS_PER_SHORT >= EMULONG_SIZE
214#define EMULONG short
215#else
216#if HOST_BITS_PER_INT >= EMULONG_SIZE
217#define EMULONG int
218#else
219#if HOST_BITS_PER_LONG >= EMULONG_SIZE
220#define EMULONG long
221#else
222#if HOST_BITS_PER_LONG_LONG >= EMULONG_SIZE
223#define EMULONG long long int
224#else
225/* You will have to modify this program to have a smaller unit size. */
226#define EMU_NON_COMPILE
227#endif
228#endif
229#endif
230#endif
231
232
233/* The host interface doesn't work if no 16-bit size exists. */
234#if EMUSHORT_SIZE != 16
235#define EMU_NON_COMPILE
236#endif
237
238/* OK to continue compilation. */
239#ifndef EMU_NON_COMPILE
240
241/* Construct macros to translate between REAL_VALUE_TYPE and e type.
242 In GET_REAL and PUT_REAL, r and e are pointers.
243 A REAL_VALUE_TYPE is guaranteed to occupy contiguous locations
244 in memory, with no holes. */
245
246#if LONG_DOUBLE_TYPE_SIZE == 96
247#define GET_REAL(r,e) bcopy (r, e, 2*NE)
248#define PUT_REAL(e,r) bcopy (e, r, 2*NE)
249#else /* no XFmode */
250
251#ifdef REAL_ARITHMETIC
252/* Emulator uses target format internally
253 but host stores it in host endian-ness. */
254
255#if defined (HOST_WORDS_BIG_ENDIAN) == WORDS_BIG_ENDIAN
256#define GET_REAL(r,e) e53toe ((r), (e))
257#define PUT_REAL(e,r) etoe53 ((e), (r))
258
259#else /* endian-ness differs */
260/* emulator uses target endian-ness internally */
261#define GET_REAL(r,e) \
262do { EMUSHORT w[4]; \
263 w[3] = ((EMUSHORT *) r)[0]; \
264 w[2] = ((EMUSHORT *) r)[1]; \
265 w[1] = ((EMUSHORT *) r)[2]; \
266 w[0] = ((EMUSHORT *) r)[3]; \
267 e53toe (w, (e)); } while (0)
268
269#define PUT_REAL(e,r) \
270do { EMUSHORT w[4]; \
271 etoe53 ((e), w); \
272 *((EMUSHORT *) r) = w[3]; \
273 *((EMUSHORT *) r + 1) = w[2]; \
274 *((EMUSHORT *) r + 2) = w[1]; \
275 *((EMUSHORT *) r + 3) = w[0]; } while (0)
276
277#endif /* endian-ness differs */
278
279#else /* not REAL_ARITHMETIC */
280
281/* emulator uses host format */
282#define GET_REAL(r,e) e53toe ((r), (e))
283#define PUT_REAL(e,r) etoe53 ((e), (r))
284
285#endif /* not REAL_ARITHMETIC */
286#endif /* no XFmode */
287
288void warning ();
289extern int extra_warnings;
290int ecmp (), enormlz (), eshift ();
291int eisneg (), eisinf (), eisnan (), eiisinf (), eiisnan ();
292void eadd (), esub (), emul (), ediv ();
293void eshup1 (), eshup8 (), eshup6 (), eshdn1 (), eshdn8 (), eshdn6 ();
294void eabs (), eneg (), emov (), eclear (), einfin (), efloor ();
295void eldexp (), efrexp (), eifrac (), euifrac (), ltoe (), ultoe ();
296void eround (), ereal_to_decimal (), eiinfin (), einan ();
297void esqrt (), elog (), eexp (), etanh (), epow ();
298void asctoe (), asctoe24 (), asctoe53 (), asctoe64 ();
299void etoasc (), e24toasc (), e53toasc (), e64toasc ();
300void etoe64 (), etoe53 (), etoe24 (), e64toe (), e53toe (), e24toe ();
301void mtherr (), make_nan ();
302void enan ();
303extern unsigned EMUSHORT ezero[], ehalf[], eone[], etwo[];
304extern unsigned EMUSHORT elog2[], esqrt2[];
305
306/* Pack output array with 32-bit numbers obtained from
307 array containing 16-bit numbers, swapping ends if required. */
308void
309endian (e, x, mode)
310 unsigned EMUSHORT e[];
311 long x[];
312 enum machine_mode mode;
313{
314 unsigned long th, t;
315
316#if WORDS_BIG_ENDIAN
317 switch (mode)
318 {
319
320 case XFmode:
321
322 /* Swap halfwords in the third long. */
323 th = (unsigned long) e[4] & 0xffff;
324 t = (unsigned long) e[5] & 0xffff;
325 t |= th << 16;
326 x[2] = (long) t;
327 /* fall into the double case */
328
329 case DFmode:
330
331 /* swap halfwords in the second word */
332 th = (unsigned long) e[2] & 0xffff;
333 t = (unsigned long) e[3] & 0xffff;
334 t |= th << 16;
335 x[1] = (long) t;
336 /* fall into the float case */
337
338 case SFmode:
339
340 /* swap halfwords in the first word */
341 th = (unsigned long) e[0] & 0xffff;
342 t = (unsigned long) e[1] & 0xffff;
343 t |= th << 16;
344 x[0] = t;
345 break;
346
347 default:
348 abort ();
349 }
350
351#else
352
353 /* Pack the output array without swapping. */
354
355 switch (mode)
356 {
357
358 case XFmode:
359
360 /* Pack the third long.
361 Each element of the input REAL_VALUE_TYPE array has 16 bit useful bits
362 in it. */
363 th = (unsigned long) e[5] & 0xffff;
364 t = (unsigned long) e[4] & 0xffff;
365 t |= th << 16;
366 x[2] = (long) t;
367 /* fall into the double case */
368
369 case DFmode:
370
371 /* pack the second long */
372 th = (unsigned long) e[3] & 0xffff;
373 t = (unsigned long) e[2] & 0xffff;
374 t |= th << 16;
375 x[1] = (long) t;
376 /* fall into the float case */
377
378 case SFmode:
379
380 /* pack the first long */
381 th = (unsigned long) e[1] & 0xffff;
382 t = (unsigned long) e[0] & 0xffff;
383 t |= th << 16;
384 x[0] = t;
385 break;
386
387 default:
388 abort ();
389 }
390
391#endif
392}
393
394
395/* This is the implementation of the REAL_ARITHMETIC macro.
396 */
397void
398earith (value, icode, r1, r2)
399 REAL_VALUE_TYPE *value;
400 int icode;
401 REAL_VALUE_TYPE *r1;
402 REAL_VALUE_TYPE *r2;
403{
404 unsigned EMUSHORT d1[NE], d2[NE], v[NE];
405 enum tree_code code;
406
407 GET_REAL (r1, d1);
408 GET_REAL (r2, d2);
409#ifdef NANS
410/* Return NaN input back to the caller. */
411 if (eisnan (d1))
412 {
413 PUT_REAL (d1, value);
414 return;
415 }
416 if (eisnan (d2))
417 {
418 PUT_REAL (d2, value);
419 return;
420 }
421#endif
422 code = (enum tree_code) icode;
423 switch (code)
424 {
425 case PLUS_EXPR:
426 eadd (d2, d1, v);
427 break;
428
429 case MINUS_EXPR:
430 esub (d2, d1, v); /* d1 - d2 */
431 break;
432
433 case MULT_EXPR:
434 emul (d2, d1, v);
435 break;
436
437 case RDIV_EXPR:
438#ifndef REAL_INFINITY
439 if (ecmp (d2, ezero) == 0)
440 {
441#ifdef NANS
442 enan (v);
443 break;
444#else
445 abort ();
446#endif
447 }
448#endif
449 ediv (d2, d1, v); /* d1/d2 */
450 break;
451
452 case MIN_EXPR: /* min (d1,d2) */
453 if (ecmp (d1, d2) < 0)
454 emov (d1, v);
455 else
456 emov (d2, v);
457 break;
458
459 case MAX_EXPR: /* max (d1,d2) */
460 if (ecmp (d1, d2) > 0)
461 emov (d1, v);
462 else
463 emov (d2, v);
464 break;
465 default:
466 emov (ezero, v);
467 break;
468 }
469PUT_REAL (v, value);
470}
471
472
473/* Truncate REAL_VALUE_TYPE toward zero to signed HOST_WIDE_INT
474 * implements REAL_VALUE_RNDZINT (x) (etrunci (x))
475 */
476REAL_VALUE_TYPE
477etrunci (x)
478 REAL_VALUE_TYPE x;
479{
480 unsigned EMUSHORT f[NE], g[NE];
481 REAL_VALUE_TYPE r;
482 long l;
483
484 GET_REAL (&x, g);
485#ifdef NANS
486 if (eisnan (g))
487 return (x);
488#endif
489 eifrac (g, &l, f);
490 ltoe (&l, g);
491 PUT_REAL (g, &r);
492 return (r);
493}
494
495
496/* Truncate REAL_VALUE_TYPE toward zero to unsigned HOST_WIDE_INT
497 * implements REAL_VALUE_UNSIGNED_RNDZINT (x) (etruncui (x))
498 */
499REAL_VALUE_TYPE
500etruncui (x)
501 REAL_VALUE_TYPE x;
502{
503 unsigned EMUSHORT f[NE], g[NE];
504 REAL_VALUE_TYPE r;
505 unsigned long l;
506
507 GET_REAL (&x, g);
508#ifdef NANS
509 if (eisnan (g))
510 return (x);
511#endif
512 euifrac (g, &l, f);
513 ultoe (&l, g);
514 PUT_REAL (g, &r);
515 return (r);
516}
517
518
519/* This is the REAL_VALUE_ATOF function.
520 * It converts a decimal string to binary, rounding off
521 * as indicated by the machine_mode argument. Then it
522 * promotes the rounded value to REAL_VALUE_TYPE.
523 */
524REAL_VALUE_TYPE
525ereal_atof (s, t)
526 char *s;
527 enum machine_mode t;
528{
529 unsigned EMUSHORT tem[NE], e[NE];
530 REAL_VALUE_TYPE r;
531
532 switch (t)
533 {
534 case SFmode:
535 asctoe24 (s, tem);
536 e24toe (tem, e);
537 break;
538 case DFmode:
539 asctoe53 (s, tem);
540 e53toe (tem, e);
541 break;
542 case XFmode:
543 asctoe64 (s, tem);
544 e64toe (tem, e);
545 break;
546 default:
547 asctoe (s, e);
548 }
549 PUT_REAL (e, &r);
550 return (r);
551}
552
553
554/* Expansion of REAL_NEGATE.
555 */
556REAL_VALUE_TYPE
557ereal_negate (x)
558 REAL_VALUE_TYPE x;
559{
560 unsigned EMUSHORT e[NE];
561 REAL_VALUE_TYPE r;
562
563 GET_REAL (&x, e);
564#ifdef NANS
565 if (eisnan (e))
566 return (x);
567#endif
568 eneg (e);
569 PUT_REAL (e, &r);
570 return (r);
571}
572
573
574/* Round real to int
575 * implements REAL_VALUE_FIX (x) (eroundi (x))
576 * The type of rounding is left unspecified by real.h.
577 * It is implemented here as round to nearest (add .5 and chop).
578 */
579int
580eroundi (x)
581 REAL_VALUE_TYPE x;
582{
583 unsigned EMUSHORT f[NE], g[NE];
584 EMULONG l;
585
586 GET_REAL (&x, f);
587#ifdef NANS
588 if (eisnan (f))
589 {
590 warning ("conversion from NaN to int");
591 return (-1);
592 }
593#endif
594 eround (f, g);
595 eifrac (g, &l, f);
596 return ((int) l);
597}
598
599/* Round real to nearest unsigned int
600 * implements REAL_VALUE_UNSIGNED_FIX (x) ((unsigned int) eroundi (x))
601 * Negative input returns zero.
602 * The type of rounding is left unspecified by real.h.
603 * It is implemented here as round to nearest (add .5 and chop).
604 */
605unsigned int
606eroundui (x)
607 REAL_VALUE_TYPE x;
608{
609 unsigned EMUSHORT f[NE], g[NE];
610 unsigned EMULONG l;
611
612 GET_REAL (&x, f);
613#ifdef NANS
614 if (eisnan (f))
615 {
616 warning ("conversion from NaN to unsigned int");
617 return (-1);
618 }
619#endif
620 eround (f, g);
621 euifrac (g, &l, f);
622 return ((unsigned int)l);
623}
624
625
626/* REAL_VALUE_FROM_INT macro.
627 */
628void
629ereal_from_int (d, i, j)
630 REAL_VALUE_TYPE *d;
631 long i, j;
632{
633 unsigned EMUSHORT df[NE], dg[NE];
634 long low, high;
635 int sign;
636
637 sign = 0;
638 low = i;
639 if ((high = j) < 0)
640 {
641 sign = 1;
642 /* complement and add 1 */
643 high = ~high;
644 if (low)
645 low = -low;
646 else
647 high += 1;
648 }
649 eldexp (eone, HOST_BITS_PER_LONG, df);
650 ultoe (&high, dg);
651 emul (dg, df, dg);
652 ultoe (&low, df);
653 eadd (df, dg, dg);
654 if (sign)
655 eneg (dg);
656 PUT_REAL (dg, d);
657}
658
659
660/* REAL_VALUE_FROM_UNSIGNED_INT macro.
661 */
662void
663ereal_from_uint (d, i, j)
664 REAL_VALUE_TYPE *d;
665 unsigned long i, j;
666{
667 unsigned EMUSHORT df[NE], dg[NE];
668 unsigned long low, high;
669
670 low = i;
671 high = j;
672 eldexp (eone, HOST_BITS_PER_LONG, df);
673 ultoe (&high, dg);
674 emul (dg, df, dg);
675 ultoe (&low, df);
676 eadd (df, dg, dg);
677 PUT_REAL (dg, d);
678}
679
680
681/* REAL_VALUE_TO_INT macro
682 */
683void
684ereal_to_int (low, high, rr)
685 long *low, *high;
686 REAL_VALUE_TYPE rr;
687{
688 unsigned EMUSHORT d[NE], df[NE], dg[NE], dh[NE];
689 int s;
690
691 GET_REAL (&rr, d);
692#ifdef NANS
693 if (eisnan (d))
694 {
695 warning ("conversion from NaN to int");
696 *low = -1;
697 *high = -1;
698 return;
699 }
700#endif
701 /* convert positive value */
702 s = 0;
703 if (eisneg (d))
704 {
705 eneg (d);
706 s = 1;
707 }
708 eldexp (eone, HOST_BITS_PER_LONG, df);
709 ediv (df, d, dg); /* dg = d / 2^32 is the high word */
710 euifrac (dg, high, dh);
711 emul (df, dh, dg); /* fractional part is the low word */
712 euifrac (dg, low, dh);
713 if (s)
714 {
715 /* complement and add 1 */
716 *high = ~(*high);
717 if (*low)
718 *low = -(*low);
719 else
720 *high += 1;
721 }
722}
723
724
725/* REAL_VALUE_LDEXP macro.
726 */
727REAL_VALUE_TYPE
728ereal_ldexp (x, n)
729 REAL_VALUE_TYPE x;
730 int n;
731{
732 unsigned EMUSHORT e[NE], y[NE];
733 REAL_VALUE_TYPE r;
734
735 GET_REAL (&x, e);
736#ifdef NANS
737 if (eisnan (e))
738 return (x);
739#endif
740 eldexp (e, n, y);
741 PUT_REAL (y, &r);
742 return (r);
743}
744
745/* These routines are conditionally compiled because functions
746 * of the same names may be defined in fold-const.c. */
747#ifdef REAL_ARITHMETIC
748
749/* Check for infinity in a REAL_VALUE_TYPE. */
750int
751target_isinf (x)
752 REAL_VALUE_TYPE x;
753{
754 unsigned EMUSHORT e[NE];
755
756#ifdef INFINITY
757 GET_REAL (&x, e);
758 return (eisinf (e));
759#else
760 return 0;
761#endif
762}
763
764
765/* Check whether a REAL_VALUE_TYPE item is a NaN. */
766
767int
768target_isnan (x)
769 REAL_VALUE_TYPE x;
770{
771 unsigned EMUSHORT e[NE];
772
773#ifdef NANS
774 GET_REAL (&x, e);
775 return (eisnan (e));
776#else
777 return (0);
778#endif
779}
780
781
782/* Check for a negative REAL_VALUE_TYPE number.
783 * this means strictly less than zero, not -0.
784 */
785
786int
787target_negative (x)
788 REAL_VALUE_TYPE x;
789{
790 unsigned EMUSHORT e[NE];
791
792 GET_REAL (&x, e);
793 if (ecmp (e, ezero) == -1)
794 return (1);
795 return (0);
796}
797
798/* Expansion of REAL_VALUE_TRUNCATE.
799 * The result is in floating point, rounded to nearest or even.
800 */
801REAL_VALUE_TYPE
802real_value_truncate (mode, arg)
803 enum machine_mode mode;
804 REAL_VALUE_TYPE arg;
805{
806 unsigned EMUSHORT e[NE], t[NE];
807 REAL_VALUE_TYPE r;
808
809 GET_REAL (&arg, e);
810#ifdef NANS
811 if (eisnan (e))
812 return (arg);
813#endif
814 eclear (t);
815 switch (mode)
816 {
817 case XFmode:
818 etoe64 (e, t);
819 e64toe (t, t);
820 break;
821
822 case DFmode:
823 etoe53 (e, t);
824 e53toe (t, t);
825 break;
826
827 case SFmode:
828 etoe24 (e, t);
829 e24toe (t, t);
830 break;
831
832 case SImode:
833 r = etrunci (e);
834 return (r);
835
836 default:
837 abort ();
838 }
839 PUT_REAL (t, &r);
840 return (r);
841}
842
843#endif /* REAL_ARITHMETIC defined */
844
845/* Target values are arrays of host longs. A long is guaranteed
846 to be at least 32 bits wide. */
847void
848etarldouble (r, l)
849 REAL_VALUE_TYPE r;
850 long l[];
851{
852 unsigned EMUSHORT e[NE];
853
854 GET_REAL (&r, e);
855 etoe64 (e, e);
856 endian (e, l, XFmode);
857}
858
859void
860etardouble (r, l)
861 REAL_VALUE_TYPE r;
862 long l[];
863{
864 unsigned EMUSHORT e[NE];
865
866 GET_REAL (&r, e);
867 etoe53 (e, e);
868 endian (e, l, DFmode);
869}
870
871long
872etarsingle (r)
873 REAL_VALUE_TYPE r;
874{
875 unsigned EMUSHORT e[NE];
876 unsigned long l;
877
878 GET_REAL (&r, e);
879 etoe24 (e, e);
880 endian (e, &l, SFmode);
881 return ((long) l);
882}
883
884void
885ereal_to_decimal (x, s)
886 REAL_VALUE_TYPE x;
887 char *s;
888{
889 unsigned EMUSHORT e[NE];
890
891 GET_REAL (&x, e);
892 etoasc (e, s, 20);
893}
894
895int
896ereal_cmp (x, y)
897 REAL_VALUE_TYPE x, y;
898{
899 unsigned EMUSHORT ex[NE], ey[NE];
900
901 GET_REAL (&x, ex);
902 GET_REAL (&y, ey);
903 return (ecmp (ex, ey));
904}
905
906int
907ereal_isneg (x)
908 REAL_VALUE_TYPE x;
909{
910 unsigned EMUSHORT ex[NE];
911
912 GET_REAL (&x, ex);
913 return (eisneg (ex));
914}
915
916/* End of REAL_ARITHMETIC interface */
917
918/* ieee.c
919 *
920 * Extended precision IEEE binary floating point arithmetic routines
921 *
922 * Numbers are stored in C language as arrays of 16-bit unsigned
923 * short integers. The arguments of the routines are pointers to
924 * the arrays.
925 *
926 *
927 * External e type data structure, simulates Intel 8087 chip
928 * temporary real format but possibly with a larger significand:
929 *
930 * NE-1 significand words (least significant word first,
931 * most significant bit is normally set)
932 * exponent (value = EXONE for 1.0,
933 * top bit is the sign)
934 *
935 *
936 * Internal data structure of a number (a "word" is 16 bits):
937 *
938 * ei[0] sign word (0 for positive, 0xffff for negative)
939 * ei[1] biased exponent (value = EXONE for the number 1.0)
940 * ei[2] high guard word (always zero after normalization)
941 * ei[3]
942 * to ei[NI-2] significand (NI-4 significand words,
943 * most significant word first,
944 * most significant bit is set)
945 * ei[NI-1] low guard word (0x8000 bit is rounding place)
946 *
947 *
948 *
949 * Routines for external format numbers
950 *
951 * asctoe (string, e) ASCII string to extended double e type
952 * asctoe64 (string, &d) ASCII string to long double
953 * asctoe53 (string, &d) ASCII string to double
954 * asctoe24 (string, &f) ASCII string to single
955 * asctoeg (string, e, prec) ASCII string to specified precision
956 * e24toe (&f, e) IEEE single precision to e type
957 * e53toe (&d, e) IEEE double precision to e type
958 * e64toe (&d, e) IEEE long double precision to e type
959 * eabs (e) absolute value
960 * eadd (a, b, c) c = b + a
961 * eclear (e) e = 0
962 * ecmp (a, b) Returns 1 if a > b, 0 if a == b,
963 * -1 if a < b, -2 if either a or b is a NaN.
964 * ediv (a, b, c) c = b / a
965 * efloor (a, b) truncate to integer, toward -infinity
966 * efrexp (a, exp, s) extract exponent and significand
967 * eifrac (e, &l, frac) e to long integer and e type fraction
968 * euifrac (e, &l, frac) e to unsigned long integer and e type fraction
969 * einfin (e) set e to infinity, leaving its sign alone
970 * eldexp (a, n, b) multiply by 2**n
971 * emov (a, b) b = a
972 * emul (a, b, c) c = b * a
973 * eneg (e) e = -e
974 * eround (a, b) b = nearest integer value to a
975 * esub (a, b, c) c = b - a
976 * e24toasc (&f, str, n) single to ASCII string, n digits after decimal
977 * e53toasc (&d, str, n) double to ASCII string, n digits after decimal
978 * e64toasc (&d, str, n) long double to ASCII string
979 * etoasc (e, str, n) e to ASCII string, n digits after decimal
980 * etoe24 (e, &f) convert e type to IEEE single precision
981 * etoe53 (e, &d) convert e type to IEEE double precision
982 * etoe64 (e, &d) convert e type to IEEE long double precision
983 * ltoe (&l, e) long (32 bit) integer to e type
984 * ultoe (&l, e) unsigned long (32 bit) integer to e type
985 * eisneg (e) 1 if sign bit of e != 0, else 0
986 * eisinf (e) 1 if e has maximum exponent (non-IEEE)
987 * or is infinite (IEEE)
988 * eisnan (e) 1 if e is a NaN
989 *
990 *
991 * Routines for internal format numbers
992 *
993 * eaddm (ai, bi) add significands, bi = bi + ai
994 * ecleaz (ei) ei = 0
995 * ecleazs (ei) set ei = 0 but leave its sign alone
996 * ecmpm (ai, bi) compare significands, return 1, 0, or -1
997 * edivm (ai, bi) divide significands, bi = bi / ai
998 * emdnorm (ai,l,s,exp) normalize and round off
999 * emovi (a, ai) convert external a to internal ai
1000 * emovo (ai, a) convert internal ai to external a
1001 * emovz (ai, bi) bi = ai, low guard word of bi = 0
1002 * emulm (ai, bi) multiply significands, bi = bi * ai
1003 * enormlz (ei) left-justify the significand
1004 * eshdn1 (ai) shift significand and guards down 1 bit
1005 * eshdn8 (ai) shift down 8 bits
1006 * eshdn6 (ai) shift down 16 bits
1007 * eshift (ai, n) shift ai n bits up (or down if n < 0)
1008 * eshup1 (ai) shift significand and guards up 1 bit
1009 * eshup8 (ai) shift up 8 bits
1010 * eshup6 (ai) shift up 16 bits
1011 * esubm (ai, bi) subtract significands, bi = bi - ai
1012 * eiisinf (ai) 1 if infinite
1013 * eiisnan (ai) 1 if a NaN
1014 * einan (ai) set ai = NaN
1015 * eiinfin (ai) set ai = infinity
1016 *
1017 *
1018 * The result is always normalized and rounded to NI-4 word precision
1019 * after each arithmetic operation.
1020 *
1021 * Exception flags are NOT fully supported.
1022 *
1023 * Signaling NaN's are NOT supported; they are treated the same
1024 * as quiet NaN's.
1025 *
1026 * Define INFINITY for support of infinity; otherwise a
1027 * saturation arithmetic is implemented.
1028 *
1029 * Define NANS for support of Not-a-Number items; otherwise the
1030 * arithmetic will never produce a NaN output, and might be confused
1031 * by a NaN input.
1032 * If NaN's are supported, the output of `ecmp (a,b)' is -2 if
1033 * either a or b is a NaN. This means asking `if (ecmp (a,b) < 0)'
1034 * may not be legitimate. Use `if (ecmp (a,b) == -1)' for `less than'
1035 * if in doubt.
1036 *
1037 * Denormals are always supported here where appropriate (e.g., not
1038 * for conversion to DEC numbers).
1039 *
1040 */
1041
1042
1043/* mconf.h
1044 *
1045 * Common include file for math routines
1046 *
1047 *
1048 *
1049 * SYNOPSIS:
1050 *
1051 * #include "mconf.h"
1052 *
1053 *
1054 *
1055 * DESCRIPTION:
1056 *
1057 * This file contains definitions for error codes that are
1058 * passed to the common error handling routine mtherr
1059 * (which see).
1060 *
1061 * The file also includes a conditional assembly definition
1062 * for the type of computer arithmetic (Intel IEEE, DEC, Motorola
1063 * IEEE, or UNKnown).
1064 *
1065 * For Digital Equipment PDP-11 and VAX computers, certain
1066 * IBM systems, and others that use numbers with a 56-bit
1067 * significand, the symbol DEC should be defined. In this
1068 * mode, most floating point constants are given as arrays
1069 * of octal integers to eliminate decimal to binary conversion
1070 * errors that might be introduced by the compiler.
1071 *
1072 * For computers, such as IBM PC, that follow the IEEE
1073 * Standard for Binary Floating Point Arithmetic (ANSI/IEEE
1074 * Std 754-1985), the symbol IBMPC or MIEEE should be defined.
1075 * These numbers have 53-bit significands. In this mode, constants
1076 * are provided as arrays of hexadecimal 16 bit integers.
1077 *
1078 * To accommodate other types of computer arithmetic, all
1079 * constants are also provided in a normal decimal radix
1080 * which one can hope are correctly converted to a suitable
1081 * format by the available C language compiler. To invoke
1082 * this mode, the symbol UNK is defined.
1083 *
1084 * An important difference among these modes is a predefined
1085 * set of machine arithmetic constants for each. The numbers
1086 * MACHEP (the machine roundoff error), MAXNUM (largest number
1087 * represented), and several other parameters are preset by
1088 * the configuration symbol. Check the file const.c to
1089 * ensure that these values are correct for your computer.
1090 *
1091 * For ANSI C compatibility, define ANSIC equal to 1. Currently
1092 * this affects only the atan2 function and others that use it.
1093 */
1094\f
1095/* Constant definitions for math error conditions. */
1096
1097#define DOMAIN 1 /* argument domain error */
1098#define SING 2 /* argument singularity */
1099#define OVERFLOW 3 /* overflow range error */
1100#define UNDERFLOW 4 /* underflow range error */
1101#define TLOSS 5 /* total loss of precision */
1102#define PLOSS 6 /* partial loss of precision */
1103#define INVALID 7 /* NaN-producing operation */
1104
1105/* e type constants used by high precision check routines */
1106
1107/*include "ehead.h"*/
1108/* 0.0 */
1109unsigned EMUSHORT ezero[NE] =
1110{
1111 0, 0000000, 0000000, 0000000, 0000000, 0000000,};
1112extern unsigned EMUSHORT ezero[];
1113
1114/* 5.0E-1 */
1115unsigned EMUSHORT ehalf[NE] =
1116{
1117 0, 0000000, 0000000, 0000000, 0100000, 0x3ffe,};
1118extern unsigned EMUSHORT ehalf[];
1119
1120/* 1.0E0 */
1121unsigned EMUSHORT eone[NE] =
1122{
1123 0, 0000000, 0000000, 0000000, 0100000, 0x3fff,};
1124extern unsigned EMUSHORT eone[];
1125
1126/* 2.0E0 */
1127unsigned EMUSHORT etwo[NE] =
1128{
1129 0, 0000000, 0000000, 0000000, 0100000, 0040000,};
1130extern unsigned EMUSHORT etwo[];
1131
1132/* 3.2E1 */
1133unsigned EMUSHORT e32[NE] =
1134{
1135 0, 0000000, 0000000, 0000000, 0100000, 0040004,};
1136extern unsigned EMUSHORT e32[];
1137
1138/* 6.93147180559945309417232121458176568075500134360255E-1 */
1139unsigned EMUSHORT elog2[NE] =
1140{
1141 0xc9e4, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1142extern unsigned EMUSHORT elog2[];
1143
1144/* 1.41421356237309504880168872420969807856967187537695E0 */
1145unsigned EMUSHORT esqrt2[NE] =
1146{
1147 0x597e, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1148extern unsigned EMUSHORT esqrt2[];
1149
1150/* 2/sqrt (PI) =
1151 * 1.12837916709551257389615890312154517168810125865800E0 */
1152unsigned EMUSHORT eoneopi[NE] =
1153{
1154 0x71d5, 0x688d, 0012333, 0135202, 0110156, 0x3fff,};
1155extern unsigned EMUSHORT eoneopi[];
1156
1157/* 3.14159265358979323846264338327950288419716939937511E0 */
1158unsigned EMUSHORT epi[NE] =
1159{
1160 0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1161extern unsigned EMUSHORT epi[];
1162
1163/* 5.7721566490153286060651209008240243104215933593992E-1 */
1164unsigned EMUSHORT eeul[NE] =
1165{
1166 0xd1be, 0xc7a4, 0076660, 0063743, 0111704, 0x3ffe,};
1167extern unsigned EMUSHORT eeul[];
1168
1169/*
1170include "ehead.h"
1171include "mconf.h"
1172*/
1173
1174
1175
1176/* Control register for rounding precision.
1177 * This can be set to 80 (if NE=6), 64, 56, 53, or 24 bits.
1178 */
1179int rndprc = NBITS;
1180extern int rndprc;
1181
1182void eaddm (), esubm (), emdnorm (), asctoeg ();
1183static void toe24 (), toe53 (), toe64 ();
1184void eremain (), einit (), eiremain ();
1185int ecmpm (), edivm (), emulm ();
1186void emovi (), emovo (), emovz (), ecleaz (), ecleazs (), eadd1 ();
1187void etodec (), todec (), dectoe ();
1188
1189
1190
1191
1192void
1193einit ()
1194{
1195}
1196
1197/*
1198; Clear out entire external format number.
1199;
1200; unsigned EMUSHORT x[];
1201; eclear (x);
1202*/
1203
1204void
1205eclear (x)
1206 register unsigned EMUSHORT *x;
1207{
1208 register int i;
1209
1210 for (i = 0; i < NE; i++)
1211 *x++ = 0;
1212}
1213
1214
1215
1216/* Move external format number from a to b.
1217 *
1218 * emov (a, b);
1219 */
1220
1221void
1222emov (a, b)
1223 register unsigned EMUSHORT *a, *b;
1224{
1225 register int i;
1226
1227 for (i = 0; i < NE; i++)
1228 *b++ = *a++;
1229}
1230
1231
1232/*
1233; Absolute value of external format number
1234;
1235; EMUSHORT x[NE];
1236; eabs (x);
1237*/
1238
1239void
1240eabs (x)
1241 unsigned EMUSHORT x[]; /* x is the memory address of a short */
1242{
1243
1244 x[NE - 1] &= 0x7fff; /* sign is top bit of last word of external format */
1245}
1246
1247
1248
1249
1250/*
1251; Negate external format number
1252;
1253; unsigned EMUSHORT x[NE];
1254; eneg (x);
1255*/
1256
1257void
1258eneg (x)
1259 unsigned EMUSHORT x[];
1260{
1261
1262#ifdef NANS
1263 if (eisnan (x))
1264 return;
1265#endif
1266 x[NE - 1] ^= 0x8000; /* Toggle the sign bit */
1267}
1268
1269
1270
1271/* Return 1 if external format number is negative,
1272 * else return zero, including when it is a NaN.
1273 */
1274int
1275eisneg (x)
1276 unsigned EMUSHORT x[];
1277{
1278
1279#ifdef NANS
1280 if (eisnan (x))
1281 return (0);
1282#endif
1283 if (x[NE - 1] & 0x8000)
1284 return (1);
1285 else
1286 return (0);
1287}
1288
1289
1290/* Return 1 if external format number is infinity.
1291 * else return zero.
1292 */
1293int
1294eisinf (x)
1295 unsigned EMUSHORT x[];
1296{
1297
1298#ifdef NANS
1299 if (eisnan (x))
1300 return (0);
1301#endif
1302 if ((x[NE - 1] & 0x7fff) == 0x7fff)
1303 return (1);
1304 else
1305 return (0);
1306}
1307
1308
1309/* Check if e-type number is not a number.
1310 The bit pattern is one that we defined, so we know for sure how to
1311 detect it. */
1312
1313int
1314eisnan (x)
1315 unsigned EMUSHORT x[];
1316{
1317
1318#ifdef NANS
1319 int i;
1320/* NaN has maximum exponent */
1321 if ((x[NE - 1] & 0x7fff) != 0x7fff)
1322 return (0);
1323/* ... and non-zero significand field. */
1324 for (i = 0; i < NE - 1; i++)
1325 {
1326 if (*x++ != 0)
1327 return (1);
1328 }
1329#endif
1330 return (0);
1331}
1332
1333/* Fill external format number with infinity pattern (IEEE)
1334 or largest possible number (non-IEEE).
1335 Before calling einfin, you should either call eclear
1336 or set up the sign bit by hand. */
1337
1338void
1339einfin (x)
1340 register unsigned EMUSHORT *x;
1341{
1342 register int i;
1343
1344#ifdef INFINITY
1345 for (i = 0; i < NE - 1; i++)
1346 *x++ = 0;
1347 *x |= 32767;
1348#else
1349 for (i = 0; i < NE - 1; i++)
1350 *x++ = 0xffff;
1351 *x |= 32766;
1352 if (rndprc < NBITS)
1353 {
1354 if (rndprc == 64)
1355 {
1356 *(x - 5) = 0;
1357 }
1358 if (rndprc == 53)
1359 {
1360 *(x - 4) = 0xf800;
1361 }
1362 else
1363 {
1364 *(x - 4) = 0;
1365 *(x - 3) = 0;
1366 *(x - 2) = 0xff00;
1367 }
1368 }
1369#endif
1370}
1371
1372
1373/* Output an e-type NaN.
1374 This generates Intel's quiet NaN pattern for extended real.
1375 The exponent is 7fff, the leading mantissa word is c000. */
1376
1377void
1378enan (x)
1379 register unsigned EMUSHORT *x;
1380{
1381 register int i;
1382
1383 for (i = 0; i < NE - 2; i++)
1384 *x++ = 0;
1385 *x++ = 0xc000;
1386 *x = 0x7fff;
1387}
1388
1389
1390/* Move in external format number,
1391 * converting it to internal format.
1392 */
1393void
1394emovi (a, b)
1395 unsigned EMUSHORT *a, *b;
1396{
1397 register unsigned EMUSHORT *p, *q;
1398 int i;
1399
1400 q = b;
1401 p = a + (NE - 1); /* point to last word of external number */
1402 /* get the sign bit */
1403 if (*p & 0x8000)
1404 *q++ = 0xffff;
1405 else
1406 *q++ = 0;
1407 /* get the exponent */
1408 *q = *p--;
1409 *q++ &= 0x7fff; /* delete the sign bit */
1410#ifdef INFINITY
1411 if ((*(q - 1) & 0x7fff) == 0x7fff)
1412 {
1413#ifdef NANS
1414 if (eisnan (a))
1415 {
1416 *q++ = 0;
1417 for (i = 3; i < NI; i++)
1418 *q++ = *p--;
1419 return;
1420 }
1421#endif
1422 for (i = 2; i < NI; i++)
1423 *q++ = 0;
1424 return;
1425 }
1426#endif
1427 /* clear high guard word */
1428 *q++ = 0;
1429 /* move in the significand */
1430 for (i = 0; i < NE - 1; i++)
1431 *q++ = *p--;
1432 /* clear low guard word */
1433 *q = 0;
1434}
1435
1436
1437/* Move internal format number out,
1438 * converting it to external format.
1439 */
1440void
1441emovo (a, b)
1442 unsigned EMUSHORT *a, *b;
1443{
1444 register unsigned EMUSHORT *p, *q;
1445 unsigned EMUSHORT i;
1446
1447 p = a;
1448 q = b + (NE - 1); /* point to output exponent */
1449 /* combine sign and exponent */
1450 i = *p++;
1451 if (i)
1452 *q-- = *p++ | 0x8000;
1453 else
1454 *q-- = *p++;
1455#ifdef INFINITY
1456 if (*(p - 1) == 0x7fff)
1457 {
1458#ifdef NANS
1459 if (eiisnan (a))
1460 {
1461 enan (b);
1462 return;
1463 }
1464#endif
1465 einfin (b);
1466 return;
1467 }
1468#endif
1469 /* skip over guard word */
1470 ++p;
1471 /* move the significand */
1472 for (i = 0; i < NE - 1; i++)
1473 *q-- = *p++;
1474}
1475
1476
1477
1478
1479/* Clear out internal format number.
1480 */
1481
1482void
1483ecleaz (xi)
1484 register unsigned EMUSHORT *xi;
1485{
1486 register int i;
1487
1488 for (i = 0; i < NI; i++)
1489 *xi++ = 0;
1490}
1491
1492
1493/* same, but don't touch the sign. */
1494
1495void
1496ecleazs (xi)
1497 register unsigned EMUSHORT *xi;
1498{
1499 register int i;
1500
1501 ++xi;
1502 for (i = 0; i < NI - 1; i++)
1503 *xi++ = 0;
1504}
1505
1506
1507
1508/* Move internal format number from a to b.
1509 */
1510void
1511emovz (a, b)
1512 register unsigned EMUSHORT *a, *b;
1513{
1514 register int i;
1515
1516 for (i = 0; i < NI - 1; i++)
1517 *b++ = *a++;
1518 /* clear low guard word */
1519 *b = 0;
1520}
1521
1522/* Generate internal format NaN.
1523 The explicit pattern for this is maximum exponent and
1524 top two significand bits set. */
1525
1526void
1527einan (x)
1528 unsigned EMUSHORT x[];
1529{
1530
1531 ecleaz (x);
1532 x[E] = 0x7fff;
1533 x[M + 1] = 0xc000;
1534}
1535
1536/* Return nonzero if internal format number is a NaN. */
1537
1538int
1539eiisnan (x)
1540 unsigned EMUSHORT x[];
1541{
1542 int i;
1543
1544 if ((x[E] & 0x7fff) == 0x7fff)
1545 {
1546 for (i = M + 1; i < NI; i++)
1547 {
1548 if (x[i] != 0)
1549 return (1);
1550 }
1551 }
1552 return (0);
1553}
1554
1555/* Fill internal format number with infinity pattern.
1556 This has maximum exponent and significand all zeros. */
1557
1558void
1559eiinfin (x)
1560 unsigned EMUSHORT x[];
1561{
1562
1563 ecleaz (x);
1564 x[E] = 0x7fff;
1565}
1566
1567/* Return nonzero if internal format number is infinite. */
1568
1569int
1570eiisinf (x)
1571 unsigned EMUSHORT x[];
1572{
1573
1574#ifdef NANS
1575 if (eiisnan (x))
1576 return (0);
1577#endif
1578 if ((x[E] & 0x7fff) == 0x7fff)
1579 return (1);
1580 return (0);
1581}
1582
1583
1584/*
1585; Compare significands of numbers in internal format.
1586; Guard words are included in the comparison.
1587;
1588; unsigned EMUSHORT a[NI], b[NI];
1589; cmpm (a, b);
1590;
1591; for the significands:
1592; returns +1 if a > b
1593; 0 if a == b
1594; -1 if a < b
1595*/
1596int
1597ecmpm (a, b)
1598 register unsigned EMUSHORT *a, *b;
1599{
1600 int i;
1601
1602 a += M; /* skip up to significand area */
1603 b += M;
1604 for (i = M; i < NI; i++)
1605 {
1606 if (*a++ != *b++)
1607 goto difrnt;
1608 }
1609 return (0);
1610
1611 difrnt:
1612 if (*(--a) > *(--b))
1613 return (1);
1614 else
1615 return (-1);
1616}
1617
1618
1619/*
1620; Shift significand down by 1 bit
1621*/
1622
1623void
1624eshdn1 (x)
1625 register unsigned EMUSHORT *x;
1626{
1627 register unsigned EMUSHORT bits;
1628 int i;
1629
1630 x += M; /* point to significand area */
1631
1632 bits = 0;
1633 for (i = M; i < NI; i++)
1634 {
1635 if (*x & 1)
1636 bits |= 1;
1637 *x >>= 1;
1638 if (bits & 2)
1639 *x |= 0x8000;
1640 bits <<= 1;
1641 ++x;
1642 }
1643}
1644
1645
1646
1647/*
1648; Shift significand up by 1 bit
1649*/
1650
1651void
1652eshup1 (x)
1653 register unsigned EMUSHORT *x;
1654{
1655 register unsigned EMUSHORT bits;
1656 int i;
1657
1658 x += NI - 1;
1659 bits = 0;
1660
1661 for (i = M; i < NI; i++)
1662 {
1663 if (*x & 0x8000)
1664 bits |= 1;
1665 *x <<= 1;
1666 if (bits & 2)
1667 *x |= 1;
1668 bits <<= 1;
1669 --x;
1670 }
1671}
1672
1673
1674
1675/*
1676; Shift significand down by 8 bits
1677*/
1678
1679void
1680eshdn8 (x)
1681 register unsigned EMUSHORT *x;
1682{
1683 register unsigned EMUSHORT newbyt, oldbyt;
1684 int i;
1685
1686 x += M;
1687 oldbyt = 0;
1688 for (i = M; i < NI; i++)
1689 {
1690 newbyt = *x << 8;
1691 *x >>= 8;
1692 *x |= oldbyt;
1693 oldbyt = newbyt;
1694 ++x;
1695 }
1696}
1697
1698/*
1699; Shift significand up by 8 bits
1700*/
1701
1702void
1703eshup8 (x)
1704 register unsigned EMUSHORT *x;
1705{
1706 int i;
1707 register unsigned EMUSHORT newbyt, oldbyt;
1708
1709 x += NI - 1;
1710 oldbyt = 0;
1711
1712 for (i = M; i < NI; i++)
1713 {
1714 newbyt = *x >> 8;
1715 *x <<= 8;
1716 *x |= oldbyt;
1717 oldbyt = newbyt;
1718 --x;
1719 }
1720}
1721
1722/*
1723; Shift significand up by 16 bits
1724*/
1725
1726void
1727eshup6 (x)
1728 register unsigned EMUSHORT *x;
1729{
1730 int i;
1731 register unsigned EMUSHORT *p;
1732
1733 p = x + M;
1734 x += M + 1;
1735
1736 for (i = M; i < NI - 1; i++)
1737 *p++ = *x++;
1738
1739 *p = 0;
1740}
1741
1742/*
1743; Shift significand down by 16 bits
1744*/
1745
1746void
1747eshdn6 (x)
1748 register unsigned EMUSHORT *x;
1749{
1750 int i;
1751 register unsigned EMUSHORT *p;
1752
1753 x += NI - 1;
1754 p = x + 1;
1755
1756 for (i = M; i < NI - 1; i++)
1757 *(--p) = *(--x);
1758
1759 *(--p) = 0;
1760}
1761\f
1762/*
1763; Add significands
1764; x + y replaces y
1765*/
1766
1767void
1768eaddm (x, y)
1769 unsigned EMUSHORT *x, *y;
1770{
1771 register unsigned EMULONG a;
1772 int i;
1773 unsigned int carry;
1774
1775 x += NI - 1;
1776 y += NI - 1;
1777 carry = 0;
1778 for (i = M; i < NI; i++)
1779 {
1780 a = (unsigned EMULONG) (*x) + (unsigned EMULONG) (*y) + carry;
1781 if (a & 0x10000)
1782 carry = 1;
1783 else
1784 carry = 0;
1785 *y = (unsigned EMUSHORT) a;
1786 --x;
1787 --y;
1788 }
1789}
1790
1791/*
1792; Subtract significands
1793; y - x replaces y
1794*/
1795
1796void
1797esubm (x, y)
1798 unsigned EMUSHORT *x, *y;
1799{
1800 unsigned EMULONG a;
1801 int i;
1802 unsigned int carry;
1803
1804 x += NI - 1;
1805 y += NI - 1;
1806 carry = 0;
1807 for (i = M; i < NI; i++)
1808 {
1809 a = (unsigned EMULONG) (*y) - (unsigned EMULONG) (*x) - carry;
1810 if (a & 0x10000)
1811 carry = 1;
1812 else
1813 carry = 0;
1814 *y = (unsigned EMUSHORT) a;
1815 --x;
1816 --y;
1817 }
1818}
1819
1820
1821/* Divide significands */
1822
1823static unsigned EMUSHORT equot[NI];
1824
1825int
1826edivm (den, num)
1827 unsigned EMUSHORT den[], num[];
1828{
1829 int i;
1830 register unsigned EMUSHORT *p, *q;
1831 unsigned EMUSHORT j;
1832
1833 p = &equot[0];
1834 *p++ = num[0];
1835 *p++ = num[1];
1836
1837 for (i = M; i < NI; i++)
1838 {
1839 *p++ = 0;
1840 }
1841
1842 /* Use faster compare and subtraction if denominator
1843 * has only 15 bits of significance.
1844 */
1845 p = &den[M + 2];
1846 if (*p++ == 0)
1847 {
1848 for (i = M + 3; i < NI; i++)
1849 {
1850 if (*p++ != 0)
1851 goto fulldiv;
1852 }
1853 if ((den[M + 1] & 1) != 0)
1854 goto fulldiv;
1855 eshdn1 (num);
1856 eshdn1 (den);
1857
1858 p = &den[M + 1];
1859 q = &num[M + 1];
1860
1861 for (i = 0; i < NBITS + 2; i++)
1862 {
1863 if (*p <= *q)
1864 {
1865 *q -= *p;
1866 j = 1;
1867 }
1868 else
1869 {
1870 j = 0;
1871 }
1872 eshup1 (equot);
1873 equot[NI - 2] |= j;
1874 eshup1 (num);
1875 }
1876 goto divdon;
1877 }
1878
1879 /* The number of quotient bits to calculate is
1880 * NBITS + 1 scaling guard bit + 1 roundoff bit.
1881 */
1882 fulldiv:
1883
1884 p = &equot[NI - 2];
1885 for (i = 0; i < NBITS + 2; i++)
1886 {
1887 if (ecmpm (den, num) <= 0)
1888 {
1889 esubm (den, num);
1890 j = 1; /* quotient bit = 1 */
1891 }
1892 else
1893 j = 0;
1894 eshup1 (equot);
1895 *p |= j;
1896 eshup1 (num);
1897 }
1898
1899 divdon:
1900
1901 eshdn1 (equot);
1902 eshdn1 (equot);
1903
1904 /* test for nonzero remainder after roundoff bit */
1905 p = &num[M];
1906 j = 0;
1907 for (i = M; i < NI; i++)
1908 {
1909 j |= *p++;
1910 }
1911 if (j)
1912 j = 1;
1913
1914
1915 for (i = 0; i < NI; i++)
1916 num[i] = equot[i];
1917 return ((int) j);
1918}
1919
1920
1921/* Multiply significands */
1922int
1923emulm (a, b)
1924 unsigned EMUSHORT a[], b[];
1925{
1926 unsigned EMUSHORT *p, *q;
1927 int i, j, k;
1928
1929 equot[0] = b[0];
1930 equot[1] = b[1];
1931 for (i = M; i < NI; i++)
1932 equot[i] = 0;
1933
1934 p = &a[NI - 2];
1935 k = NBITS;
1936 while (*p == 0) /* significand is not supposed to be all zero */
1937 {
1938 eshdn6 (a);
1939 k -= 16;
1940 }
1941 if ((*p & 0xff) == 0)
1942 {
1943 eshdn8 (a);
1944 k -= 8;
1945 }
1946
1947 q = &equot[NI - 1];
1948 j = 0;
1949 for (i = 0; i < k; i++)
1950 {
1951 if (*p & 1)
1952 eaddm (b, equot);
1953 /* remember if there were any nonzero bits shifted out */
1954 if (*q & 1)
1955 j |= 1;
1956 eshdn1 (a);
1957 eshdn1 (equot);
1958 }
1959
1960 for (i = 0; i < NI; i++)
1961 b[i] = equot[i];
1962
1963 /* return flag for lost nonzero bits */
1964 return (j);
1965}
1966
1967
1968
1969/*
1970 * Normalize and round off.
1971 *
1972 * The internal format number to be rounded is "s".
1973 * Input "lost" indicates whether or not the number is exact.
1974 * This is the so-called sticky bit.
1975 *
1976 * Input "subflg" indicates whether the number was obtained
1977 * by a subtraction operation. In that case if lost is nonzero
1978 * then the number is slightly smaller than indicated.
1979 *
1980 * Input "exp" is the biased exponent, which may be negative.
1981 * the exponent field of "s" is ignored but is replaced by
1982 * "exp" as adjusted by normalization and rounding.
1983 *
1984 * Input "rcntrl" is the rounding control.
1985 */
1986
1987static int rlast = -1;
1988static int rw = 0;
1989static unsigned EMUSHORT rmsk = 0;
1990static unsigned EMUSHORT rmbit = 0;
1991static unsigned EMUSHORT rebit = 0;
1992static int re = 0;
1993static unsigned EMUSHORT rbit[NI];
1994
1995void
1996emdnorm (s, lost, subflg, exp, rcntrl)
1997 unsigned EMUSHORT s[];
1998 int lost;
1999 int subflg;
2000 EMULONG exp;
2001 int rcntrl;
2002{
2003 int i, j;
2004 unsigned EMUSHORT r;
2005
2006 /* Normalize */
2007 j = enormlz (s);
2008
2009 /* a blank significand could mean either zero or infinity. */
2010#ifndef INFINITY
2011 if (j > NBITS)
2012 {
2013 ecleazs (s);
2014 return;
2015 }
2016#endif
2017 exp -= j;
2018#ifndef INFINITY
2019 if (exp >= 32767L)
2020 goto overf;
2021#else
2022 if ((j > NBITS) && (exp < 32767))
2023 {
2024 ecleazs (s);
2025 return;
2026 }
2027#endif
2028 if (exp < 0L)
2029 {
2030 if (exp > (EMULONG) (-NBITS - 1))
2031 {
2032 j = (int) exp;
2033 i = eshift (s, j);
2034 if (i)
2035 lost = 1;
2036 }
2037 else
2038 {
2039 ecleazs (s);
2040 return;
2041 }
2042 }
2043 /* Round off, unless told not to by rcntrl. */
2044 if (rcntrl == 0)
2045 goto mdfin;
2046 /* Set up rounding parameters if the control register changed. */
2047 if (rndprc != rlast)
2048 {
2049 ecleaz (rbit);
2050 switch (rndprc)
2051 {
2052 default:
2053 case NBITS:
2054 rw = NI - 1; /* low guard word */
2055 rmsk = 0xffff;
2056 rmbit = 0x8000;
2057 rbit[rw - 1] = 1;
2058 re = NI - 2;
2059 rebit = 1;
2060 break;
2061 case 64:
2062 rw = 7;
2063 rmsk = 0xffff;
2064 rmbit = 0x8000;
2065 rbit[rw - 1] = 1;
2066 re = rw - 1;
2067 rebit = 1;
2068 break;
2069 /* For DEC arithmetic */
2070 case 56:
2071 rw = 6;
2072 rmsk = 0xff;
2073 rmbit = 0x80;
2074 rbit[rw] = 0x100;
2075 re = rw;
2076 rebit = 0x100;
2077 break;
2078 case 53:
2079 rw = 6;
2080 rmsk = 0x7ff;
2081 rmbit = 0x0400;
2082 rbit[rw] = 0x800;
2083 re = rw;
2084 rebit = 0x800;
2085 break;
2086 case 24:
2087 rw = 4;
2088 rmsk = 0xff;
2089 rmbit = 0x80;
2090 rbit[rw] = 0x100;
2091 re = rw;
2092 rebit = 0x100;
2093 break;
2094 }
2095 rlast = rndprc;
2096 }
2097
2098 if (rndprc >= 64)
2099 {
2100 r = s[rw] & rmsk;
2101 if (rndprc == 64)
2102 {
2103 i = rw + 1;
2104 while (i < NI)
2105 {
2106 if (s[i])
2107 r |= 1;
2108 s[i] = 0;
2109 ++i;
2110 }
2111 }
2112 }
2113 else
2114 {
2115 if (exp <= 0)
2116 eshdn1 (s);
2117 r = s[rw] & rmsk;
2118 /* These tests assume NI = 8 */
2119 i = rw + 1;
2120 while (i < NI)
2121 {
2122 if (s[i])
2123 r |= 1;
2124 s[i] = 0;
2125 ++i;
2126 }
2127 /*
2128 if (rndprc == 24)
2129 {
2130 if (s[5] || s[6])
2131 r |= 1;
2132 s[5] = 0;
2133 s[6] = 0;
2134 }
2135 */
2136 }
2137 s[rw] &= ~rmsk;
2138 if ((r & rmbit) != 0)
2139 {
2140 if (r == rmbit)
2141 {
2142 if (lost == 0)
2143 { /* round to even */
2144 if ((s[re] & rebit) == 0)
2145 goto mddone;
2146 }
2147 else
2148 {
2149 if (subflg != 0)
2150 goto mddone;
2151 }
2152 }
2153 eaddm (rbit, s);
2154 }
2155 mddone:
2156 if ((rndprc < 64) && (exp <= 0))
2157 {
2158 eshup1 (s);
2159 }
2160 if (s[2] != 0)
2161 { /* overflow on roundoff */
2162 eshdn1 (s);
2163 exp += 1;
2164 }
2165 mdfin:
2166 s[NI - 1] = 0;
2167 if (exp >= 32767L)
2168 {
2169#ifndef INFINITY
2170 overf:
2171#endif
2172#ifdef INFINITY
2173 s[1] = 32767;
2174 for (i = 2; i < NI - 1; i++)
2175 s[i] = 0;
2176 if (extra_warnings)
2177 warning ("floating point overflow");
2178#else
2179 s[1] = 32766;
2180 s[2] = 0;
2181 for (i = M + 1; i < NI - 1; i++)
2182 s[i] = 0xffff;
2183 s[NI - 1] = 0;
2184 if (rndprc < 64)
2185 {
2186 s[rw] &= ~rmsk;
2187 if (rndprc == 24)
2188 {
2189 s[5] = 0;
2190 s[6] = 0;
2191 }
2192 }
2193#endif
2194 return;
2195 }
2196 if (exp < 0)
2197 s[1] = 0;
2198 else
2199 s[1] = (unsigned EMUSHORT) exp;
2200}
2201
2202
2203
2204/*
2205; Subtract external format numbers.
2206;
2207; unsigned EMUSHORT a[NE], b[NE], c[NE];
2208; esub (a, b, c); c = b - a
2209*/
2210
2211static int subflg = 0;
2212
2213void
2214esub (a, b, c)
2215 unsigned EMUSHORT *a, *b, *c;
2216{
2217
2218#ifdef NANS
2219 if (eisnan (a))
2220 {
2221 emov (a, c);
2222 return;
2223 }
2224 if (eisnan (b))
2225 {
2226 emov (b, c);
2227 return;
2228 }
2229/* Infinity minus infinity is a NaN.
2230 Test for subtracting infinities of the same sign. */
2231 if (eisinf (a) && eisinf (b)
2232 && ((eisneg (a) ^ eisneg (b)) == 0))
2233 {
2234 mtherr ("esub", INVALID);
2235 enan (c);
2236 return;
2237 }
2238#endif
2239 subflg = 1;
2240 eadd1 (a, b, c);
2241}
2242
2243
2244/*
2245; Add.
2246;
2247; unsigned EMUSHORT a[NE], b[NE], c[NE];
2248; eadd (a, b, c); c = b + a
2249*/
2250void
2251eadd (a, b, c)
2252 unsigned EMUSHORT *a, *b, *c;
2253{
2254
2255#ifdef NANS
2256/* NaN plus anything is a NaN. */
2257 if (eisnan (a))
2258 {
2259 emov (a, c);
2260 return;
2261 }
2262 if (eisnan (b))
2263 {
2264 emov (b, c);
2265 return;
2266 }
2267/* Infinity minus infinity is a NaN.
2268 Test for adding infinities of opposite signs. */
2269 if (eisinf (a) && eisinf (b)
2270 && ((eisneg (a) ^ eisneg (b)) != 0))
2271 {
2272 mtherr ("esub", INVALID);
2273 enan (c);
2274 return;
2275 }
2276#endif
2277 subflg = 0;
2278 eadd1 (a, b, c);
2279}
2280
2281void
2282eadd1 (a, b, c)
2283 unsigned EMUSHORT *a, *b, *c;
2284{
2285 unsigned EMUSHORT ai[NI], bi[NI], ci[NI];
2286 int i, lost, j, k;
2287 EMULONG lt, lta, ltb;
2288
2289#ifdef INFINITY
2290 if (eisinf (a))
2291 {
2292 emov (a, c);
2293 if (subflg)
2294 eneg (c);
2295 return;
2296 }
2297 if (eisinf (b))
2298 {
2299 emov (b, c);
2300 return;
2301 }
2302#endif
2303 emovi (a, ai);
2304 emovi (b, bi);
2305 if (subflg)
2306 ai[0] = ~ai[0];
2307
2308 /* compare exponents */
2309 lta = ai[E];
2310 ltb = bi[E];
2311 lt = lta - ltb;
2312 if (lt > 0L)
2313 { /* put the larger number in bi */
2314 emovz (bi, ci);
2315 emovz (ai, bi);
2316 emovz (ci, ai);
2317 ltb = bi[E];
2318 lt = -lt;
2319 }
2320 lost = 0;
2321 if (lt != 0L)
2322 {
2323 if (lt < (EMULONG) (-NBITS - 1))
2324 goto done; /* answer same as larger addend */
2325 k = (int) lt;
2326 lost = eshift (ai, k); /* shift the smaller number down */
2327 }
2328 else
2329 {
2330 /* exponents were the same, so must compare significands */
2331 i = ecmpm (ai, bi);
2332 if (i == 0)
2333 { /* the numbers are identical in magnitude */
2334 /* if different signs, result is zero */
2335 if (ai[0] != bi[0])
2336 {
2337 eclear (c);
2338 return;
2339 }
2340 /* if same sign, result is double */
2341 /* double denomalized tiny number */
2342 if ((bi[E] == 0) && ((bi[3] & 0x8000) == 0))
2343 {
2344 eshup1 (bi);
2345 goto done;
2346 }
2347 /* add 1 to exponent unless both are zero! */
2348 for (j = 1; j < NI - 1; j++)
2349 {
2350 if (bi[j] != 0)
2351 {
2352 /* This could overflow, but let emovo take care of that. */
2353 ltb += 1;
2354 break;
2355 }
2356 }
2357 bi[E] = (unsigned EMUSHORT) ltb;
2358 goto done;
2359 }
2360 if (i > 0)
2361 { /* put the larger number in bi */
2362 emovz (bi, ci);
2363 emovz (ai, bi);
2364 emovz (ci, ai);
2365 }
2366 }
2367 if (ai[0] == bi[0])
2368 {
2369 eaddm (ai, bi);
2370 subflg = 0;
2371 }
2372 else
2373 {
2374 esubm (ai, bi);
2375 subflg = 1;
2376 }
2377 emdnorm (bi, lost, subflg, ltb, 64);
2378
2379 done:
2380 emovo (bi, c);
2381}
2382
2383
2384
2385/*
2386; Divide.
2387;
2388; unsigned EMUSHORT a[NE], b[NE], c[NE];
2389; ediv (a, b, c); c = b / a
2390*/
2391void
2392ediv (a, b, c)
2393 unsigned EMUSHORT *a, *b, *c;
2394{
2395 unsigned EMUSHORT ai[NI], bi[NI];
2396 int i;
2397 EMULONG lt, lta, ltb;
2398
2399#ifdef NANS
2400/* Return any NaN input. */
2401 if (eisnan (a))
2402 {
2403 emov (a, c);
2404 return;
2405 }
2406 if (eisnan (b))
2407 {
2408 emov (b, c);
2409 return;
2410 }
2411/* Zero over zero, or infinity over infinity, is a NaN. */
2412 if (((ecmp (a, ezero) == 0) && (ecmp (b, ezero) == 0))
2413 || (eisinf (a) && eisinf (b)))
2414 {
2415 mtherr ("ediv", INVALID);
2416 enan (c);
2417 return;
2418 }
2419#endif
2420/* Infinity over anything else is infinity. */
2421#ifdef INFINITY
2422 if (eisinf (b))
2423 {
2424 if (eisneg (a) ^ eisneg (b))
2425 *(c + (NE - 1)) = 0x8000;
2426 else
2427 *(c + (NE - 1)) = 0;
2428 einfin (c);
2429 return;
2430 }
2431/* Anything else over infinity is zero. */
2432 if (eisinf (a))
2433 {
2434 eclear (c);
2435 return;
2436 }
2437#endif
2438 emovi (a, ai);
2439 emovi (b, bi);
2440 lta = ai[E];
2441 ltb = bi[E];
2442 if (bi[E] == 0)
2443 { /* See if numerator is zero. */
2444 for (i = 1; i < NI - 1; i++)
2445 {
2446 if (bi[i] != 0)
2447 {
2448 ltb -= enormlz (bi);
2449 goto dnzro1;
2450 }
2451 }
2452 eclear (c);
2453 return;
2454 }
2455 dnzro1:
2456
2457 if (ai[E] == 0)
2458 { /* possible divide by zero */
2459 for (i = 1; i < NI - 1; i++)
2460 {
2461 if (ai[i] != 0)
2462 {
2463 lta -= enormlz (ai);
2464 goto dnzro2;
2465 }
2466 }
2467 if (ai[0] == bi[0])
2468 *(c + (NE - 1)) = 0;
2469 else
2470 *(c + (NE - 1)) = 0x8000;
2471/* Divide by zero is not an invalid operation.
2472 It is a divide-by-zero operation! */
2473 einfin (c);
2474 mtherr ("ediv", SING);
2475 return;
2476 }
2477 dnzro2:
2478
2479 i = edivm (ai, bi);
2480 /* calculate exponent */
2481 lt = ltb - lta + EXONE;
2482 emdnorm (bi, i, 0, lt, 64);
2483 /* set the sign */
2484 if (ai[0] == bi[0])
2485 bi[0] = 0;
2486 else
2487 bi[0] = 0Xffff;
2488 emovo (bi, c);
2489}
2490
2491
2492
2493/*
2494; Multiply.
2495;
2496; unsigned EMUSHORT a[NE], b[NE], c[NE];
2497; emul (a, b, c); c = b * a
2498*/
2499void
2500emul (a, b, c)
2501 unsigned EMUSHORT *a, *b, *c;
2502{
2503 unsigned EMUSHORT ai[NI], bi[NI];
2504 int i, j;
2505 EMULONG lt, lta, ltb;
2506
2507#ifdef NANS
2508/* NaN times anything is the same NaN. */
2509 if (eisnan (a))
2510 {
2511 emov (a, c);
2512 return;
2513 }
2514 if (eisnan (b))
2515 {
2516 emov (b, c);
2517 return;
2518 }
2519/* Zero times infinity is a NaN. */
2520 if ((eisinf (a) && (ecmp (b, ezero) == 0))
2521 || (eisinf (b) && (ecmp (a, ezero) == 0)))
2522 {
2523 mtherr ("emul", INVALID);
2524 enan (c);
2525 return;
2526 }
2527#endif
2528/* Infinity times anything else is infinity. */
2529#ifdef INFINITY
2530 if (eisinf (a) || eisinf (b))
2531 {
2532 if (eisneg (a) ^ eisneg (b))
2533 *(c + (NE - 1)) = 0x8000;
2534 else
2535 *(c + (NE - 1)) = 0;
2536 einfin (c);
2537 return;
2538 }
2539#endif
2540 emovi (a, ai);
2541 emovi (b, bi);
2542 lta = ai[E];
2543 ltb = bi[E];
2544 if (ai[E] == 0)
2545 {
2546 for (i = 1; i < NI - 1; i++)
2547 {
2548 if (ai[i] != 0)
2549 {
2550 lta -= enormlz (ai);
2551 goto mnzer1;
2552 }
2553 }
2554 eclear (c);
2555 return;
2556 }
2557 mnzer1:
2558
2559 if (bi[E] == 0)
2560 {
2561 for (i = 1; i < NI - 1; i++)
2562 {
2563 if (bi[i] != 0)
2564 {
2565 ltb -= enormlz (bi);
2566 goto mnzer2;
2567 }
2568 }
2569 eclear (c);
2570 return;
2571 }
2572 mnzer2:
2573
2574 /* Multiply significands */
2575 j = emulm (ai, bi);
2576 /* calculate exponent */
2577 lt = lta + ltb - (EXONE - 1);
2578 emdnorm (bi, j, 0, lt, 64);
2579 /* calculate sign of product */
2580 if (ai[0] == bi[0])
2581 bi[0] = 0;
2582 else
2583 bi[0] = 0xffff;
2584 emovo (bi, c);
2585}
2586
2587
2588
2589
2590/*
2591; Convert IEEE double precision to e type
2592; double d;
2593; unsigned EMUSHORT x[N+2];
2594; e53toe (&d, x);
2595*/
2596void
2597e53toe (pe, y)
2598 unsigned EMUSHORT *pe, *y;
2599{
2600#ifdef DEC
2601
2602 dectoe (pe, y); /* see etodec.c */
2603
2604#else
2605
2606 register unsigned EMUSHORT r;
2607 register unsigned EMUSHORT *e, *p;
2608 unsigned EMUSHORT yy[NI];
2609 int denorm, k;
2610
2611 e = pe;
2612 denorm = 0; /* flag if denormalized number */
2613 ecleaz (yy);
2614#ifdef IBMPC
2615 e += 3;
2616#endif
2617 r = *e;
2618 yy[0] = 0;
2619 if (r & 0x8000)
2620 yy[0] = 0xffff;
2621 yy[M] = (r & 0x0f) | 0x10;
2622 r &= ~0x800f; /* strip sign and 4 significand bits */
2623#ifdef INFINITY
2624 if (r == 0x7ff0)
2625 {
2626#ifdef NANS
2627#ifdef IBMPC
2628 if (((pe[3] & 0xf) != 0) || (pe[2] != 0)
2629 || (pe[1] != 0) || (pe[0] != 0))
2630 {
2631 enan (y);
2632 return;
2633 }
2634#else
2635 if (((pe[0] & 0xf) != 0) || (pe[1] != 0)
2636 || (pe[2] != 0) || (pe[3] != 0))
2637 {
2638 enan (y);
2639 return;
2640 }
2641#endif
2642#endif /* NANS */
2643 eclear (y);
2644 einfin (y);
2645 if (yy[0])
2646 eneg (y);
2647 return;
2648 }
2649#endif /* INFINITY */
2650 r >>= 4;
2651 /* If zero exponent, then the significand is denormalized.
2652 * So, take back the understood high significand bit. */
2653 if (r == 0)
2654 {
2655 denorm = 1;
2656 yy[M] &= ~0x10;
2657 }
2658 r += EXONE - 01777;
2659 yy[E] = r;
2660 p = &yy[M + 1];
2661#ifdef IBMPC
2662 *p++ = *(--e);
2663 *p++ = *(--e);
2664 *p++ = *(--e);
2665#endif
2666#ifdef MIEEE
2667 ++e;
2668 *p++ = *e++;
2669 *p++ = *e++;
2670 *p++ = *e++;
2671#endif
2672 eshift (yy, -5);
2673 if (denorm)
2674 { /* if zero exponent, then normalize the significand */
2675 if ((k = enormlz (yy)) > NBITS)
2676 ecleazs (yy);
2677 else
2678 yy[E] -= (unsigned EMUSHORT) (k - 1);
2679 }
2680 emovo (yy, y);
2681#endif /* not DEC */
2682}
2683
2684void
2685e64toe (pe, y)
2686 unsigned EMUSHORT *pe, *y;
2687{
2688 unsigned EMUSHORT yy[NI];
2689 unsigned EMUSHORT *e, *p, *q;
2690 int i;
2691
2692 e = pe;
2693 p = yy;
2694 for (i = 0; i < NE - 5; i++)
2695 *p++ = 0;
2696#ifdef IBMPC
2697 for (i = 0; i < 5; i++)
2698 *p++ = *e++;
2699#endif
2700#ifdef DEC
2701 for (i = 0; i < 5; i++)
2702 *p++ = *e++;
2703#endif
2704#ifdef MIEEE
2705 p = &yy[0] + (NE - 1);
2706 *p-- = *e++;
2707 ++e;
2708 for (i = 0; i < 4; i++)
2709 *p-- = *e++;
2710#endif
2711 p = yy;
2712 q = y;
2713#ifdef INFINITY
2714 if (*p == 0x7fff)
2715 {
2716#ifdef NANS
2717#ifdef IBMPC
2718 for (i = 0; i < 4; i++)
2719 {
2720 if (pe[i] != 0)
2721 {
2722 enan (y);
2723 return;
2724 }
2725 }
2726#else
2727 for (i = 1; i <= 4; i++)
2728 {
2729 if (pe[i] != 0)
2730 {
2731 enan (y);
2732 return;
2733 }
2734 }
2735#endif
2736#endif /* NANS */
2737 eclear (y);
2738 einfin (y);
2739 if (*p & 0x8000)
2740 eneg (y);
2741 return;
2742 }
2743#endif /* INFINITY */
2744 for (i = 0; i < NE; i++)
2745 *q++ = *p++;
2746}
2747
2748
2749/*
2750; Convert IEEE single precision to e type
2751; float d;
2752; unsigned EMUSHORT x[N+2];
2753; dtox (&d, x);
2754*/
2755void
2756e24toe (pe, y)
2757 unsigned EMUSHORT *pe, *y;
2758{
2759 register unsigned EMUSHORT r;
2760 register unsigned EMUSHORT *e, *p;
2761 unsigned EMUSHORT yy[NI];
2762 int denorm, k;
2763
2764 e = pe;
2765 denorm = 0; /* flag if denormalized number */
2766 ecleaz (yy);
2767#ifdef IBMPC
2768 e += 1;
2769#endif
2770#ifdef DEC
2771 e += 1;
2772#endif
2773 r = *e;
2774 yy[0] = 0;
2775 if (r & 0x8000)
2776 yy[0] = 0xffff;
2777 yy[M] = (r & 0x7f) | 0200;
2778 r &= ~0x807f; /* strip sign and 7 significand bits */
2779#ifdef INFINITY
2780 if (r == 0x7f80)
2781 {
2782#ifdef NANS
2783#ifdef MIEEE
2784 if (((pe[0] & 0x7f) != 0) || (pe[1] != 0))
2785 {
2786 enan (y);
2787 return;
2788 }
2789#else
2790 if (((pe[1] & 0x7f) != 0) || (pe[0] != 0))
2791 {
2792 enan (y);
2793 return;
2794 }
2795#endif
2796#endif /* NANS */
2797 eclear (y);
2798 einfin (y);
2799 if (yy[0])
2800 eneg (y);
2801 return;
2802 }
2803#endif /* INFINITY */
2804 r >>= 7;
2805 /* If zero exponent, then the significand is denormalized.
2806 * So, take back the understood high significand bit. */
2807 if (r == 0)
2808 {
2809 denorm = 1;
2810 yy[M] &= ~0200;
2811 }
2812 r += EXONE - 0177;
2813 yy[E] = r;
2814 p = &yy[M + 1];
2815#ifdef IBMPC
2816 *p++ = *(--e);
2817#endif
2818#ifdef DEC
2819 *p++ = *(--e);
2820#endif
2821#ifdef MIEEE
2822 ++e;
2823 *p++ = *e++;
2824#endif
2825 eshift (yy, -8);
2826 if (denorm)
2827 { /* if zero exponent, then normalize the significand */
2828 if ((k = enormlz (yy)) > NBITS)
2829 ecleazs (yy);
2830 else
2831 yy[E] -= (unsigned EMUSHORT) (k - 1);
2832 }
2833 emovo (yy, y);
2834}
2835
2836
2837void
2838etoe64 (x, e)
2839 unsigned EMUSHORT *x, *e;
2840{
2841 unsigned EMUSHORT xi[NI];
2842 EMULONG exp;
2843 int rndsav;
2844
2845#ifdef NANS
2846 if (eisnan (x))
2847 {
2848 make_nan (e, XFmode);
2849 return;
2850 }
2851#endif
2852 emovi (x, xi);
2853 /* adjust exponent for offset */
2854 exp = (EMULONG) xi[E];
2855#ifdef INFINITY
2856 if (eisinf (x))
2857 goto nonorm;
2858#endif
2859 /* round off to nearest or even */
2860 rndsav = rndprc;
2861 rndprc = 64;
2862 emdnorm (xi, 0, 0, exp, 64);
2863 rndprc = rndsav;
2864 nonorm:
2865 toe64 (xi, e);
2866}
2867
2868/* move out internal format to ieee long double */
2869static void
2870toe64 (a, b)
2871 unsigned EMUSHORT *a, *b;
2872{
2873 register unsigned EMUSHORT *p, *q;
2874 unsigned EMUSHORT i;
2875
2876#ifdef NANS
2877 if (eiisnan (a))
2878 {
2879 make_nan (b, XFmode);
2880 return;
2881 }
2882#endif
2883 p = a;
2884#ifdef MIEEE
2885 q = b;
2886#else
2887 q = b + 4; /* point to output exponent */
2888#if LONG_DOUBLE_TYPE_SIZE == 96
2889 /* Clear the last two bytes of 12-byte Intel format */
2890 *(q+1) = 0;
2891#endif
2892#endif
2893
2894 /* combine sign and exponent */
2895 i = *p++;
2896#ifdef MIEEE
2897 if (i)
2898 *q++ = *p++ | 0x8000;
2899 else
2900 *q++ = *p++;
2901 *q++ = 0;
2902#else
2903 if (i)
2904 *q-- = *p++ | 0x8000;
2905 else
2906 *q-- = *p++;
2907#endif
2908 /* skip over guard word */
2909 ++p;
2910 /* move the significand */
2911#ifdef MIEEE
2912 for (i = 0; i < 4; i++)
2913 *q++ = *p++;
2914#else
2915 for (i = 0; i < 4; i++)
2916 *q-- = *p++;
2917#endif
2918}
2919
2920
2921/*
2922; e type to IEEE double precision
2923; double d;
2924; unsigned EMUSHORT x[NE];
2925; etoe53 (x, &d);
2926*/
2927
2928#ifdef DEC
2929
2930void
2931etoe53 (x, e)
2932 unsigned EMUSHORT *x, *e;
2933{
2934 etodec (x, e); /* see etodec.c */
2935}
2936
2937static void
2938toe53 (x, y)
2939 unsigned EMUSHORT *x, *y;
2940{
2941 todec (x, y);
2942}
2943
2944#else
2945
2946void
2947etoe53 (x, e)
2948 unsigned EMUSHORT *x, *e;
2949{
2950 unsigned EMUSHORT xi[NI];
2951 EMULONG exp;
2952 int rndsav;
2953
2954#ifdef NANS
2955 if (eisnan (x))
2956 {
2957 make_nan (e, DFmode);
2958 return;
2959 }
2960#endif
2961 emovi (x, xi);
2962 /* adjust exponent for offsets */
2963 exp = (EMULONG) xi[E] - (EXONE - 0x3ff);
2964#ifdef INFINITY
2965 if (eisinf (x))
2966 goto nonorm;
2967#endif
2968 /* round off to nearest or even */
2969 rndsav = rndprc;
2970 rndprc = 53;
2971 emdnorm (xi, 0, 0, exp, 64);
2972 rndprc = rndsav;
2973 nonorm:
2974 toe53 (xi, e);
2975}
2976
2977
2978static void
2979toe53 (x, y)
2980 unsigned EMUSHORT *x, *y;
2981{
2982 unsigned EMUSHORT i;
2983 unsigned EMUSHORT *p;
2984
2985#ifdef NANS
2986 if (eiisnan (x))
2987 {
2988 make_nan (y, DFmode);
2989 return;
2990 }
2991#endif
2992 p = &x[0];
2993#ifdef IBMPC
2994 y += 3;
2995#endif
2996 *y = 0; /* output high order */
2997 if (*p++)
2998 *y = 0x8000; /* output sign bit */
2999
3000 i = *p++;
3001 if (i >= (unsigned int) 2047)
3002 { /* Saturate at largest number less than infinity. */
3003#ifdef INFINITY
3004 *y |= 0x7ff0;
3005#ifdef IBMPC
3006 *(--y) = 0;
3007 *(--y) = 0;
3008 *(--y) = 0;
3009#endif
3010#ifdef MIEEE
3011 ++y;
3012 *y++ = 0;
3013 *y++ = 0;
3014 *y++ = 0;
3015#endif
3016#else
3017 *y |= (unsigned EMUSHORT) 0x7fef;
3018#ifdef IBMPC
3019 *(--y) = 0xffff;
3020 *(--y) = 0xffff;
3021 *(--y) = 0xffff;
3022#endif
3023#ifdef MIEEE
3024 ++y;
3025 *y++ = 0xffff;
3026 *y++ = 0xffff;
3027 *y++ = 0xffff;
3028#endif
3029#endif
3030 return;
3031 }
3032 if (i == 0)
3033 {
3034 eshift (x, 4);
3035 }
3036 else
3037 {
3038 i <<= 4;
3039 eshift (x, 5);
3040 }
3041 i |= *p++ & (unsigned EMUSHORT) 0x0f; /* *p = xi[M] */
3042 *y |= (unsigned EMUSHORT) i; /* high order output already has sign bit set */
3043#ifdef IBMPC
3044 *(--y) = *p++;
3045 *(--y) = *p++;
3046 *(--y) = *p;
3047#endif
3048#ifdef MIEEE
3049 ++y;
3050 *y++ = *p++;
3051 *y++ = *p++;
3052 *y++ = *p++;
3053#endif
3054}
3055
3056#endif /* not DEC */
3057
3058
3059
3060/*
3061; e type to IEEE single precision
3062; float d;
3063; unsigned EMUSHORT x[N+2];
3064; xtod (x, &d);
3065*/
3066void
3067etoe24 (x, e)
3068 unsigned EMUSHORT *x, *e;
3069{
3070 EMULONG exp;
3071 unsigned EMUSHORT xi[NI];
3072 int rndsav;
3073
3074#ifdef NANS
3075 if (eisnan (x))
3076 {
3077 make_nan (e, SFmode);
3078 return;
3079 }
3080#endif
3081 emovi (x, xi);
3082 /* adjust exponent for offsets */
3083 exp = (EMULONG) xi[E] - (EXONE - 0177);
3084#ifdef INFINITY
3085 if (eisinf (x))
3086 goto nonorm;
3087#endif
3088 /* round off to nearest or even */
3089 rndsav = rndprc;
3090 rndprc = 24;
3091 emdnorm (xi, 0, 0, exp, 64);
3092 rndprc = rndsav;
3093 nonorm:
3094 toe24 (xi, e);
3095}
3096
3097static void
3098toe24 (x, y)
3099 unsigned EMUSHORT *x, *y;
3100{
3101 unsigned EMUSHORT i;
3102 unsigned EMUSHORT *p;
3103
3104#ifdef NANS
3105 if (eiisnan (x))
3106 {
3107 make_nan (y, SFmode);
3108 return;
3109 }
3110#endif
3111 p = &x[0];
3112#ifdef IBMPC
3113 y += 1;
3114#endif
3115#ifdef DEC
3116 y += 1;
3117#endif
3118 *y = 0; /* output high order */
3119 if (*p++)
3120 *y = 0x8000; /* output sign bit */
3121
3122 i = *p++;
3123/* Handle overflow cases. */
3124 if (i >= 255)
3125 {
3126#ifdef INFINITY
3127 *y |= (unsigned EMUSHORT) 0x7f80;
3128#ifdef IBMPC
3129 *(--y) = 0;
3130#endif
3131#ifdef DEC
3132 *(--y) = 0;
3133#endif
3134#ifdef MIEEE
3135 ++y;
3136 *y = 0;
3137#endif
3138#else /* no INFINITY */
3139 *y |= (unsigned EMUSHORT) 0x7f7f;
3140#ifdef IBMPC
3141 *(--y) = 0xffff;
3142#endif
3143#ifdef DEC
3144 *(--y) = 0xffff;
3145#endif
3146#ifdef MIEEE
3147 ++y;
3148 *y = 0xffff;
3149#endif
3150#ifdef ERANGE
3151 errno = ERANGE;
3152#endif
3153#endif /* no INFINITY */
3154 return;
3155 }
3156 if (i == 0)
3157 {
3158 eshift (x, 7);
3159 }
3160 else
3161 {
3162 i <<= 7;
3163 eshift (x, 8);
3164 }
3165 i |= *p++ & (unsigned EMUSHORT) 0x7f; /* *p = xi[M] */
3166 *y |= i; /* high order output already has sign bit set */
3167#ifdef IBMPC
3168 *(--y) = *p;
3169#endif
3170#ifdef DEC
3171 *(--y) = *p;
3172#endif
3173#ifdef MIEEE
3174 ++y;
3175 *y = *p;
3176#endif
3177}
3178
3179
3180/* Compare two e type numbers.
3181 *
3182 * unsigned EMUSHORT a[NE], b[NE];
3183 * ecmp (a, b);
3184 *
3185 * returns +1 if a > b
3186 * 0 if a == b
3187 * -1 if a < b
3188 * -2 if either a or b is a NaN.
3189 */
3190int
3191ecmp (a, b)
3192 unsigned EMUSHORT *a, *b;
3193{
3194 unsigned EMUSHORT ai[NI], bi[NI];
3195 register unsigned EMUSHORT *p, *q;
3196 register int i;
3197 int msign;
3198
3199#ifdef NANS
3200 if (eisnan (a) || eisnan (b))
3201 return (-2);
3202#endif
3203 emovi (a, ai);
3204 p = ai;
3205 emovi (b, bi);
3206 q = bi;
3207
3208 if (*p != *q)
3209 { /* the signs are different */
3210 /* -0 equals + 0 */
3211 for (i = 1; i < NI - 1; i++)
3212 {
3213 if (ai[i] != 0)
3214 goto nzro;
3215 if (bi[i] != 0)
3216 goto nzro;
3217 }
3218 return (0);
3219 nzro:
3220 if (*p == 0)
3221 return (1);
3222 else
3223 return (-1);
3224 }
3225 /* both are the same sign */
3226 if (*p == 0)
3227 msign = 1;
3228 else
3229 msign = -1;
3230 i = NI - 1;
3231 do
3232 {
3233 if (*p++ != *q++)
3234 {
3235 goto diff;
3236 }
3237 }
3238 while (--i > 0);
3239
3240 return (0); /* equality */
3241
3242
3243
3244 diff:
3245
3246 if (*(--p) > *(--q))
3247 return (msign); /* p is bigger */
3248 else
3249 return (-msign); /* p is littler */
3250}
3251
3252
3253
3254
3255/* Find nearest integer to x = floor (x + 0.5)
3256 *
3257 * unsigned EMUSHORT x[NE], y[NE]
3258 * eround (x, y);
3259 */
3260void
3261eround (x, y)
3262 unsigned EMUSHORT *x, *y;
3263{
3264 eadd (ehalf, x, y);
3265 efloor (y, y);
3266}
3267
3268
3269
3270
3271/*
3272; convert long integer to e type
3273;
3274; long l;
3275; unsigned EMUSHORT x[NE];
3276; ltoe (&l, x);
3277; note &l is the memory address of l
3278*/
3279void
3280ltoe (lp, y)
3281 long *lp; /* lp is the memory address of a long integer */
3282 unsigned EMUSHORT *y; /* y is the address of a short */
3283{
3284 unsigned EMUSHORT yi[NI];
3285 unsigned long ll;
3286 int k;
3287
3288 ecleaz (yi);
3289 if (*lp < 0)
3290 {
3291 /* make it positive */
3292 ll = (unsigned long) (-(*lp));
3293 yi[0] = 0xffff; /* put correct sign in the e type number */
3294 }
3295 else
3296 {
3297 ll = (unsigned long) (*lp);
3298 }
3299 /* move the long integer to yi significand area */
3300#if HOST_BITS_PER_LONG == 64
3301 yi[M] = (unsigned EMUSHORT) (ll >> 48);
3302 yi[M + 1] = (unsigned EMUSHORT) (ll >> 32);
3303 yi[M + 2] = (unsigned EMUSHORT) (ll >> 16);
3304 yi[M + 3] = (unsigned EMUSHORT) ll;
3305 yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
3306#else
3307 yi[M] = (unsigned EMUSHORT) (ll >> 16);
3308 yi[M + 1] = (unsigned EMUSHORT) ll;
3309 yi[E] = EXONE + 15; /* exponent if normalize shift count were 0 */
3310#endif
3311
3312 if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
3313 ecleaz (yi); /* it was zero */
3314 else
3315 yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
3316 emovo (yi, y); /* output the answer */
3317}
3318
3319/*
3320; convert unsigned long integer to e type
3321;
3322; unsigned long l;
3323; unsigned EMUSHORT x[NE];
3324; ltox (&l, x);
3325; note &l is the memory address of l
3326*/
3327void
3328ultoe (lp, y)
3329 unsigned long *lp; /* lp is the memory address of a long integer */
3330 unsigned EMUSHORT *y; /* y is the address of a short */
3331{
3332 unsigned EMUSHORT yi[NI];
3333 unsigned long ll;
3334 int k;
3335
3336 ecleaz (yi);
3337 ll = *lp;
3338
3339 /* move the long integer to ayi significand area */
3340#if HOST_BITS_PER_LONG == 64
3341 yi[M] = (unsigned EMUSHORT) (ll >> 48);
3342 yi[M + 1] = (unsigned EMUSHORT) (ll >> 32);
3343 yi[M + 2] = (unsigned EMUSHORT) (ll >> 16);
3344 yi[M + 3] = (unsigned EMUSHORT) ll;
3345 yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
3346#else
3347 yi[M] = (unsigned EMUSHORT) (ll >> 16);
3348 yi[M + 1] = (unsigned EMUSHORT) ll;
3349 yi[E] = EXONE + 15; /* exponent if normalize shift count were 0 */
3350#endif
3351
3352 if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
3353 ecleaz (yi); /* it was zero */
3354 else
3355 yi[E] -= (unsigned EMUSHORT) k; /* subtract shift count from exponent */
3356 emovo (yi, y); /* output the answer */
3357}
3358
3359
3360/*
3361; Find long integer and fractional parts
3362
3363; long i;
3364; unsigned EMUSHORT x[NE], frac[NE];
3365; xifrac (x, &i, frac);
3366
3367 The integer output has the sign of the input. The fraction is
3368the positive fractional part of abs (x).
3369*/
3370void
3371eifrac (x, i, frac)
3372 unsigned EMUSHORT *x;
3373 long *i;
3374 unsigned EMUSHORT *frac;
3375{
3376 unsigned EMUSHORT xi[NI];
3377 int j, k;
3378 unsigned long ll;
3379
3380 emovi (x, xi);
3381 k = (int) xi[E] - (EXONE - 1);
3382 if (k <= 0)
3383 {
3384 /* if exponent <= 0, integer = 0 and real output is fraction */
3385 *i = 0L;
3386 emovo (xi, frac);
3387 return;
3388 }
3389 if (k > (HOST_BITS_PER_LONG - 1))
3390 {
3391 /* long integer overflow: output large integer
3392 and correct fraction */
3393 if (xi[0])
3394 *i = ((unsigned long) 1) << (HOST_BITS_PER_LONG - 1);
3395 else
3396 *i = (((unsigned long) 1) << (HOST_BITS_PER_LONG - 1)) - 1;
3397 eshift (xi, k);
3398 if (extra_warnings)
3399 warning ("overflow on truncation to integer");
3400 }
3401 else if (k > 16)
3402 {
3403 /* Shift more than 16 bits: first shift up k-16 mod 16,
3404 then shift up by 16's. */
3405 j = k - ((k >> 4) << 4);
3406 eshift (xi, j);
3407 ll = xi[M];
3408 k -= j;
3409 do
3410 {
3411 eshup6 (xi);
3412 ll = (ll << 16) | xi[M];
3413 }
3414 while ((k -= 16) > 0);
3415 *i = ll;
3416 if (xi[0])
3417 *i = -(*i);
3418 }
3419 else
3420 {
3421 /* shift not more than 16 bits */
3422 eshift (xi, k);
3423 *i = (long) xi[M] & 0xffff;
3424 if (xi[0])
3425 *i = -(*i);
3426 }
3427 xi[0] = 0;
3428 xi[E] = EXONE - 1;
3429 xi[M] = 0;
3430 if ((k = enormlz (xi)) > NBITS)
3431 ecleaz (xi);
3432 else
3433 xi[E] -= (unsigned EMUSHORT) k;
3434
3435 emovo (xi, frac);
3436}
3437
3438
3439/* Find unsigned long integer and fractional parts.
3440 A negative e type input yields integer output = 0
3441 but correct fraction. */
3442
3443void
3444euifrac (x, i, frac)
3445 unsigned EMUSHORT *x;
3446 unsigned long *i;
3447 unsigned EMUSHORT *frac;
3448{
3449 unsigned long ll;
3450 unsigned EMUSHORT xi[NI];
3451 int j, k;
3452
3453 emovi (x, xi);
3454 k = (int) xi[E] - (EXONE - 1);
3455 if (k <= 0)
3456 {
3457 /* if exponent <= 0, integer = 0 and argument is fraction */
3458 *i = 0L;
3459 emovo (xi, frac);
3460 return;
3461 }
3462 if (k > HOST_BITS_PER_LONG)
3463 {
3464 /* Long integer overflow: output large integer
3465 and correct fraction.
3466 Note, the BSD microvax compiler says that ~(0UL)
3467 is a syntax error. */
3468 *i = ~(0L);
3469 eshift (xi, k);
3470 if (extra_warnings)
3471 warning ("overflow on truncation to unsigned integer");
3472 }
3473 else if (k > 16)
3474 {
3475 /* Shift more than 16 bits: first shift up k-16 mod 16,
3476 then shift up by 16's. */
3477 j = k - ((k >> 4) << 4);
3478 eshift (xi, j);
3479 ll = xi[M];
3480 k -= j;
3481 do
3482 {
3483 eshup6 (xi);
3484 ll = (ll << 16) | xi[M];
3485 }
3486 while ((k -= 16) > 0);
3487 *i = ll;
3488 }
3489 else
3490 {
3491 /* shift not more than 16 bits */
3492 eshift (xi, k);
3493 *i = (long) xi[M] & 0xffff;
3494 }
3495
3496 if (xi[0]) /* A negative value yields unsigned integer 0. */
3497 *i = 0L;
3498 xi[0] = 0;
3499 xi[E] = EXONE - 1;
3500 xi[M] = 0;
3501 if ((k = enormlz (xi)) > NBITS)
3502 ecleaz (xi);
3503 else
3504 xi[E] -= (unsigned EMUSHORT) k;
3505
3506 emovo (xi, frac);
3507}
3508
3509
3510
3511/*
3512; Shift significand
3513;
3514; Shifts significand area up or down by the number of bits
3515; given by the variable sc.
3516*/
3517int
3518eshift (x, sc)
3519 unsigned EMUSHORT *x;
3520 int sc;
3521{
3522 unsigned EMUSHORT lost;
3523 unsigned EMUSHORT *p;
3524
3525 if (sc == 0)
3526 return (0);
3527
3528 lost = 0;
3529 p = x + NI - 1;
3530
3531 if (sc < 0)
3532 {
3533 sc = -sc;
3534 while (sc >= 16)
3535 {
3536 lost |= *p; /* remember lost bits */
3537 eshdn6 (x);
3538 sc -= 16;
3539 }
3540
3541 while (sc >= 8)
3542 {
3543 lost |= *p & 0xff;
3544 eshdn8 (x);
3545 sc -= 8;
3546 }
3547
3548 while (sc > 0)
3549 {
3550 lost |= *p & 1;
3551 eshdn1 (x);
3552 sc -= 1;
3553 }
3554 }
3555 else
3556 {
3557 while (sc >= 16)
3558 {
3559 eshup6 (x);
3560 sc -= 16;
3561 }
3562
3563 while (sc >= 8)
3564 {
3565 eshup8 (x);
3566 sc -= 8;
3567 }
3568
3569 while (sc > 0)
3570 {
3571 eshup1 (x);
3572 sc -= 1;
3573 }
3574 }
3575 if (lost)
3576 lost = 1;
3577 return ((int) lost);
3578}
3579
3580
3581
3582/*
3583; normalize
3584;
3585; Shift normalizes the significand area pointed to by argument
3586; shift count (up = positive) is returned.
3587*/
3588int
3589enormlz (x)
3590 unsigned EMUSHORT x[];
3591{
3592 register unsigned EMUSHORT *p;
3593 int sc;
3594
3595 sc = 0;
3596 p = &x[M];
3597 if (*p != 0)
3598 goto normdn;
3599 ++p;
3600 if (*p & 0x8000)
3601 return (0); /* already normalized */
3602 while (*p == 0)
3603 {
3604 eshup6 (x);
3605 sc += 16;
3606 /* With guard word, there are NBITS+16 bits available.
3607 * return true if all are zero.
3608 */
3609 if (sc > NBITS)
3610 return (sc);
3611 }
3612 /* see if high byte is zero */
3613 while ((*p & 0xff00) == 0)
3614 {
3615 eshup8 (x);
3616 sc += 8;
3617 }
3618 /* now shift 1 bit at a time */
3619 while ((*p & 0x8000) == 0)
3620 {
3621 eshup1 (x);
3622 sc += 1;
3623 if (sc > NBITS)
3624 {
3625 mtherr ("enormlz", UNDERFLOW);
3626 return (sc);
3627 }
3628 }
3629 return (sc);
3630
3631 /* Normalize by shifting down out of the high guard word
3632 of the significand */
3633 normdn:
3634
3635 if (*p & 0xff00)
3636 {
3637 eshdn8 (x);
3638 sc -= 8;
3639 }
3640 while (*p != 0)
3641 {
3642 eshdn1 (x);
3643 sc -= 1;
3644
3645 if (sc < -NBITS)
3646 {
3647 mtherr ("enormlz", OVERFLOW);
3648 return (sc);
3649 }
3650 }
3651 return (sc);
3652}
3653
3654
3655
3656
3657/* Convert e type number to decimal format ASCII string.
3658 * The constants are for 64 bit precision.
3659 */
3660
3661#define NTEN 12
3662#define MAXP 4096
3663
3664static unsigned EMUSHORT etens[NTEN + 1][NE] =
3665{
3666 {0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,}, /* 10**4096 */
3667 {0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,}, /* 10**2048 */
3668 {0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
3669 {0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
3670 {0xddbc, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
3671 {0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
3672 {0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
3673 {0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
3674 {0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
3675 {0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
3676 {0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
3677 {0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
3678 {0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,}, /* 10**1 */
3679};
3680
3681static unsigned EMUSHORT emtens[NTEN + 1][NE] =
3682{
3683 {0x2de4, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,}, /* 10**-4096 */
3684 {0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,}, /* 10**-2048 */
3685 {0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
3686 {0x7133, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
3687 {0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
3688 {0xac7d, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
3689 {0x3f24, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
3690 {0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
3691 {0x4c2f, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
3692 {0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
3693 {0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
3694 {0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
3695 {0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,}, /* 10**-1 */
3696};
3697
3698void
3699e24toasc (x, string, ndigs)
3700 unsigned EMUSHORT x[];
3701 char *string;
3702 int ndigs;
3703{
3704 unsigned EMUSHORT w[NI];
3705
3706 e24toe (x, w);
3707 etoasc (w, string, ndigs);
3708}
3709
3710
3711void
3712e53toasc (x, string, ndigs)
3713 unsigned EMUSHORT x[];
3714 char *string;
3715 int ndigs;
3716{
3717 unsigned EMUSHORT w[NI];
3718
3719 e53toe (x, w);
3720 etoasc (w, string, ndigs);
3721}
3722
3723
3724void
3725e64toasc (x, string, ndigs)
3726 unsigned EMUSHORT x[];
3727 char *string;
3728 int ndigs;
3729{
3730 unsigned EMUSHORT w[NI];
3731
3732 e64toe (x, w);
3733 etoasc (w, string, ndigs);
3734}
3735
3736
3737static char wstring[80]; /* working storage for ASCII output */
3738
3739void
3740etoasc (x, string, ndigs)
3741 unsigned EMUSHORT x[];
3742 char *string;
3743 int ndigs;
3744{
3745 EMUSHORT digit;
3746 unsigned EMUSHORT y[NI], t[NI], u[NI], w[NI];
3747 unsigned EMUSHORT *p, *r, *ten;
3748 unsigned EMUSHORT sign;
3749 int i, j, k, expon, rndsav;
3750 char *s, *ss;
3751 unsigned EMUSHORT m;
3752
3753
3754 rndsav = rndprc;
3755 ss = string;
3756 s = wstring;
3757 *ss = '\0';
3758 *s = '\0';
3759#ifdef NANS
3760 if (eisnan (x))
3761 {
3762 sprintf (wstring, " NaN ");
3763 goto bxit;
3764 }
3765#endif
3766 rndprc = NBITS; /* set to full precision */
3767 emov (x, y); /* retain external format */
3768 if (y[NE - 1] & 0x8000)
3769 {
3770 sign = 0xffff;
3771 y[NE - 1] &= 0x7fff;
3772 }
3773 else
3774 {
3775 sign = 0;
3776 }
3777 expon = 0;
3778 ten = &etens[NTEN][0];
3779 emov (eone, t);
3780 /* Test for zero exponent */
3781 if (y[NE - 1] == 0)
3782 {
3783 for (k = 0; k < NE - 1; k++)
3784 {
3785 if (y[k] != 0)
3786 goto tnzro; /* denormalized number */
3787 }
3788 goto isone; /* legal all zeros */
3789 }
3790 tnzro:
3791
3792 /* Test for infinity. */
3793 if (y[NE - 1] == 0x7fff)
3794 {
3795 if (sign)
3796 sprintf (wstring, " -Infinity ");
3797 else
3798 sprintf (wstring, " Infinity ");
3799 goto bxit;
3800 }
3801
3802 /* Test for exponent nonzero but significand denormalized.
3803 * This is an error condition.
3804 */
3805 if ((y[NE - 1] != 0) && ((y[NE - 2] & 0x8000) == 0))
3806 {
3807 mtherr ("etoasc", DOMAIN);
3808 sprintf (wstring, "NaN");
3809 goto bxit;
3810 }
3811
3812 /* Compare to 1.0 */
3813 i = ecmp (eone, y);
3814 if (i == 0)
3815 goto isone;
3816
3817 if (i == -2)
3818 abort ();
3819
3820 if (i < 0)
3821 { /* Number is greater than 1 */
3822 /* Convert significand to an integer and strip trailing decimal zeros. */
3823 emov (y, u);
3824 u[NE - 1] = EXONE + NBITS - 1;
3825
3826 p = &etens[NTEN - 4][0];
3827 m = 16;
3828 do
3829 {
3830 ediv (p, u, t);
3831 efloor (t, w);
3832 for (j = 0; j < NE - 1; j++)
3833 {
3834 if (t[j] != w[j])
3835 goto noint;
3836 }
3837 emov (t, u);
3838 expon += (int) m;
3839 noint:
3840 p += NE;
3841 m >>= 1;
3842 }
3843 while (m != 0);
3844
3845 /* Rescale from integer significand */
3846 u[NE - 1] += y[NE - 1] - (unsigned int) (EXONE + NBITS - 1);
3847 emov (u, y);
3848 /* Find power of 10 */
3849 emov (eone, t);
3850 m = MAXP;
3851 p = &etens[0][0];
3852 /* An unordered compare result shouldn't happen here. */
3853 while (ecmp (ten, u) <= 0)
3854 {
3855 if (ecmp (p, u) <= 0)
3856 {
3857 ediv (p, u, u);
3858 emul (p, t, t);
3859 expon += (int) m;
3860 }
3861 m >>= 1;
3862 if (m == 0)
3863 break;
3864 p += NE;
3865 }
3866 }
3867 else
3868 { /* Number is less than 1.0 */
3869 /* Pad significand with trailing decimal zeros. */
3870 if (y[NE - 1] == 0)
3871 {
3872 while ((y[NE - 2] & 0x8000) == 0)
3873 {
3874 emul (ten, y, y);
3875 expon -= 1;
3876 }
3877 }
3878 else
3879 {
3880 emovi (y, w);
3881 for (i = 0; i < NDEC + 1; i++)
3882 {
3883 if ((w[NI - 1] & 0x7) != 0)
3884 break;
3885 /* multiply by 10 */
3886 emovz (w, u);
3887 eshdn1 (u);
3888 eshdn1 (u);
3889 eaddm (w, u);
3890 u[1] += 3;
3891 while (u[2] != 0)
3892 {
3893 eshdn1 (u);
3894 u[1] += 1;
3895 }
3896 if (u[NI - 1] != 0)
3897 break;
3898 if (eone[NE - 1] <= u[1])
3899 break;
3900 emovz (u, w);
3901 expon -= 1;
3902 }
3903 emovo (w, y);
3904 }
3905 k = -MAXP;
3906 p = &emtens[0][0];
3907 r = &etens[0][0];
3908 emov (y, w);
3909 emov (eone, t);
3910 while (ecmp (eone, w) > 0)
3911 {
3912 if (ecmp (p, w) >= 0)
3913 {
3914 emul (r, w, w);
3915 emul (r, t, t);
3916 expon += k;
3917 }
3918 k /= 2;
3919 if (k == 0)
3920 break;
3921 p += NE;
3922 r += NE;
3923 }
3924 ediv (t, eone, t);
3925 }
3926 isone:
3927 /* Find the first (leading) digit. */
3928 emovi (t, w);
3929 emovz (w, t);
3930 emovi (y, w);
3931 emovz (w, y);
3932 eiremain (t, y);
3933 digit = equot[NI - 1];
3934 while ((digit == 0) && (ecmp (y, ezero) != 0))
3935 {
3936 eshup1 (y);
3937 emovz (y, u);
3938 eshup1 (u);
3939 eshup1 (u);
3940 eaddm (u, y);
3941 eiremain (t, y);
3942 digit = equot[NI - 1];
3943 expon -= 1;
3944 }
3945 s = wstring;
3946 if (sign)
3947 *s++ = '-';
3948 else
3949 *s++ = ' ';
3950 /* Examine number of digits requested by caller. */
3951 if (ndigs < 0)
3952 ndigs = 0;
3953 if (ndigs > NDEC)
3954 ndigs = NDEC;
3955 if (digit == 10)
3956 {
3957 *s++ = '1';
3958 *s++ = '.';
3959 if (ndigs > 0)
3960 {
3961 *s++ = '0';
3962 ndigs -= 1;
3963 }
3964 expon += 1;
3965 }
3966 else
3967 {
3968 *s++ = (char )digit + '0';
3969 *s++ = '.';
3970 }
3971 /* Generate digits after the decimal point. */
3972 for (k = 0; k <= ndigs; k++)
3973 {
3974 /* multiply current number by 10, without normalizing */
3975 eshup1 (y);
3976 emovz (y, u);
3977 eshup1 (u);
3978 eshup1 (u);
3979 eaddm (u, y);
3980 eiremain (t, y);
3981 *s++ = (char) equot[NI - 1] + '0';
3982 }
3983 digit = equot[NI - 1];
3984 --s;
3985 ss = s;
3986 /* round off the ASCII string */
3987 if (digit > 4)
3988 {
3989 /* Test for critical rounding case in ASCII output. */
3990 if (digit == 5)
3991 {
3992 emovo (y, t);
3993 if (ecmp (t, ezero) != 0)
3994 goto roun; /* round to nearest */
3995 if ((*(s - 1) & 1) == 0)
3996 goto doexp; /* round to even */
3997 }
3998 /* Round up and propagate carry-outs */
3999 roun:
4000 --s;
4001 k = *s & 0x7f;
4002 /* Carry out to most significant digit? */
4003 if (k == '.')
4004 {
4005 --s;
4006 k = *s;
4007 k += 1;
4008 *s = (char) k;
4009 /* Most significant digit carries to 10? */
4010 if (k > '9')
4011 {
4012 expon += 1;
4013 *s = '1';
4014 }
4015 goto doexp;
4016 }
4017 /* Round up and carry out from less significant digits */
4018 k += 1;
4019 *s = (char) k;
4020 if (k > '9')
4021 {
4022 *s = '0';
4023 goto roun;
4024 }
4025 }
4026 doexp:
4027 /*
4028 if (expon >= 0)
4029 sprintf (ss, "e+%d", expon);
4030 else
4031 sprintf (ss, "e%d", expon);
4032 */
4033 sprintf (ss, "e%d", expon);
4034 bxit:
4035 rndprc = rndsav;
4036 /* copy out the working string */
4037 s = string;
4038 ss = wstring;
4039 while (*ss == ' ') /* strip possible leading space */
4040 ++ss;
4041 while ((*s++ = *ss++) != '\0')
4042 ;
4043}
4044
4045
4046
4047
4048/*
4049; ASCTOQ
4050; ASCTOQ.MAC LATEST REV: 11 JAN 84
4051; SLM, 3 JAN 78
4052;
4053; Convert ASCII string to quadruple precision floating point
4054;
4055; Numeric input is free field decimal number
4056; with max of 15 digits with or without
4057; decimal point entered as ASCII from teletype.
4058; Entering E after the number followed by a second
4059; number causes the second number to be interpreted
4060; as a power of 10 to be multiplied by the first number
4061; (i.e., "scientific" notation).
4062;
4063; Usage:
4064; asctoq (string, q);
4065*/
4066
4067/* ASCII to single */
4068void
4069asctoe24 (s, y)
4070 char *s;
4071 unsigned EMUSHORT *y;
4072{
4073 asctoeg (s, y, 24);
4074}
4075
4076
4077/* ASCII to double */
4078void
4079asctoe53 (s, y)
4080 char *s;
4081 unsigned EMUSHORT *y;
4082{
4083#ifdef DEC
4084 asctoeg (s, y, 56);
4085#else
4086 asctoeg (s, y, 53);
4087#endif
4088}
4089
4090
4091/* ASCII to long double */
4092void
4093asctoe64 (s, y)
4094 char *s;
4095 unsigned EMUSHORT *y;
4096{
4097 asctoeg (s, y, 64);
4098}
4099
4100/* ASCII to super double */
4101void
4102asctoe (s, y)
4103 char *s;
4104 unsigned EMUSHORT *y;
4105{
4106 asctoeg (s, y, NBITS);
4107}
4108
4109/* Space to make a copy of the input string: */
4110static char lstr[82];
4111
4112void
4113asctoeg (ss, y, oprec)
4114 char *ss;
4115 unsigned EMUSHORT *y;
4116 int oprec;
4117{
4118 unsigned EMUSHORT yy[NI], xt[NI], tt[NI];
4119 int esign, decflg, sgnflg, nexp, exp, prec, lost;
4120 int k, trail, c, rndsav;
4121 EMULONG lexp;
4122 unsigned EMUSHORT nsign, *p;
4123 char *sp, *s;
4124
4125 /* Copy the input string. */
4126 s = ss;
4127 while (*s == ' ') /* skip leading spaces */
4128 ++s;
4129 sp = lstr;
4130 for (k = 0; k < 79; k++)
4131 {
4132 if ((*sp++ = *s++) == '\0')
4133 break;
4134 }
4135 *sp = '\0';
4136 s = lstr;
4137
4138 rndsav = rndprc;
4139 rndprc = NBITS; /* Set to full precision */
4140 lost = 0;
4141 nsign = 0;
4142 decflg = 0;
4143 sgnflg = 0;
4144 nexp = 0;
4145 exp = 0;
4146 prec = 0;
4147 ecleaz (yy);
4148 trail = 0;
4149
4150 nxtcom:
4151 k = *s - '0';
4152 if ((k >= 0) && (k <= 9))
4153 {
4154 /* Ignore leading zeros */
4155 if ((prec == 0) && (decflg == 0) && (k == 0))
4156 goto donchr;
4157 /* Identify and strip trailing zeros after the decimal point. */
4158 if ((trail == 0) && (decflg != 0))
4159 {
4160 sp = s;
4161 while ((*sp >= '0') && (*sp <= '9'))
4162 ++sp;
4163 /* Check for syntax error */
4164 c = *sp & 0x7f;
4165 if ((c != 'e') && (c != 'E') && (c != '\0')
4166 && (c != '\n') && (c != '\r') && (c != ' ')
4167 && (c != ','))
4168 goto error;
4169 --sp;
4170 while (*sp == '0')
4171 *sp-- = 'z';
4172 trail = 1;
4173 if (*s == 'z')
4174 goto donchr;
4175 }
4176 /* If enough digits were given to more than fill up the yy register,
4177 * continuing until overflow into the high guard word yy[2]
4178 * guarantees that there will be a roundoff bit at the top
4179 * of the low guard word after normalization.
4180 */
4181 if (yy[2] == 0)
4182 {
4183 if (decflg)
4184 nexp += 1; /* count digits after decimal point */
4185 eshup1 (yy); /* multiply current number by 10 */
4186 emovz (yy, xt);
4187 eshup1 (xt);
4188 eshup1 (xt);
4189 eaddm (xt, yy);
4190 ecleaz (xt);
4191 xt[NI - 2] = (unsigned EMUSHORT) k;
4192 eaddm (xt, yy);
4193 }
4194 else
4195 {
4196 lost |= k;
4197 }
4198 prec += 1;
4199 goto donchr;
4200 }
4201
4202 switch (*s)
4203 {
4204 case 'z':
4205 break;
4206 case 'E':
4207 case 'e':
4208 goto expnt;
4209 case '.': /* decimal point */
4210 if (decflg)
4211 goto error;
4212 ++decflg;
4213 break;
4214 case '-':
4215 nsign = 0xffff;
4216 if (sgnflg)
4217 goto error;
4218 ++sgnflg;
4219 break;
4220 case '+':
4221 if (sgnflg)
4222 goto error;
4223 ++sgnflg;
4224 break;
4225 case ',':
4226 case ' ':
4227 case '\0':
4228 case '\n':
4229 case '\r':
4230 goto daldone;
4231 case 'i':
4232 case 'I':
4233 goto infinite;
4234 default:
4235 error:
4236#ifdef NANS
4237 einan (yy);
4238#else
4239 mtherr ("asctoe", DOMAIN);
4240 eclear (yy);
4241#endif
4242 goto aexit;
4243 }
4244 donchr:
4245 ++s;
4246 goto nxtcom;
4247
4248 /* Exponent interpretation */
4249 expnt:
4250
4251 esign = 1;
4252 exp = 0;
4253 ++s;
4254 /* check for + or - */
4255 if (*s == '-')
4256 {
4257 esign = -1;
4258 ++s;
4259 }
4260 if (*s == '+')
4261 ++s;
4262 while ((*s >= '0') && (*s <= '9'))
4263 {
4264 exp *= 10;
4265 exp += *s++ - '0';
4266 if (exp > 4956)
4267 {
4268 if (esign < 0)
4269 goto zero;
4270 else
4271 goto infinite;
4272 }
4273 }
4274 if (esign < 0)
4275 exp = -exp;
4276 if (exp > 4932)
4277 {
4278 infinite:
4279 ecleaz (yy);
4280 yy[E] = 0x7fff; /* infinity */
4281 goto aexit;
4282 }
4283 if (exp < -4956)
4284 {
4285 zero:
4286 ecleaz (yy);
4287 goto aexit;
4288 }
4289
4290 daldone:
4291 nexp = exp - nexp;
4292 /* Pad trailing zeros to minimize power of 10, per IEEE spec. */
4293 while ((nexp > 0) && (yy[2] == 0))
4294 {
4295 emovz (yy, xt);
4296 eshup1 (xt);
4297 eshup1 (xt);
4298 eaddm (yy, xt);
4299 eshup1 (xt);
4300 if (xt[2] != 0)
4301 break;
4302 nexp -= 1;
4303 emovz (xt, yy);
4304 }
4305 if ((k = enormlz (yy)) > NBITS)
4306 {
4307 ecleaz (yy);
4308 goto aexit;
4309 }
4310 lexp = (EXONE - 1 + NBITS) - k;
4311 emdnorm (yy, lost, 0, lexp, 64);
4312 /* convert to external format */
4313
4314
4315 /* Multiply by 10**nexp. If precision is 64 bits,
4316 * the maximum relative error incurred in forming 10**n
4317 * for 0 <= n <= 324 is 8.2e-20, at 10**180.
4318 * For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
4319 * For 0 >= n >= -999, it is -1.55e-19 at 10**-435.
4320 */
4321 lexp = yy[E];
4322 if (nexp == 0)
4323 {
4324 k = 0;
4325 goto expdon;
4326 }
4327 esign = 1;
4328 if (nexp < 0)
4329 {
4330 nexp = -nexp;
4331 esign = -1;
4332 if (nexp > 4096)
4333 { /* Punt. Can't handle this without 2 divides. */
4334 emovi (etens[0], tt);
4335 lexp -= tt[E];
4336 k = edivm (tt, yy);
4337 lexp += EXONE;
4338 nexp -= 4096;
4339 }
4340 }
4341 p = &etens[NTEN][0];
4342 emov (eone, xt);
4343 exp = 1;
4344 do
4345 {
4346 if (exp & nexp)
4347 emul (p, xt, xt);
4348 p -= NE;
4349 exp = exp + exp;
4350 }
4351 while (exp <= MAXP);
4352
4353 emovi (xt, tt);
4354 if (esign < 0)
4355 {
4356 lexp -= tt[E];
4357 k = edivm (tt, yy);
4358 lexp += EXONE;
4359 }
4360 else
4361 {
4362 lexp += tt[E];
4363 k = emulm (tt, yy);
4364 lexp -= EXONE - 1;
4365 }
4366
4367 expdon:
4368
4369 /* Round and convert directly to the destination type */
4370 if (oprec == 53)
4371 lexp -= EXONE - 0x3ff;
4372 else if (oprec == 24)
4373 lexp -= EXONE - 0177;
4374#ifdef DEC
4375 else if (oprec == 56)
4376 lexp -= EXONE - 0201;
4377#endif
4378 rndprc = oprec;
4379 emdnorm (yy, k, 0, lexp, 64);
4380
4381 aexit:
4382
4383 rndprc = rndsav;
4384 yy[0] = nsign;
4385 switch (oprec)
4386 {
4387#ifdef DEC
4388 case 56:
4389 todec (yy, y); /* see etodec.c */
4390 break;
4391#endif
4392 case 53:
4393 toe53 (yy, y);
4394 break;
4395 case 24:
4396 toe24 (yy, y);
4397 break;
4398 case 64:
4399 toe64 (yy, y);
4400 break;
4401 case NBITS:
4402 emovo (yy, y);
4403 break;
4404 }
4405}
4406
4407
4408
4409/* y = largest integer not greater than x
4410 * (truncated toward minus infinity)
4411 *
4412 * unsigned EMUSHORT x[NE], y[NE]
4413 *
4414 * efloor (x, y);
4415 */
4416static unsigned EMUSHORT bmask[] =
4417{
4418 0xffff,
4419 0xfffe,
4420 0xfffc,
4421 0xfff8,
4422 0xfff0,
4423 0xffe0,
4424 0xffc0,
4425 0xff80,
4426 0xff00,
4427 0xfe00,
4428 0xfc00,
4429 0xf800,
4430 0xf000,
4431 0xe000,
4432 0xc000,
4433 0x8000,
4434 0x0000,
4435};
4436
4437void
4438efloor (x, y)
4439 unsigned EMUSHORT x[], y[];
4440{
4441 register unsigned EMUSHORT *p;
4442 int e, expon, i;
4443 unsigned EMUSHORT f[NE];
4444
4445 emov (x, f); /* leave in external format */
4446 expon = (int) f[NE - 1];
4447 e = (expon & 0x7fff) - (EXONE - 1);
4448 if (e <= 0)
4449 {
4450 eclear (y);
4451 goto isitneg;
4452 }
4453 /* number of bits to clear out */
4454 e = NBITS - e;
4455 emov (f, y);
4456 if (e <= 0)
4457 return;
4458
4459 p = &y[0];
4460 while (e >= 16)
4461 {
4462 *p++ = 0;
4463 e -= 16;
4464 }
4465 /* clear the remaining bits */
4466 *p &= bmask[e];
4467 /* truncate negatives toward minus infinity */
4468 isitneg:
4469
4470 if ((unsigned EMUSHORT) expon & (unsigned EMUSHORT) 0x8000)
4471 {
4472 for (i = 0; i < NE - 1; i++)
4473 {
4474 if (f[i] != y[i])
4475 {
4476 esub (eone, y, y);
4477 break;
4478 }
4479 }
4480 }
4481}
4482
4483
4484/* unsigned EMUSHORT x[], s[];
4485 * int *exp;
4486 *
4487 * efrexp (x, exp, s);
4488 *
4489 * Returns s and exp such that s * 2**exp = x and .5 <= s < 1.
4490 * For example, 1.1 = 0.55 * 2**1
4491 * Handles denormalized numbers properly using long integer exp.
4492 */
4493void
4494efrexp (x, exp, s)
4495 unsigned EMUSHORT x[];
4496 int *exp;
4497 unsigned EMUSHORT s[];
4498{
4499 unsigned EMUSHORT xi[NI];
4500 EMULONG li;
4501
4502 emovi (x, xi);
4503 li = (EMULONG) ((EMUSHORT) xi[1]);
4504
4505 if (li == 0)
4506 {
4507 li -= enormlz (xi);
4508 }
4509 xi[1] = 0x3ffe;
4510 emovo (xi, s);
4511 *exp = (int) (li - 0x3ffe);
4512}
4513
4514
4515
4516/* unsigned EMUSHORT x[], y[];
4517 * long pwr2;
4518 *
4519 * eldexp (x, pwr2, y);
4520 *
4521 * Returns y = x * 2**pwr2.
4522 */
4523void
4524eldexp (x, pwr2, y)
4525 unsigned EMUSHORT x[];
4526 int pwr2;
4527 unsigned EMUSHORT y[];
4528{
4529 unsigned EMUSHORT xi[NI];
4530 EMULONG li;
4531 int i;
4532
4533 emovi (x, xi);
4534 li = xi[1];
4535 li += pwr2;
4536 i = 0;
4537 emdnorm (xi, i, i, li, 64);
4538 emovo (xi, y);
4539}
4540
4541
4542/* c = remainder after dividing b by a
4543 * Least significant integer quotient bits left in equot[].
4544 */
4545void
4546eremain (a, b, c)
4547 unsigned EMUSHORT a[], b[], c[];
4548{
4549 unsigned EMUSHORT den[NI], num[NI];
4550
4551#ifdef NANS
4552 if ( eisinf (b)
4553 || (ecmp (a, ezero) == 0)
4554 || eisnan (a)
4555 || eisnan (b))
4556 {
4557 enan (c);
4558 return;
4559 }
4560#endif
4561 if (ecmp (a, ezero) == 0)
4562 {
4563 mtherr ("eremain", SING);
4564 eclear (c);
4565 return;
4566 }
4567 emovi (a, den);
4568 emovi (b, num);
4569 eiremain (den, num);
4570 /* Sign of remainder = sign of quotient */
4571 if (a[0] == b[0])
4572 num[0] = 0;
4573 else
4574 num[0] = 0xffff;
4575 emovo (num, c);
4576}
4577
4578void
4579eiremain (den, num)
4580 unsigned EMUSHORT den[], num[];
4581{
4582 EMULONG ld, ln;
4583 unsigned EMUSHORT j;
4584
4585 ld = den[E];
4586 ld -= enormlz (den);
4587 ln = num[E];
4588 ln -= enormlz (num);
4589 ecleaz (equot);
4590 while (ln >= ld)
4591 {
4592 if (ecmpm (den, num) <= 0)
4593 {
4594 esubm (den, num);
4595 j = 1;
4596 }
4597 else
4598 {
4599 j = 0;
4600 }
4601 eshup1 (equot);
4602 equot[NI - 1] |= j;
4603 eshup1 (num);
4604 ln -= 1;
4605 }
4606 emdnorm (num, 0, 0, ln, 0);
4607}
4608
4609/* mtherr.c
4610 *
4611 * Library common error handling routine
4612 *
4613 *
4614 *
4615 * SYNOPSIS:
4616 *
4617 * char *fctnam;
4618 * int code;
4619 * void mtherr ();
4620 *
4621 * mtherr (fctnam, code);
4622 *
4623 *
4624 *
4625 * DESCRIPTION:
4626 *
4627 * This routine may be called to report one of the following
4628 * error conditions (in the include file mconf.h).
4629 *
4630 * Mnemonic Value Significance
4631 *
4632 * DOMAIN 1 argument domain error
4633 * SING 2 function singularity
4634 * OVERFLOW 3 overflow range error
4635 * UNDERFLOW 4 underflow range error
4636 * TLOSS 5 total loss of precision
4637 * PLOSS 6 partial loss of precision
4638 * INVALID 7 NaN - producing operation
4639 * EDOM 33 Unix domain error code
4640 * ERANGE 34 Unix range error code
4641 *
4642 * The default version of the file prints the function name,
4643 * passed to it by the pointer fctnam, followed by the
4644 * error condition. The display is directed to the standard
4645 * output device. The routine then returns to the calling
4646 * program. Users may wish to modify the program to abort by
4647 * calling exit under severe error conditions such as domain
4648 * errors.
4649 *
4650 * Since all error conditions pass control to this function,
4651 * the display may be easily changed, eliminated, or directed
4652 * to an error logging device.
4653 *
4654 * SEE ALSO:
4655 *
4656 * mconf.h
4657 *
4658 */
4659\f
4660/*
4661Cephes Math Library Release 2.0: April, 1987
4662Copyright 1984, 1987 by Stephen L. Moshier
4663Direct inquiries to 30 Frost Street, Cambridge, MA 02140
4664*/
4665
4666/* include "mconf.h" */
4667
4668/* Notice: the order of appearance of the following
4669 * messages is bound to the error codes defined
4670 * in mconf.h.
4671 */
4672#define NMSGS 8
4673static char *ermsg[NMSGS] =
4674{
4675 "unknown", /* error code 0 */
4676 "domain", /* error code 1 */
4677 "singularity", /* et seq. */
4678 "overflow",
4679 "underflow",
4680 "total loss of precision",
4681 "partial loss of precision",
4682 "invalid operation"
4683};
4684
4685int merror = 0;
4686extern int merror;
4687
4688void
4689mtherr (name, code)
4690 char *name;
4691 int code;
4692{
4693 char errstr[80];
4694
4695 /* Display string passed by calling program,
4696 * which is supposed to be the name of the
4697 * function in which the error occurred.
4698 */
4699
4700 /* Display error message defined
4701 * by the code argument.
4702 */
4703 if ((code <= 0) || (code >= NMSGS))
4704 code = 0;
4705 sprintf (errstr, " %s %s error", name, ermsg[code]);
4706 if (extra_warnings)
4707 warning (errstr);
4708 /* Set global error message word */
4709 merror = code + 1;
4710
4711 /* Return to calling
4712 * program
4713 */
4714}
4715
4716/* Here is etodec.c .
4717 *
4718 */
4719
4720/*
4721; convert DEC double precision to e type
4722; double d;
4723; EMUSHORT e[NE];
4724; dectoe (&d, e);
4725*/
4726void
4727dectoe (d, e)
4728 unsigned EMUSHORT *d;
4729 unsigned EMUSHORT *e;
4730{
4731 unsigned EMUSHORT y[NI];
4732 register unsigned EMUSHORT r, *p;
4733
4734 ecleaz (y); /* start with a zero */
4735 p = y; /* point to our number */
4736 r = *d; /* get DEC exponent word */
4737 if (*d & (unsigned int) 0x8000)
4738 *p = 0xffff; /* fill in our sign */
4739 ++p; /* bump pointer to our exponent word */
4740 r &= 0x7fff; /* strip the sign bit */
4741 if (r == 0) /* answer = 0 if high order DEC word = 0 */
4742 goto done;
4743
4744
4745 r >>= 7; /* shift exponent word down 7 bits */
4746 r += EXONE - 0201; /* subtract DEC exponent offset */
4747 /* add our e type exponent offset */
4748 *p++ = r; /* to form our exponent */
4749
4750 r = *d++; /* now do the high order mantissa */
4751 r &= 0177; /* strip off the DEC exponent and sign bits */
4752 r |= 0200; /* the DEC understood high order mantissa bit */
4753 *p++ = r; /* put result in our high guard word */
4754
4755 *p++ = *d++; /* fill in the rest of our mantissa */
4756 *p++ = *d++;
4757 *p = *d;
4758
4759 eshdn8 (y); /* shift our mantissa down 8 bits */
4760 done:
4761 emovo (y, e);
4762}
4763
4764
4765
4766/*
4767; convert e type to DEC double precision
4768; double d;
4769; EMUSHORT e[NE];
4770; etodec (e, &d);
4771*/
4772#if 0
4773static unsigned EMUSHORT decbit[NI] = {0, 0, 0, 0, 0, 0, 0200, 0};
4774
4775void
4776etodec (x, d)
4777 unsigned EMUSHORT *x, *d;
4778{
4779 unsigned EMUSHORT xi[NI];
4780 register unsigned EMUSHORT r;
4781 int i, j;
4782
4783 emovi (x, xi);
4784 *d = 0;
4785 if (xi[0] != 0)
4786 *d = 0100000;
4787 r = xi[E];
4788 if (r < (EXONE - 128))
4789 goto zout;
4790 i = xi[M + 4];
4791 if ((i & 0200) != 0)
4792 {
4793 if ((i & 0377) == 0200)
4794 {
4795 if ((i & 0400) != 0)
4796 {
4797 /* check all less significant bits */
4798 for (j = M + 5; j < NI; j++)
4799 {
4800 if (xi[j] != 0)
4801 goto yesrnd;
4802 }
4803 }
4804 goto nornd;
4805 }
4806 yesrnd:
4807 eaddm (decbit, xi);
4808 r -= enormlz (xi);
4809 }
4810
4811 nornd:
4812
4813 r -= EXONE;
4814 r += 0201;
4815 if (r < 0)
4816 {
4817 zout:
4818 *d++ = 0;
4819 *d++ = 0;
4820 *d++ = 0;
4821 *d++ = 0;
4822 return;
4823 }
4824 if (r >= 0377)
4825 {
4826 *d++ = 077777;
4827 *d++ = -1;
4828 *d++ = -1;
4829 *d++ = -1;
4830 return;
4831 }
4832 r &= 0377;
4833 r <<= 7;
4834 eshup8 (xi);
4835 xi[M] &= 0177;
4836 r |= xi[M];
4837 *d++ |= r;
4838 *d++ = xi[M + 1];
4839 *d++ = xi[M + 2];
4840 *d++ = xi[M + 3];
4841}
4842
4843#else
4844
4845void
4846etodec (x, d)
4847 unsigned EMUSHORT *x, *d;
4848{
4849 unsigned EMUSHORT xi[NI];
4850 EMULONG exp;
4851 int rndsav;
4852
4853 emovi (x, xi);
4854 exp = (EMULONG) xi[E] - (EXONE - 0201); /* adjust exponent for offsets */
4855/* round off to nearest or even */
4856 rndsav = rndprc;
4857 rndprc = 56;
4858 emdnorm (xi, 0, 0, exp, 64);
4859 rndprc = rndsav;
4860 todec (xi, d);
4861}
4862
4863void
4864todec (x, y)
4865 unsigned EMUSHORT *x, *y;
4866{
4867 unsigned EMUSHORT i;
4868 unsigned EMUSHORT *p;
4869
4870 p = x;
4871 *y = 0;
4872 if (*p++)
4873 *y = 0100000;
4874 i = *p++;
4875 if (i == 0)
4876 {
4877 *y++ = 0;
4878 *y++ = 0;
4879 *y++ = 0;
4880 *y++ = 0;
4881 return;
4882 }
4883 if (i > 0377)
4884 {
4885 *y++ |= 077777;
4886 *y++ = 0xffff;
4887 *y++ = 0xffff;
4888 *y++ = 0xffff;
4889#ifdef ERANGE
4890 errno = ERANGE;
4891#endif
4892 return;
4893 }
4894 i &= 0377;
4895 i <<= 7;
4896 eshup8 (x);
4897 x[M] &= 0177;
4898 i |= x[M];
4899 *y++ |= i;
4900 *y++ = x[M + 1];
4901 *y++ = x[M + 2];
4902 *y++ = x[M + 3];
4903}
4904
4905#endif /* not 0 */
4906
4907
4908/* Output a binary NaN bit pattern in the target machine's format. */
4909
4910/* If special NaN bit patterns are required, define them in tm.h
4911 as arrays of unsigned 16-bit shorts. Otherwise, use the default
4912 patterns here. */
4913#ifdef TFMODE_NAN
4914TFMODE_NAN;
4915#else
4916#ifdef MIEEE
4917unsigned EMUSHORT TFnan[8] =
4918 {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
4919#endif
4920#ifdef IBMPC
4921unsigned EMUSHORT TFnan[8] = {0, 0, 0, 0, 0, 0, 0x8000, 0xffff};
4922#endif
4923#endif
4924
4925#ifdef XFMODE_NAN
4926XFMODE_NAN;
4927#else
4928#ifdef MIEEE
4929unsigned EMUSHORT XFnan[6] = {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
4930#endif
4931#ifdef IBMPC
4932unsigned EMUSHORT XFnan[6] = {0, 0, 0, 0xc000, 0xffff, 0};
4933#endif
4934#endif
4935
4936#ifdef DFMODE_NAN
4937DFMODE_NAN;
4938#else
4939#ifdef MIEEE
4940unsigned EMUSHORT DFnan[4] = {0x7fff, 0xffff, 0xffff, 0xffff};
4941#endif
4942#ifdef IBMPC
4943unsigned EMUSHORT DFnan[4] = {0, 0, 0, 0xfff8};
4944#endif
4945#endif
4946
4947#ifdef SFMODE_NAN
4948SFMODE_NAN;
4949#else
4950#ifdef MIEEE
4951unsigned EMUSHORT SFnan[2] = {0x7fff, 0xffff};
4952#endif
4953#ifdef IBMPC
4954unsigned EMUSHORT SFnan[2] = {0, 0xffc0};
4955#endif
4956#endif
4957
4958
4959void
4960make_nan (nan, mode)
4961unsigned EMUSHORT *nan;
4962enum machine_mode mode;
4963{
4964 int i, n;
4965 unsigned EMUSHORT *p;
4966
4967 switch (mode)
4968 {
4969/* Possibly the `reserved operand' patterns on a VAX can be
4970 used like NaN's, but probably not in the same way as IEEE. */
4971#ifndef DEC
4972 case TFmode:
4973 n = 8;
4974 p = TFnan;
4975 break;
4976 case XFmode:
4977 n = 6;
4978 p = XFnan;
4979 break;
4980 case DFmode:
4981 n = 4;
4982 p = DFnan;
4983 break;
4984 case SFmode:
4985 n = 2;
4986 p = SFnan;
4987 break;
4988#endif
4989 default:
4990 abort ();
4991 }
4992 for (i=0; i < n; i++)
4993 *nan++ = *p++;
4994}
4995
4996/* Convert an SFmode target `float' value to a REAL_VALUE_TYPE.
4997 This is the inverse of the function `etarsingle' invoked by
4998 REAL_VALUE_TO_TARGET_SINGLE. */
4999
5000REAL_VALUE_TYPE
5001ereal_from_float (f)
5002 unsigned long f;
5003{
5004 REAL_VALUE_TYPE r;
5005 unsigned EMUSHORT s[2];
5006 unsigned EMUSHORT e[NE];
5007
5008 /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
5009 This is the inverse operation to what the function `endian' does. */
5010#if WORDS_BIG_ENDIAN
5011 s[0] = (unsigned EMUSHORT) (f >> 16);
5012 s[1] = (unsigned EMUSHORT) f;
5013#else
5014 s[0] = (unsigned EMUSHORT) f;
5015 s[1] = (unsigned EMUSHORT) (f >> 16);
5016#endif
5017 /* Convert and promote the target float to E-type. */
5018 e24toe (s, e);
5019 /* Output E-type to REAL_VALUE_TYPE. */
5020 PUT_REAL (e, &r);
5021 return r;
5022}
5023
5024/* Convert a DFmode target `double' value to a REAL_VALUE_TYPE.
5025 This is the inverse of the function `etardouble' invoked by
5026 REAL_VALUE_TO_TARGET_DOUBLE.
5027
5028 The DFmode is stored as an array of longs (i.e., HOST_WIDE_INTs)
5029 with 32 bits of the value per each long. The first element
5030 of the input array holds the bits that would come first in the
5031 target computer's memory. */
5032
5033REAL_VALUE_TYPE
5034ereal_from_double (d)
5035 unsigned long d[];
5036{
5037 REAL_VALUE_TYPE r;
5038 unsigned EMUSHORT s[4];
5039 unsigned EMUSHORT e[NE];
5040
5041 /* Convert array of 32 bit pieces to equivalent array of 16 bit pieces.
5042 This is the inverse of `endian'. */
5043#if WORDS_BIG_ENDIAN
5044 s[0] = (unsigned EMUSHORT) (d[0] >> 16);
5045 s[1] = (unsigned EMUSHORT) d[0];
5046 s[2] = (unsigned EMUSHORT) (d[1] >> 16);
5047 s[3] = (unsigned EMUSHORT) d[1];
5048#else
5049 s[0] = (unsigned EMUSHORT) d[0];
5050 s[1] = (unsigned EMUSHORT) (d[0] >> 16);
5051 s[2] = (unsigned EMUSHORT) d[1];
5052 s[3] = (unsigned EMUSHORT) (d[1] >> 16);
5053#endif
5054 /* Convert target double to E-type. */
5055 e53toe (s, e);
5056 /* Output E-type to REAL_VALUE_TYPE. */
5057 PUT_REAL (e, &r);
5058 return r;
5059}
5060#endif /* EMU_NON_COMPILE not defined */
Continuous source code history from original PDP-7 UNIX to FreeBSD 2, the first fully unencumbered FreeBSD release.