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