* pa/pa.c (compute_movstrsi_length): Handle residuals correctly.
[official-gcc.git] / gcc / real.c
blobf7e22eae3bbee36c817fc9da86ffed201c76c1fc
1 /* real.c - implementation of REAL_ARITHMETIC, REAL_VALUE_ATOF,
2 and support for XFmode IEEE extended real floating point arithmetic.
3 Copyright (C) 1993, 1994, 1995, 1996 Free Software Foundation, Inc.
4 Contributed by Stephen L. Moshier (moshier@world.std.com).
6 This file is part of GNU CC.
8 GNU CC is free software; you can redistribute it and/or modify
9 it under the terms of the GNU General Public License as published by
10 the Free Software Foundation; either version 2, or (at your option)
11 any later version.
13 GNU CC is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
18 You should have received a copy of the GNU General Public License
19 along with GNU CC; see the file COPYING. If not, write to
20 the Free Software Foundation, 59 Temple Place - Suite 330,
21 Boston, MA 02111-1307, USA. */
23 #include <stdio.h>
24 #include <errno.h>
25 #include "config.h"
26 #include "tree.h"
28 #ifndef errno
29 extern int errno;
30 #endif
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, VAX and IBM 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
46 real.h).
48 The first part of this file interfaces gcc to a floating point
49 arithmetic suite that was not written with gcc in mind. Avoid
50 changing the low-level arithmetic routines unless you have suitable
51 test programs available. A special version of the PARANOIA floating
52 point arithmetic tester, modified for this purpose, can be found on
53 usc.edu: /pub/C-numanal/ieeetest.zoo. Other tests, and libraries of
54 XFmode and TFmode transcendental functions, can be obtained by ftp from
55 netlib.att.com: netlib/cephes. */
57 /* Type of computer arithmetic.
58 Only one of DEC, IBM, IEEE, or UNK should get defined.
60 `IEEE', when REAL_WORDS_BIG_ENDIAN is non-zero, refers generically
61 to big-endian IEEE floating-point data structure. This definition
62 should work in SFmode `float' type and DFmode `double' type on
63 virtually all big-endian IEEE machines. If LONG_DOUBLE_TYPE_SIZE
64 has been defined to be 96, then IEEE also invokes the particular
65 XFmode (`long double' type) data structure used by the Motorola
66 680x0 series processors.
68 `IEEE', when REAL_WORDS_BIG_ENDIAN is zero, refers generally to
69 little-endian IEEE machines. In this case, if LONG_DOUBLE_TYPE_SIZE
70 has been defined to be 96, then IEEE also invokes the particular
71 XFmode `long double' data structure used by the Intel 80x86 series
72 processors.
74 `DEC' refers specifically to the Digital Equipment Corp PDP-11
75 and VAX floating point data structure. This model currently
76 supports no type wider than DFmode.
78 `IBM' refers specifically to the IBM System/370 and compatible
79 floating point data structure. This model currently supports
80 no type wider than DFmode. The IBM conversions were contributed by
81 frank@atom.ansto.gov.au (Frank Crawford).
83 If LONG_DOUBLE_TYPE_SIZE = 64 (the default, unless tm.h defines it)
84 then `long double' and `double' are both implemented, but they
85 both mean DFmode. In this case, the software floating-point
86 support available here is activated by writing
87 #define REAL_ARITHMETIC
88 in tm.h.
90 The case LONG_DOUBLE_TYPE_SIZE = 128 activates TFmode support
91 and may deactivate XFmode since `long double' is used to refer
92 to both modes.
94 The macros FLOAT_WORDS_BIG_ENDIAN, HOST_FLOAT_WORDS_BIG_ENDIAN,
95 contributed by Richard Earnshaw <Richard.Earnshaw@cl.cam.ac.uk>,
96 separate the floating point unit's endian-ness from that of
97 the integer addressing. This permits one to define a big-endian
98 FPU on a little-endian machine (e.g., ARM). An extension to
99 BYTES_BIG_ENDIAN may be required for some machines in the future.
100 These optional macros may be defined in tm.h. In real.h, they
101 default to WORDS_BIG_ENDIAN, etc., so there is no need to define
102 them for any normal host or target machine on which the floats
103 and the integers have the same endian-ness. */
106 /* The following converts gcc macros into the ones used by this file. */
108 /* REAL_ARITHMETIC defined means that macros in real.h are
109 defined to call emulator functions. */
110 #ifdef REAL_ARITHMETIC
112 #if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
113 /* PDP-11, Pro350, VAX: */
114 #define DEC 1
115 #else /* it's not VAX */
116 #if TARGET_FLOAT_FORMAT == IBM_FLOAT_FORMAT
117 /* IBM System/370 style */
118 #define IBM 1
119 #else /* it's also not an IBM */
120 #if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
121 #define IEEE
122 #else /* it's not IEEE either */
123 /* UNKnown arithmetic. We don't support this and can't go on. */
124 unknown arithmetic type
125 #define UNK 1
126 #endif /* not IEEE */
127 #endif /* not IBM */
128 #endif /* not VAX */
130 #define REAL_WORDS_BIG_ENDIAN FLOAT_WORDS_BIG_ENDIAN
132 #else
133 /* REAL_ARITHMETIC not defined means that the *host's* data
134 structure will be used. It may differ by endian-ness from the
135 target machine's structure and will get its ends swapped
136 accordingly (but not here). Probably only the decimal <-> binary
137 functions in this file will actually be used in this case. */
139 #if HOST_FLOAT_FORMAT == VAX_FLOAT_FORMAT
140 #define DEC 1
141 #else /* it's not VAX */
142 #if HOST_FLOAT_FORMAT == IBM_FLOAT_FORMAT
143 /* IBM System/370 style */
144 #define IBM 1
145 #else /* it's also not an IBM */
146 #if HOST_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
147 #define IEEE
148 #else /* it's not IEEE either */
149 unknown arithmetic type
150 #define UNK 1
151 #endif /* not IEEE */
152 #endif /* not IBM */
153 #endif /* not VAX */
155 #define REAL_WORDS_BIG_ENDIAN HOST_FLOAT_WORDS_BIG_ENDIAN
157 #endif /* REAL_ARITHMETIC not defined */
159 /* Define INFINITY for support of infinity.
160 Define NANS for support of Not-a-Number's (NaN's). */
161 #if !defined(DEC) && !defined(IBM)
162 #define INFINITY
163 #define NANS
164 #endif
166 /* Support of NaNs requires support of infinity. */
167 #ifdef NANS
168 #ifndef INFINITY
169 #define INFINITY
170 #endif
171 #endif
173 /* Find a host integer type that is at least 16 bits wide,
174 and another type at least twice whatever that size is. */
176 #if HOST_BITS_PER_CHAR >= 16
177 #define EMUSHORT char
178 #define EMUSHORT_SIZE HOST_BITS_PER_CHAR
179 #define EMULONG_SIZE (2 * HOST_BITS_PER_CHAR)
180 #else
181 #if HOST_BITS_PER_SHORT >= 16
182 #define EMUSHORT short
183 #define EMUSHORT_SIZE HOST_BITS_PER_SHORT
184 #define EMULONG_SIZE (2 * HOST_BITS_PER_SHORT)
185 #else
186 #if HOST_BITS_PER_INT >= 16
187 #define EMUSHORT int
188 #define EMUSHORT_SIZE HOST_BITS_PER_INT
189 #define EMULONG_SIZE (2 * HOST_BITS_PER_INT)
190 #else
191 #if HOST_BITS_PER_LONG >= 16
192 #define EMUSHORT long
193 #define EMUSHORT_SIZE HOST_BITS_PER_LONG
194 #define EMULONG_SIZE (2 * HOST_BITS_PER_LONG)
195 #else
196 /* You will have to modify this program to have a smaller unit size. */
197 #define EMU_NON_COMPILE
198 #endif
199 #endif
200 #endif
201 #endif
203 #if HOST_BITS_PER_SHORT >= EMULONG_SIZE
204 #define EMULONG short
205 #else
206 #if HOST_BITS_PER_INT >= EMULONG_SIZE
207 #define EMULONG int
208 #else
209 #if HOST_BITS_PER_LONG >= EMULONG_SIZE
210 #define EMULONG long
211 #else
212 #if HOST_BITS_PER_LONG_LONG >= EMULONG_SIZE
213 #define EMULONG long long int
214 #else
215 /* You will have to modify this program to have a smaller unit size. */
216 #define EMU_NON_COMPILE
217 #endif
218 #endif
219 #endif
220 #endif
223 /* The host interface doesn't work if no 16-bit size exists. */
224 #if EMUSHORT_SIZE != 16
225 #define EMU_NON_COMPILE
226 #endif
228 /* OK to continue compilation. */
229 #ifndef EMU_NON_COMPILE
231 /* Construct macros to translate between REAL_VALUE_TYPE and e type.
232 In GET_REAL and PUT_REAL, r and e are pointers.
233 A REAL_VALUE_TYPE is guaranteed to occupy contiguous locations
234 in memory, with no holes. */
236 #if LONG_DOUBLE_TYPE_SIZE == 96
237 /* Number of 16 bit words in external e type format */
238 #define NE 6
239 #define MAXDECEXP 4932
240 #define MINDECEXP -4956
241 #define GET_REAL(r,e) bcopy ((char *) r, (char *) e, 2*NE)
242 #define PUT_REAL(e,r) bcopy ((char *) e, (char *) r, 2*NE)
243 #else /* no XFmode */
244 #if LONG_DOUBLE_TYPE_SIZE == 128
245 #define NE 10
246 #define MAXDECEXP 4932
247 #define MINDECEXP -4977
248 #define GET_REAL(r,e) bcopy ((char *) r, (char *) e, 2*NE)
249 #define PUT_REAL(e,r) bcopy ((char *) e, (char *) r, 2*NE)
250 #else
251 #define NE 6
252 #define MAXDECEXP 4932
253 #define MINDECEXP -4956
254 #ifdef REAL_ARITHMETIC
255 /* Emulator uses target format internally
256 but host stores it in host endian-ness. */
258 #define GET_REAL(r,e) \
259 do { \
260 if (HOST_FLOAT_WORDS_BIG_ENDIAN == REAL_WORDS_BIG_ENDIAN) \
261 e53toe ((unsigned EMUSHORT *) (r), (e)); \
262 else \
264 unsigned EMUSHORT w[4]; \
265 w[3] = ((EMUSHORT *) r)[0]; \
266 w[2] = ((EMUSHORT *) r)[1]; \
267 w[1] = ((EMUSHORT *) r)[2]; \
268 w[0] = ((EMUSHORT *) r)[3]; \
269 e53toe (w, (e)); \
271 } while (0)
273 #define PUT_REAL(e,r) \
274 do { \
275 if (HOST_FLOAT_WORDS_BIG_ENDIAN == REAL_WORDS_BIG_ENDIAN) \
276 etoe53 ((e), (unsigned EMUSHORT *) (r)); \
277 else \
279 unsigned EMUSHORT w[4]; \
280 etoe53 ((e), w); \
281 *((EMUSHORT *) r) = w[3]; \
282 *((EMUSHORT *) r + 1) = w[2]; \
283 *((EMUSHORT *) r + 2) = w[1]; \
284 *((EMUSHORT *) r + 3) = w[0]; \
286 } while (0)
288 #else /* not REAL_ARITHMETIC */
290 /* emulator uses host format */
291 #define GET_REAL(r,e) e53toe ((unsigned EMUSHORT *) (r), (e))
292 #define PUT_REAL(e,r) etoe53 ((e), (unsigned EMUSHORT *) (r))
294 #endif /* not REAL_ARITHMETIC */
295 #endif /* not TFmode */
296 #endif /* no XFmode */
299 /* Number of 16 bit words in internal format */
300 #define NI (NE+3)
302 /* Array offset to exponent */
303 #define E 1
305 /* Array offset to high guard word */
306 #define M 2
308 /* Number of bits of precision */
309 #define NBITS ((NI-4)*16)
311 /* Maximum number of decimal digits in ASCII conversion
312 * = NBITS*log10(2)
314 #define NDEC (NBITS*8/27)
316 /* The exponent of 1.0 */
317 #define EXONE (0x3fff)
319 extern int extra_warnings;
320 extern unsigned EMUSHORT ezero[], ehalf[], eone[], etwo[];
321 extern unsigned EMUSHORT elog2[], esqrt2[];
323 static void endian PROTO((unsigned EMUSHORT *, long *,
324 enum machine_mode));
325 static void eclear PROTO((unsigned EMUSHORT *));
326 static void emov PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
327 static void eabs PROTO((unsigned EMUSHORT *));
328 static void eneg PROTO((unsigned EMUSHORT *));
329 static int eisneg PROTO((unsigned EMUSHORT *));
330 static int eisinf PROTO((unsigned EMUSHORT *));
331 static int eisnan PROTO((unsigned EMUSHORT *));
332 static void einfin PROTO((unsigned EMUSHORT *));
333 static void enan PROTO((unsigned EMUSHORT *, int));
334 static void emovi PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
335 static void emovo PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
336 static void ecleaz PROTO((unsigned EMUSHORT *));
337 static void ecleazs PROTO((unsigned EMUSHORT *));
338 static void emovz PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
339 static void einan PROTO((unsigned EMUSHORT *));
340 static int eiisnan PROTO((unsigned EMUSHORT *));
341 static int eiisneg PROTO((unsigned EMUSHORT *));
342 static void eiinfin PROTO((unsigned EMUSHORT *));
343 static int eiisinf PROTO((unsigned EMUSHORT *));
344 static int ecmpm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
345 static void eshdn1 PROTO((unsigned EMUSHORT *));
346 static void eshup1 PROTO((unsigned EMUSHORT *));
347 static void eshdn8 PROTO((unsigned EMUSHORT *));
348 static void eshup8 PROTO((unsigned EMUSHORT *));
349 static void eshup6 PROTO((unsigned EMUSHORT *));
350 static void eshdn6 PROTO((unsigned EMUSHORT *));
351 static void eaddm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));\f
352 static void esubm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
353 static void m16m PROTO((unsigned int, unsigned short *,
354 unsigned short *));
355 static int edivm PROTO((unsigned short *, unsigned short *));
356 static int emulm PROTO((unsigned short *, unsigned short *));
357 static void emdnorm PROTO((unsigned EMUSHORT *, int, int, EMULONG, int));
358 static void esub PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
359 unsigned EMUSHORT *));
360 static void eadd PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
361 unsigned EMUSHORT *));
362 static void eadd1 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
363 unsigned EMUSHORT *));
364 static void ediv PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
365 unsigned EMUSHORT *));
366 static void emul PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
367 unsigned EMUSHORT *));
368 static void e53toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
369 static void e64toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
370 static void e113toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
371 static void e24toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
372 static void etoe113 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
373 static void toe113 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
374 static void etoe64 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
375 static void toe64 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
376 static void etoe53 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
377 static void toe53 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
378 static void etoe24 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
379 static void toe24 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
380 static int ecmp PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
381 static void eround PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
382 static void ltoe PROTO((HOST_WIDE_INT *, unsigned EMUSHORT *));
383 static void ultoe PROTO((unsigned HOST_WIDE_INT *, unsigned EMUSHORT *));
384 static void eifrac PROTO((unsigned EMUSHORT *, HOST_WIDE_INT *,
385 unsigned EMUSHORT *));
386 static void euifrac PROTO((unsigned EMUSHORT *, unsigned HOST_WIDE_INT *,
387 unsigned EMUSHORT *));
388 static int eshift PROTO((unsigned EMUSHORT *, int));
389 static int enormlz PROTO((unsigned EMUSHORT *));
390 static void e24toasc PROTO((unsigned EMUSHORT *, char *, int));
391 static void e53toasc PROTO((unsigned EMUSHORT *, char *, int));
392 static void e64toasc PROTO((unsigned EMUSHORT *, char *, int));
393 static void e113toasc PROTO((unsigned EMUSHORT *, char *, int));
394 static void etoasc PROTO((unsigned EMUSHORT *, char *, int));
395 static void asctoe24 PROTO((char *, unsigned EMUSHORT *));
396 static void asctoe53 PROTO((char *, unsigned EMUSHORT *));
397 static void asctoe64 PROTO((char *, unsigned EMUSHORT *));
398 static void asctoe113 PROTO((char *, unsigned EMUSHORT *));
399 static void asctoe PROTO((char *, unsigned EMUSHORT *));
400 static void asctoeg PROTO((char *, unsigned EMUSHORT *, int));
401 static void efloor PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
402 static void efrexp PROTO((unsigned EMUSHORT *, int *,
403 unsigned EMUSHORT *));
404 static void eldexp PROTO((unsigned EMUSHORT *, int, unsigned EMUSHORT *));
405 static void eremain PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
406 unsigned EMUSHORT *));
407 static void eiremain PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
408 static void mtherr PROTO((char *, int));
409 static void dectoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
410 static void etodec PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
411 static void todec PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
412 static void ibmtoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
413 enum machine_mode));
414 static void etoibm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
415 enum machine_mode));
416 static void toibm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
417 enum machine_mode));
418 static void make_nan PROTO((unsigned EMUSHORT *, int, enum machine_mode));
419 static void uditoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
420 static void ditoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
421 static void etoudi PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
422 static void etodi PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
423 static void esqrt PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
425 /* Copy 32-bit numbers obtained from array containing 16-bit numbers,
426 swapping ends if required, into output array of longs. The
427 result is normally passed to fprintf by the ASM_OUTPUT_ macros. */
429 static void
430 endian (e, x, mode)
431 unsigned EMUSHORT e[];
432 long x[];
433 enum machine_mode mode;
435 unsigned long th, t;
437 if (REAL_WORDS_BIG_ENDIAN)
439 switch (mode)
442 case TFmode:
443 /* Swap halfwords in the fourth long. */
444 th = (unsigned long) e[6] & 0xffff;
445 t = (unsigned long) e[7] & 0xffff;
446 t |= th << 16;
447 x[3] = (long) t;
449 case XFmode:
451 /* Swap halfwords in the third long. */
452 th = (unsigned long) e[4] & 0xffff;
453 t = (unsigned long) e[5] & 0xffff;
454 t |= th << 16;
455 x[2] = (long) t;
456 /* fall into the double case */
458 case DFmode:
460 /* swap halfwords in the second word */
461 th = (unsigned long) e[2] & 0xffff;
462 t = (unsigned long) e[3] & 0xffff;
463 t |= th << 16;
464 x[1] = (long) t;
465 /* fall into the float case */
467 case HFmode:
468 case SFmode:
470 /* swap halfwords in the first word */
471 th = (unsigned long) e[0] & 0xffff;
472 t = (unsigned long) e[1] & 0xffff;
473 t |= th << 16;
474 x[0] = (long) t;
475 break;
477 default:
478 abort ();
481 else
483 /* Pack the output array without swapping. */
485 switch (mode)
488 case TFmode:
490 /* Pack the fourth long. */
491 th = (unsigned long) e[7] & 0xffff;
492 t = (unsigned long) e[6] & 0xffff;
493 t |= th << 16;
494 x[3] = (long) t;
496 case XFmode:
498 /* Pack the third long.
499 Each element of the input REAL_VALUE_TYPE array has 16 useful bits
500 in it. */
501 th = (unsigned long) e[5] & 0xffff;
502 t = (unsigned long) e[4] & 0xffff;
503 t |= th << 16;
504 x[2] = (long) t;
505 /* fall into the double case */
507 case DFmode:
509 /* pack the second long */
510 th = (unsigned long) e[3] & 0xffff;
511 t = (unsigned long) e[2] & 0xffff;
512 t |= th << 16;
513 x[1] = (long) t;
514 /* fall into the float case */
516 case HFmode:
517 case SFmode:
519 /* pack the first long */
520 th = (unsigned long) e[1] & 0xffff;
521 t = (unsigned long) e[0] & 0xffff;
522 t |= th << 16;
523 x[0] = (long) t;
524 break;
526 default:
527 abort ();
533 /* This is the implementation of the REAL_ARITHMETIC macro. */
535 void
536 earith (value, icode, r1, r2)
537 REAL_VALUE_TYPE *value;
538 int icode;
539 REAL_VALUE_TYPE *r1;
540 REAL_VALUE_TYPE *r2;
542 unsigned EMUSHORT d1[NE], d2[NE], v[NE];
543 enum tree_code code;
545 GET_REAL (r1, d1);
546 GET_REAL (r2, d2);
547 #ifdef NANS
548 /* Return NaN input back to the caller. */
549 if (eisnan (d1))
551 PUT_REAL (d1, value);
552 return;
554 if (eisnan (d2))
556 PUT_REAL (d2, value);
557 return;
559 #endif
560 code = (enum tree_code) icode;
561 switch (code)
563 case PLUS_EXPR:
564 eadd (d2, d1, v);
565 break;
567 case MINUS_EXPR:
568 esub (d2, d1, v); /* d1 - d2 */
569 break;
571 case MULT_EXPR:
572 emul (d2, d1, v);
573 break;
575 case RDIV_EXPR:
576 #ifndef REAL_INFINITY
577 if (ecmp (d2, ezero) == 0)
579 #ifdef NANS
580 enan (v, eisneg (d1) ^ eisneg (d2));
581 break;
582 #else
583 abort ();
584 #endif
586 #endif
587 ediv (d2, d1, v); /* d1/d2 */
588 break;
590 case MIN_EXPR: /* min (d1,d2) */
591 if (ecmp (d1, d2) < 0)
592 emov (d1, v);
593 else
594 emov (d2, v);
595 break;
597 case MAX_EXPR: /* max (d1,d2) */
598 if (ecmp (d1, d2) > 0)
599 emov (d1, v);
600 else
601 emov (d2, v);
602 break;
603 default:
604 emov (ezero, v);
605 break;
607 PUT_REAL (v, value);
611 /* Truncate REAL_VALUE_TYPE toward zero to signed HOST_WIDE_INT.
612 implements REAL_VALUE_RNDZINT (x) (etrunci (x)). */
614 REAL_VALUE_TYPE
615 etrunci (x)
616 REAL_VALUE_TYPE x;
618 unsigned EMUSHORT f[NE], g[NE];
619 REAL_VALUE_TYPE r;
620 HOST_WIDE_INT l;
622 GET_REAL (&x, g);
623 #ifdef NANS
624 if (eisnan (g))
625 return (x);
626 #endif
627 eifrac (g, &l, f);
628 ltoe (&l, g);
629 PUT_REAL (g, &r);
630 return (r);
634 /* Truncate REAL_VALUE_TYPE toward zero to unsigned HOST_WIDE_INT;
635 implements REAL_VALUE_UNSIGNED_RNDZINT (x) (etruncui (x)). */
637 REAL_VALUE_TYPE
638 etruncui (x)
639 REAL_VALUE_TYPE x;
641 unsigned EMUSHORT f[NE], g[NE];
642 REAL_VALUE_TYPE r;
643 unsigned HOST_WIDE_INT l;
645 GET_REAL (&x, g);
646 #ifdef NANS
647 if (eisnan (g))
648 return (x);
649 #endif
650 euifrac (g, &l, f);
651 ultoe (&l, g);
652 PUT_REAL (g, &r);
653 return (r);
657 /* This is the REAL_VALUE_ATOF function. It converts a decimal string to
658 binary, rounding off as indicated by the machine_mode argument. Then it
659 promotes the rounded value to REAL_VALUE_TYPE. */
661 REAL_VALUE_TYPE
662 ereal_atof (s, t)
663 char *s;
664 enum machine_mode t;
666 unsigned EMUSHORT tem[NE], e[NE];
667 REAL_VALUE_TYPE r;
669 switch (t)
671 case HFmode:
672 case SFmode:
673 asctoe24 (s, tem);
674 e24toe (tem, e);
675 break;
676 case DFmode:
677 asctoe53 (s, tem);
678 e53toe (tem, e);
679 break;
680 case XFmode:
681 asctoe64 (s, tem);
682 e64toe (tem, e);
683 break;
684 case TFmode:
685 asctoe113 (s, tem);
686 e113toe (tem, e);
687 break;
688 default:
689 asctoe (s, e);
691 PUT_REAL (e, &r);
692 return (r);
696 /* Expansion of REAL_NEGATE. */
698 REAL_VALUE_TYPE
699 ereal_negate (x)
700 REAL_VALUE_TYPE x;
702 unsigned EMUSHORT e[NE];
703 REAL_VALUE_TYPE r;
705 GET_REAL (&x, e);
706 eneg (e);
707 PUT_REAL (e, &r);
708 return (r);
712 /* Round real toward zero to HOST_WIDE_INT;
713 implements REAL_VALUE_FIX (x). */
715 HOST_WIDE_INT
716 efixi (x)
717 REAL_VALUE_TYPE x;
719 unsigned EMUSHORT f[NE], g[NE];
720 HOST_WIDE_INT l;
722 GET_REAL (&x, f);
723 #ifdef NANS
724 if (eisnan (f))
726 warning ("conversion from NaN to int");
727 return (-1);
729 #endif
730 eifrac (f, &l, g);
731 return l;
734 /* Round real toward zero to unsigned HOST_WIDE_INT
735 implements REAL_VALUE_UNSIGNED_FIX (x).
736 Negative input returns zero. */
738 unsigned HOST_WIDE_INT
739 efixui (x)
740 REAL_VALUE_TYPE x;
742 unsigned EMUSHORT f[NE], g[NE];
743 unsigned HOST_WIDE_INT l;
745 GET_REAL (&x, f);
746 #ifdef NANS
747 if (eisnan (f))
749 warning ("conversion from NaN to unsigned int");
750 return (-1);
752 #endif
753 euifrac (f, &l, g);
754 return l;
758 /* REAL_VALUE_FROM_INT macro. */
760 void
761 ereal_from_int (d, i, j, mode)
762 REAL_VALUE_TYPE *d;
763 HOST_WIDE_INT i, j;
764 enum machine_mode mode;
766 unsigned EMUSHORT df[NE], dg[NE];
767 HOST_WIDE_INT low, high;
768 int sign;
770 if (GET_MODE_CLASS (mode) != MODE_FLOAT)
771 abort ();
772 sign = 0;
773 low = i;
774 if ((high = j) < 0)
776 sign = 1;
777 /* complement and add 1 */
778 high = ~high;
779 if (low)
780 low = -low;
781 else
782 high += 1;
784 eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
785 ultoe ((unsigned HOST_WIDE_INT *) &high, dg);
786 emul (dg, df, dg);
787 ultoe ((unsigned HOST_WIDE_INT *) &low, df);
788 eadd (df, dg, dg);
789 if (sign)
790 eneg (dg);
792 /* A REAL_VALUE_TYPE may not be wide enough to hold the two HOST_WIDE_INTS.
793 Avoid double-rounding errors later by rounding off now from the
794 extra-wide internal format to the requested precision. */
795 switch (GET_MODE_BITSIZE (mode))
797 case 32:
798 etoe24 (dg, df);
799 e24toe (df, dg);
800 break;
802 case 64:
803 etoe53 (dg, df);
804 e53toe (df, dg);
805 break;
807 case 96:
808 etoe64 (dg, df);
809 e64toe (df, dg);
810 break;
812 case 128:
813 etoe113 (dg, df);
814 e113toe (df, dg);
815 break;
817 default:
818 abort ();
821 PUT_REAL (dg, d);
825 /* REAL_VALUE_FROM_UNSIGNED_INT macro. */
827 void
828 ereal_from_uint (d, i, j, mode)
829 REAL_VALUE_TYPE *d;
830 unsigned HOST_WIDE_INT i, j;
831 enum machine_mode mode;
833 unsigned EMUSHORT df[NE], dg[NE];
834 unsigned HOST_WIDE_INT low, high;
836 if (GET_MODE_CLASS (mode) != MODE_FLOAT)
837 abort ();
838 low = i;
839 high = j;
840 eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
841 ultoe (&high, dg);
842 emul (dg, df, dg);
843 ultoe (&low, df);
844 eadd (df, dg, dg);
846 /* A REAL_VALUE_TYPE may not be wide enough to hold the two HOST_WIDE_INTS.
847 Avoid double-rounding errors later by rounding off now from the
848 extra-wide internal format to the requested precision. */
849 switch (GET_MODE_BITSIZE (mode))
851 case 32:
852 etoe24 (dg, df);
853 e24toe (df, dg);
854 break;
856 case 64:
857 etoe53 (dg, df);
858 e53toe (df, dg);
859 break;
861 case 96:
862 etoe64 (dg, df);
863 e64toe (df, dg);
864 break;
866 case 128:
867 etoe113 (dg, df);
868 e113toe (df, dg);
869 break;
871 default:
872 abort ();
875 PUT_REAL (dg, d);
879 /* REAL_VALUE_TO_INT macro. */
881 void
882 ereal_to_int (low, high, rr)
883 HOST_WIDE_INT *low, *high;
884 REAL_VALUE_TYPE rr;
886 unsigned EMUSHORT d[NE], df[NE], dg[NE], dh[NE];
887 int s;
889 GET_REAL (&rr, d);
890 #ifdef NANS
891 if (eisnan (d))
893 warning ("conversion from NaN to int");
894 *low = -1;
895 *high = -1;
896 return;
898 #endif
899 /* convert positive value */
900 s = 0;
901 if (eisneg (d))
903 eneg (d);
904 s = 1;
906 eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
907 ediv (df, d, dg); /* dg = d / 2^32 is the high word */
908 euifrac (dg, (unsigned HOST_WIDE_INT *) high, dh);
909 emul (df, dh, dg); /* fractional part is the low word */
910 euifrac (dg, (unsigned HOST_WIDE_INT *)low, dh);
911 if (s)
913 /* complement and add 1 */
914 *high = ~(*high);
915 if (*low)
916 *low = -(*low);
917 else
918 *high += 1;
923 /* REAL_VALUE_LDEXP macro. */
925 REAL_VALUE_TYPE
926 ereal_ldexp (x, n)
927 REAL_VALUE_TYPE x;
928 int n;
930 unsigned EMUSHORT e[NE], y[NE];
931 REAL_VALUE_TYPE r;
933 GET_REAL (&x, e);
934 #ifdef NANS
935 if (eisnan (e))
936 return (x);
937 #endif
938 eldexp (e, n, y);
939 PUT_REAL (y, &r);
940 return (r);
943 /* These routines are conditionally compiled because functions
944 of the same names may be defined in fold-const.c. */
946 #ifdef REAL_ARITHMETIC
948 /* Check for infinity in a REAL_VALUE_TYPE. */
951 target_isinf (x)
952 REAL_VALUE_TYPE x;
954 unsigned EMUSHORT e[NE];
956 #ifdef INFINITY
957 GET_REAL (&x, e);
958 return (eisinf (e));
959 #else
960 return 0;
961 #endif
964 /* Check whether a REAL_VALUE_TYPE item is a NaN. */
967 target_isnan (x)
968 REAL_VALUE_TYPE x;
970 unsigned EMUSHORT e[NE];
972 #ifdef NANS
973 GET_REAL (&x, e);
974 return (eisnan (e));
975 #else
976 return (0);
977 #endif
981 /* Check for a negative REAL_VALUE_TYPE number.
982 This just checks the sign bit, so that -0 counts as negative. */
985 target_negative (x)
986 REAL_VALUE_TYPE x;
988 return ereal_isneg (x);
991 /* Expansion of REAL_VALUE_TRUNCATE.
992 The result is in floating point, rounded to nearest or even. */
994 REAL_VALUE_TYPE
995 real_value_truncate (mode, arg)
996 enum machine_mode mode;
997 REAL_VALUE_TYPE arg;
999 unsigned EMUSHORT e[NE], t[NE];
1000 REAL_VALUE_TYPE r;
1002 GET_REAL (&arg, e);
1003 #ifdef NANS
1004 if (eisnan (e))
1005 return (arg);
1006 #endif
1007 eclear (t);
1008 switch (mode)
1010 case TFmode:
1011 etoe113 (e, t);
1012 e113toe (t, t);
1013 break;
1015 case XFmode:
1016 etoe64 (e, t);
1017 e64toe (t, t);
1018 break;
1020 case DFmode:
1021 etoe53 (e, t);
1022 e53toe (t, t);
1023 break;
1025 case HFmode:
1026 case SFmode:
1027 etoe24 (e, t);
1028 e24toe (t, t);
1029 break;
1031 case SImode:
1032 r = etrunci (arg);
1033 return (r);
1035 /* If an unsupported type was requested, presume that
1036 the machine files know something useful to do with
1037 the unmodified value. */
1039 default:
1040 return (arg);
1042 PUT_REAL (t, &r);
1043 return (r);
1046 /* Try to change R into its exact multiplicative inverse in machine mode
1047 MODE. Return nonzero function value if successful. */
1050 exact_real_inverse (mode, r)
1051 enum machine_mode mode;
1052 REAL_VALUE_TYPE *r;
1054 unsigned EMUSHORT e[NE], einv[NE];
1055 REAL_VALUE_TYPE rinv;
1056 int i;
1058 GET_REAL (r, e);
1060 /* Test for input in range. Don't transform IEEE special values. */
1061 if (eisinf (e) || eisnan (e) || (ecmp (e, ezero) == 0))
1062 return 0;
1064 /* Test for a power of 2: all significand bits zero except the MSB.
1065 We are assuming the target has binary (or hex) arithmetic. */
1066 if (e[NE - 2] != 0x8000)
1067 return 0;
1069 for (i = 0; i < NE - 2; i++)
1071 if (e[i] != 0)
1072 return 0;
1075 /* Compute the inverse and truncate it to the required mode. */
1076 ediv (e, eone, einv);
1077 PUT_REAL (einv, &rinv);
1078 rinv = real_value_truncate (mode, rinv);
1080 #ifdef CHECK_FLOAT_VALUE
1081 /* This check is not redundant. It may, for example, flush
1082 a supposedly IEEE denormal value to zero. */
1083 i = 0;
1084 if (CHECK_FLOAT_VALUE (mode, rinv, i))
1085 return 0;
1086 #endif
1087 GET_REAL (&rinv, einv);
1089 /* Check the bits again, because the truncation might have
1090 generated an arbitrary saturation value on overflow. */
1091 if (einv[NE - 2] != 0x8000)
1092 return 0;
1094 for (i = 0; i < NE - 2; i++)
1096 if (einv[i] != 0)
1097 return 0;
1100 /* Fail if the computed inverse is out of range. */
1101 if (eisinf (einv) || eisnan (einv) || (ecmp (einv, ezero) == 0))
1102 return 0;
1104 /* Output the reciprocal and return success flag. */
1105 PUT_REAL (einv, r);
1106 return 1;
1108 #endif /* REAL_ARITHMETIC defined */
1110 /* Used for debugging--print the value of R in human-readable format
1111 on stderr. */
1113 void
1114 debug_real (r)
1115 REAL_VALUE_TYPE r;
1117 char dstr[30];
1119 REAL_VALUE_TO_DECIMAL (r, "%.20g", dstr);
1120 fprintf (stderr, "%s", dstr);
1124 /* The following routines convert REAL_VALUE_TYPE to the various floating
1125 point formats that are meaningful to supported computers.
1127 The results are returned in 32-bit pieces, each piece stored in a `long'.
1128 This is so they can be printed by statements like
1130 fprintf (file, "%lx, %lx", L[0], L[1]);
1132 that will work on both narrow- and wide-word host computers. */
1134 /* Convert R to a 128-bit long double precision value. The output array L
1135 contains four 32-bit pieces of the result, in the order they would appear
1136 in memory. */
1138 void
1139 etartdouble (r, l)
1140 REAL_VALUE_TYPE r;
1141 long l[];
1143 unsigned EMUSHORT e[NE];
1145 GET_REAL (&r, e);
1146 etoe113 (e, e);
1147 endian (e, l, TFmode);
1150 /* Convert R to a double extended precision value. The output array L
1151 contains three 32-bit pieces of the result, in the order they would
1152 appear in memory. */
1154 void
1155 etarldouble (r, l)
1156 REAL_VALUE_TYPE r;
1157 long l[];
1159 unsigned EMUSHORT e[NE];
1161 GET_REAL (&r, e);
1162 etoe64 (e, e);
1163 endian (e, l, XFmode);
1166 /* Convert R to a double precision value. The output array L contains two
1167 32-bit pieces of the result, in the order they would appear in memory. */
1169 void
1170 etardouble (r, l)
1171 REAL_VALUE_TYPE r;
1172 long l[];
1174 unsigned EMUSHORT e[NE];
1176 GET_REAL (&r, e);
1177 etoe53 (e, e);
1178 endian (e, l, DFmode);
1181 /* Convert R to a single precision float value stored in the least-significant
1182 bits of a `long'. */
1184 long
1185 etarsingle (r)
1186 REAL_VALUE_TYPE r;
1188 unsigned EMUSHORT e[NE];
1189 long l;
1191 GET_REAL (&r, e);
1192 etoe24 (e, e);
1193 endian (e, &l, SFmode);
1194 return ((long) l);
1197 /* Convert X to a decimal ASCII string S for output to an assembly
1198 language file. Note, there is no standard way to spell infinity or
1199 a NaN, so these values may require special treatment in the tm.h
1200 macros. */
1202 void
1203 ereal_to_decimal (x, s)
1204 REAL_VALUE_TYPE x;
1205 char *s;
1207 unsigned EMUSHORT e[NE];
1209 GET_REAL (&x, e);
1210 etoasc (e, s, 20);
1213 /* Compare X and Y. Return 1 if X > Y, 0 if X == Y, -1 if X < Y,
1214 or -2 if either is a NaN. */
1217 ereal_cmp (x, y)
1218 REAL_VALUE_TYPE x, y;
1220 unsigned EMUSHORT ex[NE], ey[NE];
1222 GET_REAL (&x, ex);
1223 GET_REAL (&y, ey);
1224 return (ecmp (ex, ey));
1227 /* Return 1 if the sign bit of X is set, else return 0. */
1230 ereal_isneg (x)
1231 REAL_VALUE_TYPE x;
1233 unsigned EMUSHORT ex[NE];
1235 GET_REAL (&x, ex);
1236 return (eisneg (ex));
1239 /* End of REAL_ARITHMETIC interface */
1242 Extended precision IEEE binary floating point arithmetic routines
1244 Numbers are stored in C language as arrays of 16-bit unsigned
1245 short integers. The arguments of the routines are pointers to
1246 the arrays.
1248 External e type data structure, similar to Intel 8087 chip
1249 temporary real format but possibly with a larger significand:
1251 NE-1 significand words (least significant word first,
1252 most significant bit is normally set)
1253 exponent (value = EXONE for 1.0,
1254 top bit is the sign)
1257 Internal exploded e-type data structure of a number (a "word" is 16 bits):
1259 ei[0] sign word (0 for positive, 0xffff for negative)
1260 ei[1] biased exponent (value = EXONE for the number 1.0)
1261 ei[2] high guard word (always zero after normalization)
1262 ei[3]
1263 to ei[NI-2] significand (NI-4 significand words,
1264 most significant word first,
1265 most significant bit is set)
1266 ei[NI-1] low guard word (0x8000 bit is rounding place)
1270 Routines for external format e-type numbers
1272 asctoe (string, e) ASCII string to extended double e type
1273 asctoe64 (string, &d) ASCII string to long double
1274 asctoe53 (string, &d) ASCII string to double
1275 asctoe24 (string, &f) ASCII string to single
1276 asctoeg (string, e, prec) ASCII string to specified precision
1277 e24toe (&f, e) IEEE single precision to e type
1278 e53toe (&d, e) IEEE double precision to e type
1279 e64toe (&d, e) IEEE long double precision to e type
1280 e113toe (&d, e) 128-bit long double precision to e type
1281 eabs (e) absolute value
1282 eadd (a, b, c) c = b + a
1283 eclear (e) e = 0
1284 ecmp (a, b) Returns 1 if a > b, 0 if a == b,
1285 -1 if a < b, -2 if either a or b is a NaN.
1286 ediv (a, b, c) c = b / a
1287 efloor (a, b) truncate to integer, toward -infinity
1288 efrexp (a, exp, s) extract exponent and significand
1289 eifrac (e, &l, frac) e to HOST_WIDE_INT and e type fraction
1290 euifrac (e, &l, frac) e to unsigned HOST_WIDE_INT and e type fraction
1291 einfin (e) set e to infinity, leaving its sign alone
1292 eldexp (a, n, b) multiply by 2**n
1293 emov (a, b) b = a
1294 emul (a, b, c) c = b * a
1295 eneg (e) e = -e
1296 eround (a, b) b = nearest integer value to a
1297 esub (a, b, c) c = b - a
1298 e24toasc (&f, str, n) single to ASCII string, n digits after decimal
1299 e53toasc (&d, str, n) double to ASCII string, n digits after decimal
1300 e64toasc (&d, str, n) 80-bit long double to ASCII string
1301 e113toasc (&d, str, n) 128-bit long double to ASCII string
1302 etoasc (e, str, n) e to ASCII string, n digits after decimal
1303 etoe24 (e, &f) convert e type to IEEE single precision
1304 etoe53 (e, &d) convert e type to IEEE double precision
1305 etoe64 (e, &d) convert e type to IEEE long double precision
1306 ltoe (&l, e) HOST_WIDE_INT to e type
1307 ultoe (&l, e) unsigned HOST_WIDE_INT to e type
1308 eisneg (e) 1 if sign bit of e != 0, else 0
1309 eisinf (e) 1 if e has maximum exponent (non-IEEE)
1310 or is infinite (IEEE)
1311 eisnan (e) 1 if e is a NaN
1314 Routines for internal format exploded e-type numbers
1316 eaddm (ai, bi) add significands, bi = bi + ai
1317 ecleaz (ei) ei = 0
1318 ecleazs (ei) set ei = 0 but leave its sign alone
1319 ecmpm (ai, bi) compare significands, return 1, 0, or -1
1320 edivm (ai, bi) divide significands, bi = bi / ai
1321 emdnorm (ai,l,s,exp) normalize and round off
1322 emovi (a, ai) convert external a to internal ai
1323 emovo (ai, a) convert internal ai to external a
1324 emovz (ai, bi) bi = ai, low guard word of bi = 0
1325 emulm (ai, bi) multiply significands, bi = bi * ai
1326 enormlz (ei) left-justify the significand
1327 eshdn1 (ai) shift significand and guards down 1 bit
1328 eshdn8 (ai) shift down 8 bits
1329 eshdn6 (ai) shift down 16 bits
1330 eshift (ai, n) shift ai n bits up (or down if n < 0)
1331 eshup1 (ai) shift significand and guards up 1 bit
1332 eshup8 (ai) shift up 8 bits
1333 eshup6 (ai) shift up 16 bits
1334 esubm (ai, bi) subtract significands, bi = bi - ai
1335 eiisinf (ai) 1 if infinite
1336 eiisnan (ai) 1 if a NaN
1337 eiisneg (ai) 1 if sign bit of ai != 0, else 0
1338 einan (ai) set ai = NaN
1339 eiinfin (ai) set ai = infinity
1341 The result is always normalized and rounded to NI-4 word precision
1342 after each arithmetic operation.
1344 Exception flags are NOT fully supported.
1346 Signaling NaN's are NOT supported; they are treated the same
1347 as quiet NaN's.
1349 Define INFINITY for support of infinity; otherwise a
1350 saturation arithmetic is implemented.
1352 Define NANS for support of Not-a-Number items; otherwise the
1353 arithmetic will never produce a NaN output, and might be confused
1354 by a NaN input.
1355 If NaN's are supported, the output of `ecmp (a,b)' is -2 if
1356 either a or b is a NaN. This means asking `if (ecmp (a,b) < 0)'
1357 may not be legitimate. Use `if (ecmp (a,b) == -1)' for `less than'
1358 if in doubt.
1360 Denormals are always supported here where appropriate (e.g., not
1361 for conversion to DEC numbers). */
1363 /* Definitions for error codes that are passed to the common error handling
1364 routine mtherr.
1366 For Digital Equipment PDP-11 and VAX computers, certain
1367 IBM systems, and others that use numbers with a 56-bit
1368 significand, the symbol DEC should be defined. In this
1369 mode, most floating point constants are given as arrays
1370 of octal integers to eliminate decimal to binary conversion
1371 errors that might be introduced by the compiler.
1373 For computers, such as IBM PC, that follow the IEEE
1374 Standard for Binary Floating Point Arithmetic (ANSI/IEEE
1375 Std 754-1985), the symbol IEEE should be defined.
1376 These numbers have 53-bit significands. In this mode, constants
1377 are provided as arrays of hexadecimal 16 bit integers.
1378 The endian-ness of generated values is controlled by
1379 REAL_WORDS_BIG_ENDIAN.
1381 To accommodate other types of computer arithmetic, all
1382 constants are also provided in a normal decimal radix
1383 which one can hope are correctly converted to a suitable
1384 format by the available C language compiler. To invoke
1385 this mode, the symbol UNK is defined.
1387 An important difference among these modes is a predefined
1388 set of machine arithmetic constants for each. The numbers
1389 MACHEP (the machine roundoff error), MAXNUM (largest number
1390 represented), and several other parameters are preset by
1391 the configuration symbol. Check the file const.c to
1392 ensure that these values are correct for your computer.
1394 For ANSI C compatibility, define ANSIC equal to 1. Currently
1395 this affects only the atan2 function and others that use it. */
1397 /* Constant definitions for math error conditions. */
1399 #define DOMAIN 1 /* argument domain error */
1400 #define SING 2 /* argument singularity */
1401 #define OVERFLOW 3 /* overflow range error */
1402 #define UNDERFLOW 4 /* underflow range error */
1403 #define TLOSS 5 /* total loss of precision */
1404 #define PLOSS 6 /* partial loss of precision */
1405 #define INVALID 7 /* NaN-producing operation */
1407 /* e type constants used by high precision check routines */
1409 #if LONG_DOUBLE_TYPE_SIZE == 128
1410 /* 0.0 */
1411 unsigned EMUSHORT ezero[NE] =
1412 {0x0000, 0x0000, 0x0000, 0x0000,
1413 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,};
1414 extern unsigned EMUSHORT ezero[];
1416 /* 5.0E-1 */
1417 unsigned EMUSHORT ehalf[NE] =
1418 {0x0000, 0x0000, 0x0000, 0x0000,
1419 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3ffe,};
1420 extern unsigned EMUSHORT ehalf[];
1422 /* 1.0E0 */
1423 unsigned EMUSHORT eone[NE] =
1424 {0x0000, 0x0000, 0x0000, 0x0000,
1425 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,};
1426 extern unsigned EMUSHORT eone[];
1428 /* 2.0E0 */
1429 unsigned EMUSHORT etwo[NE] =
1430 {0x0000, 0x0000, 0x0000, 0x0000,
1431 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4000,};
1432 extern unsigned EMUSHORT etwo[];
1434 /* 3.2E1 */
1435 unsigned EMUSHORT e32[NE] =
1436 {0x0000, 0x0000, 0x0000, 0x0000,
1437 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4004,};
1438 extern unsigned EMUSHORT e32[];
1440 /* 6.93147180559945309417232121458176568075500134360255E-1 */
1441 unsigned EMUSHORT elog2[NE] =
1442 {0x40f3, 0xf6af, 0x03f2, 0xb398,
1443 0xc9e3, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1444 extern unsigned EMUSHORT elog2[];
1446 /* 1.41421356237309504880168872420969807856967187537695E0 */
1447 unsigned EMUSHORT esqrt2[NE] =
1448 {0x1d6f, 0xbe9f, 0x754a, 0x89b3,
1449 0x597d, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1450 extern unsigned EMUSHORT esqrt2[];
1452 /* 3.14159265358979323846264338327950288419716939937511E0 */
1453 unsigned EMUSHORT epi[NE] =
1454 {0x2902, 0x1cd1, 0x80dc, 0x628b,
1455 0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1456 extern unsigned EMUSHORT epi[];
1458 #else
1459 /* LONG_DOUBLE_TYPE_SIZE is other than 128 */
1460 unsigned EMUSHORT ezero[NE] =
1461 {0, 0000000, 0000000, 0000000, 0000000, 0000000,};
1462 unsigned EMUSHORT ehalf[NE] =
1463 {0, 0000000, 0000000, 0000000, 0100000, 0x3ffe,};
1464 unsigned EMUSHORT eone[NE] =
1465 {0, 0000000, 0000000, 0000000, 0100000, 0x3fff,};
1466 unsigned EMUSHORT etwo[NE] =
1467 {0, 0000000, 0000000, 0000000, 0100000, 0040000,};
1468 unsigned EMUSHORT e32[NE] =
1469 {0, 0000000, 0000000, 0000000, 0100000, 0040004,};
1470 unsigned EMUSHORT elog2[NE] =
1471 {0xc9e4, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1472 unsigned EMUSHORT esqrt2[NE] =
1473 {0x597e, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1474 unsigned EMUSHORT epi[NE] =
1475 {0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1476 #endif
1478 /* Control register for rounding precision.
1479 This can be set to 113 (if NE=10), 80 (if NE=6), 64, 56, 53, or 24 bits. */
1481 int rndprc = NBITS;
1482 extern int rndprc;
1484 /* Clear out entire e-type number X. */
1486 static void
1487 eclear (x)
1488 register unsigned EMUSHORT *x;
1490 register int i;
1492 for (i = 0; i < NE; i++)
1493 *x++ = 0;
1496 /* Move e-type number from A to B. */
1498 static void
1499 emov (a, b)
1500 register unsigned EMUSHORT *a, *b;
1502 register int i;
1504 for (i = 0; i < NE; i++)
1505 *b++ = *a++;
1509 /* Absolute value of e-type X. */
1511 static void
1512 eabs (x)
1513 unsigned EMUSHORT x[];
1515 /* sign is top bit of last word of external format */
1516 x[NE - 1] &= 0x7fff;
1519 /* Negate the e-type number X. */
1521 static void
1522 eneg (x)
1523 unsigned EMUSHORT x[];
1526 x[NE - 1] ^= 0x8000; /* Toggle the sign bit */
1529 /* Return 1 if sign bit of e-type number X is nonzero, else zero. */
1531 static int
1532 eisneg (x)
1533 unsigned EMUSHORT x[];
1536 if (x[NE - 1] & 0x8000)
1537 return (1);
1538 else
1539 return (0);
1542 /* Return 1 if e-type number X is infinity, else return zero. */
1544 static int
1545 eisinf (x)
1546 unsigned EMUSHORT x[];
1549 #ifdef NANS
1550 if (eisnan (x))
1551 return (0);
1552 #endif
1553 if ((x[NE - 1] & 0x7fff) == 0x7fff)
1554 return (1);
1555 else
1556 return (0);
1559 /* Check if e-type number is not a number. The bit pattern is one that we
1560 defined, so we know for sure how to detect it. */
1562 static int
1563 eisnan (x)
1564 unsigned EMUSHORT x[];
1566 #ifdef NANS
1567 int i;
1569 /* NaN has maximum exponent */
1570 if ((x[NE - 1] & 0x7fff) != 0x7fff)
1571 return (0);
1572 /* ... and non-zero significand field. */
1573 for (i = 0; i < NE - 1; i++)
1575 if (*x++ != 0)
1576 return (1);
1578 #endif
1580 return (0);
1583 /* Fill e-type number X with infinity pattern (IEEE)
1584 or largest possible number (non-IEEE). */
1586 static void
1587 einfin (x)
1588 register unsigned EMUSHORT *x;
1590 register int i;
1592 #ifdef INFINITY
1593 for (i = 0; i < NE - 1; i++)
1594 *x++ = 0;
1595 *x |= 32767;
1596 #else
1597 for (i = 0; i < NE - 1; i++)
1598 *x++ = 0xffff;
1599 *x |= 32766;
1600 if (rndprc < NBITS)
1602 if (rndprc == 113)
1604 *(x - 9) = 0;
1605 *(x - 8) = 0;
1607 if (rndprc == 64)
1609 *(x - 5) = 0;
1611 if (rndprc == 53)
1613 *(x - 4) = 0xf800;
1615 else
1617 *(x - 4) = 0;
1618 *(x - 3) = 0;
1619 *(x - 2) = 0xff00;
1622 #endif
1625 /* Output an e-type NaN.
1626 This generates Intel's quiet NaN pattern for extended real.
1627 The exponent is 7fff, the leading mantissa word is c000. */
1629 static void
1630 enan (x, sign)
1631 register unsigned EMUSHORT *x;
1632 int sign;
1634 register int i;
1636 for (i = 0; i < NE - 2; i++)
1637 *x++ = 0;
1638 *x++ = 0xc000;
1639 *x = (sign << 15) | 0x7fff;
1642 /* Move in an e-type number A, converting it to exploded e-type B. */
1644 static void
1645 emovi (a, b)
1646 unsigned EMUSHORT *a, *b;
1648 register unsigned EMUSHORT *p, *q;
1649 int i;
1651 q = b;
1652 p = a + (NE - 1); /* point to last word of external number */
1653 /* get the sign bit */
1654 if (*p & 0x8000)
1655 *q++ = 0xffff;
1656 else
1657 *q++ = 0;
1658 /* get the exponent */
1659 *q = *p--;
1660 *q++ &= 0x7fff; /* delete the sign bit */
1661 #ifdef INFINITY
1662 if ((*(q - 1) & 0x7fff) == 0x7fff)
1664 #ifdef NANS
1665 if (eisnan (a))
1667 *q++ = 0;
1668 for (i = 3; i < NI; i++)
1669 *q++ = *p--;
1670 return;
1672 #endif
1674 for (i = 2; i < NI; i++)
1675 *q++ = 0;
1676 return;
1678 #endif
1680 /* clear high guard word */
1681 *q++ = 0;
1682 /* move in the significand */
1683 for (i = 0; i < NE - 1; i++)
1684 *q++ = *p--;
1685 /* clear low guard word */
1686 *q = 0;
1689 /* Move out exploded e-type number A, converting it to e type B. */
1691 static void
1692 emovo (a, b)
1693 unsigned EMUSHORT *a, *b;
1695 register unsigned EMUSHORT *p, *q;
1696 unsigned EMUSHORT i;
1697 int j;
1699 p = a;
1700 q = b + (NE - 1); /* point to output exponent */
1701 /* combine sign and exponent */
1702 i = *p++;
1703 if (i)
1704 *q-- = *p++ | 0x8000;
1705 else
1706 *q-- = *p++;
1707 #ifdef INFINITY
1708 if (*(p - 1) == 0x7fff)
1710 #ifdef NANS
1711 if (eiisnan (a))
1713 enan (b, eiisneg (a));
1714 return;
1716 #endif
1717 einfin (b);
1718 return;
1720 #endif
1721 /* skip over guard word */
1722 ++p;
1723 /* move the significand */
1724 for (j = 0; j < NE - 1; j++)
1725 *q-- = *p++;
1728 /* Clear out exploded e-type number XI. */
1730 static void
1731 ecleaz (xi)
1732 register unsigned EMUSHORT *xi;
1734 register int i;
1736 for (i = 0; i < NI; i++)
1737 *xi++ = 0;
1740 /* Clear out exploded e-type XI, but don't touch the sign. */
1742 static void
1743 ecleazs (xi)
1744 register unsigned EMUSHORT *xi;
1746 register int i;
1748 ++xi;
1749 for (i = 0; i < NI - 1; i++)
1750 *xi++ = 0;
1753 /* Move exploded e-type number from A to B. */
1755 static void
1756 emovz (a, b)
1757 register unsigned EMUSHORT *a, *b;
1759 register int i;
1761 for (i = 0; i < NI - 1; i++)
1762 *b++ = *a++;
1763 /* clear low guard word */
1764 *b = 0;
1767 /* Generate exploded e-type NaN.
1768 The explicit pattern for this is maximum exponent and
1769 top two significant bits set. */
1771 static void
1772 einan (x)
1773 unsigned EMUSHORT x[];
1776 ecleaz (x);
1777 x[E] = 0x7fff;
1778 x[M + 1] = 0xc000;
1781 /* Return nonzero if exploded e-type X is a NaN. */
1783 static int
1784 eiisnan (x)
1785 unsigned EMUSHORT x[];
1787 int i;
1789 if ((x[E] & 0x7fff) == 0x7fff)
1791 for (i = M + 1; i < NI; i++)
1793 if (x[i] != 0)
1794 return (1);
1797 return (0);
1800 /* Return nonzero if sign of exploded e-type X is nonzero. */
1802 static int
1803 eiisneg (x)
1804 unsigned EMUSHORT x[];
1807 return x[0] != 0;
1810 /* Fill exploded e-type X with infinity pattern.
1811 This has maximum exponent and significand all zeros. */
1813 static void
1814 eiinfin (x)
1815 unsigned EMUSHORT x[];
1818 ecleaz (x);
1819 x[E] = 0x7fff;
1822 /* Return nonzero if exploded e-type X is infinite. */
1824 static int
1825 eiisinf (x)
1826 unsigned EMUSHORT x[];
1829 #ifdef NANS
1830 if (eiisnan (x))
1831 return (0);
1832 #endif
1833 if ((x[E] & 0x7fff) == 0x7fff)
1834 return (1);
1835 return (0);
1839 /* Compare significands of numbers in internal exploded e-type format.
1840 Guard words are included in the comparison.
1842 Returns +1 if a > b
1843 0 if a == b
1844 -1 if a < b */
1846 static int
1847 ecmpm (a, b)
1848 register unsigned EMUSHORT *a, *b;
1850 int i;
1852 a += M; /* skip up to significand area */
1853 b += M;
1854 for (i = M; i < NI; i++)
1856 if (*a++ != *b++)
1857 goto difrnt;
1859 return (0);
1861 difrnt:
1862 if (*(--a) > *(--b))
1863 return (1);
1864 else
1865 return (-1);
1868 /* Shift significand of exploded e-type X down by 1 bit. */
1870 static void
1871 eshdn1 (x)
1872 register unsigned EMUSHORT *x;
1874 register unsigned EMUSHORT bits;
1875 int i;
1877 x += M; /* point to significand area */
1879 bits = 0;
1880 for (i = M; i < NI; i++)
1882 if (*x & 1)
1883 bits |= 1;
1884 *x >>= 1;
1885 if (bits & 2)
1886 *x |= 0x8000;
1887 bits <<= 1;
1888 ++x;
1892 /* Shift significand of exploded e-type X up by 1 bit. */
1894 static void
1895 eshup1 (x)
1896 register unsigned EMUSHORT *x;
1898 register unsigned EMUSHORT bits;
1899 int i;
1901 x += NI - 1;
1902 bits = 0;
1904 for (i = M; i < NI; i++)
1906 if (*x & 0x8000)
1907 bits |= 1;
1908 *x <<= 1;
1909 if (bits & 2)
1910 *x |= 1;
1911 bits <<= 1;
1912 --x;
1917 /* Shift significand of exploded e-type X down by 8 bits. */
1919 static void
1920 eshdn8 (x)
1921 register unsigned EMUSHORT *x;
1923 register unsigned EMUSHORT newbyt, oldbyt;
1924 int i;
1926 x += M;
1927 oldbyt = 0;
1928 for (i = M; i < NI; i++)
1930 newbyt = *x << 8;
1931 *x >>= 8;
1932 *x |= oldbyt;
1933 oldbyt = newbyt;
1934 ++x;
1938 /* Shift significand of exploded e-type X up by 8 bits. */
1940 static void
1941 eshup8 (x)
1942 register unsigned EMUSHORT *x;
1944 int i;
1945 register unsigned EMUSHORT newbyt, oldbyt;
1947 x += NI - 1;
1948 oldbyt = 0;
1950 for (i = M; i < NI; i++)
1952 newbyt = *x >> 8;
1953 *x <<= 8;
1954 *x |= oldbyt;
1955 oldbyt = newbyt;
1956 --x;
1960 /* Shift significand of exploded e-type X up by 16 bits. */
1962 static void
1963 eshup6 (x)
1964 register unsigned EMUSHORT *x;
1966 int i;
1967 register unsigned EMUSHORT *p;
1969 p = x + M;
1970 x += M + 1;
1972 for (i = M; i < NI - 1; i++)
1973 *p++ = *x++;
1975 *p = 0;
1978 /* Shift significand of exploded e-type X down by 16 bits. */
1980 static void
1981 eshdn6 (x)
1982 register unsigned EMUSHORT *x;
1984 int i;
1985 register unsigned EMUSHORT *p;
1987 x += NI - 1;
1988 p = x + 1;
1990 for (i = M; i < NI - 1; i++)
1991 *(--p) = *(--x);
1993 *(--p) = 0;
1996 /* Add significands of exploded e-type X and Y. X + Y replaces Y. */
1998 static void
1999 eaddm (x, y)
2000 unsigned EMUSHORT *x, *y;
2002 register unsigned EMULONG a;
2003 int i;
2004 unsigned int carry;
2006 x += NI - 1;
2007 y += NI - 1;
2008 carry = 0;
2009 for (i = M; i < NI; i++)
2011 a = (unsigned EMULONG) (*x) + (unsigned EMULONG) (*y) + carry;
2012 if (a & 0x10000)
2013 carry = 1;
2014 else
2015 carry = 0;
2016 *y = (unsigned EMUSHORT) a;
2017 --x;
2018 --y;
2022 /* Subtract significands of exploded e-type X and Y. Y - X replaces Y. */
2024 static void
2025 esubm (x, y)
2026 unsigned EMUSHORT *x, *y;
2028 unsigned EMULONG a;
2029 int i;
2030 unsigned int carry;
2032 x += NI - 1;
2033 y += NI - 1;
2034 carry = 0;
2035 for (i = M; i < NI; i++)
2037 a = (unsigned EMULONG) (*y) - (unsigned EMULONG) (*x) - carry;
2038 if (a & 0x10000)
2039 carry = 1;
2040 else
2041 carry = 0;
2042 *y = (unsigned EMUSHORT) a;
2043 --x;
2044 --y;
2049 static unsigned EMUSHORT equot[NI];
2052 #if 0
2053 /* Radix 2 shift-and-add versions of multiply and divide */
2056 /* Divide significands */
2058 int
2059 edivm (den, num)
2060 unsigned EMUSHORT den[], num[];
2062 int i;
2063 register unsigned EMUSHORT *p, *q;
2064 unsigned EMUSHORT j;
2066 p = &equot[0];
2067 *p++ = num[0];
2068 *p++ = num[1];
2070 for (i = M; i < NI; i++)
2072 *p++ = 0;
2075 /* Use faster compare and subtraction if denominator has only 15 bits of
2076 significance. */
2078 p = &den[M + 2];
2079 if (*p++ == 0)
2081 for (i = M + 3; i < NI; i++)
2083 if (*p++ != 0)
2084 goto fulldiv;
2086 if ((den[M + 1] & 1) != 0)
2087 goto fulldiv;
2088 eshdn1 (num);
2089 eshdn1 (den);
2091 p = &den[M + 1];
2092 q = &num[M + 1];
2094 for (i = 0; i < NBITS + 2; i++)
2096 if (*p <= *q)
2098 *q -= *p;
2099 j = 1;
2101 else
2103 j = 0;
2105 eshup1 (equot);
2106 equot[NI - 2] |= j;
2107 eshup1 (num);
2109 goto divdon;
2112 /* The number of quotient bits to calculate is NBITS + 1 scaling guard
2113 bit + 1 roundoff bit. */
2115 fulldiv:
2117 p = &equot[NI - 2];
2118 for (i = 0; i < NBITS + 2; i++)
2120 if (ecmpm (den, num) <= 0)
2122 esubm (den, num);
2123 j = 1; /* quotient bit = 1 */
2125 else
2126 j = 0;
2127 eshup1 (equot);
2128 *p |= j;
2129 eshup1 (num);
2132 divdon:
2134 eshdn1 (equot);
2135 eshdn1 (equot);
2137 /* test for nonzero remainder after roundoff bit */
2138 p = &num[M];
2139 j = 0;
2140 for (i = M; i < NI; i++)
2142 j |= *p++;
2144 if (j)
2145 j = 1;
2148 for (i = 0; i < NI; i++)
2149 num[i] = equot[i];
2150 return ((int) j);
2154 /* Multiply significands */
2156 int
2157 emulm (a, b)
2158 unsigned EMUSHORT a[], b[];
2160 unsigned EMUSHORT *p, *q;
2161 int i, j, k;
2163 equot[0] = b[0];
2164 equot[1] = b[1];
2165 for (i = M; i < NI; i++)
2166 equot[i] = 0;
2168 p = &a[NI - 2];
2169 k = NBITS;
2170 while (*p == 0) /* significand is not supposed to be zero */
2172 eshdn6 (a);
2173 k -= 16;
2175 if ((*p & 0xff) == 0)
2177 eshdn8 (a);
2178 k -= 8;
2181 q = &equot[NI - 1];
2182 j = 0;
2183 for (i = 0; i < k; i++)
2185 if (*p & 1)
2186 eaddm (b, equot);
2187 /* remember if there were any nonzero bits shifted out */
2188 if (*q & 1)
2189 j |= 1;
2190 eshdn1 (a);
2191 eshdn1 (equot);
2194 for (i = 0; i < NI; i++)
2195 b[i] = equot[i];
2197 /* return flag for lost nonzero bits */
2198 return (j);
2201 #else
2203 /* Radix 65536 versions of multiply and divide. */
2205 /* Multiply significand of e-type number B
2206 by 16-bit quantity A, return e-type result to C. */
2208 static void
2209 m16m (a, b, c)
2210 unsigned int a;
2211 unsigned EMUSHORT b[], c[];
2213 register unsigned EMUSHORT *pp;
2214 register unsigned EMULONG carry;
2215 unsigned EMUSHORT *ps;
2216 unsigned EMUSHORT p[NI];
2217 unsigned EMULONG aa, m;
2218 int i;
2220 aa = a;
2221 pp = &p[NI-2];
2222 *pp++ = 0;
2223 *pp = 0;
2224 ps = &b[NI-1];
2226 for (i=M+1; i<NI; i++)
2228 if (*ps == 0)
2230 --ps;
2231 --pp;
2232 *(pp-1) = 0;
2234 else
2236 m = (unsigned EMULONG) aa * *ps--;
2237 carry = (m & 0xffff) + *pp;
2238 *pp-- = (unsigned EMUSHORT)carry;
2239 carry = (carry >> 16) + (m >> 16) + *pp;
2240 *pp = (unsigned EMUSHORT)carry;
2241 *(pp-1) = carry >> 16;
2244 for (i=M; i<NI; i++)
2245 c[i] = p[i];
2248 /* Divide significands of exploded e-types NUM / DEN. Neither the
2249 numerator NUM nor the denominator DEN is permitted to have its high guard
2250 word nonzero. */
2252 static int
2253 edivm (den, num)
2254 unsigned EMUSHORT den[], num[];
2256 int i;
2257 register unsigned EMUSHORT *p;
2258 unsigned EMULONG tnum;
2259 unsigned EMUSHORT j, tdenm, tquot;
2260 unsigned EMUSHORT tprod[NI+1];
2262 p = &equot[0];
2263 *p++ = num[0];
2264 *p++ = num[1];
2266 for (i=M; i<NI; i++)
2268 *p++ = 0;
2270 eshdn1 (num);
2271 tdenm = den[M+1];
2272 for (i=M; i<NI; i++)
2274 /* Find trial quotient digit (the radix is 65536). */
2275 tnum = (((unsigned EMULONG) num[M]) << 16) + num[M+1];
2277 /* Do not execute the divide instruction if it will overflow. */
2278 if ((tdenm * 0xffffL) < tnum)
2279 tquot = 0xffff;
2280 else
2281 tquot = tnum / tdenm;
2282 /* Multiply denominator by trial quotient digit. */
2283 m16m ((unsigned int)tquot, den, tprod);
2284 /* The quotient digit may have been overestimated. */
2285 if (ecmpm (tprod, num) > 0)
2287 tquot -= 1;
2288 esubm (den, tprod);
2289 if (ecmpm (tprod, num) > 0)
2291 tquot -= 1;
2292 esubm (den, tprod);
2295 esubm (tprod, num);
2296 equot[i] = tquot;
2297 eshup6(num);
2299 /* test for nonzero remainder after roundoff bit */
2300 p = &num[M];
2301 j = 0;
2302 for (i=M; i<NI; i++)
2304 j |= *p++;
2306 if (j)
2307 j = 1;
2309 for (i=0; i<NI; i++)
2310 num[i] = equot[i];
2312 return ((int)j);
2315 /* Multiply significands of exploded e-type A and B, result in B. */
2317 static int
2318 emulm (a, b)
2319 unsigned EMUSHORT a[], b[];
2321 unsigned EMUSHORT *p, *q;
2322 unsigned EMUSHORT pprod[NI];
2323 unsigned EMUSHORT j;
2324 int i;
2326 equot[0] = b[0];
2327 equot[1] = b[1];
2328 for (i=M; i<NI; i++)
2329 equot[i] = 0;
2331 j = 0;
2332 p = &a[NI-1];
2333 q = &equot[NI-1];
2334 for (i=M+1; i<NI; i++)
2336 if (*p == 0)
2338 --p;
2340 else
2342 m16m ((unsigned int) *p--, b, pprod);
2343 eaddm(pprod, equot);
2345 j |= *q;
2346 eshdn6(equot);
2349 for (i=0; i<NI; i++)
2350 b[i] = equot[i];
2352 /* return flag for lost nonzero bits */
2353 return ((int)j);
2355 #endif
2358 /* Normalize and round off.
2360 The internal format number to be rounded is S.
2361 Input LOST is 0 if the value is exact. This is the so-called sticky bit.
2363 Input SUBFLG indicates whether the number was obtained
2364 by a subtraction operation. In that case if LOST is nonzero
2365 then the number is slightly smaller than indicated.
2367 Input EXP is the biased exponent, which may be negative.
2368 the exponent field of S is ignored but is replaced by
2369 EXP as adjusted by normalization and rounding.
2371 Input RCNTRL is the rounding control. If it is nonzero, the
2372 returned value will be rounded to RNDPRC bits.
2374 For future reference: In order for emdnorm to round off denormal
2375 significands at the right point, the input exponent must be
2376 adjusted to be the actual value it would have after conversion to
2377 the final floating point type. This adjustment has been
2378 implemented for all type conversions (etoe53, etc.) and decimal
2379 conversions, but not for the arithmetic functions (eadd, etc.).
2380 Data types having standard 15-bit exponents are not affected by
2381 this, but SFmode and DFmode are affected. For example, ediv with
2382 rndprc = 24 will not round correctly to 24-bit precision if the
2383 result is denormal. */
2385 static int rlast = -1;
2386 static int rw = 0;
2387 static unsigned EMUSHORT rmsk = 0;
2388 static unsigned EMUSHORT rmbit = 0;
2389 static unsigned EMUSHORT rebit = 0;
2390 static int re = 0;
2391 static unsigned EMUSHORT rbit[NI];
2393 static void
2394 emdnorm (s, lost, subflg, exp, rcntrl)
2395 unsigned EMUSHORT s[];
2396 int lost;
2397 int subflg;
2398 EMULONG exp;
2399 int rcntrl;
2401 int i, j;
2402 unsigned EMUSHORT r;
2404 /* Normalize */
2405 j = enormlz (s);
2407 /* a blank significand could mean either zero or infinity. */
2408 #ifndef INFINITY
2409 if (j > NBITS)
2411 ecleazs (s);
2412 return;
2414 #endif
2415 exp -= j;
2416 #ifndef INFINITY
2417 if (exp >= 32767L)
2418 goto overf;
2419 #else
2420 if ((j > NBITS) && (exp < 32767))
2422 ecleazs (s);
2423 return;
2425 #endif
2426 if (exp < 0L)
2428 if (exp > (EMULONG) (-NBITS - 1))
2430 j = (int) exp;
2431 i = eshift (s, j);
2432 if (i)
2433 lost = 1;
2435 else
2437 ecleazs (s);
2438 return;
2441 /* Round off, unless told not to by rcntrl. */
2442 if (rcntrl == 0)
2443 goto mdfin;
2444 /* Set up rounding parameters if the control register changed. */
2445 if (rndprc != rlast)
2447 ecleaz (rbit);
2448 switch (rndprc)
2450 default:
2451 case NBITS:
2452 rw = NI - 1; /* low guard word */
2453 rmsk = 0xffff;
2454 rmbit = 0x8000;
2455 re = rw - 1;
2456 rebit = 1;
2457 break;
2458 case 113:
2459 rw = 10;
2460 rmsk = 0x7fff;
2461 rmbit = 0x4000;
2462 rebit = 0x8000;
2463 re = rw;
2464 break;
2465 case 64:
2466 rw = 7;
2467 rmsk = 0xffff;
2468 rmbit = 0x8000;
2469 re = rw - 1;
2470 rebit = 1;
2471 break;
2472 /* For DEC or IBM arithmetic */
2473 case 56:
2474 rw = 6;
2475 rmsk = 0xff;
2476 rmbit = 0x80;
2477 rebit = 0x100;
2478 re = rw;
2479 break;
2480 case 53:
2481 rw = 6;
2482 rmsk = 0x7ff;
2483 rmbit = 0x0400;
2484 rebit = 0x800;
2485 re = rw;
2486 break;
2487 case 24:
2488 rw = 4;
2489 rmsk = 0xff;
2490 rmbit = 0x80;
2491 rebit = 0x100;
2492 re = rw;
2493 break;
2495 rbit[re] = rebit;
2496 rlast = rndprc;
2499 /* Shift down 1 temporarily if the data structure has an implied
2500 most significant bit and the number is denormal.
2501 Intel long double denormals also lose one bit of precision. */
2502 if ((exp <= 0) && (rndprc != NBITS)
2503 && ((rndprc != 64) || ((rndprc == 64) && ! REAL_WORDS_BIG_ENDIAN)))
2505 lost |= s[NI - 1] & 1;
2506 eshdn1 (s);
2508 /* Clear out all bits below the rounding bit,
2509 remembering in r if any were nonzero. */
2510 r = s[rw] & rmsk;
2511 if (rndprc < NBITS)
2513 i = rw + 1;
2514 while (i < NI)
2516 if (s[i])
2517 r |= 1;
2518 s[i] = 0;
2519 ++i;
2522 s[rw] &= ~rmsk;
2523 if ((r & rmbit) != 0)
2525 if (r == rmbit)
2527 if (lost == 0)
2528 { /* round to even */
2529 if ((s[re] & rebit) == 0)
2530 goto mddone;
2532 else
2534 if (subflg != 0)
2535 goto mddone;
2538 eaddm (rbit, s);
2540 mddone:
2541 /* Undo the temporary shift for denormal values. */
2542 if ((exp <= 0) && (rndprc != NBITS)
2543 && ((rndprc != 64) || ((rndprc == 64) && ! REAL_WORDS_BIG_ENDIAN)))
2545 eshup1 (s);
2547 if (s[2] != 0)
2548 { /* overflow on roundoff */
2549 eshdn1 (s);
2550 exp += 1;
2552 mdfin:
2553 s[NI - 1] = 0;
2554 if (exp >= 32767L)
2556 #ifndef INFINITY
2557 overf:
2558 #endif
2559 #ifdef INFINITY
2560 s[1] = 32767;
2561 for (i = 2; i < NI - 1; i++)
2562 s[i] = 0;
2563 if (extra_warnings)
2564 warning ("floating point overflow");
2565 #else
2566 s[1] = 32766;
2567 s[2] = 0;
2568 for (i = M + 1; i < NI - 1; i++)
2569 s[i] = 0xffff;
2570 s[NI - 1] = 0;
2571 if ((rndprc < 64) || (rndprc == 113))
2573 s[rw] &= ~rmsk;
2574 if (rndprc == 24)
2576 s[5] = 0;
2577 s[6] = 0;
2580 #endif
2581 return;
2583 if (exp < 0)
2584 s[1] = 0;
2585 else
2586 s[1] = (unsigned EMUSHORT) exp;
2589 /* Subtract. C = B - A, all e type numbers. */
2591 static int subflg = 0;
2593 static void
2594 esub (a, b, c)
2595 unsigned EMUSHORT *a, *b, *c;
2598 #ifdef NANS
2599 if (eisnan (a))
2601 emov (a, c);
2602 return;
2604 if (eisnan (b))
2606 emov (b, c);
2607 return;
2609 /* Infinity minus infinity is a NaN.
2610 Test for subtracting infinities of the same sign. */
2611 if (eisinf (a) && eisinf (b)
2612 && ((eisneg (a) ^ eisneg (b)) == 0))
2614 mtherr ("esub", INVALID);
2615 enan (c, 0);
2616 return;
2618 #endif
2619 subflg = 1;
2620 eadd1 (a, b, c);
2623 /* Add. C = A + B, all e type. */
2625 static void
2626 eadd (a, b, c)
2627 unsigned EMUSHORT *a, *b, *c;
2630 #ifdef NANS
2631 /* NaN plus anything is a NaN. */
2632 if (eisnan (a))
2634 emov (a, c);
2635 return;
2637 if (eisnan (b))
2639 emov (b, c);
2640 return;
2642 /* Infinity minus infinity is a NaN.
2643 Test for adding infinities of opposite signs. */
2644 if (eisinf (a) && eisinf (b)
2645 && ((eisneg (a) ^ eisneg (b)) != 0))
2647 mtherr ("esub", INVALID);
2648 enan (c, 0);
2649 return;
2651 #endif
2652 subflg = 0;
2653 eadd1 (a, b, c);
2656 /* Arithmetic common to both addition and subtraction. */
2658 static void
2659 eadd1 (a, b, c)
2660 unsigned EMUSHORT *a, *b, *c;
2662 unsigned EMUSHORT ai[NI], bi[NI], ci[NI];
2663 int i, lost, j, k;
2664 EMULONG lt, lta, ltb;
2666 #ifdef INFINITY
2667 if (eisinf (a))
2669 emov (a, c);
2670 if (subflg)
2671 eneg (c);
2672 return;
2674 if (eisinf (b))
2676 emov (b, c);
2677 return;
2679 #endif
2680 emovi (a, ai);
2681 emovi (b, bi);
2682 if (subflg)
2683 ai[0] = ~ai[0];
2685 /* compare exponents */
2686 lta = ai[E];
2687 ltb = bi[E];
2688 lt = lta - ltb;
2689 if (lt > 0L)
2690 { /* put the larger number in bi */
2691 emovz (bi, ci);
2692 emovz (ai, bi);
2693 emovz (ci, ai);
2694 ltb = bi[E];
2695 lt = -lt;
2697 lost = 0;
2698 if (lt != 0L)
2700 if (lt < (EMULONG) (-NBITS - 1))
2701 goto done; /* answer same as larger addend */
2702 k = (int) lt;
2703 lost = eshift (ai, k); /* shift the smaller number down */
2705 else
2707 /* exponents were the same, so must compare significands */
2708 i = ecmpm (ai, bi);
2709 if (i == 0)
2710 { /* the numbers are identical in magnitude */
2711 /* if different signs, result is zero */
2712 if (ai[0] != bi[0])
2714 eclear (c);
2715 return;
2717 /* if same sign, result is double */
2718 /* double denormalized tiny number */
2719 if ((bi[E] == 0) && ((bi[3] & 0x8000) == 0))
2721 eshup1 (bi);
2722 goto done;
2724 /* add 1 to exponent unless both are zero! */
2725 for (j = 1; j < NI - 1; j++)
2727 if (bi[j] != 0)
2729 ltb += 1;
2730 if (ltb >= 0x7fff)
2732 eclear (c);
2733 if (ai[0] != 0)
2734 eneg (c);
2735 einfin (c);
2736 return;
2738 break;
2741 bi[E] = (unsigned EMUSHORT) ltb;
2742 goto done;
2744 if (i > 0)
2745 { /* put the larger number in bi */
2746 emovz (bi, ci);
2747 emovz (ai, bi);
2748 emovz (ci, ai);
2751 if (ai[0] == bi[0])
2753 eaddm (ai, bi);
2754 subflg = 0;
2756 else
2758 esubm (ai, bi);
2759 subflg = 1;
2761 emdnorm (bi, lost, subflg, ltb, 64);
2763 done:
2764 emovo (bi, c);
2767 /* Divide: C = B/A, all e type. */
2769 static void
2770 ediv (a, b, c)
2771 unsigned EMUSHORT *a, *b, *c;
2773 unsigned EMUSHORT ai[NI], bi[NI];
2774 int i, sign;
2775 EMULONG lt, lta, ltb;
2777 /* IEEE says if result is not a NaN, the sign is "-" if and only if
2778 operands have opposite signs -- but flush -0 to 0 later if not IEEE. */
2779 sign = eisneg(a) ^ eisneg(b);
2781 #ifdef NANS
2782 /* Return any NaN input. */
2783 if (eisnan (a))
2785 emov (a, c);
2786 return;
2788 if (eisnan (b))
2790 emov (b, c);
2791 return;
2793 /* Zero over zero, or infinity over infinity, is a NaN. */
2794 if (((ecmp (a, ezero) == 0) && (ecmp (b, ezero) == 0))
2795 || (eisinf (a) && eisinf (b)))
2797 mtherr ("ediv", INVALID);
2798 enan (c, sign);
2799 return;
2801 #endif
2802 /* Infinity over anything else is infinity. */
2803 #ifdef INFINITY
2804 if (eisinf (b))
2806 einfin (c);
2807 goto divsign;
2809 /* Anything else over infinity is zero. */
2810 if (eisinf (a))
2812 eclear (c);
2813 goto divsign;
2815 #endif
2816 emovi (a, ai);
2817 emovi (b, bi);
2818 lta = ai[E];
2819 ltb = bi[E];
2820 if (bi[E] == 0)
2821 { /* See if numerator is zero. */
2822 for (i = 1; i < NI - 1; i++)
2824 if (bi[i] != 0)
2826 ltb -= enormlz (bi);
2827 goto dnzro1;
2830 eclear (c);
2831 goto divsign;
2833 dnzro1:
2835 if (ai[E] == 0)
2836 { /* possible divide by zero */
2837 for (i = 1; i < NI - 1; i++)
2839 if (ai[i] != 0)
2841 lta -= enormlz (ai);
2842 goto dnzro2;
2845 /* Divide by zero is not an invalid operation.
2846 It is a divide-by-zero operation! */
2847 einfin (c);
2848 mtherr ("ediv", SING);
2849 goto divsign;
2851 dnzro2:
2853 i = edivm (ai, bi);
2854 /* calculate exponent */
2855 lt = ltb - lta + EXONE;
2856 emdnorm (bi, i, 0, lt, 64);
2857 emovo (bi, c);
2859 divsign:
2861 if (sign
2862 #ifndef IEEE
2863 && (ecmp (c, ezero) != 0)
2864 #endif
2866 *(c+(NE-1)) |= 0x8000;
2867 else
2868 *(c+(NE-1)) &= ~0x8000;
2871 /* Multiply e-types A and B, return e-type product C. */
2873 static void
2874 emul (a, b, c)
2875 unsigned EMUSHORT *a, *b, *c;
2877 unsigned EMUSHORT ai[NI], bi[NI];
2878 int i, j, sign;
2879 EMULONG lt, lta, ltb;
2881 /* IEEE says if result is not a NaN, the sign is "-" if and only if
2882 operands have opposite signs -- but flush -0 to 0 later if not IEEE. */
2883 sign = eisneg(a) ^ eisneg(b);
2885 #ifdef NANS
2886 /* NaN times anything is the same NaN. */
2887 if (eisnan (a))
2889 emov (a, c);
2890 return;
2892 if (eisnan (b))
2894 emov (b, c);
2895 return;
2897 /* Zero times infinity is a NaN. */
2898 if ((eisinf (a) && (ecmp (b, ezero) == 0))
2899 || (eisinf (b) && (ecmp (a, ezero) == 0)))
2901 mtherr ("emul", INVALID);
2902 enan (c, sign);
2903 return;
2905 #endif
2906 /* Infinity times anything else is infinity. */
2907 #ifdef INFINITY
2908 if (eisinf (a) || eisinf (b))
2910 einfin (c);
2911 goto mulsign;
2913 #endif
2914 emovi (a, ai);
2915 emovi (b, bi);
2916 lta = ai[E];
2917 ltb = bi[E];
2918 if (ai[E] == 0)
2920 for (i = 1; i < NI - 1; i++)
2922 if (ai[i] != 0)
2924 lta -= enormlz (ai);
2925 goto mnzer1;
2928 eclear (c);
2929 goto mulsign;
2931 mnzer1:
2933 if (bi[E] == 0)
2935 for (i = 1; i < NI - 1; i++)
2937 if (bi[i] != 0)
2939 ltb -= enormlz (bi);
2940 goto mnzer2;
2943 eclear (c);
2944 goto mulsign;
2946 mnzer2:
2948 /* Multiply significands */
2949 j = emulm (ai, bi);
2950 /* calculate exponent */
2951 lt = lta + ltb - (EXONE - 1);
2952 emdnorm (bi, j, 0, lt, 64);
2953 emovo (bi, c);
2955 mulsign:
2957 if (sign
2958 #ifndef IEEE
2959 && (ecmp (c, ezero) != 0)
2960 #endif
2962 *(c+(NE-1)) |= 0x8000;
2963 else
2964 *(c+(NE-1)) &= ~0x8000;
2967 /* Convert double precision PE to e-type Y. */
2969 static void
2970 e53toe (pe, y)
2971 unsigned EMUSHORT *pe, *y;
2973 #ifdef DEC
2975 dectoe (pe, y);
2977 #else
2978 #ifdef IBM
2980 ibmtoe (pe, y, DFmode);
2982 #else
2983 register unsigned EMUSHORT r;
2984 register unsigned EMUSHORT *e, *p;
2985 unsigned EMUSHORT yy[NI];
2986 int denorm, k;
2988 e = pe;
2989 denorm = 0; /* flag if denormalized number */
2990 ecleaz (yy);
2991 if (! REAL_WORDS_BIG_ENDIAN)
2992 e += 3;
2993 r = *e;
2994 yy[0] = 0;
2995 if (r & 0x8000)
2996 yy[0] = 0xffff;
2997 yy[M] = (r & 0x0f) | 0x10;
2998 r &= ~0x800f; /* strip sign and 4 significand bits */
2999 #ifdef INFINITY
3000 if (r == 0x7ff0)
3002 #ifdef NANS
3003 if (! REAL_WORDS_BIG_ENDIAN)
3005 if (((pe[3] & 0xf) != 0) || (pe[2] != 0)
3006 || (pe[1] != 0) || (pe[0] != 0))
3008 enan (y, yy[0] != 0);
3009 return;
3012 else
3014 if (((pe[0] & 0xf) != 0) || (pe[1] != 0)
3015 || (pe[2] != 0) || (pe[3] != 0))
3017 enan (y, yy[0] != 0);
3018 return;
3021 #endif /* NANS */
3022 eclear (y);
3023 einfin (y);
3024 if (yy[0])
3025 eneg (y);
3026 return;
3028 #endif /* INFINITY */
3029 r >>= 4;
3030 /* If zero exponent, then the significand is denormalized.
3031 So take back the understood high significand bit. */
3033 if (r == 0)
3035 denorm = 1;
3036 yy[M] &= ~0x10;
3038 r += EXONE - 01777;
3039 yy[E] = r;
3040 p = &yy[M + 1];
3041 #ifdef IEEE
3042 if (! REAL_WORDS_BIG_ENDIAN)
3044 *p++ = *(--e);
3045 *p++ = *(--e);
3046 *p++ = *(--e);
3048 else
3050 ++e;
3051 *p++ = *e++;
3052 *p++ = *e++;
3053 *p++ = *e++;
3055 #endif
3056 eshift (yy, -5);
3057 if (denorm)
3058 { /* if zero exponent, then normalize the significand */
3059 if ((k = enormlz (yy)) > NBITS)
3060 ecleazs (yy);
3061 else
3062 yy[E] -= (unsigned EMUSHORT) (k - 1);
3064 emovo (yy, y);
3065 #endif /* not IBM */
3066 #endif /* not DEC */
3069 /* Convert double extended precision float PE to e type Y. */
3071 static void
3072 e64toe (pe, y)
3073 unsigned EMUSHORT *pe, *y;
3075 unsigned EMUSHORT yy[NI];
3076 unsigned EMUSHORT *e, *p, *q;
3077 int i;
3079 e = pe;
3080 p = yy;
3081 for (i = 0; i < NE - 5; i++)
3082 *p++ = 0;
3083 /* This precision is not ordinarily supported on DEC or IBM. */
3084 #ifdef DEC
3085 for (i = 0; i < 5; i++)
3086 *p++ = *e++;
3087 #endif
3088 #ifdef IBM
3089 p = &yy[0] + (NE - 1);
3090 *p-- = *e++;
3091 ++e;
3092 for (i = 0; i < 5; i++)
3093 *p-- = *e++;
3094 #endif
3095 #ifdef IEEE
3096 if (! REAL_WORDS_BIG_ENDIAN)
3098 for (i = 0; i < 5; i++)
3099 *p++ = *e++;
3101 /* For denormal long double Intel format, shift significand up one
3102 -- but only if the top significand bit is zero. A top bit of 1
3103 is "pseudodenormal" when the exponent is zero. */
3104 if((yy[NE-1] & 0x7fff) == 0 && (yy[NE-2] & 0x8000) == 0)
3106 unsigned EMUSHORT temp[NI];
3108 emovi(yy, temp);
3109 eshup1(temp);
3110 emovo(temp,y);
3111 return;
3114 else
3116 p = &yy[0] + (NE - 1);
3117 #ifdef ARM_EXTENDED_IEEE_FORMAT
3118 /* For ARMs, the exponent is in the lowest 15 bits of the word. */
3119 *p-- = (e[0] & 0x8000) | (e[1] & 0x7ffff);
3120 e += 2;
3121 #else
3122 *p-- = *e++;
3123 ++e;
3124 #endif
3125 for (i = 0; i < 4; i++)
3126 *p-- = *e++;
3128 #endif
3129 #ifdef INFINITY
3130 /* Point to the exponent field and check max exponent cases. */
3131 p = &yy[NE - 1];
3132 if ((*p & 0x7fff) == 0x7fff)
3134 #ifdef NANS
3135 if (! REAL_WORDS_BIG_ENDIAN)
3137 for (i = 0; i < 4; i++)
3139 if ((i != 3 && pe[i] != 0)
3140 /* Anything but 0x8000 here, including 0, is a NaN. */
3141 || (i == 3 && pe[i] != 0x8000))
3143 enan (y, (*p & 0x8000) != 0);
3144 return;
3148 else
3150 #ifdef ARM_EXTENDED_IEEE_FORMAT
3151 for (i = 2; i <= 5; i++)
3153 if (pe[i] != 0)
3155 enan (y, (*p & 0x8000) != 0);
3156 return;
3159 #else /* not ARM */
3160 /* In Motorola extended precision format, the most significant
3161 bit of an infinity mantissa could be either 1 or 0. It is
3162 the lower order bits that tell whether the value is a NaN. */
3163 if ((pe[2] & 0x7fff) != 0)
3164 goto bigend_nan;
3166 for (i = 3; i <= 5; i++)
3168 if (pe[i] != 0)
3170 bigend_nan:
3171 enan (y, (*p & 0x8000) != 0);
3172 return;
3175 #endif /* not ARM */
3177 #endif /* NANS */
3178 eclear (y);
3179 einfin (y);
3180 if (*p & 0x8000)
3181 eneg (y);
3182 return;
3184 #endif /* INFINITY */
3185 p = yy;
3186 q = y;
3187 for (i = 0; i < NE; i++)
3188 *q++ = *p++;
3191 /* Convert 128-bit long double precision float PE to e type Y. */
3193 static void
3194 e113toe (pe, y)
3195 unsigned EMUSHORT *pe, *y;
3197 register unsigned EMUSHORT r;
3198 unsigned EMUSHORT *e, *p;
3199 unsigned EMUSHORT yy[NI];
3200 int denorm, i;
3202 e = pe;
3203 denorm = 0;
3204 ecleaz (yy);
3205 #ifdef IEEE
3206 if (! REAL_WORDS_BIG_ENDIAN)
3207 e += 7;
3208 #endif
3209 r = *e;
3210 yy[0] = 0;
3211 if (r & 0x8000)
3212 yy[0] = 0xffff;
3213 r &= 0x7fff;
3214 #ifdef INFINITY
3215 if (r == 0x7fff)
3217 #ifdef NANS
3218 if (! REAL_WORDS_BIG_ENDIAN)
3220 for (i = 0; i < 7; i++)
3222 if (pe[i] != 0)
3224 enan (y, yy[0] != 0);
3225 return;
3229 else
3231 for (i = 1; i < 8; i++)
3233 if (pe[i] != 0)
3235 enan (y, yy[0] != 0);
3236 return;
3240 #endif /* NANS */
3241 eclear (y);
3242 einfin (y);
3243 if (yy[0])
3244 eneg (y);
3245 return;
3247 #endif /* INFINITY */
3248 yy[E] = r;
3249 p = &yy[M + 1];
3250 #ifdef IEEE
3251 if (! REAL_WORDS_BIG_ENDIAN)
3253 for (i = 0; i < 7; i++)
3254 *p++ = *(--e);
3256 else
3258 ++e;
3259 for (i = 0; i < 7; i++)
3260 *p++ = *e++;
3262 #endif
3263 /* If denormal, remove the implied bit; else shift down 1. */
3264 if (r == 0)
3266 yy[M] = 0;
3268 else
3270 yy[M] = 1;
3271 eshift (yy, -1);
3273 emovo (yy, y);
3276 /* Convert single precision float PE to e type Y. */
3278 static void
3279 e24toe (pe, y)
3280 unsigned EMUSHORT *pe, *y;
3282 #ifdef IBM
3284 ibmtoe (pe, y, SFmode);
3286 #else
3287 register unsigned EMUSHORT r;
3288 register unsigned EMUSHORT *e, *p;
3289 unsigned EMUSHORT yy[NI];
3290 int denorm, k;
3292 e = pe;
3293 denorm = 0; /* flag if denormalized number */
3294 ecleaz (yy);
3295 #ifdef IEEE
3296 if (! REAL_WORDS_BIG_ENDIAN)
3297 e += 1;
3298 #endif
3299 #ifdef DEC
3300 e += 1;
3301 #endif
3302 r = *e;
3303 yy[0] = 0;
3304 if (r & 0x8000)
3305 yy[0] = 0xffff;
3306 yy[M] = (r & 0x7f) | 0200;
3307 r &= ~0x807f; /* strip sign and 7 significand bits */
3308 #ifdef INFINITY
3309 if (r == 0x7f80)
3311 #ifdef NANS
3312 if (REAL_WORDS_BIG_ENDIAN)
3314 if (((pe[0] & 0x7f) != 0) || (pe[1] != 0))
3316 enan (y, yy[0] != 0);
3317 return;
3320 else
3322 if (((pe[1] & 0x7f) != 0) || (pe[0] != 0))
3324 enan (y, yy[0] != 0);
3325 return;
3328 #endif /* NANS */
3329 eclear (y);
3330 einfin (y);
3331 if (yy[0])
3332 eneg (y);
3333 return;
3335 #endif /* INFINITY */
3336 r >>= 7;
3337 /* If zero exponent, then the significand is denormalized.
3338 So take back the understood high significand bit. */
3339 if (r == 0)
3341 denorm = 1;
3342 yy[M] &= ~0200;
3344 r += EXONE - 0177;
3345 yy[E] = r;
3346 p = &yy[M + 1];
3347 #ifdef DEC
3348 *p++ = *(--e);
3349 #endif
3350 #ifdef IEEE
3351 if (! REAL_WORDS_BIG_ENDIAN)
3352 *p++ = *(--e);
3353 else
3355 ++e;
3356 *p++ = *e++;
3358 #endif
3359 eshift (yy, -8);
3360 if (denorm)
3361 { /* if zero exponent, then normalize the significand */
3362 if ((k = enormlz (yy)) > NBITS)
3363 ecleazs (yy);
3364 else
3365 yy[E] -= (unsigned EMUSHORT) (k - 1);
3367 emovo (yy, y);
3368 #endif /* not IBM */
3371 /* Convert e-type X to IEEE 128-bit long double format E. */
3373 static void
3374 etoe113 (x, e)
3375 unsigned EMUSHORT *x, *e;
3377 unsigned EMUSHORT xi[NI];
3378 EMULONG exp;
3379 int rndsav;
3381 #ifdef NANS
3382 if (eisnan (x))
3384 make_nan (e, eisneg (x), TFmode);
3385 return;
3387 #endif
3388 emovi (x, xi);
3389 exp = (EMULONG) xi[E];
3390 #ifdef INFINITY
3391 if (eisinf (x))
3392 goto nonorm;
3393 #endif
3394 /* round off to nearest or even */
3395 rndsav = rndprc;
3396 rndprc = 113;
3397 emdnorm (xi, 0, 0, exp, 64);
3398 rndprc = rndsav;
3399 nonorm:
3400 toe113 (xi, e);
3403 /* Convert exploded e-type X, that has already been rounded to
3404 113-bit precision, to IEEE 128-bit long double format Y. */
3406 static void
3407 toe113 (a, b)
3408 unsigned EMUSHORT *a, *b;
3410 register unsigned EMUSHORT *p, *q;
3411 unsigned EMUSHORT i;
3413 #ifdef NANS
3414 if (eiisnan (a))
3416 make_nan (b, eiisneg (a), TFmode);
3417 return;
3419 #endif
3420 p = a;
3421 if (REAL_WORDS_BIG_ENDIAN)
3422 q = b;
3423 else
3424 q = b + 7; /* point to output exponent */
3426 /* If not denormal, delete the implied bit. */
3427 if (a[E] != 0)
3429 eshup1 (a);
3431 /* combine sign and exponent */
3432 i = *p++;
3433 if (REAL_WORDS_BIG_ENDIAN)
3435 if (i)
3436 *q++ = *p++ | 0x8000;
3437 else
3438 *q++ = *p++;
3440 else
3442 if (i)
3443 *q-- = *p++ | 0x8000;
3444 else
3445 *q-- = *p++;
3447 /* skip over guard word */
3448 ++p;
3449 /* move the significand */
3450 if (REAL_WORDS_BIG_ENDIAN)
3452 for (i = 0; i < 7; i++)
3453 *q++ = *p++;
3455 else
3457 for (i = 0; i < 7; i++)
3458 *q-- = *p++;
3462 /* Convert e-type X to IEEE double extended format E. */
3464 static void
3465 etoe64 (x, e)
3466 unsigned EMUSHORT *x, *e;
3468 unsigned EMUSHORT xi[NI];
3469 EMULONG exp;
3470 int rndsav;
3472 #ifdef NANS
3473 if (eisnan (x))
3475 make_nan (e, eisneg (x), XFmode);
3476 return;
3478 #endif
3479 emovi (x, xi);
3480 /* adjust exponent for offset */
3481 exp = (EMULONG) xi[E];
3482 #ifdef INFINITY
3483 if (eisinf (x))
3484 goto nonorm;
3485 #endif
3486 /* round off to nearest or even */
3487 rndsav = rndprc;
3488 rndprc = 64;
3489 emdnorm (xi, 0, 0, exp, 64);
3490 rndprc = rndsav;
3491 nonorm:
3492 toe64 (xi, e);
3495 /* Convert exploded e-type X, that has already been rounded to
3496 64-bit precision, to IEEE double extended format Y. */
3498 static void
3499 toe64 (a, b)
3500 unsigned EMUSHORT *a, *b;
3502 register unsigned EMUSHORT *p, *q;
3503 unsigned EMUSHORT i;
3505 #ifdef NANS
3506 if (eiisnan (a))
3508 make_nan (b, eiisneg (a), XFmode);
3509 return;
3511 #endif
3512 /* Shift denormal long double Intel format significand down one bit. */
3513 if ((a[E] == 0) && ! REAL_WORDS_BIG_ENDIAN)
3514 eshdn1 (a);
3515 p = a;
3516 #ifdef IBM
3517 q = b;
3518 #endif
3519 #ifdef DEC
3520 q = b + 4;
3521 #endif
3522 #ifdef IEEE
3523 if (REAL_WORDS_BIG_ENDIAN)
3524 q = b;
3525 else
3527 q = b + 4; /* point to output exponent */
3528 #if LONG_DOUBLE_TYPE_SIZE == 96
3529 /* Clear the last two bytes of 12-byte Intel format */
3530 *(q+1) = 0;
3531 #endif
3533 #endif
3535 /* combine sign and exponent */
3536 i = *p++;
3537 #ifdef IBM
3538 if (i)
3539 *q++ = *p++ | 0x8000;
3540 else
3541 *q++ = *p++;
3542 *q++ = 0;
3543 #endif
3544 #ifdef DEC
3545 if (i)
3546 *q-- = *p++ | 0x8000;
3547 else
3548 *q-- = *p++;
3549 #endif
3550 #ifdef IEEE
3551 if (REAL_WORDS_BIG_ENDIAN)
3553 #ifdef ARM_EXTENDED_IEEE_FORMAT
3554 /* The exponent is in the lowest 15 bits of the first word. */
3555 *q++ = i ? 0x8000 : 0;
3556 *q++ = *p++;
3557 #else
3558 if (i)
3559 *q++ = *p++ | 0x8000;
3560 else
3561 *q++ = *p++;
3562 *q++ = 0;
3563 #endif
3565 else
3567 if (i)
3568 *q-- = *p++ | 0x8000;
3569 else
3570 *q-- = *p++;
3572 #endif
3573 /* skip over guard word */
3574 ++p;
3575 /* move the significand */
3576 #ifdef IBM
3577 for (i = 0; i < 4; i++)
3578 *q++ = *p++;
3579 #endif
3580 #ifdef DEC
3581 for (i = 0; i < 4; i++)
3582 *q-- = *p++;
3583 #endif
3584 #ifdef IEEE
3585 if (REAL_WORDS_BIG_ENDIAN)
3587 for (i = 0; i < 4; i++)
3588 *q++ = *p++;
3590 else
3592 #ifdef INFINITY
3593 if (eiisinf (a))
3595 /* Intel long double infinity significand. */
3596 *q-- = 0x8000;
3597 *q-- = 0;
3598 *q-- = 0;
3599 *q = 0;
3600 return;
3602 #endif
3603 for (i = 0; i < 4; i++)
3604 *q-- = *p++;
3606 #endif
3609 /* e type to double precision. */
3611 #ifdef DEC
3612 /* Convert e-type X to DEC-format double E. */
3614 static void
3615 etoe53 (x, e)
3616 unsigned EMUSHORT *x, *e;
3618 etodec (x, e); /* see etodec.c */
3621 /* Convert exploded e-type X, that has already been rounded to
3622 56-bit double precision, to DEC double Y. */
3624 static void
3625 toe53 (x, y)
3626 unsigned EMUSHORT *x, *y;
3628 todec (x, y);
3631 #else
3632 #ifdef IBM
3633 /* Convert e-type X to IBM 370-format double E. */
3635 static void
3636 etoe53 (x, e)
3637 unsigned EMUSHORT *x, *e;
3639 etoibm (x, e, DFmode);
3642 /* Convert exploded e-type X, that has already been rounded to
3643 56-bit precision, to IBM 370 double Y. */
3645 static void
3646 toe53 (x, y)
3647 unsigned EMUSHORT *x, *y;
3649 toibm (x, y, DFmode);
3652 #else /* it's neither DEC nor IBM */
3654 /* Convert e-type X to IEEE double E. */
3656 static void
3657 etoe53 (x, e)
3658 unsigned EMUSHORT *x, *e;
3660 unsigned EMUSHORT xi[NI];
3661 EMULONG exp;
3662 int rndsav;
3664 #ifdef NANS
3665 if (eisnan (x))
3667 make_nan (e, eisneg (x), DFmode);
3668 return;
3670 #endif
3671 emovi (x, xi);
3672 /* adjust exponent for offsets */
3673 exp = (EMULONG) xi[E] - (EXONE - 0x3ff);
3674 #ifdef INFINITY
3675 if (eisinf (x))
3676 goto nonorm;
3677 #endif
3678 /* round off to nearest or even */
3679 rndsav = rndprc;
3680 rndprc = 53;
3681 emdnorm (xi, 0, 0, exp, 64);
3682 rndprc = rndsav;
3683 nonorm:
3684 toe53 (xi, e);
3687 /* Convert exploded e-type X, that has already been rounded to
3688 53-bit precision, to IEEE double Y. */
3690 static void
3691 toe53 (x, y)
3692 unsigned EMUSHORT *x, *y;
3694 unsigned EMUSHORT i;
3695 unsigned EMUSHORT *p;
3697 #ifdef NANS
3698 if (eiisnan (x))
3700 make_nan (y, eiisneg (x), DFmode);
3701 return;
3703 #endif
3704 p = &x[0];
3705 #ifdef IEEE
3706 if (! REAL_WORDS_BIG_ENDIAN)
3707 y += 3;
3708 #endif
3709 *y = 0; /* output high order */
3710 if (*p++)
3711 *y = 0x8000; /* output sign bit */
3713 i = *p++;
3714 if (i >= (unsigned int) 2047)
3716 /* Saturate at largest number less than infinity. */
3717 #ifdef INFINITY
3718 *y |= 0x7ff0;
3719 if (! REAL_WORDS_BIG_ENDIAN)
3721 *(--y) = 0;
3722 *(--y) = 0;
3723 *(--y) = 0;
3725 else
3727 ++y;
3728 *y++ = 0;
3729 *y++ = 0;
3730 *y++ = 0;
3732 #else
3733 *y |= (unsigned EMUSHORT) 0x7fef;
3734 if (! REAL_WORDS_BIG_ENDIAN)
3736 *(--y) = 0xffff;
3737 *(--y) = 0xffff;
3738 *(--y) = 0xffff;
3740 else
3742 ++y;
3743 *y++ = 0xffff;
3744 *y++ = 0xffff;
3745 *y++ = 0xffff;
3747 #endif
3748 return;
3750 if (i == 0)
3752 eshift (x, 4);
3754 else
3756 i <<= 4;
3757 eshift (x, 5);
3759 i |= *p++ & (unsigned EMUSHORT) 0x0f; /* *p = xi[M] */
3760 *y |= (unsigned EMUSHORT) i; /* high order output already has sign bit set */
3761 if (! REAL_WORDS_BIG_ENDIAN)
3763 *(--y) = *p++;
3764 *(--y) = *p++;
3765 *(--y) = *p;
3767 else
3769 ++y;
3770 *y++ = *p++;
3771 *y++ = *p++;
3772 *y++ = *p++;
3776 #endif /* not IBM */
3777 #endif /* not DEC */
3781 /* e type to single precision. */
3783 #ifdef IBM
3784 /* Convert e-type X to IBM 370 float E. */
3786 static void
3787 etoe24 (x, e)
3788 unsigned EMUSHORT *x, *e;
3790 etoibm (x, e, SFmode);
3793 /* Convert exploded e-type X, that has already been rounded to
3794 float precision, to IBM 370 float Y. */
3796 static void
3797 toe24 (x, y)
3798 unsigned EMUSHORT *x, *y;
3800 toibm (x, y, SFmode);
3803 #else
3804 /* Convert e-type X to IEEE float E. DEC float is the same as IEEE float. */
3806 static void
3807 etoe24 (x, e)
3808 unsigned EMUSHORT *x, *e;
3810 EMULONG exp;
3811 unsigned EMUSHORT xi[NI];
3812 int rndsav;
3814 #ifdef NANS
3815 if (eisnan (x))
3817 make_nan (e, eisneg (x), SFmode);
3818 return;
3820 #endif
3821 emovi (x, xi);
3822 /* adjust exponent for offsets */
3823 exp = (EMULONG) xi[E] - (EXONE - 0177);
3824 #ifdef INFINITY
3825 if (eisinf (x))
3826 goto nonorm;
3827 #endif
3828 /* round off to nearest or even */
3829 rndsav = rndprc;
3830 rndprc = 24;
3831 emdnorm (xi, 0, 0, exp, 64);
3832 rndprc = rndsav;
3833 nonorm:
3834 toe24 (xi, e);
3837 /* Convert exploded e-type X, that has already been rounded to
3838 float precision, to IEEE float Y. */
3840 static void
3841 toe24 (x, y)
3842 unsigned EMUSHORT *x, *y;
3844 unsigned EMUSHORT i;
3845 unsigned EMUSHORT *p;
3847 #ifdef NANS
3848 if (eiisnan (x))
3850 make_nan (y, eiisneg (x), SFmode);
3851 return;
3853 #endif
3854 p = &x[0];
3855 #ifdef IEEE
3856 if (! REAL_WORDS_BIG_ENDIAN)
3857 y += 1;
3858 #endif
3859 #ifdef DEC
3860 y += 1;
3861 #endif
3862 *y = 0; /* output high order */
3863 if (*p++)
3864 *y = 0x8000; /* output sign bit */
3866 i = *p++;
3867 /* Handle overflow cases. */
3868 if (i >= 255)
3870 #ifdef INFINITY
3871 *y |= (unsigned EMUSHORT) 0x7f80;
3872 #ifdef DEC
3873 *(--y) = 0;
3874 #endif
3875 #ifdef IEEE
3876 if (! REAL_WORDS_BIG_ENDIAN)
3877 *(--y) = 0;
3878 else
3880 ++y;
3881 *y = 0;
3883 #endif
3884 #else /* no INFINITY */
3885 *y |= (unsigned EMUSHORT) 0x7f7f;
3886 #ifdef DEC
3887 *(--y) = 0xffff;
3888 #endif
3889 #ifdef IEEE
3890 if (! REAL_WORDS_BIG_ENDIAN)
3891 *(--y) = 0xffff;
3892 else
3894 ++y;
3895 *y = 0xffff;
3897 #endif
3898 #ifdef ERANGE
3899 errno = ERANGE;
3900 #endif
3901 #endif /* no INFINITY */
3902 return;
3904 if (i == 0)
3906 eshift (x, 7);
3908 else
3910 i <<= 7;
3911 eshift (x, 8);
3913 i |= *p++ & (unsigned EMUSHORT) 0x7f; /* *p = xi[M] */
3914 /* High order output already has sign bit set. */
3915 *y |= i;
3916 #ifdef DEC
3917 *(--y) = *p;
3918 #endif
3919 #ifdef IEEE
3920 if (! REAL_WORDS_BIG_ENDIAN)
3921 *(--y) = *p;
3922 else
3924 ++y;
3925 *y = *p;
3927 #endif
3929 #endif /* not IBM */
3931 /* Compare two e type numbers.
3932 Return +1 if a > b
3933 0 if a == b
3934 -1 if a < b
3935 -2 if either a or b is a NaN. */
3937 static int
3938 ecmp (a, b)
3939 unsigned EMUSHORT *a, *b;
3941 unsigned EMUSHORT ai[NI], bi[NI];
3942 register unsigned EMUSHORT *p, *q;
3943 register int i;
3944 int msign;
3946 #ifdef NANS
3947 if (eisnan (a) || eisnan (b))
3948 return (-2);
3949 #endif
3950 emovi (a, ai);
3951 p = ai;
3952 emovi (b, bi);
3953 q = bi;
3955 if (*p != *q)
3956 { /* the signs are different */
3957 /* -0 equals + 0 */
3958 for (i = 1; i < NI - 1; i++)
3960 if (ai[i] != 0)
3961 goto nzro;
3962 if (bi[i] != 0)
3963 goto nzro;
3965 return (0);
3966 nzro:
3967 if (*p == 0)
3968 return (1);
3969 else
3970 return (-1);
3972 /* both are the same sign */
3973 if (*p == 0)
3974 msign = 1;
3975 else
3976 msign = -1;
3977 i = NI - 1;
3980 if (*p++ != *q++)
3982 goto diff;
3985 while (--i > 0);
3987 return (0); /* equality */
3989 diff:
3991 if (*(--p) > *(--q))
3992 return (msign); /* p is bigger */
3993 else
3994 return (-msign); /* p is littler */
3997 /* Find e-type nearest integer to X, as floor (X + 0.5). */
3999 static void
4000 eround (x, y)
4001 unsigned EMUSHORT *x, *y;
4003 eadd (ehalf, x, y);
4004 efloor (y, y);
4007 /* Convert HOST_WIDE_INT LP to e type Y. */
4009 static void
4010 ltoe (lp, y)
4011 HOST_WIDE_INT *lp;
4012 unsigned EMUSHORT *y;
4014 unsigned EMUSHORT yi[NI];
4015 unsigned HOST_WIDE_INT ll;
4016 int k;
4018 ecleaz (yi);
4019 if (*lp < 0)
4021 /* make it positive */
4022 ll = (unsigned HOST_WIDE_INT) (-(*lp));
4023 yi[0] = 0xffff; /* put correct sign in the e type number */
4025 else
4027 ll = (unsigned HOST_WIDE_INT) (*lp);
4029 /* move the long integer to yi significand area */
4030 #if HOST_BITS_PER_WIDE_INT == 64
4031 yi[M] = (unsigned EMUSHORT) (ll >> 48);
4032 yi[M + 1] = (unsigned EMUSHORT) (ll >> 32);
4033 yi[M + 2] = (unsigned EMUSHORT) (ll >> 16);
4034 yi[M + 3] = (unsigned EMUSHORT) ll;
4035 yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
4036 #else
4037 yi[M] = (unsigned EMUSHORT) (ll >> 16);
4038 yi[M + 1] = (unsigned EMUSHORT) ll;
4039 yi[E] = EXONE + 15; /* exponent if normalize shift count were 0 */
4040 #endif
4042 if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
4043 ecleaz (yi); /* it was zero */
4044 else
4045 yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
4046 emovo (yi, y); /* output the answer */
4049 /* Convert unsigned HOST_WIDE_INT LP to e type Y. */
4051 static void
4052 ultoe (lp, y)
4053 unsigned HOST_WIDE_INT *lp;
4054 unsigned EMUSHORT *y;
4056 unsigned EMUSHORT yi[NI];
4057 unsigned HOST_WIDE_INT ll;
4058 int k;
4060 ecleaz (yi);
4061 ll = *lp;
4063 /* move the long integer to ayi significand area */
4064 #if HOST_BITS_PER_WIDE_INT == 64
4065 yi[M] = (unsigned EMUSHORT) (ll >> 48);
4066 yi[M + 1] = (unsigned EMUSHORT) (ll >> 32);
4067 yi[M + 2] = (unsigned EMUSHORT) (ll >> 16);
4068 yi[M + 3] = (unsigned EMUSHORT) ll;
4069 yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
4070 #else
4071 yi[M] = (unsigned EMUSHORT) (ll >> 16);
4072 yi[M + 1] = (unsigned EMUSHORT) ll;
4073 yi[E] = EXONE + 15; /* exponent if normalize shift count were 0 */
4074 #endif
4076 if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
4077 ecleaz (yi); /* it was zero */
4078 else
4079 yi[E] -= (unsigned EMUSHORT) k; /* subtract shift count from exponent */
4080 emovo (yi, y); /* output the answer */
4084 /* Find signed HOST_WIDE_INT integer I and floating point fractional
4085 part FRAC of e-type (packed internal format) floating point input X.
4086 The integer output I has the sign of the input, except that
4087 positive overflow is permitted if FIXUNS_TRUNC_LIKE_FIX_TRUNC.
4088 The output e-type fraction FRAC is the positive fractional
4089 part of abs (X). */
4091 static void
4092 eifrac (x, i, frac)
4093 unsigned EMUSHORT *x;
4094 HOST_WIDE_INT *i;
4095 unsigned EMUSHORT *frac;
4097 unsigned EMUSHORT xi[NI];
4098 int j, k;
4099 unsigned HOST_WIDE_INT ll;
4101 emovi (x, xi);
4102 k = (int) xi[E] - (EXONE - 1);
4103 if (k <= 0)
4105 /* if exponent <= 0, integer = 0 and real output is fraction */
4106 *i = 0L;
4107 emovo (xi, frac);
4108 return;
4110 if (k > (HOST_BITS_PER_WIDE_INT - 1))
4112 /* long integer overflow: output large integer
4113 and correct fraction */
4114 if (xi[0])
4115 *i = ((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1);
4116 else
4118 #ifdef FIXUNS_TRUNC_LIKE_FIX_TRUNC
4119 /* In this case, let it overflow and convert as if unsigned. */
4120 euifrac (x, &ll, frac);
4121 *i = (HOST_WIDE_INT) ll;
4122 return;
4123 #else
4124 /* In other cases, return the largest positive integer. */
4125 *i = (((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1)) - 1;
4126 #endif
4128 eshift (xi, k);
4129 if (extra_warnings)
4130 warning ("overflow on truncation to integer");
4132 else if (k > 16)
4134 /* Shift more than 16 bits: first shift up k-16 mod 16,
4135 then shift up by 16's. */
4136 j = k - ((k >> 4) << 4);
4137 eshift (xi, j);
4138 ll = xi[M];
4139 k -= j;
4142 eshup6 (xi);
4143 ll = (ll << 16) | xi[M];
4145 while ((k -= 16) > 0);
4146 *i = ll;
4147 if (xi[0])
4148 *i = -(*i);
4150 else
4152 /* shift not more than 16 bits */
4153 eshift (xi, k);
4154 *i = (HOST_WIDE_INT) xi[M] & 0xffff;
4155 if (xi[0])
4156 *i = -(*i);
4158 xi[0] = 0;
4159 xi[E] = EXONE - 1;
4160 xi[M] = 0;
4161 if ((k = enormlz (xi)) > NBITS)
4162 ecleaz (xi);
4163 else
4164 xi[E] -= (unsigned EMUSHORT) k;
4166 emovo (xi, frac);
4170 /* Find unsigned HOST_WIDE_INT integer I and floating point fractional part
4171 FRAC of e-type X. A negative input yields integer output = 0 but
4172 correct fraction. */
4174 static void
4175 euifrac (x, i, frac)
4176 unsigned EMUSHORT *x;
4177 unsigned HOST_WIDE_INT *i;
4178 unsigned EMUSHORT *frac;
4180 unsigned HOST_WIDE_INT ll;
4181 unsigned EMUSHORT xi[NI];
4182 int j, k;
4184 emovi (x, xi);
4185 k = (int) xi[E] - (EXONE - 1);
4186 if (k <= 0)
4188 /* if exponent <= 0, integer = 0 and argument is fraction */
4189 *i = 0L;
4190 emovo (xi, frac);
4191 return;
4193 if (k > HOST_BITS_PER_WIDE_INT)
4195 /* Long integer overflow: output large integer
4196 and correct fraction.
4197 Note, the BSD microvax compiler says that ~(0UL)
4198 is a syntax error. */
4199 *i = ~(0L);
4200 eshift (xi, k);
4201 if (extra_warnings)
4202 warning ("overflow on truncation to unsigned integer");
4204 else if (k > 16)
4206 /* Shift more than 16 bits: first shift up k-16 mod 16,
4207 then shift up by 16's. */
4208 j = k - ((k >> 4) << 4);
4209 eshift (xi, j);
4210 ll = xi[M];
4211 k -= j;
4214 eshup6 (xi);
4215 ll = (ll << 16) | xi[M];
4217 while ((k -= 16) > 0);
4218 *i = ll;
4220 else
4222 /* shift not more than 16 bits */
4223 eshift (xi, k);
4224 *i = (HOST_WIDE_INT) xi[M] & 0xffff;
4227 if (xi[0]) /* A negative value yields unsigned integer 0. */
4228 *i = 0L;
4230 xi[0] = 0;
4231 xi[E] = EXONE - 1;
4232 xi[M] = 0;
4233 if ((k = enormlz (xi)) > NBITS)
4234 ecleaz (xi);
4235 else
4236 xi[E] -= (unsigned EMUSHORT) k;
4238 emovo (xi, frac);
4241 /* Shift the significand of exploded e-type X up or down by SC bits. */
4243 static int
4244 eshift (x, sc)
4245 unsigned EMUSHORT *x;
4246 int sc;
4248 unsigned EMUSHORT lost;
4249 unsigned EMUSHORT *p;
4251 if (sc == 0)
4252 return (0);
4254 lost = 0;
4255 p = x + NI - 1;
4257 if (sc < 0)
4259 sc = -sc;
4260 while (sc >= 16)
4262 lost |= *p; /* remember lost bits */
4263 eshdn6 (x);
4264 sc -= 16;
4267 while (sc >= 8)
4269 lost |= *p & 0xff;
4270 eshdn8 (x);
4271 sc -= 8;
4274 while (sc > 0)
4276 lost |= *p & 1;
4277 eshdn1 (x);
4278 sc -= 1;
4281 else
4283 while (sc >= 16)
4285 eshup6 (x);
4286 sc -= 16;
4289 while (sc >= 8)
4291 eshup8 (x);
4292 sc -= 8;
4295 while (sc > 0)
4297 eshup1 (x);
4298 sc -= 1;
4301 if (lost)
4302 lost = 1;
4303 return ((int) lost);
4306 /* Shift normalize the significand area of exploded e-type X.
4307 Return the shift count (up = positive). */
4309 static int
4310 enormlz (x)
4311 unsigned EMUSHORT x[];
4313 register unsigned EMUSHORT *p;
4314 int sc;
4316 sc = 0;
4317 p = &x[M];
4318 if (*p != 0)
4319 goto normdn;
4320 ++p;
4321 if (*p & 0x8000)
4322 return (0); /* already normalized */
4323 while (*p == 0)
4325 eshup6 (x);
4326 sc += 16;
4328 /* With guard word, there are NBITS+16 bits available.
4329 Return true if all are zero. */
4330 if (sc > NBITS)
4331 return (sc);
4333 /* see if high byte is zero */
4334 while ((*p & 0xff00) == 0)
4336 eshup8 (x);
4337 sc += 8;
4339 /* now shift 1 bit at a time */
4340 while ((*p & 0x8000) == 0)
4342 eshup1 (x);
4343 sc += 1;
4344 if (sc > NBITS)
4346 mtherr ("enormlz", UNDERFLOW);
4347 return (sc);
4350 return (sc);
4352 /* Normalize by shifting down out of the high guard word
4353 of the significand */
4354 normdn:
4356 if (*p & 0xff00)
4358 eshdn8 (x);
4359 sc -= 8;
4361 while (*p != 0)
4363 eshdn1 (x);
4364 sc -= 1;
4366 if (sc < -NBITS)
4368 mtherr ("enormlz", OVERFLOW);
4369 return (sc);
4372 return (sc);
4375 /* Powers of ten used in decimal <-> binary conversions. */
4377 #define NTEN 12
4378 #define MAXP 4096
4380 #if LONG_DOUBLE_TYPE_SIZE == 128
4381 static unsigned EMUSHORT etens[NTEN + 1][NE] =
4383 {0x6576, 0x4a92, 0x804a, 0x153f,
4384 0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,}, /* 10**4096 */
4385 {0x6a32, 0xce52, 0x329a, 0x28ce,
4386 0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,}, /* 10**2048 */
4387 {0x526c, 0x50ce, 0xf18b, 0x3d28,
4388 0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
4389 {0x9c66, 0x58f8, 0xbc50, 0x5c54,
4390 0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
4391 {0x851e, 0xeab7, 0x98fe, 0x901b,
4392 0xddbb, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
4393 {0x0235, 0x0137, 0x36b1, 0x336c,
4394 0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
4395 {0x50f8, 0x25fb, 0xc76b, 0x6b71,
4396 0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
4397 {0x0000, 0x0000, 0x0000, 0x0000,
4398 0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
4399 {0x0000, 0x0000, 0x0000, 0x0000,
4400 0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
4401 {0x0000, 0x0000, 0x0000, 0x0000,
4402 0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
4403 {0x0000, 0x0000, 0x0000, 0x0000,
4404 0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
4405 {0x0000, 0x0000, 0x0000, 0x0000,
4406 0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
4407 {0x0000, 0x0000, 0x0000, 0x0000,
4408 0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,}, /* 10**1 */
4411 static unsigned EMUSHORT emtens[NTEN + 1][NE] =
4413 {0x2030, 0xcffc, 0xa1c3, 0x8123,
4414 0x2de3, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,}, /* 10**-4096 */
4415 {0x8264, 0xd2cb, 0xf2ea, 0x12d4,
4416 0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,}, /* 10**-2048 */
4417 {0xf53f, 0xf698, 0x6bd3, 0x0158,
4418 0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
4419 {0xe731, 0x04d4, 0xe3f2, 0xd332,
4420 0x7132, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
4421 {0xa23e, 0x5308, 0xfefb, 0x1155,
4422 0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
4423 {0xe26d, 0xdbde, 0xd05d, 0xb3f6,
4424 0xac7c, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
4425 {0x2a20, 0x6224, 0x47b3, 0x98d7,
4426 0x3f23, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
4427 {0x0b5b, 0x4af2, 0xa581, 0x18ed,
4428 0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
4429 {0xbf71, 0xa9b3, 0x7989, 0xbe68,
4430 0x4c2e, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
4431 {0x3d4d, 0x7c3d, 0x36ba, 0x0d2b,
4432 0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
4433 {0xc155, 0xa4a8, 0x404e, 0x6113,
4434 0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
4435 {0xd70a, 0x70a3, 0x0a3d, 0xa3d7,
4436 0x3d70, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
4437 {0xcccd, 0xcccc, 0xcccc, 0xcccc,
4438 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,}, /* 10**-1 */
4440 #else
4441 /* LONG_DOUBLE_TYPE_SIZE is other than 128 */
4442 static unsigned EMUSHORT etens[NTEN + 1][NE] =
4444 {0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,}, /* 10**4096 */
4445 {0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,}, /* 10**2048 */
4446 {0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
4447 {0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
4448 {0xddbc, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
4449 {0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
4450 {0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
4451 {0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
4452 {0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
4453 {0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
4454 {0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
4455 {0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
4456 {0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,}, /* 10**1 */
4459 static unsigned EMUSHORT emtens[NTEN + 1][NE] =
4461 {0x2de4, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,}, /* 10**-4096 */
4462 {0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,}, /* 10**-2048 */
4463 {0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
4464 {0x7133, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
4465 {0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
4466 {0xac7d, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
4467 {0x3f24, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
4468 {0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
4469 {0x4c2f, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
4470 {0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
4471 {0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
4472 {0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
4473 {0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,}, /* 10**-1 */
4475 #endif
4477 /* Convert float value X to ASCII string STRING with NDIG digits after
4478 the decimal point. */
4480 static void
4481 e24toasc (x, string, ndigs)
4482 unsigned EMUSHORT x[];
4483 char *string;
4484 int ndigs;
4486 unsigned EMUSHORT w[NI];
4488 e24toe (x, w);
4489 etoasc (w, string, ndigs);
4492 /* Convert double value X to ASCII string STRING with NDIG digits after
4493 the decimal point. */
4495 static void
4496 e53toasc (x, string, ndigs)
4497 unsigned EMUSHORT x[];
4498 char *string;
4499 int ndigs;
4501 unsigned EMUSHORT w[NI];
4503 e53toe (x, w);
4504 etoasc (w, string, ndigs);
4507 /* Convert double extended value X to ASCII string STRING with NDIG digits
4508 after the decimal point. */
4510 static void
4511 e64toasc (x, string, ndigs)
4512 unsigned EMUSHORT x[];
4513 char *string;
4514 int ndigs;
4516 unsigned EMUSHORT w[NI];
4518 e64toe (x, w);
4519 etoasc (w, string, ndigs);
4522 /* Convert 128-bit long double value X to ASCII string STRING with NDIG digits
4523 after the decimal point. */
4525 static void
4526 e113toasc (x, string, ndigs)
4527 unsigned EMUSHORT x[];
4528 char *string;
4529 int ndigs;
4531 unsigned EMUSHORT w[NI];
4533 e113toe (x, w);
4534 etoasc (w, string, ndigs);
4537 /* Convert e-type X to ASCII string STRING with NDIGS digits after
4538 the decimal point. */
4540 static char wstring[80]; /* working storage for ASCII output */
4542 static void
4543 etoasc (x, string, ndigs)
4544 unsigned EMUSHORT x[];
4545 char *string;
4546 int ndigs;
4548 EMUSHORT digit;
4549 unsigned EMUSHORT y[NI], t[NI], u[NI], w[NI];
4550 unsigned EMUSHORT *p, *r, *ten;
4551 unsigned EMUSHORT sign;
4552 int i, j, k, expon, rndsav;
4553 char *s, *ss;
4554 unsigned EMUSHORT m;
4557 rndsav = rndprc;
4558 ss = string;
4559 s = wstring;
4560 *ss = '\0';
4561 *s = '\0';
4562 #ifdef NANS
4563 if (eisnan (x))
4565 sprintf (wstring, " NaN ");
4566 goto bxit;
4568 #endif
4569 rndprc = NBITS; /* set to full precision */
4570 emov (x, y); /* retain external format */
4571 if (y[NE - 1] & 0x8000)
4573 sign = 0xffff;
4574 y[NE - 1] &= 0x7fff;
4576 else
4578 sign = 0;
4580 expon = 0;
4581 ten = &etens[NTEN][0];
4582 emov (eone, t);
4583 /* Test for zero exponent */
4584 if (y[NE - 1] == 0)
4586 for (k = 0; k < NE - 1; k++)
4588 if (y[k] != 0)
4589 goto tnzro; /* denormalized number */
4591 goto isone; /* valid all zeros */
4593 tnzro:
4595 /* Test for infinity. */
4596 if (y[NE - 1] == 0x7fff)
4598 if (sign)
4599 sprintf (wstring, " -Infinity ");
4600 else
4601 sprintf (wstring, " Infinity ");
4602 goto bxit;
4605 /* Test for exponent nonzero but significand denormalized.
4606 * This is an error condition.
4608 if ((y[NE - 1] != 0) && ((y[NE - 2] & 0x8000) == 0))
4610 mtherr ("etoasc", DOMAIN);
4611 sprintf (wstring, "NaN");
4612 goto bxit;
4615 /* Compare to 1.0 */
4616 i = ecmp (eone, y);
4617 if (i == 0)
4618 goto isone;
4620 if (i == -2)
4621 abort ();
4623 if (i < 0)
4624 { /* Number is greater than 1 */
4625 /* Convert significand to an integer and strip trailing decimal zeros. */
4626 emov (y, u);
4627 u[NE - 1] = EXONE + NBITS - 1;
4629 p = &etens[NTEN - 4][0];
4630 m = 16;
4633 ediv (p, u, t);
4634 efloor (t, w);
4635 for (j = 0; j < NE - 1; j++)
4637 if (t[j] != w[j])
4638 goto noint;
4640 emov (t, u);
4641 expon += (int) m;
4642 noint:
4643 p += NE;
4644 m >>= 1;
4646 while (m != 0);
4648 /* Rescale from integer significand */
4649 u[NE - 1] += y[NE - 1] - (unsigned int) (EXONE + NBITS - 1);
4650 emov (u, y);
4651 /* Find power of 10 */
4652 emov (eone, t);
4653 m = MAXP;
4654 p = &etens[0][0];
4655 /* An unordered compare result shouldn't happen here. */
4656 while (ecmp (ten, u) <= 0)
4658 if (ecmp (p, u) <= 0)
4660 ediv (p, u, u);
4661 emul (p, t, t);
4662 expon += (int) m;
4664 m >>= 1;
4665 if (m == 0)
4666 break;
4667 p += NE;
4670 else
4671 { /* Number is less than 1.0 */
4672 /* Pad significand with trailing decimal zeros. */
4673 if (y[NE - 1] == 0)
4675 while ((y[NE - 2] & 0x8000) == 0)
4677 emul (ten, y, y);
4678 expon -= 1;
4681 else
4683 emovi (y, w);
4684 for (i = 0; i < NDEC + 1; i++)
4686 if ((w[NI - 1] & 0x7) != 0)
4687 break;
4688 /* multiply by 10 */
4689 emovz (w, u);
4690 eshdn1 (u);
4691 eshdn1 (u);
4692 eaddm (w, u);
4693 u[1] += 3;
4694 while (u[2] != 0)
4696 eshdn1 (u);
4697 u[1] += 1;
4699 if (u[NI - 1] != 0)
4700 break;
4701 if (eone[NE - 1] <= u[1])
4702 break;
4703 emovz (u, w);
4704 expon -= 1;
4706 emovo (w, y);
4708 k = -MAXP;
4709 p = &emtens[0][0];
4710 r = &etens[0][0];
4711 emov (y, w);
4712 emov (eone, t);
4713 while (ecmp (eone, w) > 0)
4715 if (ecmp (p, w) >= 0)
4717 emul (r, w, w);
4718 emul (r, t, t);
4719 expon += k;
4721 k /= 2;
4722 if (k == 0)
4723 break;
4724 p += NE;
4725 r += NE;
4727 ediv (t, eone, t);
4729 isone:
4730 /* Find the first (leading) digit. */
4731 emovi (t, w);
4732 emovz (w, t);
4733 emovi (y, w);
4734 emovz (w, y);
4735 eiremain (t, y);
4736 digit = equot[NI - 1];
4737 while ((digit == 0) && (ecmp (y, ezero) != 0))
4739 eshup1 (y);
4740 emovz (y, u);
4741 eshup1 (u);
4742 eshup1 (u);
4743 eaddm (u, y);
4744 eiremain (t, y);
4745 digit = equot[NI - 1];
4746 expon -= 1;
4748 s = wstring;
4749 if (sign)
4750 *s++ = '-';
4751 else
4752 *s++ = ' ';
4753 /* Examine number of digits requested by caller. */
4754 if (ndigs < 0)
4755 ndigs = 0;
4756 if (ndigs > NDEC)
4757 ndigs = NDEC;
4758 if (digit == 10)
4760 *s++ = '1';
4761 *s++ = '.';
4762 if (ndigs > 0)
4764 *s++ = '0';
4765 ndigs -= 1;
4767 expon += 1;
4769 else
4771 *s++ = (char)digit + '0';
4772 *s++ = '.';
4774 /* Generate digits after the decimal point. */
4775 for (k = 0; k <= ndigs; k++)
4777 /* multiply current number by 10, without normalizing */
4778 eshup1 (y);
4779 emovz (y, u);
4780 eshup1 (u);
4781 eshup1 (u);
4782 eaddm (u, y);
4783 eiremain (t, y);
4784 *s++ = (char) equot[NI - 1] + '0';
4786 digit = equot[NI - 1];
4787 --s;
4788 ss = s;
4789 /* round off the ASCII string */
4790 if (digit > 4)
4792 /* Test for critical rounding case in ASCII output. */
4793 if (digit == 5)
4795 emovo (y, t);
4796 if (ecmp (t, ezero) != 0)
4797 goto roun; /* round to nearest */
4798 if ((*(s - 1) & 1) == 0)
4799 goto doexp; /* round to even */
4801 /* Round up and propagate carry-outs */
4802 roun:
4803 --s;
4804 k = *s & 0x7f;
4805 /* Carry out to most significant digit? */
4806 if (k == '.')
4808 --s;
4809 k = *s;
4810 k += 1;
4811 *s = (char) k;
4812 /* Most significant digit carries to 10? */
4813 if (k > '9')
4815 expon += 1;
4816 *s = '1';
4818 goto doexp;
4820 /* Round up and carry out from less significant digits */
4821 k += 1;
4822 *s = (char) k;
4823 if (k > '9')
4825 *s = '0';
4826 goto roun;
4829 doexp:
4831 if (expon >= 0)
4832 sprintf (ss, "e+%d", expon);
4833 else
4834 sprintf (ss, "e%d", expon);
4836 sprintf (ss, "e%d", expon);
4837 bxit:
4838 rndprc = rndsav;
4839 /* copy out the working string */
4840 s = string;
4841 ss = wstring;
4842 while (*ss == ' ') /* strip possible leading space */
4843 ++ss;
4844 while ((*s++ = *ss++) != '\0')
4849 /* Convert ASCII string to floating point.
4851 Numeric input is a free format decimal number of any length, with
4852 or without decimal point. Entering E after the number followed by an
4853 integer number causes the second number to be interpreted as a power of
4854 10 to be multiplied by the first number (i.e., "scientific" notation). */
4856 /* Convert ASCII string S to single precision float value Y. */
4858 static void
4859 asctoe24 (s, y)
4860 char *s;
4861 unsigned EMUSHORT *y;
4863 asctoeg (s, y, 24);
4867 /* Convert ASCII string S to double precision value Y. */
4869 static void
4870 asctoe53 (s, y)
4871 char *s;
4872 unsigned EMUSHORT *y;
4874 #if defined(DEC) || defined(IBM)
4875 asctoeg (s, y, 56);
4876 #else
4877 asctoeg (s, y, 53);
4878 #endif
4882 /* Convert ASCII string S to double extended value Y. */
4884 static void
4885 asctoe64 (s, y)
4886 char *s;
4887 unsigned EMUSHORT *y;
4889 asctoeg (s, y, 64);
4892 /* Convert ASCII string S to 128-bit long double Y. */
4894 static void
4895 asctoe113 (s, y)
4896 char *s;
4897 unsigned EMUSHORT *y;
4899 asctoeg (s, y, 113);
4902 /* Convert ASCII string S to e type Y. */
4904 static void
4905 asctoe (s, y)
4906 char *s;
4907 unsigned EMUSHORT *y;
4909 asctoeg (s, y, NBITS);
4912 /* Convert ASCII string SS to e type Y, with a specified rounding precision
4913 of OPREC bits. */
4915 static void
4916 asctoeg (ss, y, oprec)
4917 char *ss;
4918 unsigned EMUSHORT *y;
4919 int oprec;
4921 unsigned EMUSHORT yy[NI], xt[NI], tt[NI];
4922 int esign, decflg, sgnflg, nexp, exp, prec, lost;
4923 int k, trail, c, rndsav;
4924 EMULONG lexp;
4925 unsigned EMUSHORT nsign, *p;
4926 char *sp, *s, *lstr;
4928 /* Copy the input string. */
4929 lstr = (char *) alloca (strlen (ss) + 1);
4930 s = ss;
4931 while (*s == ' ') /* skip leading spaces */
4932 ++s;
4933 sp = lstr;
4934 while ((*sp++ = *s++) != '\0')
4936 s = lstr;
4938 rndsav = rndprc;
4939 rndprc = NBITS; /* Set to full precision */
4940 lost = 0;
4941 nsign = 0;
4942 decflg = 0;
4943 sgnflg = 0;
4944 nexp = 0;
4945 exp = 0;
4946 prec = 0;
4947 ecleaz (yy);
4948 trail = 0;
4950 nxtcom:
4951 k = *s - '0';
4952 if ((k >= 0) && (k <= 9))
4954 /* Ignore leading zeros */
4955 if ((prec == 0) && (decflg == 0) && (k == 0))
4956 goto donchr;
4957 /* Identify and strip trailing zeros after the decimal point. */
4958 if ((trail == 0) && (decflg != 0))
4960 sp = s;
4961 while ((*sp >= '0') && (*sp <= '9'))
4962 ++sp;
4963 /* Check for syntax error */
4964 c = *sp & 0x7f;
4965 if ((c != 'e') && (c != 'E') && (c != '\0')
4966 && (c != '\n') && (c != '\r') && (c != ' ')
4967 && (c != ','))
4968 goto error;
4969 --sp;
4970 while (*sp == '0')
4971 *sp-- = 'z';
4972 trail = 1;
4973 if (*s == 'z')
4974 goto donchr;
4977 /* If enough digits were given to more than fill up the yy register,
4978 continuing until overflow into the high guard word yy[2]
4979 guarantees that there will be a roundoff bit at the top
4980 of the low guard word after normalization. */
4982 if (yy[2] == 0)
4984 if (decflg)
4985 nexp += 1; /* count digits after decimal point */
4986 eshup1 (yy); /* multiply current number by 10 */
4987 emovz (yy, xt);
4988 eshup1 (xt);
4989 eshup1 (xt);
4990 eaddm (xt, yy);
4991 ecleaz (xt);
4992 xt[NI - 2] = (unsigned EMUSHORT) k;
4993 eaddm (xt, yy);
4995 else
4997 /* Mark any lost non-zero digit. */
4998 lost |= k;
4999 /* Count lost digits before the decimal point. */
5000 if (decflg == 0)
5001 nexp -= 1;
5003 prec += 1;
5004 goto donchr;
5007 switch (*s)
5009 case 'z':
5010 break;
5011 case 'E':
5012 case 'e':
5013 goto expnt;
5014 case '.': /* decimal point */
5015 if (decflg)
5016 goto error;
5017 ++decflg;
5018 break;
5019 case '-':
5020 nsign = 0xffff;
5021 if (sgnflg)
5022 goto error;
5023 ++sgnflg;
5024 break;
5025 case '+':
5026 if (sgnflg)
5027 goto error;
5028 ++sgnflg;
5029 break;
5030 case ',':
5031 case ' ':
5032 case '\0':
5033 case '\n':
5034 case '\r':
5035 goto daldone;
5036 case 'i':
5037 case 'I':
5038 goto infinite;
5039 default:
5040 error:
5041 #ifdef NANS
5042 einan (yy);
5043 #else
5044 mtherr ("asctoe", DOMAIN);
5045 eclear (yy);
5046 #endif
5047 goto aexit;
5049 donchr:
5050 ++s;
5051 goto nxtcom;
5053 /* Exponent interpretation */
5054 expnt:
5055 /* 0.0eXXX is zero, regardless of XXX. Check for the 0.0. */
5056 for (k = 0; k < NI; k++)
5058 if (yy[k] != 0)
5059 goto read_expnt;
5061 goto aexit;
5063 read_expnt:
5064 esign = 1;
5065 exp = 0;
5066 ++s;
5067 /* check for + or - */
5068 if (*s == '-')
5070 esign = -1;
5071 ++s;
5073 if (*s == '+')
5074 ++s;
5075 while ((*s >= '0') && (*s <= '9'))
5077 exp *= 10;
5078 exp += *s++ - '0';
5079 if (exp > -(MINDECEXP))
5081 if (esign < 0)
5082 goto zero;
5083 else
5084 goto infinite;
5087 if (esign < 0)
5088 exp = -exp;
5089 if (exp > MAXDECEXP)
5091 infinite:
5092 ecleaz (yy);
5093 yy[E] = 0x7fff; /* infinity */
5094 goto aexit;
5096 if (exp < MINDECEXP)
5098 zero:
5099 ecleaz (yy);
5100 goto aexit;
5103 daldone:
5104 nexp = exp - nexp;
5105 /* Pad trailing zeros to minimize power of 10, per IEEE spec. */
5106 while ((nexp > 0) && (yy[2] == 0))
5108 emovz (yy, xt);
5109 eshup1 (xt);
5110 eshup1 (xt);
5111 eaddm (yy, xt);
5112 eshup1 (xt);
5113 if (xt[2] != 0)
5114 break;
5115 nexp -= 1;
5116 emovz (xt, yy);
5118 if ((k = enormlz (yy)) > NBITS)
5120 ecleaz (yy);
5121 goto aexit;
5123 lexp = (EXONE - 1 + NBITS) - k;
5124 emdnorm (yy, lost, 0, lexp, 64);
5126 /* Convert to external format:
5128 Multiply by 10**nexp. If precision is 64 bits,
5129 the maximum relative error incurred in forming 10**n
5130 for 0 <= n <= 324 is 8.2e-20, at 10**180.
5131 For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
5132 For 0 >= n >= -999, it is -1.55e-19 at 10**-435. */
5134 lexp = yy[E];
5135 if (nexp == 0)
5137 k = 0;
5138 goto expdon;
5140 esign = 1;
5141 if (nexp < 0)
5143 nexp = -nexp;
5144 esign = -1;
5145 if (nexp > 4096)
5147 /* Punt. Can't handle this without 2 divides. */
5148 emovi (etens[0], tt);
5149 lexp -= tt[E];
5150 k = edivm (tt, yy);
5151 lexp += EXONE;
5152 nexp -= 4096;
5155 p = &etens[NTEN][0];
5156 emov (eone, xt);
5157 exp = 1;
5160 if (exp & nexp)
5161 emul (p, xt, xt);
5162 p -= NE;
5163 exp = exp + exp;
5165 while (exp <= MAXP);
5167 emovi (xt, tt);
5168 if (esign < 0)
5170 lexp -= tt[E];
5171 k = edivm (tt, yy);
5172 lexp += EXONE;
5174 else
5176 lexp += tt[E];
5177 k = emulm (tt, yy);
5178 lexp -= EXONE - 1;
5181 expdon:
5183 /* Round and convert directly to the destination type */
5184 if (oprec == 53)
5185 lexp -= EXONE - 0x3ff;
5186 #ifdef IBM
5187 else if (oprec == 24 || oprec == 56)
5188 lexp -= EXONE - (0x41 << 2);
5189 #else
5190 else if (oprec == 24)
5191 lexp -= EXONE - 0177;
5192 #endif
5193 #ifdef DEC
5194 else if (oprec == 56)
5195 lexp -= EXONE - 0201;
5196 #endif
5197 rndprc = oprec;
5198 emdnorm (yy, k, 0, lexp, 64);
5200 aexit:
5202 rndprc = rndsav;
5203 yy[0] = nsign;
5204 switch (oprec)
5206 #ifdef DEC
5207 case 56:
5208 todec (yy, y); /* see etodec.c */
5209 break;
5210 #endif
5211 #ifdef IBM
5212 case 56:
5213 toibm (yy, y, DFmode);
5214 break;
5215 #endif
5216 case 53:
5217 toe53 (yy, y);
5218 break;
5219 case 24:
5220 toe24 (yy, y);
5221 break;
5222 case 64:
5223 toe64 (yy, y);
5224 break;
5225 case 113:
5226 toe113 (yy, y);
5227 break;
5228 case NBITS:
5229 emovo (yy, y);
5230 break;
5236 /* Return Y = largest integer not greater than X (truncated toward minus
5237 infinity). */
5239 static unsigned EMUSHORT bmask[] =
5241 0xffff,
5242 0xfffe,
5243 0xfffc,
5244 0xfff8,
5245 0xfff0,
5246 0xffe0,
5247 0xffc0,
5248 0xff80,
5249 0xff00,
5250 0xfe00,
5251 0xfc00,
5252 0xf800,
5253 0xf000,
5254 0xe000,
5255 0xc000,
5256 0x8000,
5257 0x0000,
5260 static void
5261 efloor (x, y)
5262 unsigned EMUSHORT x[], y[];
5264 register unsigned EMUSHORT *p;
5265 int e, expon, i;
5266 unsigned EMUSHORT f[NE];
5268 emov (x, f); /* leave in external format */
5269 expon = (int) f[NE - 1];
5270 e = (expon & 0x7fff) - (EXONE - 1);
5271 if (e <= 0)
5273 eclear (y);
5274 goto isitneg;
5276 /* number of bits to clear out */
5277 e = NBITS - e;
5278 emov (f, y);
5279 if (e <= 0)
5280 return;
5282 p = &y[0];
5283 while (e >= 16)
5285 *p++ = 0;
5286 e -= 16;
5288 /* clear the remaining bits */
5289 *p &= bmask[e];
5290 /* truncate negatives toward minus infinity */
5291 isitneg:
5293 if ((unsigned EMUSHORT) expon & (unsigned EMUSHORT) 0x8000)
5295 for (i = 0; i < NE - 1; i++)
5297 if (f[i] != y[i])
5299 esub (eone, y, y);
5300 break;
5307 /* Return S and EXP such that S * 2^EXP = X and .5 <= S < 1.
5308 For example, 1.1 = 0.55 * 2^1. */
5310 static void
5311 efrexp (x, exp, s)
5312 unsigned EMUSHORT x[];
5313 int *exp;
5314 unsigned EMUSHORT s[];
5316 unsigned EMUSHORT xi[NI];
5317 EMULONG li;
5319 emovi (x, xi);
5320 /* Handle denormalized numbers properly using long integer exponent. */
5321 li = (EMULONG) ((EMUSHORT) xi[1]);
5323 if (li == 0)
5325 li -= enormlz (xi);
5327 xi[1] = 0x3ffe;
5328 emovo (xi, s);
5329 *exp = (int) (li - 0x3ffe);
5332 /* Return e type Y = X * 2^PWR2. */
5334 static void
5335 eldexp (x, pwr2, y)
5336 unsigned EMUSHORT x[];
5337 int pwr2;
5338 unsigned EMUSHORT y[];
5340 unsigned EMUSHORT xi[NI];
5341 EMULONG li;
5342 int i;
5344 emovi (x, xi);
5345 li = xi[1];
5346 li += pwr2;
5347 i = 0;
5348 emdnorm (xi, i, i, li, 64);
5349 emovo (xi, y);
5353 /* C = remainder after dividing B by A, all e type values.
5354 Least significant integer quotient bits left in EQUOT. */
5356 static void
5357 eremain (a, b, c)
5358 unsigned EMUSHORT a[], b[], c[];
5360 unsigned EMUSHORT den[NI], num[NI];
5362 #ifdef NANS
5363 if (eisinf (b)
5364 || (ecmp (a, ezero) == 0)
5365 || eisnan (a)
5366 || eisnan (b))
5368 enan (c, 0);
5369 return;
5371 #endif
5372 if (ecmp (a, ezero) == 0)
5374 mtherr ("eremain", SING);
5375 eclear (c);
5376 return;
5378 emovi (a, den);
5379 emovi (b, num);
5380 eiremain (den, num);
5381 /* Sign of remainder = sign of quotient */
5382 if (a[0] == b[0])
5383 num[0] = 0;
5384 else
5385 num[0] = 0xffff;
5386 emovo (num, c);
5389 /* Return quotient of exploded e-types NUM / DEN in EQUOT,
5390 remainder in NUM. */
5392 static void
5393 eiremain (den, num)
5394 unsigned EMUSHORT den[], num[];
5396 EMULONG ld, ln;
5397 unsigned EMUSHORT j;
5399 ld = den[E];
5400 ld -= enormlz (den);
5401 ln = num[E];
5402 ln -= enormlz (num);
5403 ecleaz (equot);
5404 while (ln >= ld)
5406 if (ecmpm (den, num) <= 0)
5408 esubm (den, num);
5409 j = 1;
5411 else
5412 j = 0;
5413 eshup1 (equot);
5414 equot[NI - 1] |= j;
5415 eshup1 (num);
5416 ln -= 1;
5418 emdnorm (num, 0, 0, ln, 0);
5421 /* Report an error condition CODE encountered in function NAME.
5422 CODE is one of the following:
5424 Mnemonic Value Significance
5426 DOMAIN 1 argument domain error
5427 SING 2 function singularity
5428 OVERFLOW 3 overflow range error
5429 UNDERFLOW 4 underflow range error
5430 TLOSS 5 total loss of precision
5431 PLOSS 6 partial loss of precision
5432 INVALID 7 NaN - producing operation
5433 EDOM 33 Unix domain error code
5434 ERANGE 34 Unix range error code
5436 The order of appearance of the following messages is bound to the
5437 error codes defined above. */
5439 #define NMSGS 8
5440 static char *ermsg[NMSGS] =
5442 "unknown", /* error code 0 */
5443 "domain", /* error code 1 */
5444 "singularity", /* et seq. */
5445 "overflow",
5446 "underflow",
5447 "total loss of precision",
5448 "partial loss of precision",
5449 "invalid operation"
5452 int merror = 0;
5453 extern int merror;
5455 static void
5456 mtherr (name, code)
5457 char *name;
5458 int code;
5460 char errstr[80];
5462 /* The string passed by the calling program is supposed to be the
5463 name of the function in which the error occurred.
5464 The code argument selects which error message string will be printed. */
5466 if ((code <= 0) || (code >= NMSGS))
5467 code = 0;
5468 sprintf (errstr, " %s %s error", name, ermsg[code]);
5469 if (extra_warnings)
5470 warning (errstr);
5471 /* Set global error message word */
5472 merror = code + 1;
5475 #ifdef DEC
5476 /* Convert DEC double precision D to e type E. */
5478 static void
5479 dectoe (d, e)
5480 unsigned EMUSHORT *d;
5481 unsigned EMUSHORT *e;
5483 unsigned EMUSHORT y[NI];
5484 register unsigned EMUSHORT r, *p;
5486 ecleaz (y); /* start with a zero */
5487 p = y; /* point to our number */
5488 r = *d; /* get DEC exponent word */
5489 if (*d & (unsigned int) 0x8000)
5490 *p = 0xffff; /* fill in our sign */
5491 ++p; /* bump pointer to our exponent word */
5492 r &= 0x7fff; /* strip the sign bit */
5493 if (r == 0) /* answer = 0 if high order DEC word = 0 */
5494 goto done;
5497 r >>= 7; /* shift exponent word down 7 bits */
5498 r += EXONE - 0201; /* subtract DEC exponent offset */
5499 /* add our e type exponent offset */
5500 *p++ = r; /* to form our exponent */
5502 r = *d++; /* now do the high order mantissa */
5503 r &= 0177; /* strip off the DEC exponent and sign bits */
5504 r |= 0200; /* the DEC understood high order mantissa bit */
5505 *p++ = r; /* put result in our high guard word */
5507 *p++ = *d++; /* fill in the rest of our mantissa */
5508 *p++ = *d++;
5509 *p = *d;
5511 eshdn8 (y); /* shift our mantissa down 8 bits */
5512 done:
5513 emovo (y, e);
5516 /* Convert e type X to DEC double precision D. */
5518 static void
5519 etodec (x, d)
5520 unsigned EMUSHORT *x, *d;
5522 unsigned EMUSHORT xi[NI];
5523 EMULONG exp;
5524 int rndsav;
5526 emovi (x, xi);
5527 /* Adjust exponent for offsets. */
5528 exp = (EMULONG) xi[E] - (EXONE - 0201);
5529 /* Round off to nearest or even. */
5530 rndsav = rndprc;
5531 rndprc = 56;
5532 emdnorm (xi, 0, 0, exp, 64);
5533 rndprc = rndsav;
5534 todec (xi, d);
5537 /* Convert exploded e-type X, that has already been rounded to
5538 56-bit precision, to DEC format double Y. */
5540 static void
5541 todec (x, y)
5542 unsigned EMUSHORT *x, *y;
5544 unsigned EMUSHORT i;
5545 unsigned EMUSHORT *p;
5547 p = x;
5548 *y = 0;
5549 if (*p++)
5550 *y = 0100000;
5551 i = *p++;
5552 if (i == 0)
5554 *y++ = 0;
5555 *y++ = 0;
5556 *y++ = 0;
5557 *y++ = 0;
5558 return;
5560 if (i > 0377)
5562 *y++ |= 077777;
5563 *y++ = 0xffff;
5564 *y++ = 0xffff;
5565 *y++ = 0xffff;
5566 #ifdef ERANGE
5567 errno = ERANGE;
5568 #endif
5569 return;
5571 i &= 0377;
5572 i <<= 7;
5573 eshup8 (x);
5574 x[M] &= 0177;
5575 i |= x[M];
5576 *y++ |= i;
5577 *y++ = x[M + 1];
5578 *y++ = x[M + 2];
5579 *y++ = x[M + 3];
5581 #endif /* DEC */
5583 #ifdef IBM
5584 /* Convert IBM single/double precision to e type. */
5586 static void
5587 ibmtoe (d, e, mode)
5588 unsigned EMUSHORT *d;
5589 unsigned EMUSHORT *e;
5590 enum machine_mode mode;
5592 unsigned EMUSHORT y[NI];
5593 register unsigned EMUSHORT r, *p;
5594 int rndsav;
5596 ecleaz (y); /* start with a zero */
5597 p = y; /* point to our number */
5598 r = *d; /* get IBM exponent word */
5599 if (*d & (unsigned int) 0x8000)
5600 *p = 0xffff; /* fill in our sign */
5601 ++p; /* bump pointer to our exponent word */
5602 r &= 0x7f00; /* strip the sign bit */
5603 r >>= 6; /* shift exponent word down 6 bits */
5604 /* in fact shift by 8 right and 2 left */
5605 r += EXONE - (0x41 << 2); /* subtract IBM exponent offset */
5606 /* add our e type exponent offset */
5607 *p++ = r; /* to form our exponent */
5609 *p++ = *d++ & 0xff; /* now do the high order mantissa */
5610 /* strip off the IBM exponent and sign bits */
5611 if (mode != SFmode) /* there are only 2 words in SFmode */
5613 *p++ = *d++; /* fill in the rest of our mantissa */
5614 *p++ = *d++;
5616 *p = *d;
5618 if (y[M] == 0 && y[M+1] == 0 && y[M+2] == 0 && y[M+3] == 0)
5619 y[0] = y[E] = 0;
5620 else
5621 y[E] -= 5 + enormlz (y); /* now normalise the mantissa */
5622 /* handle change in RADIX */
5623 emovo (y, e);
5628 /* Convert e type to IBM single/double precision. */
5630 static void
5631 etoibm (x, d, mode)
5632 unsigned EMUSHORT *x, *d;
5633 enum machine_mode mode;
5635 unsigned EMUSHORT xi[NI];
5636 EMULONG exp;
5637 int rndsav;
5639 emovi (x, xi);
5640 exp = (EMULONG) xi[E] - (EXONE - (0x41 << 2)); /* adjust exponent for offsets */
5641 /* round off to nearest or even */
5642 rndsav = rndprc;
5643 rndprc = 56;
5644 emdnorm (xi, 0, 0, exp, 64);
5645 rndprc = rndsav;
5646 toibm (xi, d, mode);
5649 static void
5650 toibm (x, y, mode)
5651 unsigned EMUSHORT *x, *y;
5652 enum machine_mode mode;
5654 unsigned EMUSHORT i;
5655 unsigned EMUSHORT *p;
5656 int r;
5658 p = x;
5659 *y = 0;
5660 if (*p++)
5661 *y = 0x8000;
5662 i = *p++;
5663 if (i == 0)
5665 *y++ = 0;
5666 *y++ = 0;
5667 if (mode != SFmode)
5669 *y++ = 0;
5670 *y++ = 0;
5672 return;
5674 r = i & 0x3;
5675 i >>= 2;
5676 if (i > 0x7f)
5678 *y++ |= 0x7fff;
5679 *y++ = 0xffff;
5680 if (mode != SFmode)
5682 *y++ = 0xffff;
5683 *y++ = 0xffff;
5685 #ifdef ERANGE
5686 errno = ERANGE;
5687 #endif
5688 return;
5690 i &= 0x7f;
5691 *y |= (i << 8);
5692 eshift (x, r + 5);
5693 *y++ |= x[M];
5694 *y++ = x[M + 1];
5695 if (mode != SFmode)
5697 *y++ = x[M + 2];
5698 *y++ = x[M + 3];
5701 #endif /* IBM */
5703 /* Output a binary NaN bit pattern in the target machine's format. */
5705 /* If special NaN bit patterns are required, define them in tm.h
5706 as arrays of unsigned 16-bit shorts. Otherwise, use the default
5707 patterns here. */
5708 #ifdef TFMODE_NAN
5709 TFMODE_NAN;
5710 #else
5711 #ifdef IEEE
5712 unsigned EMUSHORT TFbignan[8] =
5713 {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
5714 unsigned EMUSHORT TFlittlenan[8] = {0, 0, 0, 0, 0, 0, 0x8000, 0xffff};
5715 #endif
5716 #endif
5718 #ifdef XFMODE_NAN
5719 XFMODE_NAN;
5720 #else
5721 #ifdef IEEE
5722 unsigned EMUSHORT XFbignan[6] =
5723 {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
5724 unsigned EMUSHORT XFlittlenan[6] = {0, 0, 0, 0xc000, 0xffff, 0};
5725 #endif
5726 #endif
5728 #ifdef DFMODE_NAN
5729 DFMODE_NAN;
5730 #else
5731 #ifdef IEEE
5732 unsigned EMUSHORT DFbignan[4] = {0x7fff, 0xffff, 0xffff, 0xffff};
5733 unsigned EMUSHORT DFlittlenan[4] = {0, 0, 0, 0xfff8};
5734 #endif
5735 #endif
5737 #ifdef SFMODE_NAN
5738 SFMODE_NAN;
5739 #else
5740 #ifdef IEEE
5741 unsigned EMUSHORT SFbignan[2] = {0x7fff, 0xffff};
5742 unsigned EMUSHORT SFlittlenan[2] = {0, 0xffc0};
5743 #endif
5744 #endif
5747 static void
5748 make_nan (nan, sign, mode)
5749 unsigned EMUSHORT *nan;
5750 int sign;
5751 enum machine_mode mode;
5753 int n;
5754 unsigned EMUSHORT *p;
5756 switch (mode)
5758 /* Possibly the `reserved operand' patterns on a VAX can be
5759 used like NaN's, but probably not in the same way as IEEE. */
5760 #if !defined(DEC) && !defined(IBM)
5761 case TFmode:
5762 n = 8;
5763 if (REAL_WORDS_BIG_ENDIAN)
5764 p = TFbignan;
5765 else
5766 p = TFlittlenan;
5767 break;
5768 case XFmode:
5769 n = 6;
5770 if (REAL_WORDS_BIG_ENDIAN)
5771 p = XFbignan;
5772 else
5773 p = XFlittlenan;
5774 break;
5775 case DFmode:
5776 n = 4;
5777 if (REAL_WORDS_BIG_ENDIAN)
5778 p = DFbignan;
5779 else
5780 p = DFlittlenan;
5781 break;
5782 case HFmode:
5783 case SFmode:
5784 n = 2;
5785 if (REAL_WORDS_BIG_ENDIAN)
5786 p = SFbignan;
5787 else
5788 p = SFlittlenan;
5789 break;
5790 #endif
5791 default:
5792 abort ();
5794 if (REAL_WORDS_BIG_ENDIAN)
5795 *nan++ = (sign << 15) | *p++;
5796 while (--n != 0)
5797 *nan++ = *p++;
5798 if (! REAL_WORDS_BIG_ENDIAN)
5799 *nan = (sign << 15) | *p;
5802 /* Convert an SFmode target `float' value to a REAL_VALUE_TYPE.
5803 This is the inverse of the function `etarsingle' invoked by
5804 REAL_VALUE_TO_TARGET_SINGLE. */
5806 REAL_VALUE_TYPE
5807 ereal_from_float (f)
5808 HOST_WIDE_INT f;
5810 REAL_VALUE_TYPE r;
5811 unsigned EMUSHORT s[2];
5812 unsigned EMUSHORT e[NE];
5814 /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
5815 This is the inverse operation to what the function `endian' does. */
5816 if (REAL_WORDS_BIG_ENDIAN)
5818 s[0] = (unsigned EMUSHORT) (f >> 16);
5819 s[1] = (unsigned EMUSHORT) f;
5821 else
5823 s[0] = (unsigned EMUSHORT) f;
5824 s[1] = (unsigned EMUSHORT) (f >> 16);
5826 /* Convert and promote the target float to E-type. */
5827 e24toe (s, e);
5828 /* Output E-type to REAL_VALUE_TYPE. */
5829 PUT_REAL (e, &r);
5830 return r;
5834 /* Convert a DFmode target `double' value to a REAL_VALUE_TYPE.
5835 This is the inverse of the function `etardouble' invoked by
5836 REAL_VALUE_TO_TARGET_DOUBLE.
5838 The DFmode is stored as an array of HOST_WIDE_INT in the target's
5839 data format, with no holes in the bit packing. The first element
5840 of the input array holds the bits that would come first in the
5841 target computer's memory. */
5843 REAL_VALUE_TYPE
5844 ereal_from_double (d)
5845 HOST_WIDE_INT d[];
5847 REAL_VALUE_TYPE r;
5848 unsigned EMUSHORT s[4];
5849 unsigned EMUSHORT e[NE];
5851 /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces. */
5852 if (REAL_WORDS_BIG_ENDIAN)
5854 s[0] = (unsigned EMUSHORT) (d[0] >> 16);
5855 s[1] = (unsigned EMUSHORT) d[0];
5856 #if HOST_BITS_PER_WIDE_INT == 32
5857 s[2] = (unsigned EMUSHORT) (d[1] >> 16);
5858 s[3] = (unsigned EMUSHORT) d[1];
5859 #else
5860 /* In this case the entire target double is contained in the
5861 first array element. The second element of the input is
5862 ignored. */
5863 s[2] = (unsigned EMUSHORT) (d[0] >> 48);
5864 s[3] = (unsigned EMUSHORT) (d[0] >> 32);
5865 #endif
5867 else
5869 /* Target float words are little-endian. */
5870 s[0] = (unsigned EMUSHORT) d[0];
5871 s[1] = (unsigned EMUSHORT) (d[0] >> 16);
5872 #if HOST_BITS_PER_WIDE_INT == 32
5873 s[2] = (unsigned EMUSHORT) d[1];
5874 s[3] = (unsigned EMUSHORT) (d[1] >> 16);
5875 #else
5876 s[2] = (unsigned EMUSHORT) (d[0] >> 32);
5877 s[3] = (unsigned EMUSHORT) (d[0] >> 48);
5878 #endif
5880 /* Convert target double to E-type. */
5881 e53toe (s, e);
5882 /* Output E-type to REAL_VALUE_TYPE. */
5883 PUT_REAL (e, &r);
5884 return r;
5888 /* Convert target computer unsigned 64-bit integer to e-type.
5889 The endian-ness of DImode follows the convention for integers,
5890 so we use WORDS_BIG_ENDIAN here, not REAL_WORDS_BIG_ENDIAN. */
5892 static void
5893 uditoe (di, e)
5894 unsigned EMUSHORT *di; /* Address of the 64-bit int. */
5895 unsigned EMUSHORT *e;
5897 unsigned EMUSHORT yi[NI];
5898 int k;
5900 ecleaz (yi);
5901 if (WORDS_BIG_ENDIAN)
5903 for (k = M; k < M + 4; k++)
5904 yi[k] = *di++;
5906 else
5908 for (k = M + 3; k >= M; k--)
5909 yi[k] = *di++;
5911 yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
5912 if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
5913 ecleaz (yi); /* it was zero */
5914 else
5915 yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
5916 emovo (yi, e);
5919 /* Convert target computer signed 64-bit integer to e-type. */
5921 static void
5922 ditoe (di, e)
5923 unsigned EMUSHORT *di; /* Address of the 64-bit int. */
5924 unsigned EMUSHORT *e;
5926 unsigned EMULONG acc;
5927 unsigned EMUSHORT yi[NI];
5928 unsigned EMUSHORT carry;
5929 int k, sign;
5931 ecleaz (yi);
5932 if (WORDS_BIG_ENDIAN)
5934 for (k = M; k < M + 4; k++)
5935 yi[k] = *di++;
5937 else
5939 for (k = M + 3; k >= M; k--)
5940 yi[k] = *di++;
5942 /* Take absolute value */
5943 sign = 0;
5944 if (yi[M] & 0x8000)
5946 sign = 1;
5947 carry = 0;
5948 for (k = M + 3; k >= M; k--)
5950 acc = (unsigned EMULONG) (~yi[k] & 0xffff) + carry;
5951 yi[k] = acc;
5952 carry = 0;
5953 if (acc & 0x10000)
5954 carry = 1;
5957 yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
5958 if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
5959 ecleaz (yi); /* it was zero */
5960 else
5961 yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
5962 emovo (yi, e);
5963 if (sign)
5964 eneg (e);
5968 /* Convert e-type to unsigned 64-bit int. */
5970 static void
5971 etoudi (x, i)
5972 unsigned EMUSHORT *x;
5973 unsigned EMUSHORT *i;
5975 unsigned EMUSHORT xi[NI];
5976 int j, k;
5978 emovi (x, xi);
5979 if (xi[0])
5981 xi[M] = 0;
5982 goto noshift;
5984 k = (int) xi[E] - (EXONE - 1);
5985 if (k <= 0)
5987 for (j = 0; j < 4; j++)
5988 *i++ = 0;
5989 return;
5991 if (k > 64)
5993 for (j = 0; j < 4; j++)
5994 *i++ = 0xffff;
5995 if (extra_warnings)
5996 warning ("overflow on truncation to integer");
5997 return;
5999 if (k > 16)
6001 /* Shift more than 16 bits: first shift up k-16 mod 16,
6002 then shift up by 16's. */
6003 j = k - ((k >> 4) << 4);
6004 if (j == 0)
6005 j = 16;
6006 eshift (xi, j);
6007 if (WORDS_BIG_ENDIAN)
6008 *i++ = xi[M];
6009 else
6011 i += 3;
6012 *i-- = xi[M];
6014 k -= j;
6017 eshup6 (xi);
6018 if (WORDS_BIG_ENDIAN)
6019 *i++ = xi[M];
6020 else
6021 *i-- = xi[M];
6023 while ((k -= 16) > 0);
6025 else
6027 /* shift not more than 16 bits */
6028 eshift (xi, k);
6030 noshift:
6032 if (WORDS_BIG_ENDIAN)
6034 i += 3;
6035 *i-- = xi[M];
6036 *i-- = 0;
6037 *i-- = 0;
6038 *i = 0;
6040 else
6042 *i++ = xi[M];
6043 *i++ = 0;
6044 *i++ = 0;
6045 *i = 0;
6051 /* Convert e-type to signed 64-bit int. */
6053 static void
6054 etodi (x, i)
6055 unsigned EMUSHORT *x;
6056 unsigned EMUSHORT *i;
6058 unsigned EMULONG acc;
6059 unsigned EMUSHORT xi[NI];
6060 unsigned EMUSHORT carry;
6061 unsigned EMUSHORT *isave;
6062 int j, k;
6064 emovi (x, xi);
6065 k = (int) xi[E] - (EXONE - 1);
6066 if (k <= 0)
6068 for (j = 0; j < 4; j++)
6069 *i++ = 0;
6070 return;
6072 if (k > 64)
6074 for (j = 0; j < 4; j++)
6075 *i++ = 0xffff;
6076 if (extra_warnings)
6077 warning ("overflow on truncation to integer");
6078 return;
6080 isave = i;
6081 if (k > 16)
6083 /* Shift more than 16 bits: first shift up k-16 mod 16,
6084 then shift up by 16's. */
6085 j = k - ((k >> 4) << 4);
6086 if (j == 0)
6087 j = 16;
6088 eshift (xi, j);
6089 if (WORDS_BIG_ENDIAN)
6090 *i++ = xi[M];
6091 else
6093 i += 3;
6094 *i-- = xi[M];
6096 k -= j;
6099 eshup6 (xi);
6100 if (WORDS_BIG_ENDIAN)
6101 *i++ = xi[M];
6102 else
6103 *i-- = xi[M];
6105 while ((k -= 16) > 0);
6107 else
6109 /* shift not more than 16 bits */
6110 eshift (xi, k);
6112 if (WORDS_BIG_ENDIAN)
6114 i += 3;
6115 *i = xi[M];
6116 *i-- = 0;
6117 *i-- = 0;
6118 *i = 0;
6120 else
6122 *i++ = xi[M];
6123 *i++ = 0;
6124 *i++ = 0;
6125 *i = 0;
6128 /* Negate if negative */
6129 if (xi[0])
6131 carry = 0;
6132 if (WORDS_BIG_ENDIAN)
6133 isave += 3;
6134 for (k = 0; k < 4; k++)
6136 acc = (unsigned EMULONG) (~(*isave) & 0xffff) + carry;
6137 if (WORDS_BIG_ENDIAN)
6138 *isave-- = acc;
6139 else
6140 *isave++ = acc;
6141 carry = 0;
6142 if (acc & 0x10000)
6143 carry = 1;
6149 /* Longhand square root routine. */
6152 static int esqinited = 0;
6153 static unsigned short sqrndbit[NI];
6155 static void
6156 esqrt (x, y)
6157 unsigned EMUSHORT *x, *y;
6159 unsigned EMUSHORT temp[NI], num[NI], sq[NI], xx[NI];
6160 EMULONG m, exp;
6161 int i, j, k, n, nlups;
6163 if (esqinited == 0)
6165 ecleaz (sqrndbit);
6166 sqrndbit[NI - 2] = 1;
6167 esqinited = 1;
6169 /* Check for arg <= 0 */
6170 i = ecmp (x, ezero);
6171 if (i <= 0)
6173 if (i == -1)
6175 mtherr ("esqrt", DOMAIN);
6176 eclear (y);
6178 else
6179 emov (x, y);
6180 return;
6183 #ifdef INFINITY
6184 if (eisinf (x))
6186 eclear (y);
6187 einfin (y);
6188 return;
6190 #endif
6191 /* Bring in the arg and renormalize if it is denormal. */
6192 emovi (x, xx);
6193 m = (EMULONG) xx[1]; /* local long word exponent */
6194 if (m == 0)
6195 m -= enormlz (xx);
6197 /* Divide exponent by 2 */
6198 m -= 0x3ffe;
6199 exp = (unsigned short) ((m / 2) + 0x3ffe);
6201 /* Adjust if exponent odd */
6202 if ((m & 1) != 0)
6204 if (m > 0)
6205 exp += 1;
6206 eshdn1 (xx);
6209 ecleaz (sq);
6210 ecleaz (num);
6211 n = 8; /* get 8 bits of result per inner loop */
6212 nlups = rndprc;
6213 j = 0;
6215 while (nlups > 0)
6217 /* bring in next word of arg */
6218 if (j < NE)
6219 num[NI - 1] = xx[j + 3];
6220 /* Do additional bit on last outer loop, for roundoff. */
6221 if (nlups <= 8)
6222 n = nlups + 1;
6223 for (i = 0; i < n; i++)
6225 /* Next 2 bits of arg */
6226 eshup1 (num);
6227 eshup1 (num);
6228 /* Shift up answer */
6229 eshup1 (sq);
6230 /* Make trial divisor */
6231 for (k = 0; k < NI; k++)
6232 temp[k] = sq[k];
6233 eshup1 (temp);
6234 eaddm (sqrndbit, temp);
6235 /* Subtract and insert answer bit if it goes in */
6236 if (ecmpm (temp, num) <= 0)
6238 esubm (temp, num);
6239 sq[NI - 2] |= 1;
6242 nlups -= n;
6243 j += 1;
6246 /* Adjust for extra, roundoff loop done. */
6247 exp += (NBITS - 1) - rndprc;
6249 /* Sticky bit = 1 if the remainder is nonzero. */
6250 k = 0;
6251 for (i = 3; i < NI; i++)
6252 k |= (int) num[i];
6254 /* Renormalize and round off. */
6255 emdnorm (sq, k, 0, exp, 64);
6256 emovo (sq, y);
6258 #endif /* EMU_NON_COMPILE not defined */
6260 /* Return the binary precision of the significand for a given
6261 floating point mode. The mode can hold an integer value
6262 that many bits wide, without losing any bits. */
6265 significand_size (mode)
6266 enum machine_mode mode;
6269 /* Don't test the modes, but their sizes, lest this
6270 code won't work for BITS_PER_UNIT != 8 . */
6272 switch (GET_MODE_BITSIZE (mode))
6274 case 32:
6275 return 24;
6277 case 64:
6278 #if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
6279 return 53;
6280 #else
6281 #if TARGET_FLOAT_FORMAT == IBM_FLOAT_FORMAT
6282 return 56;
6283 #else
6284 #if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
6285 return 56;
6286 #else
6287 abort ();
6288 #endif
6289 #endif
6290 #endif
6292 case 96:
6293 return 64;
6294 case 128:
6295 return 113;
6297 default:
6298 abort ();