]>
gcc.gnu.org Git - gcc.git/blob - gcc/real.c
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. */
28 /* To enable support of XFmode extended real floating point, define
29 LONG_DOUBLE_TYPE_SIZE 96 in the tm.h file (m68k.h or i386.h).
31 To support cross compilation between IEEE and VAX floating
32 point formats, define REAL_ARITHMETIC in the tm.h file.
34 In either case the machine files (tm.h) must not contain any code
35 that tries to use host floating point arithmetic to convert
36 REAL_VALUE_TYPEs from `double' to `float', pass them to fprintf,
37 etc. In cross-compile situations a REAL_VALUE_TYPE may not
38 be intelligible to the host computer's native arithmetic.
40 The emulator defaults to the host's floating point format so that
41 its decimal conversion functions can be used if desired (see
44 The first part of this file interfaces gcc to ieee.c, which is a
45 floating point arithmetic suite that was not written with gcc in
46 mind. The interface is followed by ieee.c itself and related
47 items. Avoid changing ieee.c unless you have suitable test
48 programs available. A special version of the PARANOIA floating
49 point arithmetic tester, modified for this purpose, can be found
50 on usc.edu : /pub/C-numanal/ieeetest.zoo. Some tutorial
51 information on ieee.c is given in my book: S. L. Moshier,
52 _Methods and Programs for Mathematical Functions_, Prentice-Hall
53 or Simon & Schuster Int'l, 1989. A library of XFmode elementary
54 transcendental functions can be obtained by ftp from
55 research.att.com: netlib/cephes/ldouble.shar.Z */
57 /* Type of computer arithmetic.
58 * Only one of DEC, MIEEE, IBMPC, or UNK should get defined.
59 * The following modification converts gcc macros into the ones
62 * Note: long double formats differ between IBMPC and MIEEE
63 * by more than just endian-ness.
66 /* REAL_ARITHMETIC defined means that macros in real.h are
67 defined to call emulator functions. */
68 #ifdef REAL_ARITHMETIC
70 #if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
71 /* PDP-11, Pro350, VAX: */
73 #else /* it's not VAX */
74 #if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
76 /* Motorola IEEE, high order words come first (Sun workstation): */
78 #else /* not big-endian */
79 /* Intel IEEE, low order words come first:
82 #endif /* big-endian */
83 #else /* it's not IEEE either */
84 /* UNKnown arithmetic. We don't support this and can't go on. */
85 unknown arithmetic type
91 /* REAL_ARITHMETIC not defined means that the *host's* data
92 structure will be used. It may differ by endian-ness from the
93 target machine's structure and will get its ends swapped
94 accordingly (but not here). Probably only the decimal <-> binary
95 functions in this file will actually be used in this case. */
96 #if HOST_FLOAT_FORMAT == VAX_FLOAT_FORMAT
98 #else /* it's not VAX */
99 #if HOST_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
100 #ifdef HOST_WORDS_BIG_ENDIAN
102 #else /* not big-endian */
104 #endif /* big-endian */
105 #else /* it's not IEEE either */
106 unknown arithmetic type
108 #endif /* not IEEE */
111 #endif /* REAL_ARITHMETIC not defined */
113 /* Define for support of infinity. */
121 * Include file for extended precision arithmetic programs.
124 /* Number of 16 bit words in external e type format */
127 /* Number of 16 bit words in internal format */
130 /* Array offset to exponent */
133 /* Array offset to high guard word */
136 /* Number of bits of precision */
137 #define NBITS ((NI-4)*16)
139 /* Maximum number of decimal digits in ASCII conversion
142 #define NDEC (NBITS*8/27)
144 /* The exponent of 1.0 */
145 #define EXONE (0x3fff)
147 /* Find a host integer type that is at least 16 bits wide,
148 and another type at least twice whatever that size is. */
150 #if HOST_BITS_PER_CHAR >= 16
151 #define EMUSHORT char
152 #define EMUSHORT_SIZE HOST_BITS_PER_CHAR
153 #define EMULONG_SIZE (2 * HOST_BITS_PER_CHAR)
155 #if HOST_BITS_PER_SHORT >= 16
156 #define EMUSHORT short
157 #define EMUSHORT_SIZE HOST_BITS_PER_SHORT
158 #define EMULONG_SIZE (2 * HOST_BITS_PER_SHORT)
160 #if HOST_BITS_PER_INT >= 16
162 #define EMUSHORT_SIZE HOST_BITS_PER_INT
163 #define EMULONG_SIZE (2 * HOST_BITS_PER_INT)
165 #if HOST_BITS_PER_LONG >= 16
166 #define EMUSHORT long
167 #define EMUSHORT_SIZE HOST_BITS_PER_LONG
168 #define EMULONG_SIZE (2 * HOST_BITS_PER_LONG)
170 /* You will have to modify this program to have a smaller unit size. */
171 #define EMU_NON_COMPILE
177 #if HOST_BITS_PER_SHORT >= EMULONG_SIZE
178 #define EMULONG short
180 #if HOST_BITS_PER_INT >= EMULONG_SIZE
183 #if HOST_BITS_PER_LONG >= EMULONG_SIZE
186 #if HOST_BITS_PER_LONG_LONG >= EMULONG_SIZE
187 #define EMULONG long long int
189 /* You will have to modify this program to have a smaller unit size. */
190 #define EMU_NON_COMPILE
197 /* The host interface doesn't work if no 16-bit size exists. */
198 #if EMUSHORT_SIZE != 16
199 #define EMU_NON_COMPILE
202 /* OK to continue compilation. */
203 #ifndef EMU_NON_COMPILE
205 /* Construct macros to translate between REAL_VALUE_TYPE and e type.
206 In GET_REAL and PUT_REAL, r and e are pointers.
207 A REAL_VALUE_TYPE is guaranteed to occupy contiguous locations
208 in memory, with no holes. */
210 #if LONG_DOUBLE_TYPE_SIZE == 96
211 #define GET_REAL(r,e) bcopy (r, e, 2*NE)
212 #define PUT_REAL(e,r) bcopy (e, r, 2*NE)
213 #else /* no XFmode */
215 #ifdef REAL_ARITHMETIC
216 /* Emulator uses target format internally
217 but host stores it in host endian-ness. */
219 #if defined (HOST_WORDS_BIG_ENDIAN) == WORDS_BIG_ENDIAN
220 #define GET_REAL(r,e) e53toe ((r), (e))
221 #define PUT_REAL(e,r) etoe53 ((e), (r))
223 #else /* endian-ness differs */
224 /* emulator uses target endian-ness internally */
225 #define GET_REAL(r,e) \
226 do { EMUSHORT w[4]; \
227 w[3] = ((EMUSHORT *) r)[0]; \
228 w[2] = ((EMUSHORT *) r)[1]; \
229 w[1] = ((EMUSHORT *) r)[2]; \
230 w[0] = ((EMUSHORT *) r)[3]; \
231 e53toe (w, (e)); } while (0)
233 #define PUT_REAL(e,r) \
234 do { EMUSHORT w[4]; \
236 *((EMUSHORT *) r) = w[3]; \
237 *((EMUSHORT *) r + 1) = w[2]; \
238 *((EMUSHORT *) r + 2) = w[1]; \
239 *((EMUSHORT *) r + 3) = w[0]; } while (0)
241 #endif /* endian-ness differs */
243 #else /* not REAL_ARITHMETIC */
245 /* emulator uses host format */
246 #define GET_REAL(r,e) e53toe ((r), (e))
247 #define PUT_REAL(e,r) etoe53 ((e), (r))
249 #endif /* not REAL_ARITHMETIC */
250 #endif /* no XFmode */
253 int ecmp (), enormlz (), eshift (), eisneg (), eisinf ();
254 void eadd (), esub (), emul (), ediv ();
255 void eshup1 (), eshup8 (), eshup6 (), eshdn1 (), eshdn8 (), eshdn6 ();
256 void eabs (), eneg (), emov (), eclear (), einfin (), efloor ();
257 void eldexp (), efrexp (), eifrac (), euifrac (), ltoe (), ultoe ();
258 void eround (), ereal_to_decimal ();
259 void esqrt (), elog (), eexp (), etanh (), epow ();
260 void asctoe (), asctoe24 (), asctoe53 (), asctoe64 ();
261 void etoasc (), e24toasc (), e53toasc (), e64toasc ();
262 void etoe64 (), etoe53 (), etoe24 (), e64toe (), e53toe (), e24toe ();
264 extern unsigned EMUSHORT ezero
[], ehalf
[], eone
[], etwo
[];
265 extern unsigned EMUSHORT elog2
[], esqrt2
[];
267 /* Pack output array with 32-bit numbers obtained from
268 array containing 16-bit numbers, swapping ends if required. */
271 unsigned EMUSHORT e
[];
273 enum machine_mode mode
;
283 /* Swap halfwords in the third long. */
284 th
= (unsigned long) e
[4] & 0xffff;
285 t
= (unsigned long) e
[5] & 0xffff;
288 /* fall into the double case */
292 /* swap halfwords in the second word */
293 th
= (unsigned long) e
[2] & 0xffff;
294 t
= (unsigned long) e
[3] & 0xffff;
297 /* fall into the float case */
301 /* swap halfwords in the first word */
302 th
= (unsigned long) e
[0] & 0xffff;
303 t
= (unsigned long) e
[1] & 0xffff;
314 /* Pack the output array without swapping. */
321 /* Pack the third long.
322 Each element of the input REAL_VALUE_TYPE array has 16 bit useful bits
324 th
= (unsigned long) e
[5] & 0xffff;
325 t
= (unsigned long) e
[4] & 0xffff;
328 /* fall into the double case */
332 /* pack the second long */
333 th
= (unsigned long) e
[3] & 0xffff;
334 t
= (unsigned long) e
[2] & 0xffff;
337 /* fall into the float case */
341 /* pack the first long */
342 th
= (unsigned long) e
[1] & 0xffff;
343 t
= (unsigned long) e
[0] & 0xffff;
356 /* This is the implementation of the REAL_ARITHMETIC macro.
359 earith (value
, icode
, r1
, r2
)
360 REAL_VALUE_TYPE
*value
;
365 unsigned EMUSHORT d1
[NE
], d2
[NE
], v
[NE
];
370 code
= (enum tree_code
) icode
;
378 esub (d2
, d1
, v
); /* d1 - d2 */
386 #ifndef REAL_INFINITY
387 if (ecmp (d2
, ezero
) == 0)
390 ediv (d2
, d1
, v
); /* d1/d2 */
393 case MIN_EXPR
: /* min (d1,d2) */
394 if (ecmp (d1
, d2
) < 0)
400 case MAX_EXPR
: /* max (d1,d2) */
401 if (ecmp (d1
, d2
) > 0)
414 /* Truncate REAL_VALUE_TYPE toward zero to signed HOST_WIDE_INT
415 * implements REAL_VALUE_FIX_TRUNCATE (x) (etrunci (x))
421 unsigned EMUSHORT f
[NE
], g
[NE
];
433 /* Truncate REAL_VALUE_TYPE toward zero to unsigned HOST_WIDE_INT
434 * implements REAL_VALUE_UNSIGNED_FIX_TRUNCATE (x) (etruncui (x))
440 unsigned EMUSHORT f
[NE
], g
[NE
];
452 /* This is the REAL_VALUE_ATOF function.
453 * It converts a decimal string to binary, rounding off
454 * as indicated by the machine_mode argument. Then it
455 * promotes the rounded value to REAL_VALUE_TYPE.
462 unsigned EMUSHORT tem
[NE
], e
[NE
];
487 /* Expansion of REAL_NEGATE.
493 unsigned EMUSHORT e
[NE
];
504 * implements REAL_VALUE_FIX (x) (eroundi (x))
505 * The type of rounding is left unspecified by real.h.
506 * It is implemented here as round to nearest (add .5 and chop).
512 unsigned EMUSHORT f
[NE
], g
[NE
];
521 /* Round real to nearest unsigned int
522 * implements REAL_VALUE_UNSIGNED_FIX (x) ((unsigned int) eroundi (x))
523 * Negative input returns zero.
524 * The type of rounding is left unspecified by real.h.
525 * It is implemented here as round to nearest (add .5 and chop).
531 unsigned EMUSHORT f
[NE
], g
[NE
];
537 return ((unsigned int)l
);
541 /* REAL_VALUE_FROM_INT macro.
544 ereal_from_int (d
, i
, j
)
548 unsigned EMUSHORT df
[NE
], dg
[NE
];
557 /* complement and add 1 */
564 eldexp (eone
, HOST_BITS_PER_LONG
, df
);
575 /* REAL_VALUE_FROM_UNSIGNED_INT macro.
578 ereal_from_uint (d
, i
, j
)
582 unsigned EMUSHORT df
[NE
], dg
[NE
];
583 unsigned long low
, high
;
587 eldexp (eone
, HOST_BITS_PER_LONG
, df
);
596 /* REAL_VALUE_TO_INT macro
599 ereal_to_int (low
, high
, rr
)
603 unsigned EMUSHORT d
[NE
], df
[NE
], dg
[NE
], dh
[NE
];
607 /* convert positive value */
614 eldexp (eone
, HOST_BITS_PER_LONG
, df
);
615 ediv (df
, d
, dg
); /* dg = d / 2^32 is the high word */
616 euifrac (dg
, high
, dh
);
617 emul (df
, dh
, dg
); /* fractional part is the low word */
618 euifrac (dg
, low
, dh
);
621 /* complement and add 1 */
631 /* REAL_VALUE_LDEXP macro.
638 unsigned EMUSHORT e
[NE
], y
[NE
];
647 /* These routines are conditionally compiled because functions
648 * of the same names may be defined in fold-const.c. */
649 #ifdef REAL_ARITHMETIC
651 /* Check for infinity in a REAL_VALUE_TYPE. */
656 unsigned EMUSHORT e
[NE
];
667 /* Check whether an IEEE double precision number is a NaN. */
677 /* Check for a negative IEEE double precision number.
678 * this means strictly less than zero, not -0.
685 unsigned EMUSHORT e
[NE
];
688 if (ecmp (e
, ezero
) < 0)
693 /* Expansion of REAL_VALUE_TRUNCATE.
694 * The result is in floating point, rounded to nearest or even.
697 real_value_truncate (mode
, arg
)
698 enum machine_mode mode
;
701 unsigned EMUSHORT e
[NE
], t
[NE
];
734 #endif /* REAL_ARITHMETIC defined */
736 /* Target values are arrays of host longs. A long is guaranteed
737 to be at least 32 bits wide. */
743 unsigned EMUSHORT e
[NE
];
747 endian (e
, l
, XFmode
);
755 unsigned EMUSHORT e
[NE
];
759 endian (e
, l
, DFmode
);
766 unsigned EMUSHORT e
[NE
];
771 endian (e
, &l
, SFmode
);
776 ereal_to_decimal (x
, s
)
780 unsigned EMUSHORT e
[NE
];
788 REAL_VALUE_TYPE x
, y
;
790 unsigned EMUSHORT ex
[NE
], ey
[NE
];
794 return (ecmp (ex
, ey
));
801 unsigned EMUSHORT ex
[NE
];
804 return (eisneg (ex
));
807 /* End of REAL_ARITHMETIC interface */
811 * Extended precision IEEE binary floating point arithmetic routines
813 * Numbers are stored in C language as arrays of 16-bit unsigned
814 * short integers. The arguments of the routines are pointers to
818 * External e type data structure, simulates Intel 8087 chip
819 * temporary real format but possibly with a larger significand:
821 * NE-1 significand words (least significant word first,
822 * most significant bit is normally set)
823 * exponent (value = EXONE for 1.0,
824 * top bit is the sign)
827 * Internal data structure of a number (a "word" is 16 bits):
829 * ei[0] sign word (0 for positive, 0xffff for negative)
830 * ei[1] biased exponent (value = EXONE for the number 1.0)
831 * ei[2] high guard word (always zero after normalization)
833 * to ei[NI-2] significand (NI-4 significand words,
834 * most significant word first,
835 * most significant bit is set)
836 * ei[NI-1] low guard word (0x8000 bit is rounding place)
840 * Routines for external format numbers
842 * asctoe (string, e) ASCII string to extended double e type
843 * asctoe64 (string, &d) ASCII string to long double
844 * asctoe53 (string, &d) ASCII string to double
845 * asctoe24 (string, &f) ASCII string to single
846 * asctoeg (string, e, prec) ASCII string to specified precision
847 * e24toe (&f, e) IEEE single precision to e type
848 * e53toe (&d, e) IEEE double precision to e type
849 * e64toe (&d, e) IEEE long double precision to e type
850 * eabs (e) absolute value
851 * eadd (a, b, c) c = b + a
853 * ecmp (a, b) compare a to b, return 1, 0, or -1
854 * ediv (a, b, c) c = b / a
855 * efloor (a, b) truncate to integer, toward -infinity
856 * efrexp (a, exp, s) extract exponent and significand
857 * eifrac (e, &l, frac) e to long integer and e type fraction
858 * euifrac (e, &l, frac) e to unsigned long integer and e type fraction
859 * einfin (e) set e to infinity, leaving its sign alone
860 * eldexp (a, n, b) multiply by 2**n
862 * emul (a, b, c) c = b * a
864 * eround (a, b) b = nearest integer value to a
865 * esub (a, b, c) c = b - a
866 * e24toasc (&f, str, n) single to ASCII string, n digits after decimal
867 * e53toasc (&d, str, n) double to ASCII string, n digits after decimal
868 * e64toasc (&d, str, n) long double to ASCII string
869 * etoasc (e, str, n) e to ASCII string, n digits after decimal
870 * etoe24 (e, &f) convert e type to IEEE single precision
871 * etoe53 (e, &d) convert e type to IEEE double precision
872 * etoe64 (e, &d) convert e type to IEEE long double precision
873 * ltoe (&l, e) long (32 bit) integer to e type
874 * ultoe (&l, e) unsigned long (32 bit) integer to e type
875 * eisneg (e) 1 if sign bit of e != 0, else 0
876 * eisinf (e) 1 if e has maximum exponent
879 * Routines for internal format numbers
881 * eaddm (ai, bi) add significands, bi = bi + ai
883 * ecleazs (ei) set ei = 0 but leave its sign alone
884 * ecmpm (ai, bi) compare significands, return 1, 0, or -1
885 * edivm (ai, bi) divide significands, bi = bi / ai
886 * emdnorm (ai,l,s,exp) normalize and round off
887 * emovi (a, ai) convert external a to internal ai
888 * emovo (ai, a) convert internal ai to external a
889 * emovz (ai, bi) bi = ai, low guard word of bi = 0
890 * emulm (ai, bi) multiply significands, bi = bi * ai
891 * enormlz (ei) left-justify the significand
892 * eshdn1 (ai) shift significand and guards down 1 bit
893 * eshdn8 (ai) shift down 8 bits
894 * eshdn6 (ai) shift down 16 bits
895 * eshift (ai, n) shift ai n bits up (or down if n < 0)
896 * eshup1 (ai) shift significand and guards up 1 bit
897 * eshup8 (ai) shift up 8 bits
898 * eshup6 (ai) shift up 16 bits
899 * esubm (ai, bi) subtract significands, bi = bi - ai
902 * The result is always normalized and rounded to NI-4 word precision
903 * after each arithmetic operation.
905 * Exception flags and NaNs are NOT fully supported.
906 * This arithmetic should never produce a NaN output, but it might
907 * be confused by a NaN input.
908 * Define INFINITY in mconf.h for support of infinity; otherwise a
909 * saturation arithmetic is implemented.
910 * Denormals are always supported here where appropriate (e.g., not
911 * for conversion to DEC numbers).
918 * 5 Jan 84 PDP-11 assembly language version
919 * 2 Mar 86 fixed bug in asctoq
920 * 6 Dec 86 C language version
921 * 30 Aug 88 100 digit version, improved rounding
922 * 15 May 92 80-bit long double support
924 * Author: S. L. Moshier.
930 * Common include file for math routines
942 * This file contains definitions for error codes that are
943 * passed to the common error handling routine mtherr
946 * The file also includes a conditional assembly definition
947 * for the type of computer arithmetic (Intel IEEE, DEC, Motorola
950 * For Digital Equipment PDP-11 and VAX computers, certain
951 * IBM systems, and others that use numbers with a 56-bit
952 * significand, the symbol DEC should be defined. In this
953 * mode, most floating point constants are given as arrays
954 * of octal integers to eliminate decimal to binary conversion
955 * errors that might be introduced by the compiler.
957 * For computers, such as IBM PC, that follow the IEEE
958 * Standard for Binary Floating Point Arithmetic (ANSI/IEEE
959 * Std 754-1985), the symbol IBMPC or MIEEE should be defined.
960 * These numbers have 53-bit significands. In this mode, constants
961 * are provided as arrays of hexadecimal 16 bit integers.
963 * To accommodate other types of computer arithmetic, all
964 * constants are also provided in a normal decimal radix
965 * which one can hope are correctly converted to a suitable
966 * format by the available C language compiler. To invoke
967 * this mode, the symbol UNK is defined.
969 * An important difference among these modes is a predefined
970 * set of machine arithmetic constants for each. The numbers
971 * MACHEP (the machine roundoff error), MAXNUM (largest number
972 * represented), and several other parameters are preset by
973 * the configuration symbol. Check the file const.c to
974 * ensure that these values are correct for your computer.
976 * For ANSI C compatibility, define ANSIC equal to 1. Currently
977 * this affects only the atan2 function and others that use it.
980 /* Constant definitions for math error conditions. */
982 #define DOMAIN 1 /* argument domain error */
983 #define SING 2 /* argument singularity */
984 #define OVERFLOW 3 /* overflow range error */
985 #define UNDERFLOW 4 /* underflow range error */
986 #define TLOSS 5 /* total loss of precision */
987 #define PLOSS 6 /* partial loss of precision */
992 /* e type constants used by high precision check routines */
994 /*include "ehead.h"*/
996 unsigned EMUSHORT ezero
[NE
] =
998 0, 0000000, 0000000, 0000000, 0000000, 0000000,};
999 extern unsigned EMUSHORT ezero
[];
1002 unsigned EMUSHORT ehalf
[NE
] =
1004 0, 0000000, 0000000, 0000000, 0100000, 0x3ffe,};
1005 extern unsigned EMUSHORT ehalf
[];
1008 unsigned EMUSHORT eone
[NE
] =
1010 0, 0000000, 0000000, 0000000, 0100000, 0x3fff,};
1011 extern unsigned EMUSHORT eone
[];
1014 unsigned EMUSHORT etwo
[NE
] =
1016 0, 0000000, 0000000, 0000000, 0100000, 0040000,};
1017 extern unsigned EMUSHORT etwo
[];
1020 unsigned EMUSHORT e32
[NE
] =
1022 0, 0000000, 0000000, 0000000, 0100000, 0040004,};
1023 extern unsigned EMUSHORT e32
[];
1025 /* 6.93147180559945309417232121458176568075500134360255E-1 */
1026 unsigned EMUSHORT elog2
[NE
] =
1028 0xc9e4, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1029 extern unsigned EMUSHORT elog2
[];
1031 /* 1.41421356237309504880168872420969807856967187537695E0 */
1032 unsigned EMUSHORT esqrt2
[NE
] =
1034 0x597e, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1035 extern unsigned EMUSHORT esqrt2
[];
1038 * 1.12837916709551257389615890312154517168810125865800E0 */
1039 unsigned EMUSHORT eoneopi
[NE
] =
1041 0x71d5, 0x688d, 0012333, 0135202, 0110156, 0x3fff,};
1042 extern unsigned EMUSHORT eoneopi
[];
1044 /* 3.14159265358979323846264338327950288419716939937511E0 */
1045 unsigned EMUSHORT epi
[NE
] =
1047 0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1048 extern unsigned EMUSHORT epi
[];
1050 /* 5.7721566490153286060651209008240243104215933593992E-1 */
1051 unsigned EMUSHORT eeul
[NE
] =
1053 0xd1be, 0xc7a4, 0076660, 0063743, 0111704, 0x3ffe,};
1054 extern unsigned EMUSHORT eeul
[];
1063 /* Control register for rounding precision.
1064 * This can be set to 80 (if NE=6), 64, 56, 53, or 24 bits.
1069 void eaddm (), esubm (), emdnorm (), asctoeg ();
1070 static void toe24 (), toe53 (), toe64 ();
1071 void eremain (), einit (), eiremain ();
1072 int ecmpm (), edivm (), emulm ();
1073 void emovi (), emovo (), emovz (), ecleaz (), ecleazs (), eadd1 ();
1074 void etodec (), todec (), dectoe ();
1085 ; Clear out entire external format number.
1087 ; unsigned EMUSHORT x[];
1093 register unsigned EMUSHORT
*x
;
1097 for (i
= 0; i
< NE
; i
++)
1103 /* Move external format number from a to b.
1110 register unsigned EMUSHORT
*a
, *b
;
1114 for (i
= 0; i
< NE
; i
++)
1120 ; Absolute value of external format number
1128 unsigned EMUSHORT x
[]; /* x is the memory address of a short */
1131 x
[NE
- 1] &= 0x7fff; /* sign is top bit of last word of external format */
1138 ; Negate external format number
1140 ; unsigned EMUSHORT x[NE];
1146 unsigned EMUSHORT x
[];
1149 x
[NE
- 1] ^= 0x8000; /* Toggle the sign bit */
1154 /* Return 1 if external format number is negative,
1159 unsigned EMUSHORT x
[];
1162 if (x
[NE
- 1] & 0x8000)
1169 /* Return 1 if external format number has maximum possible exponent,
1174 unsigned EMUSHORT x
[];
1177 if ((x
[NE
- 1] & 0x7fff) == 0x7fff)
1185 ; Fill entire number, including exponent and significand, with
1186 ; largest possible number. These programs implement a saturation
1187 ; value that is an ordinary, legal number. A special value
1188 ; "infinity" may also be implemented; this would require tests
1189 ; for that value and implementation of special rules for arithmetic
1190 ; operations involving inifinity.
1195 register unsigned EMUSHORT
*x
;
1200 for (i
= 0; i
< NE
- 1; i
++)
1204 for (i
= 0; i
< NE
- 1; i
++)
1229 /* Move in external format number,
1230 * converting it to internal format.
1234 unsigned EMUSHORT
*a
, *b
;
1236 register unsigned EMUSHORT
*p
, *q
;
1240 p
= a
+ (NE
- 1); /* point to last word of external number */
1241 /* get the sign bit */
1246 /* get the exponent */
1248 *q
++ &= 0x7fff; /* delete the sign bit */
1250 if ((*(q
- 1) & 0x7fff) == 0x7fff)
1252 for (i
= 2; i
< NI
; i
++)
1257 /* clear high guard word */
1259 /* move in the significand */
1260 for (i
= 0; i
< NE
- 1; i
++)
1262 /* clear low guard word */
1267 /* Move internal format number out,
1268 * converting it to external format.
1272 unsigned EMUSHORT
*a
, *b
;
1274 register unsigned EMUSHORT
*p
, *q
;
1275 unsigned EMUSHORT i
;
1278 q
= b
+ (NE
- 1); /* point to output exponent */
1279 /* combine sign and exponent */
1282 *q
-- = *p
++ | 0x8000;
1286 if (*(p
- 1) == 0x7fff)
1292 /* skip over guard word */
1294 /* move the significand */
1295 for (i
= 0; i
< NE
- 1; i
++)
1302 /* Clear out internal format number.
1307 register unsigned EMUSHORT
*xi
;
1311 for (i
= 0; i
< NI
; i
++)
1316 /* same, but don't touch the sign. */
1320 register unsigned EMUSHORT
*xi
;
1325 for (i
= 0; i
< NI
- 1; i
++)
1331 /* Move internal format number from a to b.
1335 register unsigned EMUSHORT
*a
, *b
;
1339 for (i
= 0; i
< NI
- 1; i
++)
1341 /* clear low guard word */
1347 ; Compare significands of numbers in internal format.
1348 ; Guard words are included in the comparison.
1350 ; unsigned EMUSHORT a[NI], b[NI];
1353 ; for the significands:
1354 ; returns +1 if a > b
1360 register unsigned EMUSHORT
*a
, *b
;
1364 a
+= M
; /* skip up to significand area */
1366 for (i
= M
; i
< NI
; i
++)
1374 if (*(--a
) > *(--b
))
1382 ; Shift significand down by 1 bit
1387 register unsigned EMUSHORT
*x
;
1389 register unsigned EMUSHORT bits
;
1392 x
+= M
; /* point to significand area */
1395 for (i
= M
; i
< NI
; i
++)
1410 ; Shift significand up by 1 bit
1415 register unsigned EMUSHORT
*x
;
1417 register unsigned EMUSHORT bits
;
1423 for (i
= M
; i
< NI
; i
++)
1438 ; Shift significand down by 8 bits
1443 register unsigned EMUSHORT
*x
;
1445 register unsigned EMUSHORT newbyt
, oldbyt
;
1450 for (i
= M
; i
< NI
; i
++)
1461 ; Shift significand up by 8 bits
1466 register unsigned EMUSHORT
*x
;
1469 register unsigned EMUSHORT newbyt
, oldbyt
;
1474 for (i
= M
; i
< NI
; i
++)
1485 ; Shift significand up by 16 bits
1490 register unsigned EMUSHORT
*x
;
1493 register unsigned EMUSHORT
*p
;
1498 for (i
= M
; i
< NI
- 1; i
++)
1505 ; Shift significand down by 16 bits
1510 register unsigned EMUSHORT
*x
;
1513 register unsigned EMUSHORT
*p
;
1518 for (i
= M
; i
< NI
- 1; i
++)
1531 unsigned EMUSHORT
*x
, *y
;
1533 register unsigned EMULONG a
;
1540 for (i
= M
; i
< NI
; i
++)
1542 a
= (unsigned EMULONG
) (*x
) + (unsigned EMULONG
) (*y
) + carry
;
1547 *y
= (unsigned EMUSHORT
) a
;
1554 ; Subtract significands
1560 unsigned EMUSHORT
*x
, *y
;
1569 for (i
= M
; i
< NI
; i
++)
1571 a
= (unsigned EMULONG
) (*y
) - (unsigned EMULONG
) (*x
) - carry
;
1576 *y
= (unsigned EMUSHORT
) a
;
1583 /* Divide significands */
1585 static unsigned EMUSHORT equot
[NI
];
1589 unsigned EMUSHORT den
[], num
[];
1592 register unsigned EMUSHORT
*p
, *q
;
1593 unsigned EMUSHORT j
;
1599 for (i
= M
; i
< NI
; i
++)
1604 /* Use faster compare and subtraction if denominator
1605 * has only 15 bits of significance.
1610 for (i
= M
+ 3; i
< NI
; i
++)
1615 if ((den
[M
+ 1] & 1) != 0)
1623 for (i
= 0; i
< NBITS
+ 2; i
++)
1641 /* The number of quotient bits to calculate is
1642 * NBITS + 1 scaling guard bit + 1 roundoff bit.
1647 for (i
= 0; i
< NBITS
+ 2; i
++)
1649 if (ecmpm (den
, num
) <= 0)
1652 j
= 1; /* quotient bit = 1 */
1666 /* test for nonzero remainder after roundoff bit */
1669 for (i
= M
; i
< NI
; i
++)
1677 for (i
= 0; i
< NI
; i
++)
1683 /* Multiply significands */
1686 unsigned EMUSHORT a
[], b
[];
1688 unsigned EMUSHORT
*p
, *q
;
1693 for (i
= M
; i
< NI
; i
++)
1698 while (*p
== 0) /* significand is not supposed to be all zero */
1703 if ((*p
& 0xff) == 0)
1711 for (i
= 0; i
< k
; i
++)
1715 /* remember if there were any nonzero bits shifted out */
1722 for (i
= 0; i
< NI
; i
++)
1725 /* return flag for lost nonzero bits */
1732 * Normalize and round off.
1734 * The internal format number to be rounded is "s".
1735 * Input "lost" indicates whether or not the number is exact.
1736 * This is the so-called sticky bit.
1738 * Input "subflg" indicates whether the number was obtained
1739 * by a subtraction operation. In that case if lost is nonzero
1740 * then the number is slightly smaller than indicated.
1742 * Input "exp" is the biased exponent, which may be negative.
1743 * the exponent field of "s" is ignored but is replaced by
1744 * "exp" as adjusted by normalization and rounding.
1746 * Input "rcntrl" is the rounding control.
1749 static int rlast
= -1;
1751 static unsigned EMUSHORT rmsk
= 0;
1752 static unsigned EMUSHORT rmbit
= 0;
1753 static unsigned EMUSHORT rebit
= 0;
1755 static unsigned EMUSHORT rbit
[NI
];
1758 emdnorm (s
, lost
, subflg
, exp
, rcntrl
)
1759 unsigned EMUSHORT s
[];
1766 unsigned EMUSHORT r
;
1771 /* a blank significand could mean either zero or infinity. */
1784 if ((j
> NBITS
) && (exp
< 32767))
1792 if (exp
> (EMULONG
) (-NBITS
- 1))
1805 /* Round off, unless told not to by rcntrl. */
1808 /* Set up rounding parameters if the control register changed. */
1809 if (rndprc
!= rlast
)
1816 rw
= NI
- 1; /* low guard word */
1831 /* For DEC arithmetic */
1880 /* These tests assume NI = 8 */
1900 if ((r
& rmbit
) != 0)
1905 { /* round to even */
1906 if ((s
[re
] & rebit
) == 0)
1918 if ((rndprc
< 64) && (exp
<= 0))
1923 { /* overflow on roundoff */
1936 for (i
= 2; i
< NI
- 1; i
++)
1938 pedwarn ("floating point overflow");
1942 for (i
= M
+ 1; i
< NI
- 1; i
++)
1960 s
[1] = (unsigned EMUSHORT
) exp
;
1966 ; Subtract external format numbers.
1968 ; unsigned EMUSHORT a[NE], b[NE], c[NE];
1969 ; esub (a, b, c); c = b - a
1972 static int subflg
= 0;
1976 unsigned EMUSHORT
*a
, *b
, *c
;
1987 ; unsigned EMUSHORT a[NE], b[NE], c[NE];
1988 ; eadd (a, b, c); c = b + a
1992 unsigned EMUSHORT
*a
, *b
, *c
;
2001 unsigned EMUSHORT
*a
, *b
, *c
;
2003 unsigned EMUSHORT ai
[NI
], bi
[NI
], ci
[NI
];
2005 EMULONG lt
, lta
, ltb
;
2026 /* compare exponents */
2031 { /* put the larger number in bi */
2041 if (lt
< (EMULONG
) (-NBITS
- 1))
2042 goto done
; /* answer same as larger addend */
2044 lost
= eshift (ai
, k
); /* shift the smaller number down */
2048 /* exponents were the same, so must compare significands */
2051 { /* the numbers are identical in magnitude */
2052 /* if different signs, result is zero */
2058 /* if same sign, result is double */
2059 /* double denomalized tiny number */
2060 if ((bi
[E
] == 0) && ((bi
[3] & 0x8000) == 0))
2065 /* add 1 to exponent unless both are zero! */
2066 for (j
= 1; j
< NI
- 1; j
++)
2070 /* This could overflow, but let emovo take care of that. */
2075 bi
[E
] = (unsigned EMUSHORT
) ltb
;
2079 { /* put the larger number in bi */
2095 emdnorm (bi
, lost
, subflg
, ltb
, 64);
2106 ; unsigned EMUSHORT a[NE], b[NE], c[NE];
2107 ; ediv (a, b, c); c = b / a
2111 unsigned EMUSHORT
*a
, *b
, *c
;
2113 unsigned EMUSHORT ai
[NI
], bi
[NI
];
2115 EMULONG lt
, lta
, ltb
;
2120 if (eisneg (a
) ^ eisneg (b
))
2121 *(c
+ (NE
- 1)) = 0x8000;
2123 *(c
+ (NE
- 1)) = 0;
2138 { /* See if numerator is zero. */
2139 for (i
= 1; i
< NI
- 1; i
++)
2143 ltb
-= enormlz (bi
);
2153 { /* possible divide by zero */
2154 for (i
= 1; i
< NI
- 1; i
++)
2158 lta
-= enormlz (ai
);
2163 *(c
+ (NE
- 1)) = 0;
2165 *(c
+ (NE
- 1)) = 0x8000;
2167 mtherr ("ediv", SING
);
2173 /* calculate exponent */
2174 lt
= ltb
- lta
+ EXONE
;
2175 emdnorm (bi
, i
, 0, lt
, 64);
2189 ; unsigned EMUSHORT a[NE], b[NE], c[NE];
2190 ; emul (a, b, c); c = b * a
2194 unsigned EMUSHORT
*a
, *b
, *c
;
2196 unsigned EMUSHORT ai
[NI
], bi
[NI
];
2198 EMULONG lt
, lta
, ltb
;
2201 if (eisinf (a
) || eisinf (b
))
2203 if (eisneg (a
) ^ eisneg (b
))
2204 *(c
+ (NE
- 1)) = 0x8000;
2206 *(c
+ (NE
- 1)) = 0;
2217 for (i
= 1; i
< NI
- 1; i
++)
2221 lta
-= enormlz (ai
);
2232 for (i
= 1; i
< NI
- 1; i
++)
2236 ltb
-= enormlz (bi
);
2245 /* Multiply significands */
2247 /* calculate exponent */
2248 lt
= lta
+ ltb
- (EXONE
- 1);
2249 emdnorm (bi
, j
, 0, lt
, 64);
2250 /* calculate sign of product */
2262 ; Convert IEEE double precision to e type
2264 ; unsigned EMUSHORT x[N+2];
2269 unsigned EMUSHORT
*e
, *y
;
2273 dectoe (e
, y
); /* see etodec.c */
2277 register unsigned EMUSHORT r
;
2278 register unsigned EMUSHORT
*p
;
2279 unsigned EMUSHORT yy
[NI
];
2282 denorm
= 0; /* flag if denormalized number */
2291 yy
[M
] = (r
& 0x0f) | 0x10;
2292 r
&= ~0x800f; /* strip sign and 4 significand bits */
2303 /* If zero exponent, then the significand is denormalized.
2304 * So, take back the understood high significand bit. */
2324 (void) eshift (yy
, -5);
2326 { /* if zero exponent, then normalize the significand */
2327 if ((k
= enormlz (yy
)) > NBITS
)
2330 yy
[E
] -= (unsigned EMUSHORT
) (k
- 1);
2333 #endif /* not DEC */
2338 unsigned EMUSHORT
*e
, *y
;
2340 unsigned EMUSHORT yy
[NI
];
2341 unsigned EMUSHORT
*p
, *q
;
2345 for (i
= 0; i
< NE
- 5; i
++)
2348 for (i
= 0; i
< 5; i
++)
2352 for (i
= 0; i
< 5; i
++)
2356 p
= &yy
[0] + (NE
- 1);
2359 for (i
= 0; i
< 4; i
++)
2373 for (i
= 0; i
< NE
; i
++)
2379 ; Convert IEEE single precision to e type
2381 ; unsigned EMUSHORT x[N+2];
2386 unsigned EMUSHORT
*e
, *y
;
2388 register unsigned EMUSHORT r
;
2389 register unsigned EMUSHORT
*p
;
2390 unsigned EMUSHORT yy
[NI
];
2393 denorm
= 0; /* flag if denormalized number */
2405 yy
[M
] = (r
& 0x7f) | 0200;
2406 r
&= ~0x807f; /* strip sign and 7 significand bits */
2417 /* If zero exponent, then the significand is denormalized.
2418 * So, take back the understood high significand bit. */
2437 (void) eshift (yy
, -8);
2439 { /* if zero exponent, then normalize the significand */
2440 if ((k
= enormlz (yy
)) > NBITS
)
2443 yy
[E
] -= (unsigned EMUSHORT
) (k
- 1);
2451 unsigned EMUSHORT
*x
, *e
;
2453 unsigned EMUSHORT xi
[NI
];
2458 /* adjust exponent for offset */
2459 exp
= (EMULONG
) xi
[E
];
2464 /* round off to nearest or even */
2467 emdnorm (xi
, 0, 0, exp
, 64);
2473 /* move out internal format to ieee long double */
2476 unsigned EMUSHORT
*a
, *b
;
2478 register unsigned EMUSHORT
*p
, *q
;
2479 unsigned EMUSHORT i
;
2485 q
= b
+ 4; /* point to output exponent */
2486 #if LONG_DOUBLE_TYPE_SIZE == 96
2487 /* Clear the last two bytes of 12-byte Intel format */
2492 /* combine sign and exponent */
2496 *q
++ = *p
++ | 0x8000;
2502 *q
-- = *p
++ | 0x8000;
2506 /* skip over guard word */
2508 /* move the significand */
2510 for (i
= 0; i
< 4; i
++)
2513 for (i
= 0; i
< 4; i
++)
2520 ; e type to IEEE double precision
2522 ; unsigned EMUSHORT x[NE];
2530 unsigned EMUSHORT
*x
, *e
;
2532 etodec (x
, e
); /* see etodec.c */
2537 unsigned EMUSHORT
*x
, *y
;
2546 unsigned EMUSHORT
*x
, *e
;
2548 unsigned EMUSHORT xi
[NI
];
2553 /* adjust exponent for offsets */
2554 exp
= (EMULONG
) xi
[E
] - (EXONE
- 0x3ff);
2559 /* round off to nearest or even */
2562 emdnorm (xi
, 0, 0, exp
, 64);
2571 unsigned EMUSHORT
*x
, *y
;
2573 unsigned EMUSHORT i
;
2574 unsigned EMUSHORT
*p
;
2581 *y
= 0; /* output high order */
2583 *y
= 0x8000; /* output sign bit */
2586 if (i
>= (unsigned int) 2047)
2587 { /* Saturate at largest number less than infinity. */
2602 *y
|= (unsigned EMUSHORT
) 0x7fef;
2619 (void) eshift (x
, 4);
2624 (void) eshift (x
, 5);
2626 i
|= *p
++ & (unsigned EMUSHORT
) 0x0f; /* *p = xi[M] */
2627 *y
|= (unsigned EMUSHORT
) i
; /* high order output already has sign bit set */
2641 #endif /* not DEC */
2646 ; e type to IEEE single precision
2648 ; unsigned EMUSHORT x[N+2];
2653 unsigned EMUSHORT
*x
, *e
;
2656 unsigned EMUSHORT xi
[NI
];
2660 /* adjust exponent for offsets */
2661 exp
= (EMULONG
) xi
[E
] - (EXONE
- 0177);
2666 /* round off to nearest or even */
2669 emdnorm (xi
, 0, 0, exp
, 64);
2677 unsigned EMUSHORT
*x
, *y
;
2679 unsigned EMUSHORT i
;
2680 unsigned EMUSHORT
*p
;
2689 *y
= 0; /* output high order */
2691 *y
= 0x8000; /* output sign bit */
2695 { /* Saturate at largest number less than infinity. */
2697 *y
|= (unsigned EMUSHORT
) 0x7f80;
2709 *y
|= (unsigned EMUSHORT
) 0x7f7f;
2725 (void) eshift (x
, 7);
2730 (void) eshift (x
, 8);
2732 i
|= *p
++ & (unsigned EMUSHORT
) 0x7f; /* *p = xi[M] */
2733 *y
|= i
; /* high order output already has sign bit set */
2747 /* Compare two e type numbers.
2749 * unsigned EMUSHORT a[NE], b[NE];
2752 * returns +1 if a > b
2758 unsigned EMUSHORT
*a
, *b
;
2760 unsigned EMUSHORT ai
[NI
], bi
[NI
];
2761 register unsigned EMUSHORT
*p
, *q
;
2771 { /* the signs are different */
2773 for (i
= 1; i
< NI
- 1; i
++)
2787 /* both are the same sign */
2802 return (0); /* equality */
2808 if (*(--p
) > *(--q
))
2809 return (msign
); /* p is bigger */
2811 return (-msign
); /* p is littler */
2817 /* Find nearest integer to x = floor (x + 0.5)
2819 * unsigned EMUSHORT x[NE], y[NE]
2824 unsigned EMUSHORT
*x
, *y
;
2834 ; convert long integer to e type
2837 ; unsigned EMUSHORT x[NE];
2839 ; note &l is the memory address of l
2843 long *lp
; /* lp is the memory address of a long integer */
2844 unsigned EMUSHORT
*y
; /* y is the address of a short */
2846 unsigned EMUSHORT yi
[NI
];
2853 /* make it positive */
2854 ll
= (unsigned long) (-(*lp
));
2855 yi
[0] = 0xffff; /* put correct sign in the e type number */
2859 ll
= (unsigned long) (*lp
);
2861 /* move the long integer to yi significand area */
2862 yi
[M
] = (unsigned EMUSHORT
) (ll
>> 16);
2863 yi
[M
+ 1] = (unsigned EMUSHORT
) ll
;
2865 yi
[E
] = EXONE
+ 15; /* exponent if normalize shift count were 0 */
2866 if ((k
= enormlz (yi
)) > NBITS
)/* normalize the significand */
2867 ecleaz (yi
); /* it was zero */
2869 yi
[E
] -= (unsigned EMUSHORT
) k
;/* subtract shift count from exponent */
2870 emovo (yi
, y
); /* output the answer */
2874 ; convert unsigned long integer to e type
2877 ; unsigned EMUSHORT x[NE];
2879 ; note &l is the memory address of l
2883 unsigned long *lp
; /* lp is the memory address of a long integer */
2884 unsigned EMUSHORT
*y
; /* y is the address of a short */
2886 unsigned EMUSHORT yi
[NI
];
2893 /* move the long integer to ayi significand area */
2894 yi
[M
] = (unsigned EMUSHORT
) (ll
>> 16);
2895 yi
[M
+ 1] = (unsigned EMUSHORT
) ll
;
2897 yi
[E
] = EXONE
+ 15; /* exponent if normalize shift count were 0 */
2898 if ((k
= enormlz (yi
)) > NBITS
)/* normalize the significand */
2899 ecleaz (yi
); /* it was zero */
2901 yi
[E
] -= (unsigned EMUSHORT
) k
; /* subtract shift count from exponent */
2902 emovo (yi
, y
); /* output the answer */
2907 ; Find long integer and fractional parts
2910 ; unsigned EMUSHORT x[NE], frac[NE];
2911 ; xifrac (x, &i, frac);
2913 The integer output has the sign of the input. The fraction is
2914 the positive fractional part of abs (x).
2918 unsigned EMUSHORT
*x
;
2920 unsigned EMUSHORT
*frac
;
2922 unsigned EMUSHORT xi
[NI
];
2926 k
= (int) xi
[E
] - (EXONE
- 1);
2929 /* if exponent <= 0, integer = 0 and real output is fraction */
2934 if (k
> (HOST_BITS_PER_LONG
- 1))
2937 ; long integer overflow: output large integer
2938 ; and correct fraction
2941 *i
= ((unsigned long) 1) << (HOST_BITS_PER_LONG
- 1);
2943 *i
= (((unsigned long) 1) << (HOST_BITS_PER_LONG
- 1)) - 1;
2944 (void) eshift (xi
, k
);
2945 pedwarn ("Overflow on truncation to integer");
2952 ; shift more than 16 bits: shift up k-16, output the integer,
2953 ; then complete the shift to get the fraction.
2956 (void) eshift (xi
, k
);
2958 *i
= (long) (((unsigned long) xi
[M
] << 16) | xi
[M
+ 1]);
2963 /* shift not more than 16 bits */
2964 (void) eshift (xi
, k
);
2965 *i
= (long) xi
[M
] & 0xffff;
2976 if ((k
= enormlz (xi
)) > NBITS
)
2979 xi
[E
] -= (unsigned EMUSHORT
) k
;
2986 ; Find unsigned long integer and fractional parts
2989 ; unsigned EMUSHORT x[NE], frac[NE];
2990 ; xifrac (x, &i, frac);
2992 A negative e type input yields integer output = 0
2993 but correct fraction.
2996 euifrac (x
, i
, frac
)
2997 unsigned EMUSHORT
*x
;
2999 unsigned EMUSHORT
*frac
;
3001 unsigned EMUSHORT xi
[NI
];
3005 k
= (int) xi
[E
] - (EXONE
- 1);
3008 /* if exponent <= 0, integer = 0 and argument is fraction */
3016 ; long integer overflow: output large integer
3017 ; and correct fraction
3020 (void) eshift (xi
, k
);
3021 pedwarn ("Overflow on truncation to unsigned integer");
3028 ; shift more than 16 bits: shift up k-16, output the integer,
3029 ; then complete the shift to get the fraction.
3032 (void) eshift (xi
, k
);
3034 *i
= (long) (((unsigned long) xi
[M
] << 16) | xi
[M
+ 1]);
3039 /* shift not more than 16 bits */
3040 (void) eshift (xi
, k
);
3041 *i
= (long) xi
[M
] & 0xffff;
3051 if ((k
= enormlz (xi
)) > NBITS
)
3054 xi
[E
] -= (unsigned EMUSHORT
) k
;
3064 ; Shifts significand area up or down by the number of bits
3065 ; given by the variable sc.
3069 unsigned EMUSHORT
*x
;
3072 unsigned EMUSHORT lost
;
3073 unsigned EMUSHORT
*p
;
3086 lost
|= *p
; /* remember lost bits */
3127 return ((int) lost
);
3135 ; Shift normalizes the significand area pointed to by argument
3136 ; shift count (up = positive) is returned.
3140 unsigned EMUSHORT x
[];
3142 register unsigned EMUSHORT
*p
;
3151 return (0); /* already normalized */
3156 /* With guard word, there are NBITS+16 bits available.
3157 * return true if all are zero.
3162 /* see if high byte is zero */
3163 while ((*p
& 0xff00) == 0)
3168 /* now shift 1 bit at a time */
3169 while ((*p
& 0x8000) == 0)
3175 mtherr ("enormlz", UNDERFLOW
);
3181 /* Normalize by shifting down out of the high guard word
3182 of the significand */
3197 mtherr ("enormlz", OVERFLOW
);
3207 /* Convert e type number to decimal format ASCII string.
3208 * The constants are for 64 bit precision.
3214 static unsigned EMUSHORT etens
[NTEN
+ 1][NE
] =
3216 {0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,}, /* 10**4096 */
3217 {0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,}, /* 10**2048 */
3218 {0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
3219 {0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
3220 {0xddbc, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
3221 {0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
3222 {0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
3223 {0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
3224 {0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
3225 {0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
3226 {0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
3227 {0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
3228 {0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,}, /* 10**1 */
3231 static unsigned EMUSHORT emtens
[NTEN
+ 1][NE
] =
3233 {0x2de4, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,}, /* 10**-4096 */
3234 {0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,}, /* 10**-2048 */
3235 {0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
3236 {0x7133, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
3237 {0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
3238 {0xac7d, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
3239 {0x3f24, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
3240 {0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
3241 {0x4c2f, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
3242 {0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
3243 {0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
3244 {0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
3245 {0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,}, /* 10**-1 */
3249 e24toasc (x
, string
, ndigs
)
3250 unsigned EMUSHORT x
[];
3254 unsigned EMUSHORT w
[NI
];
3258 if ((x
[1] & 0x7f80) == 0x7f80)
3260 if ((x
[0] & 0x7f80) == 0x7f80)
3268 sprintf (string
, " -Infinity ");
3270 sprintf (string
, " Infinity ");
3275 etoasc (w
, string
, ndigs
);
3280 e53toasc (x
, string
, ndigs
)
3281 unsigned EMUSHORT x
[];
3285 unsigned EMUSHORT w
[NI
];
3289 if ((x
[3] & 0x7ff0) == 0x7ff0)
3291 if ((x
[0] & 0x7ff0) == 0x7ff0)
3299 sprintf (string
, " -Infinity ");
3301 sprintf (string
, " Infinity ");
3306 etoasc (w
, string
, ndigs
);
3311 e64toasc (x
, string
, ndigs
)
3312 unsigned EMUSHORT x
[];
3316 unsigned EMUSHORT w
[NI
];
3320 if ((x
[4] & 0x7fff) == 0x7fff)
3322 if ((x
[0] & 0x7fff) == 0x7fff)
3330 sprintf (string
, " -Infinity ");
3332 sprintf (string
, " Infinity ");
3337 etoasc (w
, string
, ndigs
);
3341 static char wstring
[80]; /* working storage for ASCII output */
3344 etoasc (x
, string
, ndigs
)
3345 unsigned EMUSHORT x
[];
3350 unsigned EMUSHORT y
[NI
], t
[NI
], u
[NI
], w
[NI
];
3351 unsigned EMUSHORT
*p
, *r
, *ten
;
3352 unsigned EMUSHORT sign
;
3353 int i
, j
, k
, expon
, rndsav
;
3355 unsigned EMUSHORT m
;
3359 while ((*s
++ = *ss
++) != '\0')
3362 rndprc
= NBITS
; /* set to full precision */
3363 emov (x
, y
); /* retain external format */
3364 if (y
[NE
- 1] & 0x8000)
3367 y
[NE
- 1] &= 0x7fff;
3374 ten
= &etens
[NTEN
][0];
3376 /* Test for zero exponent */
3379 for (k
= 0; k
< NE
- 1; k
++)
3382 goto tnzro
; /* denormalized number */
3384 goto isone
; /* legal all zeros */
3388 /* Test for infinity. Don't bother with illegal infinities.
3390 if (y
[NE
- 1] == 0x7fff)
3393 sprintf (wstring
, " -Infinity ");
3395 sprintf (wstring
, " Infinity ");
3399 /* Test for exponent nonzero but significand denormalized.
3400 * This is an error condition.
3402 if ((y
[NE
- 1] != 0) && ((y
[NE
- 2] & 0x8000) == 0))
3404 mtherr ("etoasc", DOMAIN
);
3405 sprintf (wstring
, "NaN");
3409 /* Compare to 1.0 */
3415 { /* Number is greater than 1 */
3416 /* Convert significand to an integer and strip trailing decimal zeros. */
3418 u
[NE
- 1] = EXONE
+ NBITS
- 1;
3420 p
= &etens
[NTEN
- 4][0];
3426 for (j
= 0; j
< NE
- 1; j
++)
3439 /* Rescale from integer significand */
3440 u
[NE
- 1] += y
[NE
- 1] - (unsigned int) (EXONE
+ NBITS
- 1);
3442 /* Find power of 10 */
3446 while (ecmp (ten
, u
) <= 0)
3448 if (ecmp (p
, u
) <= 0)
3461 { /* Number is less than 1.0 */
3462 /* Pad significand with trailing decimal zeros. */
3465 while ((y
[NE
- 2] & 0x8000) == 0)
3474 for (i
= 0; i
< NDEC
+ 1; i
++)
3476 if ((w
[NI
- 1] & 0x7) != 0)
3478 /* multiply by 10 */
3491 if (eone
[NE
- 1] <= u
[1])
3503 while (ecmp (eone
, w
) > 0)
3505 if (ecmp (p
, w
) >= 0)
3520 /* Find the first (leading) digit. */
3526 digit
= equot
[NI
- 1];
3527 while ((digit
== 0) && (ecmp (y
, ezero
) != 0))
3535 digit
= equot
[NI
- 1];
3543 *s
++ = (char) digit
+ '0';
3545 /* Examine number of digits requested by caller. */
3550 /* Generate digits after the decimal point. */
3551 for (k
= 0; k
<= ndigs
; k
++)
3553 /* multiply current number by 10, without normalizing */
3560 *s
++ = (char) equot
[NI
- 1] + '0';
3562 digit
= equot
[NI
- 1];
3565 /* round off the ASCII string */
3568 /* Test for critical rounding case in ASCII output. */
3572 if (ecmp (t
, ezero
) != 0)
3573 goto roun
; /* round to nearest */
3574 if ((*(s
- 1) & 1) == 0)
3575 goto doexp
; /* round to even */
3577 /* Round up and propagate carry-outs */
3581 /* Carry out to most significant digit? */
3588 /* Most significant digit carries to 10? */
3596 /* Round up and carry out from less significant digits */
3608 sprintf (ss, "e+%d", expon);
3610 sprintf (ss, "e%d", expon);
3612 sprintf (ss
, "e%d", expon
);
3615 /* copy out the working string */
3618 while (*ss
== ' ') /* strip possible leading space */
3620 while ((*s
++ = *ss
++) != '\0')
3629 ; ASCTOQ.MAC LATEST REV: 11 JAN 84
3632 ; Convert ASCII string to quadruple precision floating point
3634 ; Numeric input is free field decimal number
3635 ; with max of 15 digits with or without
3636 ; decimal point entered as ASCII from teletype.
3637 ; Entering E after the number followed by a second
3638 ; number causes the second number to be interpreted
3639 ; as a power of 10 to be multiplied by the first number
3640 ; (i.e., "scientific" notation).
3643 ; asctoq (string, q);
3646 /* ASCII to single */
3650 unsigned EMUSHORT
*y
;
3656 /* ASCII to double */
3660 unsigned EMUSHORT
*y
;
3670 /* ASCII to long double */
3674 unsigned EMUSHORT
*y
;
3679 /* ASCII to super double */
3683 unsigned EMUSHORT
*y
;
3685 asctoeg (s
, y
, NBITS
);
3688 /* Space to make a copy of the input string: */
3689 static char lstr
[82];
3692 asctoeg (ss
, y
, oprec
)
3694 unsigned EMUSHORT
*y
;
3697 unsigned EMUSHORT yy
[NI
], xt
[NI
], tt
[NI
];
3698 int esign
, decflg
, sgnflg
, nexp
, exp
, prec
, lost
;
3699 int k
, trail
, c
, rndsav
;
3701 unsigned EMUSHORT nsign
, *p
;
3704 /* Copy the input string. */
3706 while (*s
== ' ') /* skip leading spaces */
3709 for (k
= 0; k
< 79; k
++)
3711 if ((*sp
++ = *s
++) == '\0')
3718 rndprc
= NBITS
; /* Set to full precision */
3731 if ((k
>= 0) && (k
<= 9))
3733 /* Ignore leading zeros */
3734 if ((prec
== 0) && (decflg
== 0) && (k
== 0))
3736 /* Identify and strip trailing zeros after the decimal point. */
3737 if ((trail
== 0) && (decflg
!= 0))
3740 while ((*sp
>= '0') && (*sp
<= '9'))
3742 /* Check for syntax error */
3744 if ((c
!= 'e') && (c
!= 'E') && (c
!= '\0')
3745 && (c
!= '\n') && (c
!= '\r') && (c
!= ' ')
3755 /* If enough digits were given to more than fill up the yy register,
3756 * continuing until overflow into the high guard word yy[2]
3757 * guarantees that there will be a roundoff bit at the top
3758 * of the low guard word after normalization.
3763 nexp
+= 1; /* count digits after decimal point */
3764 eshup1 (yy
); /* multiply current number by 10 */
3770 xt
[NI
- 2] = (unsigned EMUSHORT
) k
;
3788 case '.': /* decimal point */
3813 yy
[E
] = 0x7fff; /* infinity */
3817 mtherr ("asctoe", DOMAIN
);
3825 /* Exponent interpretation */
3831 /* check for + or - */
3839 while ((*s
>= '0') && (*s
<= '9'))
3849 /* Pad trailing zeros to minimize power of 10, per IEEE spec. */
3850 while ((nexp
> 0) && (yy
[2] == 0))
3862 if ((k
= enormlz (yy
)) > NBITS
)
3867 lexp
= (EXONE
- 1 + NBITS
) - k
;
3868 emdnorm (yy
, lost
, 0, lexp
, 64);
3869 /* convert to external format */
3872 /* Multiply by 10**nexp. If precision is 64 bits,
3873 * the maximum relative error incurred in forming 10**n
3874 * for 0 <= n <= 324 is 8.2e-20, at 10**180.
3875 * For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
3876 * For 0 >= n >= -999, it is -1.55e-19 at 10**-435.
3890 { /* Punt. Can't handle this without 2 divides. */
3891 emovi (etens
[0], tt
);
3898 p
= &etens
[NTEN
][0];
3908 while (exp
<= MAXP
);
3926 /* Round and convert directly to the destination type */
3928 lexp
-= EXONE
- 0x3ff;
3929 else if (oprec
== 24)
3930 lexp
-= EXONE
- 0177;
3932 else if (oprec
== 56)
3933 lexp
-= EXONE
- 0201;
3936 emdnorm (yy
, k
, 0, lexp
, 64);
3946 todec (yy
, y
); /* see etodec.c */
3966 /* y = largest integer not greater than x
3967 * (truncated toward minus infinity)
3969 * unsigned EMUSHORT x[NE], y[NE]
3973 static unsigned EMUSHORT bmask
[] =
3996 unsigned EMUSHORT x
[], y
[];
3998 register unsigned EMUSHORT
*p
;
4000 unsigned EMUSHORT f
[NE
];
4002 emov (x
, f
); /* leave in external format */
4003 expon
= (int) f
[NE
- 1];
4004 e
= (expon
& 0x7fff) - (EXONE
- 1);
4010 /* number of bits to clear out */
4022 /* clear the remaining bits */
4024 /* truncate negatives toward minus infinity */
4027 if ((unsigned EMUSHORT
) expon
& (unsigned EMUSHORT
) 0x8000)
4029 for (i
= 0; i
< NE
- 1; i
++)
4041 /* unsigned EMUSHORT x[], s[];
4044 * efrexp (x, exp, s);
4046 * Returns s and exp such that s * 2**exp = x and .5 <= s < 1.
4047 * For example, 1.1 = 0.55 * 2**1
4048 * Handles denormalized numbers properly using long integer exp.
4052 unsigned EMUSHORT x
[];
4054 unsigned EMUSHORT s
[];
4056 unsigned EMUSHORT xi
[NI
];
4060 li
= (EMULONG
) ((EMUSHORT
) xi
[1]);
4068 *exp
= (int) (li
- 0x3ffe);
4073 /* unsigned EMUSHORT x[], y[];
4076 * eldexp (x, pwr2, y);
4078 * Returns y = x * 2**pwr2.
4082 unsigned EMUSHORT x
[];
4084 unsigned EMUSHORT y
[];
4086 unsigned EMUSHORT xi
[NI
];
4094 emdnorm (xi
, i
, i
, li
, 64);
4099 /* c = remainder after dividing b by a
4100 * Least significant integer quotient bits left in equot[].
4104 unsigned EMUSHORT a
[], b
[], c
[];
4106 unsigned EMUSHORT den
[NI
], num
[NI
];
4108 if (ecmp (a
, ezero
) == 0)
4110 mtherr ("eremain", SING
);
4116 eiremain (den
, num
);
4117 /* Sign of remainder = sign of quotient */
4127 unsigned EMUSHORT den
[], num
[];
4130 unsigned EMUSHORT j
;
4133 ld
-= enormlz (den
);
4135 ln
-= enormlz (num
);
4139 if (ecmpm (den
, num
) <= 0)
4153 emdnorm (num
, 0, 0, ln
, 0);
4158 * Library common error handling routine
4168 * mtherr (fctnam, code);
4174 * This routine may be called to report one of the following
4175 * error conditions (in the include file mconf.h).
4177 * Mnemonic Value Significance
4179 * DOMAIN 1 argument domain error
4180 * SING 2 function singularity
4181 * OVERFLOW 3 overflow range error
4182 * UNDERFLOW 4 underflow range error
4183 * TLOSS 5 total loss of precision
4184 * PLOSS 6 partial loss of precision
4185 * EDOM 33 Unix domain error code
4186 * ERANGE 34 Unix range error code
4188 * The default version of the file prints the function name,
4189 * passed to it by the pointer fctnam, followed by the
4190 * error condition. The display is directed to the standard
4191 * output device. The routine then returns to the calling
4192 * program. Users may wish to modify the program to abort by
4193 * calling exit under severe error conditions such as domain
4196 * Since all error conditions pass control to this function,
4197 * the display may be easily changed, eliminated, or directed
4198 * to an error logging device.
4207 Cephes Math Library Release 2.0: April, 1987
4208 Copyright 1984, 1987 by Stephen L. Moshier
4209 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
4212 /* include "mconf.h" */
4214 /* Notice: the order of appearance of the following
4215 * messages is bound to the error codes defined
4218 static char *ermsg
[7] =
4220 "unknown", /* error code 0 */
4221 "domain", /* error code 1 */
4222 "singularity", /* et seq. */
4225 "total loss of precision",
4226 "partial loss of precision"
4239 /* Display string passed by calling program,
4240 * which is supposed to be the name of the
4241 * function in which the error occurred.
4244 /* Display error message defined
4245 * by the code argument.
4247 if ((code
<= 0) || (code
>= 6))
4249 sprintf (errstr
, "\n%s %s error\n", name
, ermsg
[code
]);
4251 /* Set global error message word */
4254 /* Return to calling
4259 /* Here is etodec.c .
4264 ; convert DEC double precision to e type
4271 unsigned EMUSHORT
*d
;
4272 unsigned EMUSHORT
*e
;
4274 unsigned EMUSHORT y
[NI
];
4275 register unsigned EMUSHORT r
, *p
;
4277 ecleaz (y
); /* start with a zero */
4278 p
= y
; /* point to our number */
4279 r
= *d
; /* get DEC exponent word */
4280 if (*d
& (unsigned int) 0x8000)
4281 *p
= 0xffff; /* fill in our sign */
4282 ++p
; /* bump pointer to our exponent word */
4283 r
&= 0x7fff; /* strip the sign bit */
4284 if (r
== 0) /* answer = 0 if high order DEC word = 0 */
4288 r
>>= 7; /* shift exponent word down 7 bits */
4289 r
+= EXONE
- 0201; /* subtract DEC exponent offset */
4290 /* add our e type exponent offset */
4291 *p
++ = r
; /* to form our exponent */
4293 r
= *d
++; /* now do the high order mantissa */
4294 r
&= 0177; /* strip off the DEC exponent and sign bits */
4295 r
|= 0200; /* the DEC understood high order mantissa bit */
4296 *p
++ = r
; /* put result in our high guard word */
4298 *p
++ = *d
++; /* fill in the rest of our mantissa */
4302 eshdn8 (y
); /* shift our mantissa down 8 bits */
4310 ; convert e type to DEC double precision
4316 static unsigned EMUSHORT decbit
[NI
] = {0, 0, 0, 0, 0, 0, 0200, 0};
4320 unsigned EMUSHORT
*x
, *d
;
4322 unsigned EMUSHORT xi
[NI
];
4323 register unsigned EMUSHORT r
;
4331 if (r
< (EXONE
- 128))
4334 if ((i
& 0200) != 0)
4336 if ((i
& 0377) == 0200)
4338 if ((i
& 0400) != 0)
4340 /* check all less significant bits */
4341 for (j
= M
+ 5; j
< NI
; j
++)
4390 unsigned EMUSHORT
*x
, *d
;
4392 unsigned EMUSHORT xi
[NI
];
4397 exp
= (EMULONG
) xi
[E
] - (EXONE
- 0201); /* adjust exponent for offsets */
4398 /* round off to nearest or even */
4401 emdnorm (xi
, 0, 0, exp
, 64);
4408 unsigned EMUSHORT
*x
, *y
;
4410 unsigned EMUSHORT i
;
4411 unsigned EMUSHORT
*p
;
4447 #endif /* EMU_NON_COMPILE not defined */
This page took 0.22875 seconds and 5 git commands to generate.