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 (moshier@world.std.com).
5 Copyright (C) 1993 Free Software Foundation, Inc.
7 This file is part of GNU CC.
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)
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.
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. */
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).
35 To support cross compilation between IEEE and VAX floating
36 point formats, define REAL_ARITHMETIC in the tm.h file.
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.
44 The emulator defaults to the host's floating point format so that
45 its decimal conversion functions can be used if desired (see
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 */
61 /* Type of computer arithmetic.
62 * Only one of DEC, MIEEE, IBMPC, or UNK should get defined.
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.
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.
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.
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
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. */
92 /* The following converts gcc macros into the ones used by this file. */
94 /* REAL_ARITHMETIC defined means that macros in real.h are
95 defined to call emulator functions. */
96 #ifdef REAL_ARITHMETIC
98 #if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
99 /* PDP-11, Pro350, VAX: */
101 #else /* it's not VAX */
102 #if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
104 /* Motorola IEEE, high order words come first (Sun workstation): */
106 #else /* not big-endian */
107 /* Intel IEEE, low order words come first:
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
115 #endif /* not IEEE */
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
126 #else /* it's not VAX */
127 #if HOST_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
128 #ifdef HOST_WORDS_BIG_ENDIAN
130 #else /* not big-endian */
132 #endif /* big-endian */
133 #else /* it's not IEEE either */
134 unknown arithmetic type
136 #endif /* not IEEE */
139 #endif /* REAL_ARITHMETIC not defined */
141 /* Define INFINITY for support of infinity.
142 Define NANS for support of Not-a-Number's (NaN's). */
148 /* Support of NaNs requires support of infinity. */
157 * Include file for extended precision arithmetic programs.
160 /* Number of 16 bit words in external e type format */
163 /* Number of 16 bit words in internal format */
166 /* Array offset to exponent */
169 /* Array offset to high guard word */
172 /* Number of bits of precision */
173 #define NBITS ((NI-4)*16)
175 /* Maximum number of decimal digits in ASCII conversion
178 #define NDEC (NBITS*8/27)
180 /* The exponent of 1.0 */
181 #define EXONE (0x3fff)
183 /* Find a host integer type that is at least 16 bits wide,
184 and another type at least twice whatever that size is. */
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)
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)
196 #if HOST_BITS_PER_INT >= 16
198 #define EMUSHORT_SIZE HOST_BITS_PER_INT
199 #define EMULONG_SIZE (2 * HOST_BITS_PER_INT)
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)
206 /* You will have to modify this program to have a smaller unit size. */
207 #define EMU_NON_COMPILE
213 #if HOST_BITS_PER_SHORT >= EMULONG_SIZE
214 #define EMULONG short
216 #if HOST_BITS_PER_INT >= EMULONG_SIZE
219 #if HOST_BITS_PER_LONG >= EMULONG_SIZE
222 #if HOST_BITS_PER_LONG_LONG >= EMULONG_SIZE
223 #define EMULONG long long int
225 /* You will have to modify this program to have a smaller unit size. */
226 #define EMU_NON_COMPILE
233 /* The host interface doesn't work if no 16-bit size exists. */
234 #if EMUSHORT_SIZE != 16
235 #define EMU_NON_COMPILE
238 /* OK to continue compilation. */
239 #ifndef EMU_NON_COMPILE
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. */
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 */
251 #ifdef REAL_ARITHMETIC
252 /* Emulator uses target format internally
253 but host stores it in host endian-ness. */
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))
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)
269 #define PUT_REAL(e,r) \
270 do { EMUSHORT w[4]; \
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)
277 #endif /* endian-ness differs */
279 #else /* not REAL_ARITHMETIC */
281 /* emulator uses host format */
282 #define GET_REAL(r,e) e53toe ((r), (e))
283 #define PUT_REAL(e,r) etoe53 ((e), (r))
285 #endif /* not REAL_ARITHMETIC */
286 #endif /* no XFmode */
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 ();
303 extern unsigned EMUSHORT ezero[], ehalf[], eone[], etwo[];
304 extern unsigned EMUSHORT elog2[], esqrt2[];
306 /* Pack output array with 32-bit numbers obtained from
307 array containing 16-bit numbers, swapping ends if required. */
310 unsigned EMUSHORT e[];
312 enum machine_mode mode;
322 /* Swap halfwords in the third long. */
323 th = (unsigned long) e[4] & 0xffff;
324 t = (unsigned long) e[5] & 0xffff;
327 /* fall into the double case */
331 /* swap halfwords in the second word */
332 th = (unsigned long) e[2] & 0xffff;
333 t = (unsigned long) e[3] & 0xffff;
336 /* fall into the float case */
340 /* swap halfwords in the first word */
341 th = (unsigned long) e[0] & 0xffff;
342 t = (unsigned long) e[1] & 0xffff;
353 /* Pack the output array without swapping. */
360 /* Pack the third long.
361 Each element of the input REAL_VALUE_TYPE array has 16 bit useful bits
363 th = (unsigned long) e[5] & 0xffff;
364 t = (unsigned long) e[4] & 0xffff;
367 /* fall into the double case */
371 /* pack the second long */
372 th = (unsigned long) e[3] & 0xffff;
373 t = (unsigned long) e[2] & 0xffff;
376 /* fall into the float case */
380 /* pack the first long */
381 th = (unsigned long) e[1] & 0xffff;
382 t = (unsigned long) e[0] & 0xffff;
395 /* This is the implementation of the REAL_ARITHMETIC macro.
398 earith (value, icode, r1, r2)
399 REAL_VALUE_TYPE *value;
404 unsigned EMUSHORT d1[NE], d2[NE], v[NE];
410 /* Return NaN input back to the caller. */
413 PUT_REAL (d1, value);
418 PUT_REAL (d2, value);
422 code = (enum tree_code) icode;
430 esub (d2, d1, v); /* d1 - d2 */
438 #ifndef REAL_INFINITY
439 if (ecmp (d2, ezero) == 0)
449 ediv (d2, d1, v); /* d1/d2 */
452 case MIN_EXPR: /* min (d1,d2) */
453 if (ecmp (d1, d2) < 0)
459 case MAX_EXPR: /* max (d1,d2) */
460 if (ecmp (d1, d2) > 0)
473 /* Truncate REAL_VALUE_TYPE toward zero to signed HOST_WIDE_INT
474 * implements REAL_VALUE_RNDZINT (x) (etrunci (x))
480 unsigned EMUSHORT f[NE], g[NE];
496 /* Truncate REAL_VALUE_TYPE toward zero to unsigned HOST_WIDE_INT
497 * implements REAL_VALUE_UNSIGNED_RNDZINT (x) (etruncui (x))
503 unsigned EMUSHORT f[NE], g[NE];
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.
529 unsigned EMUSHORT tem[NE], e[NE];
554 /* Expansion of REAL_NEGATE.
560 unsigned EMUSHORT e[NE];
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).
583 unsigned EMUSHORT f[NE], g[NE];
590 warning ("conversion from NaN to int");
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).
609 unsigned EMUSHORT f[NE], g[NE];
616 warning ("conversion from NaN to unsigned int");
622 return ((unsigned int)l);
626 /* REAL_VALUE_FROM_INT macro.
629 ereal_from_int (d, i, j)
633 unsigned EMUSHORT df[NE], dg[NE];
642 /* complement and add 1 */
649 eldexp (eone, HOST_BITS_PER_LONG, df);
660 /* REAL_VALUE_FROM_UNSIGNED_INT macro.
663 ereal_from_uint (d, i, j)
667 unsigned EMUSHORT df[NE], dg[NE];
668 unsigned long low, high;
672 eldexp (eone, HOST_BITS_PER_LONG, df);
681 /* REAL_VALUE_TO_INT macro
684 ereal_to_int (low, high, rr)
688 unsigned EMUSHORT d[NE], df[NE], dg[NE], dh[NE];
695 warning ("conversion from NaN to int");
701 /* convert positive value */
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);
715 /* complement and add 1 */
725 /* REAL_VALUE_LDEXP macro.
732 unsigned EMUSHORT e[NE], y[NE];
745 /* These routines are conditionally compiled because functions
746 * of the same names may be defined in fold-const.c. */
747 #ifdef REAL_ARITHMETIC
749 /* Check for infinity in a REAL_VALUE_TYPE. */
754 unsigned EMUSHORT e[NE];
765 /* Check whether a REAL_VALUE_TYPE item is a NaN. */
771 unsigned EMUSHORT e[NE];
782 /* Check for a negative REAL_VALUE_TYPE number.
783 * this means strictly less than zero, not -0.
790 unsigned EMUSHORT e[NE];
793 if (ecmp (e, ezero) == -1)
798 /* Expansion of REAL_VALUE_TRUNCATE.
799 * The result is in floating point, rounded to nearest or even.
802 real_value_truncate (mode, arg)
803 enum machine_mode mode;
806 unsigned EMUSHORT e[NE], t[NE];
843 #endif /* REAL_ARITHMETIC defined */
845 /* Target values are arrays of host longs. A long is guaranteed
846 to be at least 32 bits wide. */
852 unsigned EMUSHORT e[NE];
856 endian (e, l, XFmode);
864 unsigned EMUSHORT e[NE];
868 endian (e, l, DFmode);
875 unsigned EMUSHORT e[NE];
880 endian (e, &l, SFmode);
885 ereal_to_decimal (x, s)
889 unsigned EMUSHORT e[NE];
897 REAL_VALUE_TYPE x, y;
899 unsigned EMUSHORT ex[NE], ey[NE];
903 return (ecmp (ex, ey));
910 unsigned EMUSHORT ex[NE];
913 return (eisneg (ex));
916 /* End of REAL_ARITHMETIC interface */
920 * Extended precision IEEE binary floating point arithmetic routines
922 * Numbers are stored in C language as arrays of 16-bit unsigned
923 * short integers. The arguments of the routines are pointers to
927 * External e type data structure, simulates Intel 8087 chip
928 * temporary real format but possibly with a larger significand:
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)
936 * Internal data structure of a number (a "word" is 16 bits):
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)
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)
949 * Routines for external format numbers
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
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
972 * emul (a, b, c) c = b * a
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
991 * Routines for internal format numbers
993 * eaddm (ai, bi) add significands, bi = bi + ai
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
1018 * The result is always normalized and rounded to NI-4 word precision
1019 * after each arithmetic operation.
1021 * Exception flags are NOT fully supported.
1023 * Signaling NaN's are NOT supported; they are treated the same
1026 * Define INFINITY for support of infinity; otherwise a
1027 * saturation arithmetic is implemented.
1029 * Define NANS for support of Not-a-Number items; otherwise the
1030 * arithmetic will never produce a NaN output, and might be confused
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'
1037 * Denormals are always supported here where appropriate (e.g., not
1038 * for conversion to DEC numbers).
1045 * Common include file for math routines
1051 * #include "mconf.h"
1057 * This file contains definitions for error codes that are
1058 * passed to the common error handling routine mtherr
1061 * The file also includes a conditional assembly definition
1062 * for the type of computer arithmetic (Intel IEEE, DEC, Motorola
1063 * IEEE, or UNKnown).
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.
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.
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.
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.
1091 * For ANSI C compatibility, define ANSIC equal to 1. Currently
1092 * this affects only the atan2 function and others that use it.
1095 /* Constant definitions for math error conditions. */
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 */
1105 /* e type constants used by high precision check routines */
1107 /*include "ehead.h"*/
1109 unsigned EMUSHORT ezero[NE] =
1111 0, 0000000, 0000000, 0000000, 0000000, 0000000,};
1112 extern unsigned EMUSHORT ezero[];
1115 unsigned EMUSHORT ehalf[NE] =
1117 0, 0000000, 0000000, 0000000, 0100000, 0x3ffe,};
1118 extern unsigned EMUSHORT ehalf[];
1121 unsigned EMUSHORT eone[NE] =
1123 0, 0000000, 0000000, 0000000, 0100000, 0x3fff,};
1124 extern unsigned EMUSHORT eone[];
1127 unsigned EMUSHORT etwo[NE] =
1129 0, 0000000, 0000000, 0000000, 0100000, 0040000,};
1130 extern unsigned EMUSHORT etwo[];
1133 unsigned EMUSHORT e32[NE] =
1135 0, 0000000, 0000000, 0000000, 0100000, 0040004,};
1136 extern unsigned EMUSHORT e32[];
1138 /* 6.93147180559945309417232121458176568075500134360255E-1 */
1139 unsigned EMUSHORT elog2[NE] =
1141 0xc9e4, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1142 extern unsigned EMUSHORT elog2[];
1144 /* 1.41421356237309504880168872420969807856967187537695E0 */
1145 unsigned EMUSHORT esqrt2[NE] =
1147 0x597e, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1148 extern unsigned EMUSHORT esqrt2[];
1151 * 1.12837916709551257389615890312154517168810125865800E0 */
1152 unsigned EMUSHORT eoneopi[NE] =
1154 0x71d5, 0x688d, 0012333, 0135202, 0110156, 0x3fff,};
1155 extern unsigned EMUSHORT eoneopi[];
1157 /* 3.14159265358979323846264338327950288419716939937511E0 */
1158 unsigned EMUSHORT epi[NE] =
1160 0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1161 extern unsigned EMUSHORT epi[];
1163 /* 5.7721566490153286060651209008240243104215933593992E-1 */
1164 unsigned EMUSHORT eeul[NE] =
1166 0xd1be, 0xc7a4, 0076660, 0063743, 0111704, 0x3ffe,};
1167 extern unsigned EMUSHORT eeul[];
1176 /* Control register for rounding precision.
1177 * This can be set to 80 (if NE=6), 64, 56, 53, or 24 bits.
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 ();
1198 ; Clear out entire external format number.
1200 ; unsigned EMUSHORT x[];
1206 register unsigned EMUSHORT *x;
1210 for (i = 0; i < NE; i++)
1216 /* Move external format number from a to b.
1223 register unsigned EMUSHORT *a, *b;
1227 for (i = 0; i < NE; i++)
1233 ; Absolute value of external format number
1241 unsigned EMUSHORT x[]; /* x is the memory address of a short */
1244 x[NE - 1] &= 0x7fff; /* sign is top bit of last word of external format */
1251 ; Negate external format number
1253 ; unsigned EMUSHORT x[NE];
1259 unsigned EMUSHORT x[];
1266 x[NE - 1] ^= 0x8000; /* Toggle the sign bit */
1271 /* Return 1 if external format number is negative,
1272 * else return zero, including when it is a NaN.
1276 unsigned EMUSHORT x[];
1283 if (x[NE - 1] & 0x8000)
1290 /* Return 1 if external format number is infinity.
1295 unsigned EMUSHORT x[];
1302 if ((x[NE - 1] & 0x7fff) == 0x7fff)
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
1315 unsigned EMUSHORT x[];
1320 /* NaN has maximum exponent */
1321 if ((x[NE - 1] & 0x7fff) != 0x7fff)
1323 /* ... and non-zero significand field. */
1324 for (i = 0; i < NE - 1; i++)
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. */
1340 register unsigned EMUSHORT *x;
1345 for (i = 0; i < NE - 1; i++)
1349 for (i = 0; i < NE - 1; i++)
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. */
1379 register unsigned EMUSHORT *x;
1383 for (i = 0; i < NE - 2; i++)
1390 /* Move in external format number,
1391 * converting it to internal format.
1395 unsigned EMUSHORT *a, *b;
1397 register unsigned EMUSHORT *p, *q;
1401 p = a + (NE - 1); /* point to last word of external number */
1402 /* get the sign bit */
1407 /* get the exponent */
1409 *q++ &= 0x7fff; /* delete the sign bit */
1411 if ((*(q - 1) & 0x7fff) == 0x7fff)
1417 for (i = 3; i < NI; i++)
1422 for (i = 2; i < NI; i++)
1427 /* clear high guard word */
1429 /* move in the significand */
1430 for (i = 0; i < NE - 1; i++)
1432 /* clear low guard word */
1437 /* Move internal format number out,
1438 * converting it to external format.
1442 unsigned EMUSHORT *a, *b;
1444 register unsigned EMUSHORT *p, *q;
1445 unsigned EMUSHORT i;
1448 q = b + (NE - 1); /* point to output exponent */
1449 /* combine sign and exponent */
1452 *q-- = *p++ | 0x8000;
1456 if (*(p - 1) == 0x7fff)
1469 /* skip over guard word */
1471 /* move the significand */
1472 for (i = 0; i < NE - 1; i++)
1479 /* Clear out internal format number.
1484 register unsigned EMUSHORT *xi;
1488 for (i = 0; i < NI; i++)
1493 /* same, but don't touch the sign. */
1497 register unsigned EMUSHORT *xi;
1502 for (i = 0; i < NI - 1; i++)
1508 /* Move internal format number from a to b.
1512 register unsigned EMUSHORT *a, *b;
1516 for (i = 0; i < NI - 1; i++)
1518 /* clear low guard word */
1522 /* Generate internal format NaN.
1523 The explicit pattern for this is maximum exponent and
1524 top two significand bits set. */
1528 unsigned EMUSHORT x[];
1536 /* Return nonzero if internal format number is a NaN. */
1540 unsigned EMUSHORT x[];
1544 if ((x[E] & 0x7fff) == 0x7fff)
1546 for (i = M + 1; i < NI; i++)
1555 /* Fill internal format number with infinity pattern.
1556 This has maximum exponent and significand all zeros. */
1560 unsigned EMUSHORT x[];
1567 /* Return nonzero if internal format number is infinite. */
1571 unsigned EMUSHORT x[];
1578 if ((x[E] & 0x7fff) == 0x7fff)
1585 ; Compare significands of numbers in internal format.
1586 ; Guard words are included in the comparison.
1588 ; unsigned EMUSHORT a[NI], b[NI];
1591 ; for the significands:
1592 ; returns +1 if a > b
1598 register unsigned EMUSHORT *a, *b;
1602 a += M; /* skip up to significand area */
1604 for (i = M; i < NI; i++)
1612 if (*(--a) > *(--b))
1620 ; Shift significand down by 1 bit
1625 register unsigned EMUSHORT *x;
1627 register unsigned EMUSHORT bits;
1630 x += M; /* point to significand area */
1633 for (i = M; i < NI; i++)
1648 ; Shift significand up by 1 bit
1653 register unsigned EMUSHORT *x;
1655 register unsigned EMUSHORT bits;
1661 for (i = M; i < NI; i++)
1676 ; Shift significand down by 8 bits
1681 register unsigned EMUSHORT *x;
1683 register unsigned EMUSHORT newbyt, oldbyt;
1688 for (i = M; i < NI; i++)
1699 ; Shift significand up by 8 bits
1704 register unsigned EMUSHORT *x;
1707 register unsigned EMUSHORT newbyt, oldbyt;
1712 for (i = M; i < NI; i++)
1723 ; Shift significand up by 16 bits
1728 register unsigned EMUSHORT *x;
1731 register unsigned EMUSHORT *p;
1736 for (i = M; i < NI - 1; i++)
1743 ; Shift significand down by 16 bits
1748 register unsigned EMUSHORT *x;
1751 register unsigned EMUSHORT *p;
1756 for (i = M; i < NI - 1; i++)
1769 unsigned EMUSHORT *x, *y;
1771 register unsigned EMULONG a;
1778 for (i = M; i < NI; i++)
1780 a = (unsigned EMULONG) (*x) + (unsigned EMULONG) (*y) + carry;
1785 *y = (unsigned EMUSHORT) a;
1792 ; Subtract significands
1798 unsigned EMUSHORT *x, *y;
1807 for (i = M; i < NI; i++)
1809 a = (unsigned EMULONG) (*y) - (unsigned EMULONG) (*x) - carry;
1814 *y = (unsigned EMUSHORT) a;
1821 /* Divide significands */
1823 static unsigned EMUSHORT equot[NI];
1827 unsigned EMUSHORT den[], num[];
1830 register unsigned EMUSHORT *p, *q;
1831 unsigned EMUSHORT j;
1837 for (i = M; i < NI; i++)
1842 /* Use faster compare and subtraction if denominator
1843 * has only 15 bits of significance.
1848 for (i = M + 3; i < NI; i++)
1853 if ((den[M + 1] & 1) != 0)
1861 for (i = 0; i < NBITS + 2; i++)
1879 /* The number of quotient bits to calculate is
1880 * NBITS + 1 scaling guard bit + 1 roundoff bit.
1885 for (i = 0; i < NBITS + 2; i++)
1887 if (ecmpm (den, num) <= 0)
1890 j = 1; /* quotient bit = 1 */
1904 /* test for nonzero remainder after roundoff bit */
1907 for (i = M; i < NI; i++)
1915 for (i = 0; i < NI; i++)
1921 /* Multiply significands */
1924 unsigned EMUSHORT a[], b[];
1926 unsigned EMUSHORT *p, *q;
1931 for (i = M; i < NI; i++)
1936 while (*p == 0) /* significand is not supposed to be all zero */
1941 if ((*p & 0xff) == 0)
1949 for (i = 0; i < k; i++)
1953 /* remember if there were any nonzero bits shifted out */
1960 for (i = 0; i < NI; i++)
1963 /* return flag for lost nonzero bits */
1970 * Normalize and round off.
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.
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.
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.
1984 * Input "rcntrl" is the rounding control.
1987 static int rlast = -1;
1989 static unsigned EMUSHORT rmsk = 0;
1990 static unsigned EMUSHORT rmbit = 0;
1991 static unsigned EMUSHORT rebit = 0;
1993 static unsigned EMUSHORT rbit[NI];
1996 emdnorm (s, lost, subflg, exp, rcntrl)
1997 unsigned EMUSHORT s[];
2004 unsigned EMUSHORT r;
2009 /* a blank significand could mean either zero or infinity. */
2022 if ((j > NBITS) && (exp < 32767))
2030 if (exp > (EMULONG) (-NBITS - 1))
2043 /* Round off, unless told not to by rcntrl. */
2046 /* Set up rounding parameters if the control register changed. */
2047 if (rndprc != rlast)
2054 rw = NI - 1; /* low guard word */
2069 /* For DEC arithmetic */
2118 /* These tests assume NI = 8 */
2138 if ((r & rmbit) != 0)
2143 { /* round to even */
2144 if ((s[re] & rebit) == 0)
2156 if ((rndprc < 64) && (exp <= 0))
2161 { /* overflow on roundoff */
2174 for (i = 2; i < NI - 1; i++)
2177 warning ("floating point overflow");
2181 for (i = M + 1; i < NI - 1; i++)
2199 s[1] = (unsigned EMUSHORT) exp;
2205 ; Subtract external format numbers.
2207 ; unsigned EMUSHORT a[NE], b[NE], c[NE];
2208 ; esub (a, b, c); c = b - a
2211 static int subflg = 0;
2215 unsigned EMUSHORT *a, *b, *c;
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))
2234 mtherr ("esub", INVALID);
2247 ; unsigned EMUSHORT a[NE], b[NE], c[NE];
2248 ; eadd (a, b, c); c = b + a
2252 unsigned EMUSHORT *a, *b, *c;
2256 /* NaN plus anything is a NaN. */
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))
2272 mtherr ("esub", INVALID);
2283 unsigned EMUSHORT *a, *b, *c;
2285 unsigned EMUSHORT ai[NI], bi[NI], ci[NI];
2287 EMULONG lt, lta, ltb;
2308 /* compare exponents */
2313 { /* put the larger number in bi */
2323 if (lt < (EMULONG) (-NBITS - 1))
2324 goto done; /* answer same as larger addend */
2326 lost = eshift (ai, k); /* shift the smaller number down */
2330 /* exponents were the same, so must compare significands */
2333 { /* the numbers are identical in magnitude */
2334 /* if different signs, result is zero */
2340 /* if same sign, result is double */
2341 /* double denomalized tiny number */
2342 if ((bi[E] == 0) && ((bi[3] & 0x8000) == 0))
2347 /* add 1 to exponent unless both are zero! */
2348 for (j = 1; j < NI - 1; j++)
2352 /* This could overflow, but let emovo take care of that. */
2357 bi[E] = (unsigned EMUSHORT) ltb;
2361 { /* put the larger number in bi */
2377 emdnorm (bi, lost, subflg, ltb, 64);
2388 ; unsigned EMUSHORT a[NE], b[NE], c[NE];
2389 ; ediv (a, b, c); c = b / a
2393 unsigned EMUSHORT *a, *b, *c;
2395 unsigned EMUSHORT ai[NI], bi[NI];
2397 EMULONG lt, lta, ltb;
2400 /* Return any NaN input. */
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)))
2415 mtherr ("ediv", INVALID);
2420 /* Infinity over anything else is infinity. */
2424 if (eisneg (a) ^ eisneg (b))
2425 *(c + (NE - 1)) = 0x8000;
2427 *(c + (NE - 1)) = 0;
2431 /* Anything else over infinity is zero. */
2443 { /* See if numerator is zero. */
2444 for (i = 1; i < NI - 1; i++)
2448 ltb -= enormlz (bi);
2458 { /* possible divide by zero */
2459 for (i = 1; i < NI - 1; i++)
2463 lta -= enormlz (ai);
2468 *(c + (NE - 1)) = 0;
2470 *(c + (NE - 1)) = 0x8000;
2471 /* Divide by zero is not an invalid operation.
2472 It is a divide-by-zero operation! */
2474 mtherr ("ediv", SING);
2480 /* calculate exponent */
2481 lt = ltb - lta + EXONE;
2482 emdnorm (bi, i, 0, lt, 64);
2496 ; unsigned EMUSHORT a[NE], b[NE], c[NE];
2497 ; emul (a, b, c); c = b * a
2501 unsigned EMUSHORT *a, *b, *c;
2503 unsigned EMUSHORT ai[NI], bi[NI];
2505 EMULONG lt, lta, ltb;
2508 /* NaN times anything is the same NaN. */
2519 /* Zero times infinity is a NaN. */
2520 if ((eisinf (a) && (ecmp (b, ezero) == 0))
2521 || (eisinf (b) && (ecmp (a, ezero) == 0)))
2523 mtherr ("emul", INVALID);
2528 /* Infinity times anything else is infinity. */
2530 if (eisinf (a) || eisinf (b))
2532 if (eisneg (a) ^ eisneg (b))
2533 *(c + (NE - 1)) = 0x8000;
2535 *(c + (NE - 1)) = 0;
2546 for (i = 1; i < NI - 1; i++)
2550 lta -= enormlz (ai);
2561 for (i = 1; i < NI - 1; i++)
2565 ltb -= enormlz (bi);
2574 /* Multiply significands */
2576 /* calculate exponent */
2577 lt = lta + ltb - (EXONE - 1);
2578 emdnorm (bi, j, 0, lt, 64);
2579 /* calculate sign of product */
2591 ; Convert IEEE double precision to e type
2593 ; unsigned EMUSHORT x[N+2];
2598 unsigned EMUSHORT *pe, *y;
2602 dectoe (pe, y); /* see etodec.c */
2606 register unsigned EMUSHORT r;
2607 register unsigned EMUSHORT *e, *p;
2608 unsigned EMUSHORT yy[NI];
2612 denorm = 0; /* flag if denormalized number */
2621 yy[M] = (r & 0x0f) | 0x10;
2622 r &= ~0x800f; /* strip sign and 4 significand bits */
2628 if (((pe[3] & 0xf) != 0) || (pe[2] != 0)
2629 || (pe[1] != 0) || (pe[0] != 0))
2635 if (((pe[0] & 0xf) != 0) || (pe[1] != 0)
2636 || (pe[2] != 0) || (pe[3] != 0))
2649 #endif /* INFINITY */
2651 /* If zero exponent, then the significand is denormalized.
2652 * So, take back the understood high significand bit. */
2674 { /* if zero exponent, then normalize the significand */
2675 if ((k = enormlz (yy)) > NBITS)
2678 yy[E] -= (unsigned EMUSHORT) (k - 1);
2681 #endif /* not DEC */
2686 unsigned EMUSHORT *pe, *y;
2688 unsigned EMUSHORT yy[NI];
2689 unsigned EMUSHORT *e, *p, *q;
2694 for (i = 0; i < NE - 5; i++)
2697 for (i = 0; i < 5; i++)
2701 for (i = 0; i < 5; i++)
2705 p = &yy[0] + (NE - 1);
2708 for (i = 0; i < 4; i++)
2718 for (i = 0; i < 4; i++)
2727 for (i = 1; i <= 4; i++)
2743 #endif /* INFINITY */
2744 for (i = 0; i < NE; i++)
2750 ; Convert IEEE single precision to e type
2752 ; unsigned EMUSHORT x[N+2];
2757 unsigned EMUSHORT *pe, *y;
2759 register unsigned EMUSHORT r;
2760 register unsigned EMUSHORT *e, *p;
2761 unsigned EMUSHORT yy[NI];
2765 denorm = 0; /* flag if denormalized number */
2777 yy[M] = (r & 0x7f) | 0200;
2778 r &= ~0x807f; /* strip sign and 7 significand bits */
2784 if (((pe[0] & 0x7f) != 0) || (pe[1] != 0))
2790 if (((pe[1] & 0x7f) != 0) || (pe[0] != 0))
2803 #endif /* INFINITY */
2805 /* If zero exponent, then the significand is denormalized.
2806 * So, take back the understood high significand bit. */
2827 { /* if zero exponent, then normalize the significand */
2828 if ((k = enormlz (yy)) > NBITS)
2831 yy[E] -= (unsigned EMUSHORT) (k - 1);
2839 unsigned EMUSHORT *x, *e;
2841 unsigned EMUSHORT xi[NI];
2848 make_nan (e, XFmode);
2853 /* adjust exponent for offset */
2854 exp = (EMULONG) xi[E];
2859 /* round off to nearest or even */
2862 emdnorm (xi, 0, 0, exp, 64);
2868 /* move out internal format to ieee long double */
2871 unsigned EMUSHORT *a, *b;
2873 register unsigned EMUSHORT *p, *q;
2874 unsigned EMUSHORT i;
2879 make_nan (b, XFmode);
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 */
2894 /* combine sign and exponent */
2898 *q++ = *p++ | 0x8000;
2904 *q-- = *p++ | 0x8000;
2908 /* skip over guard word */
2910 /* move the significand */
2912 for (i = 0; i < 4; i++)
2915 for (i = 0; i < 4; i++)
2922 ; e type to IEEE double precision
2924 ; unsigned EMUSHORT x[NE];
2932 unsigned EMUSHORT *x, *e;
2934 etodec (x, e); /* see etodec.c */
2939 unsigned EMUSHORT *x, *y;
2948 unsigned EMUSHORT *x, *e;
2950 unsigned EMUSHORT xi[NI];
2957 make_nan (e, DFmode);
2962 /* adjust exponent for offsets */
2963 exp = (EMULONG) xi[E] - (EXONE - 0x3ff);
2968 /* round off to nearest or even */
2971 emdnorm (xi, 0, 0, exp, 64);
2980 unsigned EMUSHORT *x, *y;
2982 unsigned EMUSHORT i;
2983 unsigned EMUSHORT *p;
2988 make_nan (y, DFmode);
2996 *y = 0; /* output high order */
2998 *y = 0x8000; /* output sign bit */
3001 if (i >= (unsigned int) 2047)
3002 { /* Saturate at largest number less than infinity. */
3017 *y |= (unsigned EMUSHORT) 0x7fef;
3041 i |= *p++ & (unsigned EMUSHORT) 0x0f; /* *p = xi[M] */
3042 *y |= (unsigned EMUSHORT) i; /* high order output already has sign bit set */
3056 #endif /* not DEC */
3061 ; e type to IEEE single precision
3063 ; unsigned EMUSHORT x[N+2];
3068 unsigned EMUSHORT *x, *e;
3071 unsigned EMUSHORT xi[NI];
3077 make_nan (e, SFmode);
3082 /* adjust exponent for offsets */
3083 exp = (EMULONG) xi[E] - (EXONE - 0177);
3088 /* round off to nearest or even */
3091 emdnorm (xi, 0, 0, exp, 64);
3099 unsigned EMUSHORT *x, *y;
3101 unsigned EMUSHORT i;
3102 unsigned EMUSHORT *p;
3107 make_nan (y, SFmode);
3118 *y = 0; /* output high order */
3120 *y = 0x8000; /* output sign bit */
3123 /* Handle overflow cases. */
3127 *y |= (unsigned EMUSHORT) 0x7f80;
3138 #else /* no INFINITY */
3139 *y |= (unsigned EMUSHORT) 0x7f7f;
3153 #endif /* no INFINITY */
3165 i |= *p++ & (unsigned EMUSHORT) 0x7f; /* *p = xi[M] */
3166 *y |= i; /* high order output already has sign bit set */
3180 /* Compare two e type numbers.
3182 * unsigned EMUSHORT a[NE], b[NE];
3185 * returns +1 if a > b
3188 * -2 if either a or b is a NaN.
3192 unsigned EMUSHORT *a, *b;
3194 unsigned EMUSHORT ai[NI], bi[NI];
3195 register unsigned EMUSHORT *p, *q;
3200 if (eisnan (a) || eisnan (b))
3209 { /* the signs are different */
3211 for (i = 1; i < NI - 1; i++)
3225 /* both are the same sign */
3240 return (0); /* equality */
3246 if (*(--p) > *(--q))
3247 return (msign); /* p is bigger */
3249 return (-msign); /* p is littler */
3255 /* Find nearest integer to x = floor (x + 0.5)
3257 * unsigned EMUSHORT x[NE], y[NE]
3262 unsigned EMUSHORT *x, *y;
3272 ; convert long integer to e type
3275 ; unsigned EMUSHORT x[NE];
3277 ; note &l is the memory address of l
3281 long *lp; /* lp is the memory address of a long integer */
3282 unsigned EMUSHORT *y; /* y is the address of a short */
3284 unsigned EMUSHORT yi[NI];
3291 /* make it positive */
3292 ll = (unsigned long) (-(*lp));
3293 yi[0] = 0xffff; /* put correct sign in the e type number */
3297 ll = (unsigned long) (*lp);
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 */
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 */
3312 if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
3313 ecleaz (yi); /* it was zero */
3315 yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
3316 emovo (yi, y); /* output the answer */
3320 ; convert unsigned long integer to e type
3323 ; unsigned EMUSHORT x[NE];
3325 ; note &l is the memory address of l
3329 unsigned long *lp; /* lp is the memory address of a long integer */
3330 unsigned EMUSHORT *y; /* y is the address of a short */
3332 unsigned EMUSHORT yi[NI];
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 */
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 */
3352 if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
3353 ecleaz (yi); /* it was zero */
3355 yi[E] -= (unsigned EMUSHORT) k; /* subtract shift count from exponent */
3356 emovo (yi, y); /* output the answer */
3361 ; Find long integer and fractional parts
3364 ; unsigned EMUSHORT x[NE], frac[NE];
3365 ; xifrac (x, &i, frac);
3367 The integer output has the sign of the input. The fraction is
3368 the positive fractional part of abs (x).
3372 unsigned EMUSHORT *x;
3374 unsigned EMUSHORT *frac;
3376 unsigned EMUSHORT xi[NI];
3381 k = (int) xi[E] - (EXONE - 1);
3384 /* if exponent <= 0, integer = 0 and real output is fraction */
3389 if (k > (HOST_BITS_PER_LONG - 1))
3391 /* long integer overflow: output large integer
3392 and correct fraction */
3394 *i = ((unsigned long) 1) << (HOST_BITS_PER_LONG - 1);
3396 *i = (((unsigned long) 1) << (HOST_BITS_PER_LONG - 1)) - 1;
3399 warning ("overflow on truncation to integer");
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);
3412 ll = (ll << 16) | xi[M];
3414 while ((k -= 16) > 0);
3421 /* shift not more than 16 bits */
3423 *i = (long) xi[M] & 0xffff;
3430 if ((k = enormlz (xi)) > NBITS)
3433 xi[E] -= (unsigned EMUSHORT) k;
3439 /* Find unsigned long integer and fractional parts.
3440 A negative e type input yields integer output = 0
3441 but correct fraction. */
3444 euifrac (x, i, frac)
3445 unsigned EMUSHORT *x;
3447 unsigned EMUSHORT *frac;
3450 unsigned EMUSHORT xi[NI];
3454 k = (int) xi[E] - (EXONE - 1);
3457 /* if exponent <= 0, integer = 0 and argument is fraction */
3462 if (k > HOST_BITS_PER_LONG)
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. */
3471 warning ("overflow on truncation to unsigned integer");
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);
3484 ll = (ll << 16) | xi[M];
3486 while ((k -= 16) > 0);
3491 /* shift not more than 16 bits */
3493 *i = (long) xi[M] & 0xffff;
3496 if (xi[0]) /* A negative value yields unsigned integer 0. */
3501 if ((k = enormlz (xi)) > NBITS)
3504 xi[E] -= (unsigned EMUSHORT) k;
3514 ; Shifts significand area up or down by the number of bits
3515 ; given by the variable sc.
3519 unsigned EMUSHORT *x;
3522 unsigned EMUSHORT lost;
3523 unsigned EMUSHORT *p;
3536 lost |= *p; /* remember lost bits */
3577 return ((int) lost);
3585 ; Shift normalizes the significand area pointed to by argument
3586 ; shift count (up = positive) is returned.
3590 unsigned EMUSHORT x[];
3592 register unsigned EMUSHORT *p;
3601 return (0); /* already normalized */
3606 /* With guard word, there are NBITS+16 bits available.
3607 * return true if all are zero.
3612 /* see if high byte is zero */
3613 while ((*p & 0xff00) == 0)
3618 /* now shift 1 bit at a time */
3619 while ((*p & 0x8000) == 0)
3625 mtherr ("enormlz", UNDERFLOW);
3631 /* Normalize by shifting down out of the high guard word
3632 of the significand */
3647 mtherr ("enormlz", OVERFLOW);
3657 /* Convert e type number to decimal format ASCII string.
3658 * The constants are for 64 bit precision.
3664 static unsigned EMUSHORT etens[NTEN + 1][NE] =
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 */
3681 static unsigned EMUSHORT emtens[NTEN + 1][NE] =
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 */
3699 e24toasc (x, string, ndigs)
3700 unsigned EMUSHORT x[];
3704 unsigned EMUSHORT w[NI];
3707 etoasc (w, string, ndigs);
3712 e53toasc (x, string, ndigs)
3713 unsigned EMUSHORT x[];
3717 unsigned EMUSHORT w[NI];
3720 etoasc (w, string, ndigs);
3725 e64toasc (x, string, ndigs)
3726 unsigned EMUSHORT x[];
3730 unsigned EMUSHORT w[NI];
3733 etoasc (w, string, ndigs);
3737 static char wstring[80]; /* working storage for ASCII output */
3740 etoasc (x, string, ndigs)
3741 unsigned EMUSHORT x[];
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;
3751 unsigned EMUSHORT m;
3762 sprintf (wstring, " NaN ");
3766 rndprc = NBITS; /* set to full precision */
3767 emov (x, y); /* retain external format */
3768 if (y[NE - 1] & 0x8000)
3771 y[NE - 1] &= 0x7fff;
3778 ten = &etens[NTEN][0];
3780 /* Test for zero exponent */
3783 for (k = 0; k < NE - 1; k++)
3786 goto tnzro; /* denormalized number */
3788 goto isone; /* legal all zeros */
3792 /* Test for infinity. */
3793 if (y[NE - 1] == 0x7fff)
3796 sprintf (wstring, " -Infinity ");
3798 sprintf (wstring, " Infinity ");
3802 /* Test for exponent nonzero but significand denormalized.
3803 * This is an error condition.
3805 if ((y[NE - 1] != 0) && ((y[NE - 2] & 0x8000) == 0))
3807 mtherr ("etoasc", DOMAIN);
3808 sprintf (wstring, "NaN");
3812 /* Compare to 1.0 */
3821 { /* Number is greater than 1 */
3822 /* Convert significand to an integer and strip trailing decimal zeros. */
3824 u[NE - 1] = EXONE + NBITS - 1;
3826 p = &etens[NTEN - 4][0];
3832 for (j = 0; j < NE - 1; j++)
3845 /* Rescale from integer significand */
3846 u[NE - 1] += y[NE - 1] - (unsigned int) (EXONE + NBITS - 1);
3848 /* Find power of 10 */
3852 /* An unordered compare result shouldn't happen here. */
3853 while (ecmp (ten, u) <= 0)
3855 if (ecmp (p, u) <= 0)
3868 { /* Number is less than 1.0 */
3869 /* Pad significand with trailing decimal zeros. */
3872 while ((y[NE - 2] & 0x8000) == 0)
3881 for (i = 0; i < NDEC + 1; i++)
3883 if ((w[NI - 1] & 0x7) != 0)
3885 /* multiply by 10 */
3898 if (eone[NE - 1] <= u[1])
3910 while (ecmp (eone, w) > 0)
3912 if (ecmp (p, w) >= 0)
3927 /* Find the first (leading) digit. */
3933 digit = equot[NI - 1];
3934 while ((digit == 0) && (ecmp (y, ezero) != 0))
3942 digit = equot[NI - 1];
3950 /* Examine number of digits requested by caller. */
3968 *s++ = (char )digit + '0';
3971 /* Generate digits after the decimal point. */
3972 for (k = 0; k <= ndigs; k++)
3974 /* multiply current number by 10, without normalizing */
3981 *s++ = (char) equot[NI - 1] + '0';
3983 digit = equot[NI - 1];
3986 /* round off the ASCII string */
3989 /* Test for critical rounding case in ASCII output. */
3993 if (ecmp (t, ezero) != 0)
3994 goto roun; /* round to nearest */
3995 if ((*(s - 1) & 1) == 0)
3996 goto doexp; /* round to even */
3998 /* Round up and propagate carry-outs */
4002 /* Carry out to most significant digit? */
4009 /* Most significant digit carries to 10? */
4017 /* Round up and carry out from less significant digits */
4029 sprintf (ss, "e+%d", expon);
4031 sprintf (ss, "e%d", expon);
4033 sprintf (ss, "e%d", expon);
4036 /* copy out the working string */
4039 while (*ss == ' ') /* strip possible leading space */
4041 while ((*s++ = *ss++) != '\0')
4050 ; ASCTOQ.MAC LATEST REV: 11 JAN 84
4053 ; Convert ASCII string to quadruple precision floating point
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).
4064 ; asctoq (string, q);
4067 /* ASCII to single */
4071 unsigned EMUSHORT *y;
4077 /* ASCII to double */
4081 unsigned EMUSHORT *y;
4091 /* ASCII to long double */
4095 unsigned EMUSHORT *y;
4100 /* ASCII to super double */
4104 unsigned EMUSHORT *y;
4106 asctoeg (s, y, NBITS);
4109 /* Space to make a copy of the input string: */
4110 static char lstr[82];
4113 asctoeg (ss, y, oprec)
4115 unsigned EMUSHORT *y;
4118 unsigned EMUSHORT yy[NI], xt[NI], tt[NI];
4119 int esign, decflg, sgnflg, nexp, exp, prec, lost;
4120 int k, trail, c, rndsav;
4122 unsigned EMUSHORT nsign, *p;
4125 /* Copy the input string. */
4127 while (*s == ' ') /* skip leading spaces */
4130 for (k = 0; k < 79; k++)
4132 if ((*sp++ = *s++) == '\0')
4139 rndprc = NBITS; /* Set to full precision */
4152 if ((k >= 0) && (k <= 9))
4154 /* Ignore leading zeros */
4155 if ((prec == 0) && (decflg == 0) && (k == 0))
4157 /* Identify and strip trailing zeros after the decimal point. */
4158 if ((trail == 0) && (decflg != 0))
4161 while ((*sp >= '0') && (*sp <= '9'))
4163 /* Check for syntax error */
4165 if ((c != 'e') && (c != 'E') && (c != '\0')
4166 && (c != '\n') && (c != '\r') && (c != ' ')
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.
4184 nexp += 1; /* count digits after decimal point */
4185 eshup1 (yy); /* multiply current number by 10 */
4191 xt[NI - 2] = (unsigned EMUSHORT) k;
4209 case '.': /* decimal point */
4239 mtherr ("asctoe", DOMAIN);
4248 /* Exponent interpretation */
4254 /* check for + or - */
4262 while ((*s >= '0') && (*s <= '9'))
4280 yy[E] = 0x7fff; /* infinity */
4292 /* Pad trailing zeros to minimize power of 10, per IEEE spec. */
4293 while ((nexp > 0) && (yy[2] == 0))
4305 if ((k = enormlz (yy)) > NBITS)
4310 lexp = (EXONE - 1 + NBITS) - k;
4311 emdnorm (yy, lost, 0, lexp, 64);
4312 /* convert to external format */
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.
4333 { /* Punt. Can't handle this without 2 divides. */
4334 emovi (etens[0], tt);
4341 p = &etens[NTEN][0];
4351 while (exp <= MAXP);
4369 /* Round and convert directly to the destination type */
4371 lexp -= EXONE - 0x3ff;
4372 else if (oprec == 24)
4373 lexp -= EXONE - 0177;
4375 else if (oprec == 56)
4376 lexp -= EXONE - 0201;
4379 emdnorm (yy, k, 0, lexp, 64);
4389 todec (yy, y); /* see etodec.c */
4409 /* y = largest integer not greater than x
4410 * (truncated toward minus infinity)
4412 * unsigned EMUSHORT x[NE], y[NE]
4416 static unsigned EMUSHORT bmask[] =
4439 unsigned EMUSHORT x[], y[];
4441 register unsigned EMUSHORT *p;
4443 unsigned EMUSHORT f[NE];
4445 emov (x, f); /* leave in external format */
4446 expon = (int) f[NE - 1];
4447 e = (expon & 0x7fff) - (EXONE - 1);
4453 /* number of bits to clear out */
4465 /* clear the remaining bits */
4467 /* truncate negatives toward minus infinity */
4470 if ((unsigned EMUSHORT) expon & (unsigned EMUSHORT) 0x8000)
4472 for (i = 0; i < NE - 1; i++)
4484 /* unsigned EMUSHORT x[], s[];
4487 * efrexp (x, exp, s);
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.
4495 unsigned EMUSHORT x[];
4497 unsigned EMUSHORT s[];
4499 unsigned EMUSHORT xi[NI];
4503 li = (EMULONG) ((EMUSHORT) xi[1]);
4511 *exp = (int) (li - 0x3ffe);
4516 /* unsigned EMUSHORT x[], y[];
4519 * eldexp (x, pwr2, y);
4521 * Returns y = x * 2**pwr2.
4525 unsigned EMUSHORT x[];
4527 unsigned EMUSHORT y[];
4529 unsigned EMUSHORT xi[NI];
4537 emdnorm (xi, i, i, li, 64);
4542 /* c = remainder after dividing b by a
4543 * Least significant integer quotient bits left in equot[].
4547 unsigned EMUSHORT a[], b[], c[];
4549 unsigned EMUSHORT den[NI], num[NI];
4553 || (ecmp (a, ezero) == 0)
4561 if (ecmp (a, ezero) == 0)
4563 mtherr ("eremain", SING);
4569 eiremain (den, num);
4570 /* Sign of remainder = sign of quotient */
4580 unsigned EMUSHORT den[], num[];
4583 unsigned EMUSHORT j;
4586 ld -= enormlz (den);
4588 ln -= enormlz (num);
4592 if (ecmpm (den, num) <= 0)
4606 emdnorm (num, 0, 0, ln, 0);
4611 * Library common error handling routine
4621 * mtherr (fctnam, code);
4627 * This routine may be called to report one of the following
4628 * error conditions (in the include file mconf.h).
4630 * Mnemonic Value Significance
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
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
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.
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
4666 /* include "mconf.h" */
4668 /* Notice: the order of appearance of the following
4669 * messages is bound to the error codes defined
4673 static char *ermsg[NMSGS] =
4675 "unknown", /* error code 0 */
4676 "domain", /* error code 1 */
4677 "singularity", /* et seq. */
4680 "total loss of precision",
4681 "partial loss of precision",
4695 /* Display string passed by calling program,
4696 * which is supposed to be the name of the
4697 * function in which the error occurred.
4700 /* Display error message defined
4701 * by the code argument.
4703 if ((code <= 0) || (code >= NMSGS))
4705 sprintf (errstr, " %s %s error", name, ermsg[code]);
4708 /* Set global error message word */
4711 /* Return to calling
4716 /* Here is etodec.c .
4721 ; convert DEC double precision to e type
4728 unsigned EMUSHORT *d;
4729 unsigned EMUSHORT *e;
4731 unsigned EMUSHORT y[NI];
4732 register unsigned EMUSHORT r, *p;
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 */
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 */
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 */
4755 *p++ = *d++; /* fill in the rest of our mantissa */
4759 eshdn8 (y); /* shift our mantissa down 8 bits */
4767 ; convert e type to DEC double precision
4773 static unsigned EMUSHORT decbit[NI] = {0, 0, 0, 0, 0, 0, 0200, 0};
4777 unsigned EMUSHORT *x, *d;
4779 unsigned EMUSHORT xi[NI];
4780 register unsigned EMUSHORT r;
4788 if (r < (EXONE - 128))
4791 if ((i & 0200) != 0)
4793 if ((i & 0377) == 0200)
4795 if ((i & 0400) != 0)
4797 /* check all less significant bits */
4798 for (j = M + 5; j < NI; j++)
4847 unsigned EMUSHORT *x, *d;
4849 unsigned EMUSHORT xi[NI];
4854 exp = (EMULONG) xi[E] - (EXONE - 0201); /* adjust exponent for offsets */
4855 /* round off to nearest or even */
4858 emdnorm (xi, 0, 0, exp, 64);
4865 unsigned EMUSHORT *x, *y;
4867 unsigned EMUSHORT i;
4868 unsigned EMUSHORT *p;
4908 /* Output a binary NaN bit pattern in the target machine's format. */
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
4917 unsigned EMUSHORT TFnan[8] =
4918 {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
4921 unsigned EMUSHORT TFnan[8] = {0, 0, 0, 0, 0, 0, 0x8000, 0xffff};
4929 unsigned EMUSHORT XFnan[6] = {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
4932 unsigned EMUSHORT XFnan[6] = {0, 0, 0, 0xc000, 0xffff, 0};
4940 unsigned EMUSHORT DFnan[4] = {0x7fff, 0xffff, 0xffff, 0xffff};
4943 unsigned EMUSHORT DFnan[4] = {0, 0, 0, 0xfff8};
4951 unsigned EMUSHORT SFnan[2] = {0x7fff, 0xffff};
4954 unsigned EMUSHORT SFnan[2] = {0, 0xffc0};
4960 make_nan (nan, mode)
4961 unsigned EMUSHORT *nan;
4962 enum machine_mode mode;
4965 unsigned EMUSHORT *p;
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. */
4992 for (i=0; i < n; i++)
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. */
5001 ereal_from_float (f)
5005 unsigned EMUSHORT s[2];
5006 unsigned EMUSHORT e[NE];
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;
5014 s[0] = (unsigned EMUSHORT) f;
5015 s[1] = (unsigned EMUSHORT) (f >> 16);
5017 /* Convert and promote the target float to E-type. */
5019 /* Output E-type to REAL_VALUE_TYPE. */
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.
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. */
5034 ereal_from_double (d)
5038 unsigned EMUSHORT s[4];
5039 unsigned EMUSHORT e[NE];
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];
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);
5054 /* Convert target double to E-type. */
5056 /* Output E-type to REAL_VALUE_TYPE. */
5060 #endif /* EMU_NON_COMPILE not defined */