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