]> gcc.gnu.org Git - gcc.git/blame_incremental - gcc/real.c
git://gcc.gnu.org / gcc.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
90th Cygnus<->FSF quick merge
[gcc.git] / gcc / real.c
... / ...
CommitLineData
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).
5
6This file is part of GNU CC.
7
8GNU CC is free software; you can redistribute it and/or modify
9it under the terms of the GNU General Public License as published by
10the Free Software Foundation; either version 2, or (at your option)
11any later version.
12
13GNU CC is distributed in the hope that it will be useful,
14but WITHOUT ANY WARRANTY; without even the implied warranty of
15MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16GNU General Public License for more details.
17
18You should have received a copy of the GNU General Public License
19along with GNU CC; see the file COPYING. If not, write to
20the Free Software Foundation, 59 Temple Place - Suite 330,
21Boston, MA 02111-1307, USA. */
22
23#include <stdio.h>
24#include <errno.h>
25#include "config.h"
26#include "tree.h"
27
28#ifndef errno
29extern int errno;
30#endif
31
32/* To enable support of XFmode extended real floating point, define
33LONG_DOUBLE_TYPE_SIZE 96 in the tm.h file (m68k.h or i386.h).
34
35To support cross compilation between IEEE, VAX and IBM floating
36point formats, define REAL_ARITHMETIC in the tm.h file.
37
38In either case the machine files (tm.h) must not contain any code
39that tries to use host floating point arithmetic to convert
40REAL_VALUE_TYPEs from `double' to `float', pass them to fprintf,
41etc. In cross-compile situations a REAL_VALUE_TYPE may not
42be intelligible to the host computer's native arithmetic.
43
44The emulator defaults to the host's floating point format so that
45its decimal conversion functions can be used if desired (see
46real.h).
47
48The first part of this file interfaces gcc to a floating point
49arithmetic suite that was not written with gcc in mind. Avoid
50changing the low-level arithmetic routines unless you have suitable
51test programs available. A special version of the PARANOIA floating
52point arithmetic tester, modified for this purpose, can be found on
53usc.edu: /pub/C-numanal/ieeetest.zoo. Other tests, and libraries of
54XFmode and TFmode transcendental functions, can be obtained by ftp from
55netlib.att.com: netlib/cephes. */
56\f
57/* Type of computer arithmetic.
58 Only one of DEC, IBM, IEEE, or UNK should get defined.
59
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.
67
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.
73
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.
77
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).
82
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.
89
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.
93
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. */
104
105
106/* The following converts gcc macros into the ones used by this file. */
107
108/* REAL_ARITHMETIC defined means that macros in real.h are
109 defined to call emulator functions. */
110#ifdef REAL_ARITHMETIC
111
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. */
124unknown arithmetic type
125#define UNK 1
126#endif /* not IEEE */
127#endif /* not IBM */
128#endif /* not VAX */
129
130#define REAL_WORDS_BIG_ENDIAN FLOAT_WORDS_BIG_ENDIAN
131
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. */
138
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 */
149unknown arithmetic type
150#define UNK 1
151#endif /* not IEEE */
152#endif /* not IBM */
153#endif /* not VAX */
154
155#define REAL_WORDS_BIG_ENDIAN HOST_FLOAT_WORDS_BIG_ENDIAN
156
157#endif /* REAL_ARITHMETIC not defined */
158
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
165
166/* Support of NaNs requires support of infinity. */
167#ifdef NANS
168#ifndef INFINITY
169#define INFINITY
170#endif
171#endif
172\f
173/* Find a host integer type that is at least 16 bits wide,
174 and another type at least twice whatever that size is. */
175
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
202
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
221
222
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
227
228/* OK to continue compilation. */
229#ifndef EMU_NON_COMPILE
230
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. */
235
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. */
257
258#define GET_REAL(r,e) \
259do { \
260 if (HOST_FLOAT_WORDS_BIG_ENDIAN == REAL_WORDS_BIG_ENDIAN) \
261 e53toe ((unsigned EMUSHORT *) (r), (e)); \
262 else \
263 { \
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)); \
270 } \
271 } while (0)
272
273#define PUT_REAL(e,r) \
274do { \
275 if (HOST_FLOAT_WORDS_BIG_ENDIAN == REAL_WORDS_BIG_ENDIAN) \
276 etoe53 ((e), (unsigned EMUSHORT *) (r)); \
277 else \
278 { \
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]; \
285 } \
286 } while (0)
287
288#else /* not REAL_ARITHMETIC */
289
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))
293
294#endif /* not REAL_ARITHMETIC */
295#endif /* not TFmode */
296#endif /* no XFmode */
297
298
299/* Number of 16 bit words in internal format */
300#define NI (NE+3)
301
302/* Array offset to exponent */
303#define E 1
304
305/* Array offset to high guard word */
306#define M 2
307
308/* Number of bits of precision */
309#define NBITS ((NI-4)*16)
310
311/* Maximum number of decimal digits in ASCII conversion
312 * = NBITS*log10(2)
313 */
314#define NDEC (NBITS*8/27)
315
316/* The exponent of 1.0 */
317#define EXONE (0x3fff)
318
319extern int extra_warnings;
320extern unsigned EMUSHORT ezero[], ehalf[], eone[], etwo[];
321extern unsigned EMUSHORT elog2[], esqrt2[];
322
323static void endian PROTO((unsigned EMUSHORT *, long *,
324 enum machine_mode));
325static void eclear PROTO((unsigned EMUSHORT *));
326static void emov PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
327static void eabs PROTO((unsigned EMUSHORT *));
328static void eneg PROTO((unsigned EMUSHORT *));
329static int eisneg PROTO((unsigned EMUSHORT *));
330static int eisinf PROTO((unsigned EMUSHORT *));
331static int eisnan PROTO((unsigned EMUSHORT *));
332static void einfin PROTO((unsigned EMUSHORT *));
333static void enan PROTO((unsigned EMUSHORT *, int));
334static void emovi PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
335static void emovo PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
336static void ecleaz PROTO((unsigned EMUSHORT *));
337static void ecleazs PROTO((unsigned EMUSHORT *));
338static void emovz PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
339static void einan PROTO((unsigned EMUSHORT *));
340static int eiisnan PROTO((unsigned EMUSHORT *));
341static int eiisneg PROTO((unsigned EMUSHORT *));
342static void eiinfin PROTO((unsigned EMUSHORT *));
343static int eiisinf PROTO((unsigned EMUSHORT *));
344static int ecmpm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
345static void eshdn1 PROTO((unsigned EMUSHORT *));
346static void eshup1 PROTO((unsigned EMUSHORT *));
347static void eshdn8 PROTO((unsigned EMUSHORT *));
348static void eshup8 PROTO((unsigned EMUSHORT *));
349static void eshup6 PROTO((unsigned EMUSHORT *));
350static void eshdn6 PROTO((unsigned EMUSHORT *));
351static void eaddm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));\f
352static void esubm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
353static void m16m PROTO((unsigned int, unsigned short *,
354 unsigned short *));
355static int edivm PROTO((unsigned short *, unsigned short *));
356static int emulm PROTO((unsigned short *, unsigned short *));
357static void emdnorm PROTO((unsigned EMUSHORT *, int, int, EMULONG, int));
358static void esub PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
359 unsigned EMUSHORT *));
360static void eadd PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
361 unsigned EMUSHORT *));
362static void eadd1 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
363 unsigned EMUSHORT *));
364static void ediv PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
365 unsigned EMUSHORT *));
366static void emul PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
367 unsigned EMUSHORT *));
368static void e53toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
369static void e64toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
370static void e113toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
371static void e24toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
372static void etoe113 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
373static void toe113 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
374static void etoe64 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
375static void toe64 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
376static void etoe53 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
377static void toe53 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
378static void etoe24 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
379static void toe24 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
380static int ecmp PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
381static void eround PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
382static void ltoe PROTO((HOST_WIDE_INT *, unsigned EMUSHORT *));
383static void ultoe PROTO((unsigned HOST_WIDE_INT *, unsigned EMUSHORT *));
384static void eifrac PROTO((unsigned EMUSHORT *, HOST_WIDE_INT *,
385 unsigned EMUSHORT *));
386static void euifrac PROTO((unsigned EMUSHORT *, unsigned HOST_WIDE_INT *,
387 unsigned EMUSHORT *));
388static int eshift PROTO((unsigned EMUSHORT *, int));
389static int enormlz PROTO((unsigned EMUSHORT *));
390static void e24toasc PROTO((unsigned EMUSHORT *, char *, int));
391static void e53toasc PROTO((unsigned EMUSHORT *, char *, int));
392static void e64toasc PROTO((unsigned EMUSHORT *, char *, int));
393static void e113toasc PROTO((unsigned EMUSHORT *, char *, int));
394static void etoasc PROTO((unsigned EMUSHORT *, char *, int));
395static void asctoe24 PROTO((char *, unsigned EMUSHORT *));
396static void asctoe53 PROTO((char *, unsigned EMUSHORT *));
397static void asctoe64 PROTO((char *, unsigned EMUSHORT *));
398static void asctoe113 PROTO((char *, unsigned EMUSHORT *));
399static void asctoe PROTO((char *, unsigned EMUSHORT *));
400static void asctoeg PROTO((char *, unsigned EMUSHORT *, int));
401static void efloor PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
402static void efrexp PROTO((unsigned EMUSHORT *, int *,
403 unsigned EMUSHORT *));
404static void eldexp PROTO((unsigned EMUSHORT *, int, unsigned EMUSHORT *));
405static void eremain PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
406 unsigned EMUSHORT *));
407static void eiremain PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
408static void mtherr PROTO((char *, int));
409static void dectoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
410static void etodec PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
411static void todec PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
412static void ibmtoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
413 enum machine_mode));
414static void etoibm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
415 enum machine_mode));
416static void toibm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
417 enum machine_mode));
418static void make_nan PROTO((unsigned EMUSHORT *, int, enum machine_mode));
419static void uditoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
420static void ditoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
421static void etoudi PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
422static void etodi PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
423static void esqrt PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
424\f
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. */
428
429static void
430endian (e, x, mode)
431 unsigned EMUSHORT e[];
432 long x[];
433 enum machine_mode mode;
434{
435 unsigned long th, t;
436
437 if (REAL_WORDS_BIG_ENDIAN)
438 {
439 switch (mode)
440 {
441
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;
448
449 case XFmode:
450
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 */
457
458 case DFmode:
459
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 */
466
467 case HFmode:
468 case SFmode:
469
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;
476
477 default:
478 abort ();
479 }
480 }
481 else
482 {
483 /* Pack the output array without swapping. */
484
485 switch (mode)
486 {
487
488 case TFmode:
489
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;
495
496 case XFmode:
497
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 */
506
507 case DFmode:
508
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 */
515
516 case HFmode:
517 case SFmode:
518
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;
525
526 default:
527 abort ();
528 }
529 }
530}
531
532
533/* This is the implementation of the REAL_ARITHMETIC macro. */
534
535void
536earith (value, icode, r1, r2)
537 REAL_VALUE_TYPE *value;
538 int icode;
539 REAL_VALUE_TYPE *r1;
540 REAL_VALUE_TYPE *r2;
541{
542 unsigned EMUSHORT d1[NE], d2[NE], v[NE];
543 enum tree_code code;
544
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))
550 {
551 PUT_REAL (d1, value);
552 return;
553 }
554 if (eisnan (d2))
555 {
556 PUT_REAL (d2, value);
557 return;
558 }
559#endif
560 code = (enum tree_code) icode;
561 switch (code)
562 {
563 case PLUS_EXPR:
564 eadd (d2, d1, v);
565 break;
566
567 case MINUS_EXPR:
568 esub (d2, d1, v); /* d1 - d2 */
569 break;
570
571 case MULT_EXPR:
572 emul (d2, d1, v);
573 break;
574
575 case RDIV_EXPR:
576#ifndef REAL_INFINITY
577 if (ecmp (d2, ezero) == 0)
578 {
579#ifdef NANS
580 enan (v, eisneg (d1) ^ eisneg (d2));
581 break;
582#else
583 abort ();
584#endif
585 }
586#endif
587 ediv (d2, d1, v); /* d1/d2 */
588 break;
589
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;
596
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;
606 }
607PUT_REAL (v, value);
608}
609
610
611/* Truncate REAL_VALUE_TYPE toward zero to signed HOST_WIDE_INT.
612 implements REAL_VALUE_RNDZINT (x) (etrunci (x)). */
613
614REAL_VALUE_TYPE
615etrunci (x)
616 REAL_VALUE_TYPE x;
617{
618 unsigned EMUSHORT f[NE], g[NE];
619 REAL_VALUE_TYPE r;
620 HOST_WIDE_INT l;
621
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);
631}
632
633
634/* Truncate REAL_VALUE_TYPE toward zero to unsigned HOST_WIDE_INT;
635 implements REAL_VALUE_UNSIGNED_RNDZINT (x) (etruncui (x)). */
636
637REAL_VALUE_TYPE
638etruncui (x)
639 REAL_VALUE_TYPE x;
640{
641 unsigned EMUSHORT f[NE], g[NE];
642 REAL_VALUE_TYPE r;
643 unsigned HOST_WIDE_INT l;
644
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);
654}
655
656
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. */
660
661REAL_VALUE_TYPE
662ereal_atof (s, t)
663 char *s;
664 enum machine_mode t;
665{
666 unsigned EMUSHORT tem[NE], e[NE];
667 REAL_VALUE_TYPE r;
668
669 switch (t)
670 {
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);
690 }
691 PUT_REAL (e, &r);
692 return (r);
693}
694
695
696/* Expansion of REAL_NEGATE. */
697
698REAL_VALUE_TYPE
699ereal_negate (x)
700 REAL_VALUE_TYPE x;
701{
702 unsigned EMUSHORT e[NE];
703 REAL_VALUE_TYPE r;
704
705 GET_REAL (&x, e);
706 eneg (e);
707 PUT_REAL (e, &r);
708 return (r);
709}
710
711
712/* Round real toward zero to HOST_WIDE_INT;
713 implements REAL_VALUE_FIX (x). */
714
715HOST_WIDE_INT
716efixi (x)
717 REAL_VALUE_TYPE x;
718{
719 unsigned EMUSHORT f[NE], g[NE];
720 HOST_WIDE_INT l;
721
722 GET_REAL (&x, f);
723#ifdef NANS
724 if (eisnan (f))
725 {
726 warning ("conversion from NaN to int");
727 return (-1);
728 }
729#endif
730 eifrac (f, &l, g);
731 return l;
732}
733
734/* Round real toward zero to unsigned HOST_WIDE_INT
735 implements REAL_VALUE_UNSIGNED_FIX (x).
736 Negative input returns zero. */
737
738unsigned HOST_WIDE_INT
739efixui (x)
740 REAL_VALUE_TYPE x;
741{
742 unsigned EMUSHORT f[NE], g[NE];
743 unsigned HOST_WIDE_INT l;
744
745 GET_REAL (&x, f);
746#ifdef NANS
747 if (eisnan (f))
748 {
749 warning ("conversion from NaN to unsigned int");
750 return (-1);
751 }
752#endif
753 euifrac (f, &l, g);
754 return l;
755}
756
757
758/* REAL_VALUE_FROM_INT macro. */
759
760void
761ereal_from_int (d, i, j, mode)
762 REAL_VALUE_TYPE *d;
763 HOST_WIDE_INT i, j;
764 enum machine_mode mode;
765{
766 unsigned EMUSHORT df[NE], dg[NE];
767 HOST_WIDE_INT low, high;
768 int sign;
769
770 if (GET_MODE_CLASS (mode) != MODE_FLOAT)
771 abort ();
772 sign = 0;
773 low = i;
774 if ((high = j) < 0)
775 {
776 sign = 1;
777 /* complement and add 1 */
778 high = ~high;
779 if (low)
780 low = -low;
781 else
782 high += 1;
783 }
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);
791
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))
796 {
797 case 32:
798 etoe24 (dg, df);
799 e24toe (df, dg);
800 break;
801
802 case 64:
803 etoe53 (dg, df);
804 e53toe (df, dg);
805 break;
806
807 case 96:
808 etoe64 (dg, df);
809 e64toe (df, dg);
810 break;
811
812 case 128:
813 etoe113 (dg, df);
814 e113toe (df, dg);
815 break;
816
817 default:
818 abort ();
819 }
820
821 PUT_REAL (dg, d);
822}
823
824
825/* REAL_VALUE_FROM_UNSIGNED_INT macro. */
826
827void
828ereal_from_uint (d, i, j, mode)
829 REAL_VALUE_TYPE *d;
830 unsigned HOST_WIDE_INT i, j;
831 enum machine_mode mode;
832{
833 unsigned EMUSHORT df[NE], dg[NE];
834 unsigned HOST_WIDE_INT low, high;
835
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);
845
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))
850 {
851 case 32:
852 etoe24 (dg, df);
853 e24toe (df, dg);
854 break;
855
856 case 64:
857 etoe53 (dg, df);
858 e53toe (df, dg);
859 break;
860
861 case 96:
862 etoe64 (dg, df);
863 e64toe (df, dg);
864 break;
865
866 case 128:
867 etoe113 (dg, df);
868 e113toe (df, dg);
869 break;
870
871 default:
872 abort ();
873 }
874
875 PUT_REAL (dg, d);
876}
877
878
879/* REAL_VALUE_TO_INT macro. */
880
881void
882ereal_to_int (low, high, rr)
883 HOST_WIDE_INT *low, *high;
884 REAL_VALUE_TYPE rr;
885{
886 unsigned EMUSHORT d[NE], df[NE], dg[NE], dh[NE];
887 int s;
888
889 GET_REAL (&rr, d);
890#ifdef NANS
891 if (eisnan (d))
892 {
893 warning ("conversion from NaN to int");
894 *low = -1;
895 *high = -1;
896 return;
897 }
898#endif
899 /* convert positive value */
900 s = 0;
901 if (eisneg (d))
902 {
903 eneg (d);
904 s = 1;
905 }
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)
912 {
913 /* complement and add 1 */
914 *high = ~(*high);
915 if (*low)
916 *low = -(*low);
917 else
918 *high += 1;
919 }
920}
921
922
923/* REAL_VALUE_LDEXP macro. */
924
925REAL_VALUE_TYPE
926ereal_ldexp (x, n)
927 REAL_VALUE_TYPE x;
928 int n;
929{
930 unsigned EMUSHORT e[NE], y[NE];
931 REAL_VALUE_TYPE r;
932
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);
941}
942
943/* These routines are conditionally compiled because functions
944 of the same names may be defined in fold-const.c. */
945
946#ifdef REAL_ARITHMETIC
947
948/* Check for infinity in a REAL_VALUE_TYPE. */
949
950int
951target_isinf (x)
952 REAL_VALUE_TYPE x;
953{
954 unsigned EMUSHORT e[NE];
955
956#ifdef INFINITY
957 GET_REAL (&x, e);
958 return (eisinf (e));
959#else
960 return 0;
961#endif
962}
963
964/* Check whether a REAL_VALUE_TYPE item is a NaN. */
965
966int
967target_isnan (x)
968 REAL_VALUE_TYPE x;
969{
970 unsigned EMUSHORT e[NE];
971
972#ifdef NANS
973 GET_REAL (&x, e);
974 return (eisnan (e));
975#else
976 return (0);
977#endif
978}
979
980
981/* Check for a negative REAL_VALUE_TYPE number.
982 This just checks the sign bit, so that -0 counts as negative. */
983
984int
985target_negative (x)
986 REAL_VALUE_TYPE x;
987{
988 return ereal_isneg (x);
989}
990
991/* Expansion of REAL_VALUE_TRUNCATE.
992 The result is in floating point, rounded to nearest or even. */
993
994REAL_VALUE_TYPE
995real_value_truncate (mode, arg)
996 enum machine_mode mode;
997 REAL_VALUE_TYPE arg;
998{
999 unsigned EMUSHORT e[NE], t[NE];
1000 REAL_VALUE_TYPE r;
1001
1002 GET_REAL (&arg, e);
1003#ifdef NANS
1004 if (eisnan (e))
1005 return (arg);
1006#endif
1007 eclear (t);
1008 switch (mode)
1009 {
1010 case TFmode:
1011 etoe113 (e, t);
1012 e113toe (t, t);
1013 break;
1014
1015 case XFmode:
1016 etoe64 (e, t);
1017 e64toe (t, t);
1018 break;
1019
1020 case DFmode:
1021 etoe53 (e, t);
1022 e53toe (t, t);
1023 break;
1024
1025 case HFmode:
1026 case SFmode:
1027 etoe24 (e, t);
1028 e24toe (t, t);
1029 break;
1030
1031 case SImode:
1032 r = etrunci (arg);
1033 return (r);
1034
1035 /* If an unsupported type was requested, presume that
1036 the machine files know something useful to do with
1037 the unmodified value. */
1038
1039 default:
1040 return (arg);
1041 }
1042 PUT_REAL (t, &r);
1043 return (r);
1044}
1045
1046/* Try to change R into its exact multiplicative inverse in machine mode
1047 MODE. Return nonzero function value if successful. */
1048
1049int
1050exact_real_inverse (mode, r)
1051 enum machine_mode mode;
1052 REAL_VALUE_TYPE *r;
1053{
1054 unsigned EMUSHORT e[NE], einv[NE];
1055 REAL_VALUE_TYPE rinv;
1056 int i;
1057
1058 GET_REAL (r, e);
1059
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;
1063
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;
1068
1069 for (i = 0; i < NE - 2; i++)
1070 {
1071 if (e[i] != 0)
1072 return 0;
1073 }
1074
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);
1079
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);
1088
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;
1093
1094 for (i = 0; i < NE - 2; i++)
1095 {
1096 if (einv[i] != 0)
1097 return 0;
1098 }
1099
1100 /* Fail if the computed inverse is out of range. */
1101 if (eisinf (einv) || eisnan (einv) || (ecmp (einv, ezero) == 0))
1102 return 0;
1103
1104 /* Output the reciprocal and return success flag. */
1105 PUT_REAL (einv, r);
1106 return 1;
1107}
1108#endif /* REAL_ARITHMETIC defined */
1109
1110/* Used for debugging--print the value of R in human-readable format
1111 on stderr. */
1112
1113void
1114debug_real (r)
1115 REAL_VALUE_TYPE r;
1116{
1117 char dstr[30];
1118
1119 REAL_VALUE_TO_DECIMAL (r, "%.20g", dstr);
1120 fprintf (stderr, "%s", dstr);
1121}
1122
1123\f
1124/* The following routines convert REAL_VALUE_TYPE to the various floating
1125 point formats that are meaningful to supported computers.
1126
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
1129
1130 fprintf (file, "%lx, %lx", L[0], L[1]);
1131
1132 that will work on both narrow- and wide-word host computers. */
1133
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. */
1137
1138void
1139etartdouble (r, l)
1140 REAL_VALUE_TYPE r;
1141 long l[];
1142{
1143 unsigned EMUSHORT e[NE];
1144
1145 GET_REAL (&r, e);
1146 etoe113 (e, e);
1147 endian (e, l, TFmode);
1148}
1149
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. */
1153
1154void
1155etarldouble (r, l)
1156 REAL_VALUE_TYPE r;
1157 long l[];
1158{
1159 unsigned EMUSHORT e[NE];
1160
1161 GET_REAL (&r, e);
1162 etoe64 (e, e);
1163 endian (e, l, XFmode);
1164}
1165
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. */
1168
1169void
1170etardouble (r, l)
1171 REAL_VALUE_TYPE r;
1172 long l[];
1173{
1174 unsigned EMUSHORT e[NE];
1175
1176 GET_REAL (&r, e);
1177 etoe53 (e, e);
1178 endian (e, l, DFmode);
1179}
1180
1181/* Convert R to a single precision float value stored in the least-significant
1182 bits of a `long'. */
1183
1184long
1185etarsingle (r)
1186 REAL_VALUE_TYPE r;
1187{
1188 unsigned EMUSHORT e[NE];
1189 long l;
1190
1191 GET_REAL (&r, e);
1192 etoe24 (e, e);
1193 endian (e, &l, SFmode);
1194 return ((long) l);
1195}
1196
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. */
1201
1202void
1203ereal_to_decimal (x, s)
1204 REAL_VALUE_TYPE x;
1205 char *s;
1206{
1207 unsigned EMUSHORT e[NE];
1208
1209 GET_REAL (&x, e);
1210 etoasc (e, s, 20);
1211}
1212
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. */
1215
1216int
1217ereal_cmp (x, y)
1218 REAL_VALUE_TYPE x, y;
1219{
1220 unsigned EMUSHORT ex[NE], ey[NE];
1221
1222 GET_REAL (&x, ex);
1223 GET_REAL (&y, ey);
1224 return (ecmp (ex, ey));
1225}
1226
1227/* Return 1 if the sign bit of X is set, else return 0. */
1228
1229int
1230ereal_isneg (x)
1231 REAL_VALUE_TYPE x;
1232{
1233 unsigned EMUSHORT ex[NE];
1234
1235 GET_REAL (&x, ex);
1236 return (eisneg (ex));
1237}
1238
1239/* End of REAL_ARITHMETIC interface */
1240\f
1241/*
1242 Extended precision IEEE binary floating point arithmetic routines
1243
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.
1247
1248 External e type data structure, similar to Intel 8087 chip
1249 temporary real format but possibly with a larger significand:
1250
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)
1255
1256
1257 Internal exploded e-type data structure of a number (a "word" is 16 bits):
1258
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)
1267
1268
1269
1270 Routines for external format e-type numbers
1271
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
1312
1313
1314 Routines for internal format exploded e-type numbers
1315
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
1340
1341 The result is always normalized and rounded to NI-4 word precision
1342 after each arithmetic operation.
1343
1344 Exception flags are NOT fully supported.
1345
1346 Signaling NaN's are NOT supported; they are treated the same
1347 as quiet NaN's.
1348
1349 Define INFINITY for support of infinity; otherwise a
1350 saturation arithmetic is implemented.
1351
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.
1359
1360 Denormals are always supported here where appropriate (e.g., not
1361 for conversion to DEC numbers). */
1362
1363/* Definitions for error codes that are passed to the common error handling
1364 routine mtherr.
1365
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.
1372
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.
1380
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.
1386
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.
1393
1394 For ANSI C compatibility, define ANSIC equal to 1. Currently
1395 this affects only the atan2 function and others that use it. */
1396
1397/* Constant definitions for math error conditions. */
1398
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 */
1406
1407/* e type constants used by high precision check routines */
1408
1409#if LONG_DOUBLE_TYPE_SIZE == 128
1410/* 0.0 */
1411unsigned EMUSHORT ezero[NE] =
1412 {0x0000, 0x0000, 0x0000, 0x0000,
1413 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,};
1414extern unsigned EMUSHORT ezero[];
1415
1416/* 5.0E-1 */
1417unsigned EMUSHORT ehalf[NE] =
1418 {0x0000, 0x0000, 0x0000, 0x0000,
1419 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3ffe,};
1420extern unsigned EMUSHORT ehalf[];
1421
1422/* 1.0E0 */
1423unsigned EMUSHORT eone[NE] =
1424 {0x0000, 0x0000, 0x0000, 0x0000,
1425 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,};
1426extern unsigned EMUSHORT eone[];
1427
1428/* 2.0E0 */
1429unsigned EMUSHORT etwo[NE] =
1430 {0x0000, 0x0000, 0x0000, 0x0000,
1431 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4000,};
1432extern unsigned EMUSHORT etwo[];
1433
1434/* 3.2E1 */
1435unsigned EMUSHORT e32[NE] =
1436 {0x0000, 0x0000, 0x0000, 0x0000,
1437 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4004,};
1438extern unsigned EMUSHORT e32[];
1439
1440/* 6.93147180559945309417232121458176568075500134360255E-1 */
1441unsigned EMUSHORT elog2[NE] =
1442 {0x40f3, 0xf6af, 0x03f2, 0xb398,
1443 0xc9e3, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1444extern unsigned EMUSHORT elog2[];
1445
1446/* 1.41421356237309504880168872420969807856967187537695E0 */
1447unsigned EMUSHORT esqrt2[NE] =
1448 {0x1d6f, 0xbe9f, 0x754a, 0x89b3,
1449 0x597d, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1450extern unsigned EMUSHORT esqrt2[];
1451
1452/* 3.14159265358979323846264338327950288419716939937511E0 */
1453unsigned EMUSHORT epi[NE] =
1454 {0x2902, 0x1cd1, 0x80dc, 0x628b,
1455 0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1456extern unsigned EMUSHORT epi[];
1457
1458#else
1459/* LONG_DOUBLE_TYPE_SIZE is other than 128 */
1460unsigned EMUSHORT ezero[NE] =
1461 {0, 0000000, 0000000, 0000000, 0000000, 0000000,};
1462unsigned EMUSHORT ehalf[NE] =
1463 {0, 0000000, 0000000, 0000000, 0100000, 0x3ffe,};
1464unsigned EMUSHORT eone[NE] =
1465 {0, 0000000, 0000000, 0000000, 0100000, 0x3fff,};
1466unsigned EMUSHORT etwo[NE] =
1467 {0, 0000000, 0000000, 0000000, 0100000, 0040000,};
1468unsigned EMUSHORT e32[NE] =
1469 {0, 0000000, 0000000, 0000000, 0100000, 0040004,};
1470unsigned EMUSHORT elog2[NE] =
1471 {0xc9e4, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1472unsigned EMUSHORT esqrt2[NE] =
1473 {0x597e, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1474unsigned EMUSHORT epi[NE] =
1475 {0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1476#endif
1477
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. */
1480
1481int rndprc = NBITS;
1482extern int rndprc;
1483
1484/* Clear out entire e-type number X. */
1485
1486static void
1487eclear (x)
1488 register unsigned EMUSHORT *x;
1489{
1490 register int i;
1491
1492 for (i = 0; i < NE; i++)
1493 *x++ = 0;
1494}
1495
1496/* Move e-type number from A to B. */
1497
1498static void
1499emov (a, b)
1500 register unsigned EMUSHORT *a, *b;
1501{
1502 register int i;
1503
1504 for (i = 0; i < NE; i++)
1505 *b++ = *a++;
1506}
1507
1508
1509/* Absolute value of e-type X. */
1510
1511static void
1512eabs (x)
1513 unsigned EMUSHORT x[];
1514{
1515 /* sign is top bit of last word of external format */
1516 x[NE - 1] &= 0x7fff;
1517}
1518
1519/* Negate the e-type number X. */
1520
1521static void
1522eneg (x)
1523 unsigned EMUSHORT x[];
1524{
1525
1526 x[NE - 1] ^= 0x8000; /* Toggle the sign bit */
1527}
1528
1529/* Return 1 if sign bit of e-type number X is nonzero, else zero. */
1530
1531static int
1532eisneg (x)
1533 unsigned EMUSHORT x[];
1534{
1535
1536 if (x[NE - 1] & 0x8000)
1537 return (1);
1538 else
1539 return (0);
1540}
1541
1542/* Return 1 if e-type number X is infinity, else return zero. */
1543
1544static int
1545eisinf (x)
1546 unsigned EMUSHORT x[];
1547{
1548
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);
1557}
1558
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. */
1561
1562static int
1563eisnan (x)
1564 unsigned EMUSHORT x[];
1565{
1566#ifdef NANS
1567 int i;
1568
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++)
1574 {
1575 if (*x++ != 0)
1576 return (1);
1577 }
1578#endif
1579
1580 return (0);
1581}
1582
1583/* Fill e-type number X with infinity pattern (IEEE)
1584 or largest possible number (non-IEEE). */
1585
1586static void
1587einfin (x)
1588 register unsigned EMUSHORT *x;
1589{
1590 register int i;
1591
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)
1601 {
1602 if (rndprc == 113)
1603 {
1604 *(x - 9) = 0;
1605 *(x - 8) = 0;
1606 }
1607 if (rndprc == 64)
1608 {
1609 *(x - 5) = 0;
1610 }
1611 if (rndprc == 53)
1612 {
1613 *(x - 4) = 0xf800;
1614 }
1615 else
1616 {
1617 *(x - 4) = 0;
1618 *(x - 3) = 0;
1619 *(x - 2) = 0xff00;
1620 }
1621 }
1622#endif
1623}
1624
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. */
1628
1629static void
1630enan (x, sign)
1631 register unsigned EMUSHORT *x;
1632 int sign;
1633{
1634 register int i;
1635
1636 for (i = 0; i < NE - 2; i++)
1637 *x++ = 0;
1638 *x++ = 0xc000;
1639 *x = (sign << 15) | 0x7fff;
1640}
1641
1642/* Move in an e-type number A, converting it to exploded e-type B. */
1643
1644static void
1645emovi (a, b)
1646 unsigned EMUSHORT *a, *b;
1647{
1648 register unsigned EMUSHORT *p, *q;
1649 int i;
1650
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)
1663 {
1664#ifdef NANS
1665 if (eisnan (a))
1666 {
1667 *q++ = 0;
1668 for (i = 3; i < NI; i++)
1669 *q++ = *p--;
1670 return;
1671 }
1672#endif
1673
1674 for (i = 2; i < NI; i++)
1675 *q++ = 0;
1676 return;
1677 }
1678#endif
1679
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;
1687}
1688
1689/* Move out exploded e-type number A, converting it to e type B. */
1690
1691static void
1692emovo (a, b)
1693 unsigned EMUSHORT *a, *b;
1694{
1695 register unsigned EMUSHORT *p, *q;
1696 unsigned EMUSHORT i;
1697 int j;
1698
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)
1709 {
1710#ifdef NANS
1711 if (eiisnan (a))
1712 {
1713 enan (b, eiisneg (a));
1714 return;
1715 }
1716#endif
1717 einfin (b);
1718 return;
1719 }
1720#endif
1721 /* skip over guard word */
1722 ++p;
1723 /* move the significand */
1724 for (j = 0; j < NE - 1; j++)
1725 *q-- = *p++;
1726}
1727
1728/* Clear out exploded e-type number XI. */
1729
1730static void
1731ecleaz (xi)
1732 register unsigned EMUSHORT *xi;
1733{
1734 register int i;
1735
1736 for (i = 0; i < NI; i++)
1737 *xi++ = 0;
1738}
1739
1740/* Clear out exploded e-type XI, but don't touch the sign. */
1741
1742static void
1743ecleazs (xi)
1744 register unsigned EMUSHORT *xi;
1745{
1746 register int i;
1747
1748 ++xi;
1749 for (i = 0; i < NI - 1; i++)
1750 *xi++ = 0;
1751}
1752
1753/* Move exploded e-type number from A to B. */
1754
1755static void
1756emovz (a, b)
1757 register unsigned EMUSHORT *a, *b;
1758{
1759 register int i;
1760
1761 for (i = 0; i < NI - 1; i++)
1762 *b++ = *a++;
1763 /* clear low guard word */
1764 *b = 0;
1765}
1766
1767/* Generate exploded e-type NaN.
1768 The explicit pattern for this is maximum exponent and
1769 top two significant bits set. */
1770
1771static void
1772einan (x)
1773 unsigned EMUSHORT x[];
1774{
1775
1776 ecleaz (x);
1777 x[E] = 0x7fff;
1778 x[M + 1] = 0xc000;
1779}
1780
1781/* Return nonzero if exploded e-type X is a NaN. */
1782
1783static int
1784eiisnan (x)
1785 unsigned EMUSHORT x[];
1786{
1787 int i;
1788
1789 if ((x[E] & 0x7fff) == 0x7fff)
1790 {
1791 for (i = M + 1; i < NI; i++)
1792 {
1793 if (x[i] != 0)
1794 return (1);
1795 }
1796 }
1797 return (0);
1798}
1799
1800/* Return nonzero if sign of exploded e-type X is nonzero. */
1801
1802static int
1803eiisneg (x)
1804 unsigned EMUSHORT x[];
1805{
1806
1807 return x[0] != 0;
1808}
1809
1810/* Fill exploded e-type X with infinity pattern.
1811 This has maximum exponent and significand all zeros. */
1812
1813static void
1814eiinfin (x)
1815 unsigned EMUSHORT x[];
1816{
1817
1818 ecleaz (x);
1819 x[E] = 0x7fff;
1820}
1821
1822/* Return nonzero if exploded e-type X is infinite. */
1823
1824static int
1825eiisinf (x)
1826 unsigned EMUSHORT x[];
1827{
1828
1829#ifdef NANS
1830 if (eiisnan (x))
1831 return (0);
1832#endif
1833 if ((x[E] & 0x7fff) == 0x7fff)
1834 return (1);
1835 return (0);
1836}
1837
1838
1839/* Compare significands of numbers in internal exploded e-type format.
1840 Guard words are included in the comparison.
1841
1842 Returns +1 if a > b
1843 0 if a == b
1844 -1 if a < b */
1845
1846static int
1847ecmpm (a, b)
1848 register unsigned EMUSHORT *a, *b;
1849{
1850 int i;
1851
1852 a += M; /* skip up to significand area */
1853 b += M;
1854 for (i = M; i < NI; i++)
1855 {
1856 if (*a++ != *b++)
1857 goto difrnt;
1858 }
1859 return (0);
1860
1861 difrnt:
1862 if (*(--a) > *(--b))
1863 return (1);
1864 else
1865 return (-1);
1866}
1867
1868/* Shift significand of exploded e-type X down by 1 bit. */
1869
1870static void
1871eshdn1 (x)
1872 register unsigned EMUSHORT *x;
1873{
1874 register unsigned EMUSHORT bits;
1875 int i;
1876
1877 x += M; /* point to significand area */
1878
1879 bits = 0;
1880 for (i = M; i < NI; i++)
1881 {
1882 if (*x & 1)
1883 bits |= 1;
1884 *x >>= 1;
1885 if (bits & 2)
1886 *x |= 0x8000;
1887 bits <<= 1;
1888 ++x;
1889 }
1890}
1891
1892/* Shift significand of exploded e-type X up by 1 bit. */
1893
1894static void
1895eshup1 (x)
1896 register unsigned EMUSHORT *x;
1897{
1898 register unsigned EMUSHORT bits;
1899 int i;
1900
1901 x += NI - 1;
1902 bits = 0;
1903
1904 for (i = M; i < NI; i++)
1905 {
1906 if (*x & 0x8000)
1907 bits |= 1;
1908 *x <<= 1;
1909 if (bits & 2)
1910 *x |= 1;
1911 bits <<= 1;
1912 --x;
1913 }
1914}
1915
1916
1917/* Shift significand of exploded e-type X down by 8 bits. */
1918
1919static void
1920eshdn8 (x)
1921 register unsigned EMUSHORT *x;
1922{
1923 register unsigned EMUSHORT newbyt, oldbyt;
1924 int i;
1925
1926 x += M;
1927 oldbyt = 0;
1928 for (i = M; i < NI; i++)
1929 {
1930 newbyt = *x << 8;
1931 *x >>= 8;
1932 *x |= oldbyt;
1933 oldbyt = newbyt;
1934 ++x;
1935 }
1936}
1937
1938/* Shift significand of exploded e-type X up by 8 bits. */
1939
1940static void
1941eshup8 (x)
1942 register unsigned EMUSHORT *x;
1943{
1944 int i;
1945 register unsigned EMUSHORT newbyt, oldbyt;
1946
1947 x += NI - 1;
1948 oldbyt = 0;
1949
1950 for (i = M; i < NI; i++)
1951 {
1952 newbyt = *x >> 8;
1953 *x <<= 8;
1954 *x |= oldbyt;
1955 oldbyt = newbyt;
1956 --x;
1957 }
1958}
1959
1960/* Shift significand of exploded e-type X up by 16 bits. */
1961
1962static void
1963eshup6 (x)
1964 register unsigned EMUSHORT *x;
1965{
1966 int i;
1967 register unsigned EMUSHORT *p;
1968
1969 p = x + M;
1970 x += M + 1;
1971
1972 for (i = M; i < NI - 1; i++)
1973 *p++ = *x++;
1974
1975 *p = 0;
1976}
1977
1978/* Shift significand of exploded e-type X down by 16 bits. */
1979
1980static void
1981eshdn6 (x)
1982 register unsigned EMUSHORT *x;
1983{
1984 int i;
1985 register unsigned EMUSHORT *p;
1986
1987 x += NI - 1;
1988 p = x + 1;
1989
1990 for (i = M; i < NI - 1; i++)
1991 *(--p) = *(--x);
1992
1993 *(--p) = 0;
1994}
1995
1996/* Add significands of exploded e-type X and Y. X + Y replaces Y. */
1997
1998static void
1999eaddm (x, y)
2000 unsigned EMUSHORT *x, *y;
2001{
2002 register unsigned EMULONG a;
2003 int i;
2004 unsigned int carry;
2005
2006 x += NI - 1;
2007 y += NI - 1;
2008 carry = 0;
2009 for (i = M; i < NI; i++)
2010 {
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;
2019 }
2020}
2021
2022/* Subtract significands of exploded e-type X and Y. Y - X replaces Y. */
2023
2024static void
2025esubm (x, y)
2026 unsigned EMUSHORT *x, *y;
2027{
2028 unsigned EMULONG a;
2029 int i;
2030 unsigned int carry;
2031
2032 x += NI - 1;
2033 y += NI - 1;
2034 carry = 0;
2035 for (i = M; i < NI; i++)
2036 {
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;
2045 }
2046}
2047
2048
2049static unsigned EMUSHORT equot[NI];
2050
2051
2052#if 0
2053/* Radix 2 shift-and-add versions of multiply and divide */
2054
2055
2056/* Divide significands */
2057
2058int
2059edivm (den, num)
2060 unsigned EMUSHORT den[], num[];
2061{
2062 int i;
2063 register unsigned EMUSHORT *p, *q;
2064 unsigned EMUSHORT j;
2065
2066 p = &equot[0];
2067 *p++ = num[0];
2068 *p++ = num[1];
2069
2070 for (i = M; i < NI; i++)
2071 {
2072 *p++ = 0;
2073 }
2074
2075 /* Use faster compare and subtraction if denominator has only 15 bits of
2076 significance. */
2077
2078 p = &den[M + 2];
2079 if (*p++ == 0)
2080 {
2081 for (i = M + 3; i < NI; i++)
2082 {
2083 if (*p++ != 0)
2084 goto fulldiv;
2085 }
2086 if ((den[M + 1] & 1) != 0)
2087 goto fulldiv;
2088 eshdn1 (num);
2089 eshdn1 (den);
2090
2091 p = &den[M + 1];
2092 q = &num[M + 1];
2093
2094 for (i = 0; i < NBITS + 2; i++)
2095 {
2096 if (*p <= *q)
2097 {
2098 *q -= *p;
2099 j = 1;
2100 }
2101 else
2102 {
2103 j = 0;
2104 }
2105 eshup1 (equot);
2106 equot[NI - 2] |= j;
2107 eshup1 (num);
2108 }
2109 goto divdon;
2110 }
2111
2112 /* The number of quotient bits to calculate is NBITS + 1 scaling guard
2113 bit + 1 roundoff bit. */
2114
2115 fulldiv:
2116
2117 p = &equot[NI - 2];
2118 for (i = 0; i < NBITS + 2; i++)
2119 {
2120 if (ecmpm (den, num) <= 0)
2121 {
2122 esubm (den, num);
2123 j = 1; /* quotient bit = 1 */
2124 }
2125 else
2126 j = 0;
2127 eshup1 (equot);
2128 *p |= j;
2129 eshup1 (num);
2130 }
2131
2132 divdon:
2133
2134 eshdn1 (equot);
2135 eshdn1 (equot);
2136
2137 /* test for nonzero remainder after roundoff bit */
2138 p = &num[M];
2139 j = 0;
2140 for (i = M; i < NI; i++)
2141 {
2142 j |= *p++;
2143 }
2144 if (j)
2145 j = 1;
2146
2147
2148 for (i = 0; i < NI; i++)
2149 num[i] = equot[i];
2150 return ((int) j);
2151}
2152
2153
2154/* Multiply significands */
2155
2156int
2157emulm (a, b)
2158 unsigned EMUSHORT a[], b[];
2159{
2160 unsigned EMUSHORT *p, *q;
2161 int i, j, k;
2162
2163 equot[0] = b[0];
2164 equot[1] = b[1];
2165 for (i = M; i < NI; i++)
2166 equot[i] = 0;
2167
2168 p = &a[NI - 2];
2169 k = NBITS;
2170 while (*p == 0) /* significand is not supposed to be zero */
2171 {
2172 eshdn6 (a);
2173 k -= 16;
2174 }
2175 if ((*p & 0xff) == 0)
2176 {
2177 eshdn8 (a);
2178 k -= 8;
2179 }
2180
2181 q = &equot[NI - 1];
2182 j = 0;
2183 for (i = 0; i < k; i++)
2184 {
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);
2192 }
2193
2194 for (i = 0; i < NI; i++)
2195 b[i] = equot[i];
2196
2197 /* return flag for lost nonzero bits */
2198 return (j);
2199}
2200
2201#else
2202
2203/* Radix 65536 versions of multiply and divide. */
2204
2205/* Multiply significand of e-type number B
2206 by 16-bit quantity A, return e-type result to C. */
2207
2208static void
2209m16m (a, b, c)
2210 unsigned int a;
2211 unsigned EMUSHORT b[], c[];
2212{
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;
2219
2220 aa = a;
2221 pp = &p[NI-2];
2222 *pp++ = 0;
2223 *pp = 0;
2224 ps = &b[NI-1];
2225
2226 for (i=M+1; i<NI; i++)
2227 {
2228 if (*ps == 0)
2229 {
2230 --ps;
2231 --pp;
2232 *(pp-1) = 0;
2233 }
2234 else
2235 {
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;
2242 }
2243 }
2244 for (i=M; i<NI; i++)
2245 c[i] = p[i];
2246}
2247
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. */
2251
2252static int
2253edivm (den, num)
2254 unsigned EMUSHORT den[], num[];
2255{
2256 int i;
2257 register unsigned EMUSHORT *p;
2258 unsigned EMULONG tnum;
2259 unsigned EMUSHORT j, tdenm, tquot;
2260 unsigned EMUSHORT tprod[NI+1];
2261
2262 p = &equot[0];
2263 *p++ = num[0];
2264 *p++ = num[1];
2265
2266 for (i=M; i<NI; i++)
2267 {
2268 *p++ = 0;
2269 }
2270 eshdn1 (num);
2271 tdenm = den[M+1];
2272 for (i=M; i<NI; i++)
2273 {
2274 /* Find trial quotient digit (the radix is 65536). */
2275 tnum = (((unsigned EMULONG) num[M]) << 16) + num[M+1];
2276
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)
2286 {
2287 tquot -= 1;
2288 esubm (den, tprod);
2289 if (ecmpm (tprod, num) > 0)
2290 {
2291 tquot -= 1;
2292 esubm (den, tprod);
2293 }
2294 }
2295 esubm (tprod, num);
2296 equot[i] = tquot;
2297 eshup6(num);
2298 }
2299 /* test for nonzero remainder after roundoff bit */
2300 p = &num[M];
2301 j = 0;
2302 for (i=M; i<NI; i++)
2303 {
2304 j |= *p++;
2305 }
2306 if (j)
2307 j = 1;
2308
2309 for (i=0; i<NI; i++)
2310 num[i] = equot[i];
2311
2312 return ((int)j);
2313}
2314
2315/* Multiply significands of exploded e-type A and B, result in B. */
2316
2317static int
2318emulm (a, b)
2319 unsigned EMUSHORT a[], b[];
2320{
2321 unsigned EMUSHORT *p, *q;
2322 unsigned EMUSHORT pprod[NI];
2323 unsigned EMUSHORT j;
2324 int i;
2325
2326 equot[0] = b[0];
2327 equot[1] = b[1];
2328 for (i=M; i<NI; i++)
2329 equot[i] = 0;
2330
2331 j = 0;
2332 p = &a[NI-1];
2333 q = &equot[NI-1];
2334 for (i=M+1; i<NI; i++)
2335 {
2336 if (*p == 0)
2337 {
2338 --p;
2339 }
2340 else
2341 {
2342 m16m ((unsigned int) *p--, b, pprod);
2343 eaddm(pprod, equot);
2344 }
2345 j |= *q;
2346 eshdn6(equot);
2347 }
2348
2349 for (i=0; i<NI; i++)
2350 b[i] = equot[i];
2351
2352 /* return flag for lost nonzero bits */
2353 return ((int)j);
2354}
2355#endif
2356
2357
2358/* Normalize and round off.
2359
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.
2362
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.
2366
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.
2370
2371 Input RCNTRL is the rounding control. If it is nonzero, the
2372 returned value will be rounded to RNDPRC bits.
2373
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. */
2384
2385static int rlast = -1;
2386static int rw = 0;
2387static unsigned EMUSHORT rmsk = 0;
2388static unsigned EMUSHORT rmbit = 0;
2389static unsigned EMUSHORT rebit = 0;
2390static int re = 0;
2391static unsigned EMUSHORT rbit[NI];
2392
2393static void
2394emdnorm (s, lost, subflg, exp, rcntrl)
2395 unsigned EMUSHORT s[];
2396 int lost;
2397 int subflg;
2398 EMULONG exp;
2399 int rcntrl;
2400{
2401 int i, j;
2402 unsigned EMUSHORT r;
2403
2404 /* Normalize */
2405 j = enormlz (s);
2406
2407 /* a blank significand could mean either zero or infinity. */
2408#ifndef INFINITY
2409 if (j > NBITS)
2410 {
2411 ecleazs (s);
2412 return;
2413 }
2414#endif
2415 exp -= j;
2416#ifndef INFINITY
2417 if (exp >= 32767L)
2418 goto overf;
2419#else
2420 if ((j > NBITS) && (exp < 32767))
2421 {
2422 ecleazs (s);
2423 return;
2424 }
2425#endif
2426 if (exp < 0L)
2427 {
2428 if (exp > (EMULONG) (-NBITS - 1))
2429 {
2430 j = (int) exp;
2431 i = eshift (s, j);
2432 if (i)
2433 lost = 1;
2434 }
2435 else
2436 {
2437 ecleazs (s);
2438 return;
2439 }
2440 }
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)
2446 {
2447 ecleaz (rbit);
2448 switch (rndprc)
2449 {
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;
2494 }
2495 rbit[re] = rebit;
2496 rlast = rndprc;
2497 }
2498
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)))
2504 {
2505 lost |= s[NI - 1] & 1;
2506 eshdn1 (s);
2507 }
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)
2512 {
2513 i = rw + 1;
2514 while (i < NI)
2515 {
2516 if (s[i])
2517 r |= 1;
2518 s[i] = 0;
2519 ++i;
2520 }
2521 }
2522 s[rw] &= ~rmsk;
2523 if ((r & rmbit) != 0)
2524 {
2525 if (r == rmbit)
2526 {
2527 if (lost == 0)
2528 { /* round to even */
2529 if ((s[re] & rebit) == 0)
2530 goto mddone;
2531 }
2532 else
2533 {
2534 if (subflg != 0)
2535 goto mddone;
2536 }
2537 }
2538 eaddm (rbit, s);
2539 }
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)))
2544 {
2545 eshup1 (s);
2546 }
2547 if (s[2] != 0)
2548 { /* overflow on roundoff */
2549 eshdn1 (s);
2550 exp += 1;
2551 }
2552 mdfin:
2553 s[NI - 1] = 0;
2554 if (exp >= 32767L)
2555 {
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))
2572 {
2573 s[rw] &= ~rmsk;
2574 if (rndprc == 24)
2575 {
2576 s[5] = 0;
2577 s[6] = 0;
2578 }
2579 }
2580#endif
2581 return;
2582 }
2583 if (exp < 0)
2584 s[1] = 0;
2585 else
2586 s[1] = (unsigned EMUSHORT) exp;
2587}
2588
2589/* Subtract. C = B - A, all e type numbers. */
2590
2591static int subflg = 0;
2592
2593static void
2594esub (a, b, c)
2595 unsigned EMUSHORT *a, *b, *c;
2596{
2597
2598#ifdef NANS
2599 if (eisnan (a))
2600 {
2601 emov (a, c);
2602 return;
2603 }
2604 if (eisnan (b))
2605 {
2606 emov (b, c);
2607 return;
2608 }
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))
2613 {
2614 mtherr ("esub", INVALID);
2615 enan (c, 0);
2616 return;
2617 }
2618#endif
2619 subflg = 1;
2620 eadd1 (a, b, c);
2621}
2622
2623/* Add. C = A + B, all e type. */
2624
2625static void
2626eadd (a, b, c)
2627 unsigned EMUSHORT *a, *b, *c;
2628{
2629
2630#ifdef NANS
2631/* NaN plus anything is a NaN. */
2632 if (eisnan (a))
2633 {
2634 emov (a, c);
2635 return;
2636 }
2637 if (eisnan (b))
2638 {
2639 emov (b, c);
2640 return;
2641 }
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))
2646 {
2647 mtherr ("esub", INVALID);
2648 enan (c, 0);
2649 return;
2650 }
2651#endif
2652 subflg = 0;
2653 eadd1 (a, b, c);
2654}
2655
2656/* Arithmetic common to both addition and subtraction. */
2657
2658static void
2659eadd1 (a, b, c)
2660 unsigned EMUSHORT *a, *b, *c;
2661{
2662 unsigned EMUSHORT ai[NI], bi[NI], ci[NI];
2663 int i, lost, j, k;
2664 EMULONG lt, lta, ltb;
2665
2666#ifdef INFINITY
2667 if (eisinf (a))
2668 {
2669 emov (a, c);
2670 if (subflg)
2671 eneg (c);
2672 return;
2673 }
2674 if (eisinf (b))
2675 {
2676 emov (b, c);
2677 return;
2678 }
2679#endif
2680 emovi (a, ai);
2681 emovi (b, bi);
2682 if (subflg)
2683 ai[0] = ~ai[0];
2684
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;
2696 }
2697 lost = 0;
2698 if (lt != 0L)
2699 {
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 */
2704 }
2705 else
2706 {
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])
2713 {
2714 eclear (c);
2715 return;
2716 }
2717 /* if same sign, result is double */
2718 /* double denormalized tiny number */
2719 if ((bi[E] == 0) && ((bi[3] & 0x8000) == 0))
2720 {
2721 eshup1 (bi);
2722 goto done;
2723 }
2724 /* add 1 to exponent unless both are zero! */
2725 for (j = 1; j < NI - 1; j++)
2726 {
2727 if (bi[j] != 0)
2728 {
2729 ltb += 1;
2730 if (ltb >= 0x7fff)
2731 {
2732 eclear (c);
2733 if (ai[0] != 0)
2734 eneg (c);
2735 einfin (c);
2736 return;
2737 }
2738 break;
2739 }
2740 }
2741 bi[E] = (unsigned EMUSHORT) ltb;
2742 goto done;
2743 }
2744 if (i > 0)
2745 { /* put the larger number in bi */
2746 emovz (bi, ci);
2747 emovz (ai, bi);
2748 emovz (ci, ai);
2749 }
2750 }
2751 if (ai[0] == bi[0])
2752 {
2753 eaddm (ai, bi);
2754 subflg = 0;
2755 }
2756 else
2757 {
2758 esubm (ai, bi);
2759 subflg = 1;
2760 }
2761 emdnorm (bi, lost, subflg, ltb, 64);
2762
2763 done:
2764 emovo (bi, c);
2765}
2766
2767/* Divide: C = B/A, all e type. */
2768
2769static void
2770ediv (a, b, c)
2771 unsigned EMUSHORT *a, *b, *c;
2772{
2773 unsigned EMUSHORT ai[NI], bi[NI];
2774 int i, sign;
2775 EMULONG lt, lta, ltb;
2776
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);
2780
2781#ifdef NANS
2782/* Return any NaN input. */
2783 if (eisnan (a))
2784 {
2785 emov (a, c);
2786 return;
2787 }
2788 if (eisnan (b))
2789 {
2790 emov (b, c);
2791 return;
2792 }
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)))
2796 {
2797 mtherr ("ediv", INVALID);
2798 enan (c, sign);
2799 return;
2800 }
2801#endif
2802/* Infinity over anything else is infinity. */
2803#ifdef INFINITY
2804 if (eisinf (b))
2805 {
2806 einfin (c);
2807 goto divsign;
2808 }
2809/* Anything else over infinity is zero. */
2810 if (eisinf (a))
2811 {
2812 eclear (c);
2813 goto divsign;
2814 }
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++)
2823 {
2824 if (bi[i] != 0)
2825 {
2826 ltb -= enormlz (bi);
2827 goto dnzro1;
2828 }
2829 }
2830 eclear (c);
2831 goto divsign;
2832 }
2833 dnzro1:
2834
2835 if (ai[E] == 0)
2836 { /* possible divide by zero */
2837 for (i = 1; i < NI - 1; i++)
2838 {
2839 if (ai[i] != 0)
2840 {
2841 lta -= enormlz (ai);
2842 goto dnzro2;
2843 }
2844 }
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;
2850 }
2851 dnzro2:
2852
2853 i = edivm (ai, bi);
2854 /* calculate exponent */
2855 lt = ltb - lta + EXONE;
2856 emdnorm (bi, i, 0, lt, 64);
2857 emovo (bi, c);
2858
2859 divsign:
2860
2861 if (sign
2862#ifndef IEEE
2863 && (ecmp (c, ezero) != 0)
2864#endif
2865 )
2866 *(c+(NE-1)) |= 0x8000;
2867 else
2868 *(c+(NE-1)) &= ~0x8000;
2869}
2870
2871/* Multiply e-types A and B, return e-type product C. */
2872
2873static void
2874emul (a, b, c)
2875 unsigned EMUSHORT *a, *b, *c;
2876{
2877 unsigned EMUSHORT ai[NI], bi[NI];
2878 int i, j, sign;
2879 EMULONG lt, lta, ltb;
2880
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);
2884
2885#ifdef NANS
2886/* NaN times anything is the same NaN. */
2887 if (eisnan (a))
2888 {
2889 emov (a, c);
2890 return;
2891 }
2892 if (eisnan (b))
2893 {
2894 emov (b, c);
2895 return;
2896 }
2897/* Zero times infinity is a NaN. */
2898 if ((eisinf (a) && (ecmp (b, ezero) == 0))
2899 || (eisinf (b) && (ecmp (a, ezero) == 0)))
2900 {
2901 mtherr ("emul", INVALID);
2902 enan (c, sign);
2903 return;
2904 }
2905#endif
2906/* Infinity times anything else is infinity. */
2907#ifdef INFINITY
2908 if (eisinf (a) || eisinf (b))
2909 {
2910 einfin (c);
2911 goto mulsign;
2912 }
2913#endif
2914 emovi (a, ai);
2915 emovi (b, bi);
2916 lta = ai[E];
2917 ltb = bi[E];
2918 if (ai[E] == 0)
2919 {
2920 for (i = 1; i < NI - 1; i++)
2921 {
2922 if (ai[i] != 0)
2923 {
2924 lta -= enormlz (ai);
2925 goto mnzer1;
2926 }
2927 }
2928 eclear (c);
2929 goto mulsign;
2930 }
2931 mnzer1:
2932
2933 if (bi[E] == 0)
2934 {
2935 for (i = 1; i < NI - 1; i++)
2936 {
2937 if (bi[i] != 0)
2938 {
2939 ltb -= enormlz (bi);
2940 goto mnzer2;
2941 }
2942 }
2943 eclear (c);
2944 goto mulsign;
2945 }
2946 mnzer2:
2947
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);
2954
2955 mulsign:
2956
2957 if (sign
2958#ifndef IEEE
2959 && (ecmp (c, ezero) != 0)
2960#endif
2961 )
2962 *(c+(NE-1)) |= 0x8000;
2963 else
2964 *(c+(NE-1)) &= ~0x8000;
2965}
2966
2967/* Convert double precision PE to e-type Y. */
2968
2969static void
2970e53toe (pe, y)
2971 unsigned EMUSHORT *pe, *y;
2972{
2973#ifdef DEC
2974
2975 dectoe (pe, y);
2976
2977#else
2978#ifdef IBM
2979
2980 ibmtoe (pe, y, DFmode);
2981
2982#else
2983 register unsigned EMUSHORT r;
2984 register unsigned EMUSHORT *e, *p;
2985 unsigned EMUSHORT yy[NI];
2986 int denorm, k;
2987
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)
3001 {
3002#ifdef NANS
3003 if (! REAL_WORDS_BIG_ENDIAN)
3004 {
3005 if (((pe[3] & 0xf) != 0) || (pe[2] != 0)
3006 || (pe[1] != 0) || (pe[0] != 0))
3007 {
3008 enan (y, yy[0] != 0);
3009 return;
3010 }
3011 }
3012 else
3013 {
3014 if (((pe[0] & 0xf) != 0) || (pe[1] != 0)
3015 || (pe[2] != 0) || (pe[3] != 0))
3016 {
3017 enan (y, yy[0] != 0);
3018 return;
3019 }
3020 }
3021#endif /* NANS */
3022 eclear (y);
3023 einfin (y);
3024 if (yy[0])
3025 eneg (y);
3026 return;
3027 }
3028#endif /* INFINITY */
3029 r >>= 4;
3030 /* If zero exponent, then the significand is denormalized.
3031 So take back the understood high significand bit. */
3032
3033 if (r == 0)
3034 {
3035 denorm = 1;
3036 yy[M] &= ~0x10;
3037 }
3038 r += EXONE - 01777;
3039 yy[E] = r;
3040 p = &yy[M + 1];
3041#ifdef IEEE
3042 if (! REAL_WORDS_BIG_ENDIAN)
3043 {
3044 *p++ = *(--e);
3045 *p++ = *(--e);
3046 *p++ = *(--e);
3047 }
3048 else
3049 {
3050 ++e;
3051 *p++ = *e++;
3052 *p++ = *e++;
3053 *p++ = *e++;
3054 }
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);
3063 }
3064 emovo (yy, y);
3065#endif /* not IBM */
3066#endif /* not DEC */
3067}
3068
3069/* Convert double extended precision float PE to e type Y. */
3070
3071static void
3072e64toe (pe, y)
3073 unsigned EMUSHORT *pe, *y;
3074{
3075 unsigned EMUSHORT yy[NI];
3076 unsigned EMUSHORT *e, *p, *q;
3077 int i;
3078
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)
3097 {
3098 for (i = 0; i < 5; i++)
3099 *p++ = *e++;
3100
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)
3105 {
3106 unsigned EMUSHORT temp[NI];
3107
3108 emovi(yy, temp);
3109 eshup1(temp);
3110 emovo(temp,y);
3111 return;
3112 }
3113 }
3114 else
3115 {
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++;
3127 }
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)
3133 {
3134#ifdef NANS
3135 if (! REAL_WORDS_BIG_ENDIAN)
3136 {
3137 for (i = 0; i < 4; i++)
3138 {
3139 if ((i != 3 && pe[i] != 0)
3140 /* Anything but 0x8000 here, including 0, is a NaN. */
3141 || (i == 3 && pe[i] != 0x8000))
3142 {
3143 enan (y, (*p & 0x8000) != 0);
3144 return;
3145 }
3146 }
3147 }
3148 else
3149 {
3150#ifdef ARM_EXTENDED_IEEE_FORMAT
3151 for (i = 2; i <= 5; i++)
3152 {
3153 if (pe[i] != 0)
3154 {
3155 enan (y, (*p & 0x8000) != 0);
3156 return;
3157 }
3158 }
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;
3165
3166 for (i = 3; i <= 5; i++)
3167 {
3168 if (pe[i] != 0)
3169 {
3170bigend_nan:
3171 enan (y, (*p & 0x8000) != 0);
3172 return;
3173 }
3174 }
3175#endif /* not ARM */
3176 }
3177#endif /* NANS */
3178 eclear (y);
3179 einfin (y);
3180 if (*p & 0x8000)
3181 eneg (y);
3182 return;
3183 }
3184#endif /* INFINITY */
3185 p = yy;
3186 q = y;
3187 for (i = 0; i < NE; i++)
3188 *q++ = *p++;
3189}
3190
3191/* Convert 128-bit long double precision float PE to e type Y. */
3192
3193static void
3194e113toe (pe, y)
3195 unsigned EMUSHORT *pe, *y;
3196{
3197 register unsigned EMUSHORT r;
3198 unsigned EMUSHORT *e, *p;
3199 unsigned EMUSHORT yy[NI];
3200 int denorm, i;
3201
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)
3216 {
3217#ifdef NANS
3218 if (! REAL_WORDS_BIG_ENDIAN)
3219 {
3220 for (i = 0; i < 7; i++)
3221 {
3222 if (pe[i] != 0)
3223 {
3224 enan (y, yy[0] != 0);
3225 return;
3226 }
3227 }
3228 }
3229 else
3230 {
3231 for (i = 1; i < 8; i++)
3232 {
3233 if (pe[i] != 0)
3234 {
3235 enan (y, yy[0] != 0);
3236 return;
3237 }
3238 }
3239 }
3240#endif /* NANS */
3241 eclear (y);
3242 einfin (y);
3243 if (yy[0])
3244 eneg (y);
3245 return;
3246 }
3247#endif /* INFINITY */
3248 yy[E] = r;
3249 p = &yy[M + 1];
3250#ifdef IEEE
3251 if (! REAL_WORDS_BIG_ENDIAN)
3252 {
3253 for (i = 0; i < 7; i++)
3254 *p++ = *(--e);
3255 }
3256 else
3257 {
3258 ++e;
3259 for (i = 0; i < 7; i++)
3260 *p++ = *e++;
3261 }
3262#endif
3263/* If denormal, remove the implied bit; else shift down 1. */
3264 if (r == 0)
3265 {
3266 yy[M] = 0;
3267 }
3268 else
3269 {
3270 yy[M] = 1;
3271 eshift (yy, -1);
3272 }
3273 emovo (yy, y);
3274}
3275
3276/* Convert single precision float PE to e type Y. */
3277
3278static void
3279e24toe (pe, y)
3280 unsigned EMUSHORT *pe, *y;
3281{
3282#ifdef IBM
3283
3284 ibmtoe (pe, y, SFmode);
3285
3286#else
3287 register unsigned EMUSHORT r;
3288 register unsigned EMUSHORT *e, *p;
3289 unsigned EMUSHORT yy[NI];
3290 int denorm, k;
3291
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)
3310 {
3311#ifdef NANS
3312 if (REAL_WORDS_BIG_ENDIAN)
3313 {
3314 if (((pe[0] & 0x7f) != 0) || (pe[1] != 0))
3315 {
3316 enan (y, yy[0] != 0);
3317 return;
3318 }
3319 }
3320 else
3321 {
3322 if (((pe[1] & 0x7f) != 0) || (pe[0] != 0))
3323 {
3324 enan (y, yy[0] != 0);
3325 return;
3326 }
3327 }
3328#endif /* NANS */
3329 eclear (y);
3330 einfin (y);
3331 if (yy[0])
3332 eneg (y);
3333 return;
3334 }
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)
3340 {
3341 denorm = 1;
3342 yy[M] &= ~0200;
3343 }
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
3354 {
3355 ++e;
3356 *p++ = *e++;
3357 }
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);
3366 }
3367 emovo (yy, y);
3368#endif /* not IBM */
3369}
3370
3371/* Convert e-type X to IEEE 128-bit long double format E. */
3372
3373static void
3374etoe113 (x, e)
3375 unsigned EMUSHORT *x, *e;
3376{
3377 unsigned EMUSHORT xi[NI];
3378 EMULONG exp;
3379 int rndsav;
3380
3381#ifdef NANS
3382 if (eisnan (x))
3383 {
3384 make_nan (e, eisneg (x), TFmode);
3385 return;
3386 }
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);
3401}
3402
3403/* Convert exploded e-type X, that has already been rounded to
3404 113-bit precision, to IEEE 128-bit long double format Y. */
3405
3406static void
3407toe113 (a, b)
3408 unsigned EMUSHORT *a, *b;
3409{
3410 register unsigned EMUSHORT *p, *q;
3411 unsigned EMUSHORT i;
3412
3413#ifdef NANS
3414 if (eiisnan (a))
3415 {
3416 make_nan (b, eiisneg (a), TFmode);
3417 return;
3418 }
3419#endif
3420 p = a;
3421 if (REAL_WORDS_BIG_ENDIAN)
3422 q = b;
3423 else
3424 q = b + 7; /* point to output exponent */
3425
3426 /* If not denormal, delete the implied bit. */
3427 if (a[E] != 0)
3428 {
3429 eshup1 (a);
3430 }
3431 /* combine sign and exponent */
3432 i = *p++;
3433 if (REAL_WORDS_BIG_ENDIAN)
3434 {
3435 if (i)
3436 *q++ = *p++ | 0x8000;
3437 else
3438 *q++ = *p++;
3439 }
3440 else
3441 {
3442 if (i)
3443 *q-- = *p++ | 0x8000;
3444 else
3445 *q-- = *p++;
3446 }
3447 /* skip over guard word */
3448 ++p;
3449 /* move the significand */
3450 if (REAL_WORDS_BIG_ENDIAN)
3451 {
3452 for (i = 0; i < 7; i++)
3453 *q++ = *p++;
3454 }
3455 else
3456 {
3457 for (i = 0; i < 7; i++)
3458 *q-- = *p++;
3459 }
3460}
3461
3462/* Convert e-type X to IEEE double extended format E. */
3463
3464static void
3465etoe64 (x, e)
3466 unsigned EMUSHORT *x, *e;
3467{
3468 unsigned EMUSHORT xi[NI];
3469 EMULONG exp;
3470 int rndsav;
3471
3472#ifdef NANS
3473 if (eisnan (x))
3474 {
3475 make_nan (e, eisneg (x), XFmode);
3476 return;
3477 }
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);
3493}
3494
3495/* Convert exploded e-type X, that has already been rounded to
3496 64-bit precision, to IEEE double extended format Y. */
3497
3498static void
3499toe64 (a, b)
3500 unsigned EMUSHORT *a, *b;
3501{
3502 register unsigned EMUSHORT *p, *q;
3503 unsigned EMUSHORT i;
3504
3505#ifdef NANS
3506 if (eiisnan (a))
3507 {
3508 make_nan (b, eiisneg (a), XFmode);
3509 return;
3510 }
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
3526 {
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
3532 }
3533#endif
3534
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)
3552 {
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
3564 }
3565 else
3566 {
3567 if (i)
3568 *q-- = *p++ | 0x8000;
3569 else
3570 *q-- = *p++;
3571 }
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)
3586 {
3587 for (i = 0; i < 4; i++)
3588 *q++ = *p++;
3589 }
3590 else
3591 {
3592#ifdef INFINITY
3593 if (eiisinf (a))
3594 {
3595 /* Intel long double infinity significand. */
3596 *q-- = 0x8000;
3597 *q-- = 0;
3598 *q-- = 0;
3599 *q = 0;
3600 return;
3601 }
3602#endif
3603 for (i = 0; i < 4; i++)
3604 *q-- = *p++;
3605 }
3606#endif
3607}
3608
3609/* e type to double precision. */
3610
3611#ifdef DEC
3612/* Convert e-type X to DEC-format double E. */
3613
3614static void
3615etoe53 (x, e)
3616 unsigned EMUSHORT *x, *e;
3617{
3618 etodec (x, e); /* see etodec.c */
3619}
3620
3621/* Convert exploded e-type X, that has already been rounded to
3622 56-bit double precision, to DEC double Y. */
3623
3624static void
3625toe53 (x, y)
3626 unsigned EMUSHORT *x, *y;
3627{
3628 todec (x, y);
3629}
3630
3631#else
3632#ifdef IBM
3633/* Convert e-type X to IBM 370-format double E. */
3634
3635static void
3636etoe53 (x, e)
3637 unsigned EMUSHORT *x, *e;
3638{
3639 etoibm (x, e, DFmode);
3640}
3641
3642/* Convert exploded e-type X, that has already been rounded to
3643 56-bit precision, to IBM 370 double Y. */
3644
3645static void
3646toe53 (x, y)
3647 unsigned EMUSHORT *x, *y;
3648{
3649 toibm (x, y, DFmode);
3650}
3651
3652#else /* it's neither DEC nor IBM */
3653
3654/* Convert e-type X to IEEE double E. */
3655
3656static void
3657etoe53 (x, e)
3658 unsigned EMUSHORT *x, *e;
3659{
3660 unsigned EMUSHORT xi[NI];
3661 EMULONG exp;
3662 int rndsav;
3663
3664#ifdef NANS
3665 if (eisnan (x))
3666 {
3667 make_nan (e, eisneg (x), DFmode);
3668 return;
3669 }
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);
3685}
3686
3687/* Convert exploded e-type X, that has already been rounded to
3688 53-bit precision, to IEEE double Y. */
3689
3690static void
3691toe53 (x, y)
3692 unsigned EMUSHORT *x, *y;
3693{
3694 unsigned EMUSHORT i;
3695 unsigned EMUSHORT *p;
3696
3697#ifdef NANS
3698 if (eiisnan (x))
3699 {
3700 make_nan (y, eiisneg (x), DFmode);
3701 return;
3702 }
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 */
3712
3713 i = *p++;
3714 if (i >= (unsigned int) 2047)
3715 {
3716 /* Saturate at largest number less than infinity. */
3717#ifdef INFINITY
3718 *y |= 0x7ff0;
3719 if (! REAL_WORDS_BIG_ENDIAN)
3720 {
3721 *(--y) = 0;
3722 *(--y) = 0;
3723 *(--y) = 0;
3724 }
3725 else
3726 {
3727 ++y;
3728 *y++ = 0;
3729 *y++ = 0;
3730 *y++ = 0;
3731 }
3732#else
3733 *y |= (unsigned EMUSHORT) 0x7fef;
3734 if (! REAL_WORDS_BIG_ENDIAN)
3735 {
3736 *(--y) = 0xffff;
3737 *(--y) = 0xffff;
3738 *(--y) = 0xffff;
3739 }
3740 else
3741 {
3742 ++y;
3743 *y++ = 0xffff;
3744 *y++ = 0xffff;
3745 *y++ = 0xffff;
3746 }
3747#endif
3748 return;
3749 }
3750 if (i == 0)
3751 {
3752 eshift (x, 4);
3753 }
3754 else
3755 {
3756 i <<= 4;
3757 eshift (x, 5);
3758 }
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)
3762 {
3763 *(--y) = *p++;
3764 *(--y) = *p++;
3765 *(--y) = *p;
3766 }
3767 else
3768 {
3769 ++y;
3770 *y++ = *p++;
3771 *y++ = *p++;
3772 *y++ = *p++;
3773 }
3774}
3775
3776#endif /* not IBM */
3777#endif /* not DEC */
3778
3779
3780
3781/* e type to single precision. */
3782
3783#ifdef IBM
3784/* Convert e-type X to IBM 370 float E. */
3785
3786static void
3787etoe24 (x, e)
3788 unsigned EMUSHORT *x, *e;
3789{
3790 etoibm (x, e, SFmode);
3791}
3792
3793/* Convert exploded e-type X, that has already been rounded to
3794 float precision, to IBM 370 float Y. */
3795
3796static void
3797toe24 (x, y)
3798 unsigned EMUSHORT *x, *y;
3799{
3800 toibm (x, y, SFmode);
3801}
3802
3803#else
3804/* Convert e-type X to IEEE float E. DEC float is the same as IEEE float. */
3805
3806static void
3807etoe24 (x, e)
3808 unsigned EMUSHORT *x, *e;
3809{
3810 EMULONG exp;
3811 unsigned EMUSHORT xi[NI];
3812 int rndsav;
3813
3814#ifdef NANS
3815 if (eisnan (x))
3816 {
3817 make_nan (e, eisneg (x), SFmode);
3818 return;
3819 }
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);
3835}
3836
3837/* Convert exploded e-type X, that has already been rounded to
3838 float precision, to IEEE float Y. */
3839
3840static void
3841toe24 (x, y)
3842 unsigned EMUSHORT *x, *y;
3843{
3844 unsigned EMUSHORT i;
3845 unsigned EMUSHORT *p;
3846
3847#ifdef NANS
3848 if (eiisnan (x))
3849 {
3850 make_nan (y, eiisneg (x), SFmode);
3851 return;
3852 }
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 */
3865
3866 i = *p++;
3867/* Handle overflow cases. */
3868 if (i >= 255)
3869 {
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
3879 {
3880 ++y;
3881 *y = 0;
3882 }
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
3893 {
3894 ++y;
3895 *y = 0xffff;
3896 }
3897#endif
3898#ifdef ERANGE
3899 errno = ERANGE;
3900#endif
3901#endif /* no INFINITY */
3902 return;
3903 }
3904 if (i == 0)
3905 {
3906 eshift (x, 7);
3907 }
3908 else
3909 {
3910 i <<= 7;
3911 eshift (x, 8);
3912 }
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
3923 {
3924 ++y;
3925 *y = *p;
3926 }
3927#endif
3928}
3929#endif /* not IBM */
3930
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. */
3936
3937static int
3938ecmp (a, b)
3939 unsigned EMUSHORT *a, *b;
3940{
3941 unsigned EMUSHORT ai[NI], bi[NI];
3942 register unsigned EMUSHORT *p, *q;
3943 register int i;
3944 int msign;
3945
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;
3954
3955 if (*p != *q)
3956 { /* the signs are different */
3957 /* -0 equals + 0 */
3958 for (i = 1; i < NI - 1; i++)
3959 {
3960 if (ai[i] != 0)
3961 goto nzro;
3962 if (bi[i] != 0)
3963 goto nzro;
3964 }
3965 return (0);
3966 nzro:
3967 if (*p == 0)
3968 return (1);
3969 else
3970 return (-1);
3971 }
3972 /* both are the same sign */
3973 if (*p == 0)
3974 msign = 1;
3975 else
3976 msign = -1;
3977 i = NI - 1;
3978 do
3979 {
3980 if (*p++ != *q++)
3981 {
3982 goto diff;
3983 }
3984 }
3985 while (--i > 0);
3986
3987 return (0); /* equality */
3988
3989 diff:
3990
3991 if (*(--p) > *(--q))
3992 return (msign); /* p is bigger */
3993 else
3994 return (-msign); /* p is littler */
3995}
3996
3997/* Find e-type nearest integer to X, as floor (X + 0.5). */
3998
3999static void
4000eround (x, y)
4001 unsigned EMUSHORT *x, *y;
4002{
4003 eadd (ehalf, x, y);
4004 efloor (y, y);
4005}
4006
4007/* Convert HOST_WIDE_INT LP to e type Y. */
4008
4009static void
4010ltoe (lp, y)
4011 HOST_WIDE_INT *lp;
4012 unsigned EMUSHORT *y;
4013{
4014 unsigned EMUSHORT yi[NI];
4015 unsigned HOST_WIDE_INT ll;
4016 int k;
4017
4018 ecleaz (yi);
4019 if (*lp < 0)
4020 {
4021 /* make it positive */
4022 ll = (unsigned HOST_WIDE_INT) (-(*lp));
4023 yi[0] = 0xffff; /* put correct sign in the e type number */
4024 }
4025 else
4026 {
4027 ll = (unsigned HOST_WIDE_INT) (*lp);
4028 }
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
4041
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 */
4047}
4048
4049/* Convert unsigned HOST_WIDE_INT LP to e type Y. */
4050
4051static void
4052ultoe (lp, y)
4053 unsigned HOST_WIDE_INT *lp;
4054 unsigned EMUSHORT *y;
4055{
4056 unsigned EMUSHORT yi[NI];
4057 unsigned HOST_WIDE_INT ll;
4058 int k;
4059
4060 ecleaz (yi);
4061 ll = *lp;
4062
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
4075
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 */
4081}
4082
4083
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). */
4090
4091static void
4092eifrac (x, i, frac)
4093 unsigned EMUSHORT *x;
4094 HOST_WIDE_INT *i;
4095 unsigned EMUSHORT *frac;
4096{
4097 unsigned EMUSHORT xi[NI];
4098 int j, k;
4099 unsigned HOST_WIDE_INT ll;
4100
4101 emovi (x, xi);
4102 k = (int) xi[E] - (EXONE - 1);
4103 if (k <= 0)
4104 {
4105 /* if exponent <= 0, integer = 0 and real output is fraction */
4106 *i = 0L;
4107 emovo (xi, frac);
4108 return;
4109 }
4110 if (k > (HOST_BITS_PER_WIDE_INT - 1))
4111 {
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
4117 {
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
4127 }
4128 eshift (xi, k);
4129 if (extra_warnings)
4130 warning ("overflow on truncation to integer");
4131 }
4132 else if (k > 16)
4133 {
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;
4140 do
4141 {
4142 eshup6 (xi);
4143 ll = (ll << 16) | xi[M];
4144 }
4145 while ((k -= 16) > 0);
4146 *i = ll;
4147 if (xi[0])
4148 *i = -(*i);
4149 }
4150 else
4151 {
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);
4157 }
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;
4165
4166 emovo (xi, frac);
4167}
4168
4169
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. */
4173
4174static void
4175euifrac (x, i, frac)
4176 unsigned EMUSHORT *x;
4177 unsigned HOST_WIDE_INT *i;
4178 unsigned EMUSHORT *frac;
4179{
4180 unsigned HOST_WIDE_INT ll;
4181 unsigned EMUSHORT xi[NI];
4182 int j, k;
4183
4184 emovi (x, xi);
4185 k = (int) xi[E] - (EXONE - 1);
4186 if (k <= 0)
4187 {
4188 /* if exponent <= 0, integer = 0 and argument is fraction */
4189 *i = 0L;
4190 emovo (xi, frac);
4191 return;
4192 }
4193 if (k > HOST_BITS_PER_WIDE_INT)
4194 {
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");
4203 }
4204 else if (k > 16)
4205 {
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;
4212 do
4213 {
4214 eshup6 (xi);
4215 ll = (ll << 16) | xi[M];
4216 }
4217 while ((k -= 16) > 0);
4218 *i = ll;
4219 }
4220 else
4221 {
4222 /* shift not more than 16 bits */
4223 eshift (xi, k);
4224 *i = (HOST_WIDE_INT) xi[M] & 0xffff;
4225 }
4226
4227 if (xi[0]) /* A negative value yields unsigned integer 0. */
4228 *i = 0L;
4229
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;
4237
4238 emovo (xi, frac);
4239}
4240
4241/* Shift the significand of exploded e-type X up or down by SC bits. */
4242
4243static int
4244eshift (x, sc)
4245 unsigned EMUSHORT *x;
4246 int sc;
4247{
4248 unsigned EMUSHORT lost;
4249 unsigned EMUSHORT *p;
4250
4251 if (sc == 0)
4252 return (0);
4253
4254 lost = 0;
4255 p = x + NI - 1;
4256
4257 if (sc < 0)
4258 {
4259 sc = -sc;
4260 while (sc >= 16)
4261 {
4262 lost |= *p; /* remember lost bits */
4263 eshdn6 (x);
4264 sc -= 16;
4265 }
4266
4267 while (sc >= 8)
4268 {
4269 lost |= *p & 0xff;
4270 eshdn8 (x);
4271 sc -= 8;
4272 }
4273
4274 while (sc > 0)
4275 {
4276 lost |= *p & 1;
4277 eshdn1 (x);
4278 sc -= 1;
4279 }
4280 }
4281 else
4282 {
4283 while (sc >= 16)
4284 {
4285 eshup6 (x);
4286 sc -= 16;
4287 }
4288
4289 while (sc >= 8)
4290 {
4291 eshup8 (x);
4292 sc -= 8;
4293 }
4294
4295 while (sc > 0)
4296 {
4297 eshup1 (x);
4298 sc -= 1;
4299 }
4300 }
4301 if (lost)
4302 lost = 1;
4303 return ((int) lost);
4304}
4305
4306/* Shift normalize the significand area of exploded e-type X.
4307 Return the shift count (up = positive). */
4308
4309static int
4310enormlz (x)
4311 unsigned EMUSHORT x[];
4312{
4313 register unsigned EMUSHORT *p;
4314 int sc;
4315
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)
4324 {
4325 eshup6 (x);
4326 sc += 16;
4327
4328 /* With guard word, there are NBITS+16 bits available.
4329 Return true if all are zero. */
4330 if (sc > NBITS)
4331 return (sc);
4332 }
4333 /* see if high byte is zero */
4334 while ((*p & 0xff00) == 0)
4335 {
4336 eshup8 (x);
4337 sc += 8;
4338 }
4339 /* now shift 1 bit at a time */
4340 while ((*p & 0x8000) == 0)
4341 {
4342 eshup1 (x);
4343 sc += 1;
4344 if (sc > NBITS)
4345 {
4346 mtherr ("enormlz", UNDERFLOW);
4347 return (sc);
4348 }
4349 }
4350 return (sc);
4351
4352 /* Normalize by shifting down out of the high guard word
4353 of the significand */
4354 normdn:
4355
4356 if (*p & 0xff00)
4357 {
4358 eshdn8 (x);
4359 sc -= 8;
4360 }
4361 while (*p != 0)
4362 {
4363 eshdn1 (x);
4364 sc -= 1;
4365
4366 if (sc < -NBITS)
4367 {
4368 mtherr ("enormlz", OVERFLOW);
4369 return (sc);
4370 }
4371 }
4372 return (sc);
4373}
4374
4375/* Powers of ten used in decimal <-> binary conversions. */
4376
4377#define NTEN 12
4378#define MAXP 4096
4379
4380#if LONG_DOUBLE_TYPE_SIZE == 128
4381static unsigned EMUSHORT etens[NTEN + 1][NE] =
4382{
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 */
4409};
4410
4411static unsigned EMUSHORT emtens[NTEN + 1][NE] =
4412{
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 */
4439};
4440#else
4441/* LONG_DOUBLE_TYPE_SIZE is other than 128 */
4442static unsigned EMUSHORT etens[NTEN + 1][NE] =
4443{
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 */
4457};
4458
4459static unsigned EMUSHORT emtens[NTEN + 1][NE] =
4460{
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 */
4474};
4475#endif
4476
4477/* Convert float value X to ASCII string STRING with NDIG digits after
4478 the decimal point. */
4479
4480static void
4481e24toasc (x, string, ndigs)
4482 unsigned EMUSHORT x[];
4483 char *string;
4484 int ndigs;
4485{
4486 unsigned EMUSHORT w[NI];
4487
4488 e24toe (x, w);
4489 etoasc (w, string, ndigs);
4490}
4491
4492/* Convert double value X to ASCII string STRING with NDIG digits after
4493 the decimal point. */
4494
4495static void
4496e53toasc (x, string, ndigs)
4497 unsigned EMUSHORT x[];
4498 char *string;
4499 int ndigs;
4500{
4501 unsigned EMUSHORT w[NI];
4502
4503 e53toe (x, w);
4504 etoasc (w, string, ndigs);
4505}
4506
4507/* Convert double extended value X to ASCII string STRING with NDIG digits
4508 after the decimal point. */
4509
4510static void
4511e64toasc (x, string, ndigs)
4512 unsigned EMUSHORT x[];
4513 char *string;
4514 int ndigs;
4515{
4516 unsigned EMUSHORT w[NI];
4517
4518 e64toe (x, w);
4519 etoasc (w, string, ndigs);
4520}
4521
4522/* Convert 128-bit long double value X to ASCII string STRING with NDIG digits
4523 after the decimal point. */
4524
4525static void
4526e113toasc (x, string, ndigs)
4527 unsigned EMUSHORT x[];
4528 char *string;
4529 int ndigs;
4530{
4531 unsigned EMUSHORT w[NI];
4532
4533 e113toe (x, w);
4534 etoasc (w, string, ndigs);
4535}
4536
4537/* Convert e-type X to ASCII string STRING with NDIGS digits after
4538 the decimal point. */
4539
4540static char wstring[80]; /* working storage for ASCII output */
4541
4542static void
4543etoasc (x, string, ndigs)
4544 unsigned EMUSHORT x[];
4545 char *string;
4546 int ndigs;
4547{
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;
4555
4556
4557 rndsav = rndprc;
4558 ss = string;
4559 s = wstring;
4560 *ss = '\0';
4561 *s = '\0';
4562#ifdef NANS
4563 if (eisnan (x))
4564 {
4565 sprintf (wstring, " NaN ");
4566 goto bxit;
4567 }
4568#endif
4569 rndprc = NBITS; /* set to full precision */
4570 emov (x, y); /* retain external format */
4571 if (y[NE - 1] & 0x8000)
4572 {
4573 sign = 0xffff;
4574 y[NE - 1] &= 0x7fff;
4575 }
4576 else
4577 {
4578 sign = 0;
4579 }
4580 expon = 0;
4581 ten = &etens[NTEN][0];
4582 emov (eone, t);
4583 /* Test for zero exponent */
4584 if (y[NE - 1] == 0)
4585 {
4586 for (k = 0; k < NE - 1; k++)
4587 {
4588 if (y[k] != 0)
4589 goto tnzro; /* denormalized number */
4590 }
4591 goto isone; /* valid all zeros */
4592 }
4593 tnzro:
4594
4595 /* Test for infinity. */
4596 if (y[NE - 1] == 0x7fff)
4597 {
4598 if (sign)
4599 sprintf (wstring, " -Infinity ");
4600 else
4601 sprintf (wstring, " Infinity ");
4602 goto bxit;
4603 }
4604
4605 /* Test for exponent nonzero but significand denormalized.
4606 * This is an error condition.
4607 */
4608 if ((y[NE - 1] != 0) && ((y[NE - 2] & 0x8000) == 0))
4609 {
4610 mtherr ("etoasc", DOMAIN);
4611 sprintf (wstring, "NaN");
4612 goto bxit;
4613 }
4614
4615 /* Compare to 1.0 */
4616 i = ecmp (eone, y);
4617 if (i == 0)
4618 goto isone;
4619
4620 if (i == -2)
4621 abort ();
4622
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;
4628
4629 p = &etens[NTEN - 4][0];
4630 m = 16;
4631 do
4632 {
4633 ediv (p, u, t);
4634 efloor (t, w);
4635 for (j = 0; j < NE - 1; j++)
4636 {
4637 if (t[j] != w[j])
4638 goto noint;
4639 }
4640 emov (t, u);
4641 expon += (int) m;
4642 noint:
4643 p += NE;
4644 m >>= 1;
4645 }
4646 while (m != 0);
4647
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)
4657 {
4658 if (ecmp (p, u) <= 0)
4659 {
4660 ediv (p, u, u);
4661 emul (p, t, t);
4662 expon += (int) m;
4663 }
4664 m >>= 1;
4665 if (m == 0)
4666 break;
4667 p += NE;
4668 }
4669 }
4670 else
4671 { /* Number is less than 1.0 */
4672 /* Pad significand with trailing decimal zeros. */
4673 if (y[NE - 1] == 0)
4674 {
4675 while ((y[NE - 2] & 0x8000) == 0)
4676 {
4677 emul (ten, y, y);
4678 expon -= 1;
4679 }
4680 }
4681 else
4682 {
4683 emovi (y, w);
4684 for (i = 0; i < NDEC + 1; i++)
4685 {
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)
4695 {
4696 eshdn1 (u);
4697 u[1] += 1;
4698 }
4699 if (u[NI - 1] != 0)
4700 break;
4701 if (eone[NE - 1] <= u[1])
4702 break;
4703 emovz (u, w);
4704 expon -= 1;
4705 }
4706 emovo (w, y);
4707 }
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)
4714 {
4715 if (ecmp (p, w) >= 0)
4716 {
4717 emul (r, w, w);
4718 emul (r, t, t);
4719 expon += k;
4720 }
4721 k /= 2;
4722 if (k == 0)
4723 break;
4724 p += NE;
4725 r += NE;
4726 }
4727 ediv (t, eone, t);
4728 }
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))
4738 {
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;
4747 }
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)
4759 {
4760 *s++ = '1';
4761 *s++ = '.';
4762 if (ndigs > 0)
4763 {
4764 *s++ = '0';
4765 ndigs -= 1;
4766 }
4767 expon += 1;
4768 }
4769 else
4770 {
4771 *s++ = (char)digit + '0';
4772 *s++ = '.';
4773 }
4774 /* Generate digits after the decimal point. */
4775 for (k = 0; k <= ndigs; k++)
4776 {
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';
4785 }
4786 digit = equot[NI - 1];
4787 --s;
4788 ss = s;
4789 /* round off the ASCII string */
4790 if (digit > 4)
4791 {
4792 /* Test for critical rounding case in ASCII output. */
4793 if (digit == 5)
4794 {
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 */
4800 }
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 == '.')
4807 {
4808 --s;
4809 k = *s;
4810 k += 1;
4811 *s = (char) k;
4812 /* Most significant digit carries to 10? */
4813 if (k > '9')
4814 {
4815 expon += 1;
4816 *s = '1';
4817 }
4818 goto doexp;
4819 }
4820 /* Round up and carry out from less significant digits */
4821 k += 1;
4822 *s = (char) k;
4823 if (k > '9')
4824 {
4825 *s = '0';
4826 goto roun;
4827 }
4828 }
4829 doexp:
4830 /*
4831 if (expon >= 0)
4832 sprintf (ss, "e+%d", expon);
4833 else
4834 sprintf (ss, "e%d", expon);
4835 */
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')
4845 ;
4846}
4847
4848
4849/* Convert ASCII string to floating point.
4850
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). */
4855
4856/* Convert ASCII string S to single precision float value Y. */
4857
4858static void
4859asctoe24 (s, y)
4860 char *s;
4861 unsigned EMUSHORT *y;
4862{
4863 asctoeg (s, y, 24);
4864}
4865
4866
4867/* Convert ASCII string S to double precision value Y. */
4868
4869static void
4870asctoe53 (s, y)
4871 char *s;
4872 unsigned EMUSHORT *y;
4873{
4874#if defined(DEC) || defined(IBM)
4875 asctoeg (s, y, 56);
4876#else
4877 asctoeg (s, y, 53);
4878#endif
4879}
4880
4881
4882/* Convert ASCII string S to double extended value Y. */
4883
4884static void
4885asctoe64 (s, y)
4886 char *s;
4887 unsigned EMUSHORT *y;
4888{
4889 asctoeg (s, y, 64);
4890}
4891
4892/* Convert ASCII string S to 128-bit long double Y. */
4893
4894static void
4895asctoe113 (s, y)
4896 char *s;
4897 unsigned EMUSHORT *y;
4898{
4899 asctoeg (s, y, 113);
4900}
4901
4902/* Convert ASCII string S to e type Y. */
4903
4904static void
4905asctoe (s, y)
4906 char *s;
4907 unsigned EMUSHORT *y;
4908{
4909 asctoeg (s, y, NBITS);
4910}
4911
4912/* Convert ASCII string SS to e type Y, with a specified rounding precision
4913 of OPREC bits. */
4914
4915static void
4916asctoeg (ss, y, oprec)
4917 char *ss;
4918 unsigned EMUSHORT *y;
4919 int oprec;
4920{
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;
4927
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')
4935 ;
4936 s = lstr;
4937
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;
4949
4950 nxtcom:
4951 k = *s - '0';
4952 if ((k >= 0) && (k <= 9))
4953 {
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))
4959 {
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;
4975 }
4976
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. */
4981
4982 if (yy[2] == 0)
4983 {
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);
4994 }
4995 else
4996 {
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;
5002 }
5003 prec += 1;
5004 goto donchr;
5005 }
5006
5007 switch (*s)
5008 {
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;
5048 }
5049 donchr:
5050 ++s;
5051 goto nxtcom;
5052
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++)
5057 {
5058 if (yy[k] != 0)
5059 goto read_expnt;
5060 }
5061 goto aexit;
5062
5063read_expnt:
5064 esign = 1;
5065 exp = 0;
5066 ++s;
5067 /* check for + or - */
5068 if (*s == '-')
5069 {
5070 esign = -1;
5071 ++s;
5072 }
5073 if (*s == '+')
5074 ++s;
5075 while ((*s >= '0') && (*s <= '9'))
5076 {
5077 exp *= 10;
5078 exp += *s++ - '0';
5079 if (exp > -(MINDECEXP))
5080 {
5081 if (esign < 0)
5082 goto zero;
5083 else
5084 goto infinite;
5085 }
5086 }
5087 if (esign < 0)
5088 exp = -exp;
5089 if (exp > MAXDECEXP)
5090 {
5091 infinite:
5092 ecleaz (yy);
5093 yy[E] = 0x7fff; /* infinity */
5094 goto aexit;
5095 }
5096 if (exp < MINDECEXP)
5097 {
5098 zero:
5099 ecleaz (yy);
5100 goto aexit;
5101 }
5102
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))
5107 {
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);
5117 }
5118 if ((k = enormlz (yy)) > NBITS)
5119 {
5120 ecleaz (yy);
5121 goto aexit;
5122 }
5123 lexp = (EXONE - 1 + NBITS) - k;
5124 emdnorm (yy, lost, 0, lexp, 64);
5125
5126 /* Convert to external format:
5127
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. */
5133
5134 lexp = yy[E];
5135 if (nexp == 0)
5136 {
5137 k = 0;
5138 goto expdon;
5139 }
5140 esign = 1;
5141 if (nexp < 0)
5142 {
5143 nexp = -nexp;
5144 esign = -1;
5145 if (nexp > 4096)
5146 {
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;
5153 }
5154 }
5155 p = &etens[NTEN][0];
5156 emov (eone, xt);
5157 exp = 1;
5158 do
5159 {
5160 if (exp & nexp)
5161 emul (p, xt, xt);
5162 p -= NE;
5163 exp = exp + exp;
5164 }
5165 while (exp <= MAXP);
5166
5167 emovi (xt, tt);
5168 if (esign < 0)
5169 {
5170 lexp -= tt[E];
5171 k = edivm (tt, yy);
5172 lexp += EXONE;
5173 }
5174 else
5175 {
5176 lexp += tt[E];
5177 k = emulm (tt, yy);
5178 lexp -= EXONE - 1;
5179 }
5180
5181 expdon:
5182
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);
5199
5200 aexit:
5201
5202 rndprc = rndsav;
5203 yy[0] = nsign;
5204 switch (oprec)
5205 {
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;
5231 }
5232}
5233
5234
5235
5236/* Return Y = largest integer not greater than X (truncated toward minus
5237 infinity). */
5238
5239static unsigned EMUSHORT bmask[] =
5240{
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,
5258};
5259
5260static void
5261efloor (x, y)
5262 unsigned EMUSHORT x[], y[];
5263{
5264 register unsigned EMUSHORT *p;
5265 int e, expon, i;
5266 unsigned EMUSHORT f[NE];
5267
5268 emov (x, f); /* leave in external format */
5269 expon = (int) f[NE - 1];
5270 e = (expon & 0x7fff) - (EXONE - 1);
5271 if (e <= 0)
5272 {
5273 eclear (y);
5274 goto isitneg;
5275 }
5276 /* number of bits to clear out */
5277 e = NBITS - e;
5278 emov (f, y);
5279 if (e <= 0)
5280 return;
5281
5282 p = &y[0];
5283 while (e >= 16)
5284 {
5285 *p++ = 0;
5286 e -= 16;
5287 }
5288 /* clear the remaining bits */
5289 *p &= bmask[e];
5290 /* truncate negatives toward minus infinity */
5291 isitneg:
5292
5293 if ((unsigned EMUSHORT) expon & (unsigned EMUSHORT) 0x8000)
5294 {
5295 for (i = 0; i < NE - 1; i++)
5296 {
5297 if (f[i] != y[i])
5298 {
5299 esub (eone, y, y);
5300 break;
5301 }
5302 }
5303 }
5304}
5305
5306
5307/* Return S and EXP such that S * 2^EXP = X and .5 <= S < 1.
5308 For example, 1.1 = 0.55 * 2^1. */
5309
5310static void
5311efrexp (x, exp, s)
5312 unsigned EMUSHORT x[];
5313 int *exp;
5314 unsigned EMUSHORT s[];
5315{
5316 unsigned EMUSHORT xi[NI];
5317 EMULONG li;
5318
5319 emovi (x, xi);
5320 /* Handle denormalized numbers properly using long integer exponent. */
5321 li = (EMULONG) ((EMUSHORT) xi[1]);
5322
5323 if (li == 0)
5324 {
5325 li -= enormlz (xi);
5326 }
5327 xi[1] = 0x3ffe;
5328 emovo (xi, s);
5329 *exp = (int) (li - 0x3ffe);
5330}
5331
5332/* Return e type Y = X * 2^PWR2. */
5333
5334static void
5335eldexp (x, pwr2, y)
5336 unsigned EMUSHORT x[];
5337 int pwr2;
5338 unsigned EMUSHORT y[];
5339{
5340 unsigned EMUSHORT xi[NI];
5341 EMULONG li;
5342 int i;
5343
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);
5350}
5351
5352
5353/* C = remainder after dividing B by A, all e type values.
5354 Least significant integer quotient bits left in EQUOT. */
5355
5356static void
5357eremain (a, b, c)
5358 unsigned EMUSHORT a[], b[], c[];
5359{
5360 unsigned EMUSHORT den[NI], num[NI];
5361
5362#ifdef NANS
5363 if (eisinf (b)
5364 || (ecmp (a, ezero) == 0)
5365 || eisnan (a)
5366 || eisnan (b))
5367 {
5368 enan (c, 0);
5369 return;
5370 }
5371#endif
5372 if (ecmp (a, ezero) == 0)
5373 {
5374 mtherr ("eremain", SING);
5375 eclear (c);
5376 return;
5377 }
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);
5387}
5388
5389/* Return quotient of exploded e-types NUM / DEN in EQUOT,
5390 remainder in NUM. */
5391
5392static void
5393eiremain (den, num)
5394 unsigned EMUSHORT den[], num[];
5395{
5396 EMULONG ld, ln;
5397 unsigned EMUSHORT j;
5398
5399 ld = den[E];
5400 ld -= enormlz (den);
5401 ln = num[E];
5402 ln -= enormlz (num);
5403 ecleaz (equot);
5404 while (ln >= ld)
5405 {
5406 if (ecmpm (den, num) <= 0)
5407 {
5408 esubm (den, num);
5409 j = 1;
5410 }
5411 else
5412 j = 0;
5413 eshup1 (equot);
5414 equot[NI - 1] |= j;
5415 eshup1 (num);
5416 ln -= 1;
5417 }
5418 emdnorm (num, 0, 0, ln, 0);
5419}
5420
5421/* Report an error condition CODE encountered in function NAME.
5422 CODE is one of the following:
5423
5424 Mnemonic Value Significance
5425
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
5435
5436 The order of appearance of the following messages is bound to the
5437 error codes defined above. */
5438
5439#define NMSGS 8
5440static char *ermsg[NMSGS] =
5441{
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"
5450};
5451
5452int merror = 0;
5453extern int merror;
5454
5455static void
5456mtherr (name, code)
5457 char *name;
5458 int code;
5459{
5460 char errstr[80];
5461
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. */
5465
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;
5473}
5474
5475#ifdef DEC
5476/* Convert DEC double precision D to e type E. */
5477
5478static void
5479dectoe (d, e)
5480 unsigned EMUSHORT *d;
5481 unsigned EMUSHORT *e;
5482{
5483 unsigned EMUSHORT y[NI];
5484 register unsigned EMUSHORT r, *p;
5485
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;
5495
5496
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 */
5501
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 */
5506
5507 *p++ = *d++; /* fill in the rest of our mantissa */
5508 *p++ = *d++;
5509 *p = *d;
5510
5511 eshdn8 (y); /* shift our mantissa down 8 bits */
5512 done:
5513 emovo (y, e);
5514}
5515
5516/* Convert e type X to DEC double precision D. */
5517
5518static void
5519etodec (x, d)
5520 unsigned EMUSHORT *x, *d;
5521{
5522 unsigned EMUSHORT xi[NI];
5523 EMULONG exp;
5524 int rndsav;
5525
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);
5535}
5536
5537/* Convert exploded e-type X, that has already been rounded to
5538 56-bit precision, to DEC format double Y. */
5539
5540static void
5541todec (x, y)
5542 unsigned EMUSHORT *x, *y;
5543{
5544 unsigned EMUSHORT i;
5545 unsigned EMUSHORT *p;
5546
5547 p = x;
5548 *y = 0;
5549 if (*p++)
5550 *y = 0100000;
5551 i = *p++;
5552 if (i == 0)
5553 {
5554 *y++ = 0;
5555 *y++ = 0;
5556 *y++ = 0;
5557 *y++ = 0;
5558 return;
5559 }
5560 if (i > 0377)
5561 {
5562 *y++ |= 077777;
5563 *y++ = 0xffff;
5564 *y++ = 0xffff;
5565 *y++ = 0xffff;
5566#ifdef ERANGE
5567 errno = ERANGE;
5568#endif
5569 return;
5570 }
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];
5580}
5581#endif /* DEC */
5582
5583#ifdef IBM
5584/* Convert IBM single/double precision to e type. */
5585
5586static void
5587ibmtoe (d, e, mode)
5588 unsigned EMUSHORT *d;
5589 unsigned EMUSHORT *e;
5590 enum machine_mode mode;
5591{
5592 unsigned EMUSHORT y[NI];
5593 register unsigned EMUSHORT r, *p;
5594 int rndsav;
5595
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 */
5608
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 */
5612 {
5613 *p++ = *d++; /* fill in the rest of our mantissa */
5614 *p++ = *d++;
5615 }
5616 *p = *d;
5617
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);
5624}
5625
5626
5627
5628/* Convert e type to IBM single/double precision. */
5629
5630static void
5631etoibm (x, d, mode)
5632 unsigned EMUSHORT *x, *d;
5633 enum machine_mode mode;
5634{
5635 unsigned EMUSHORT xi[NI];
5636 EMULONG exp;
5637 int rndsav;
5638
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);
5647}
5648
5649static void
5650toibm (x, y, mode)
5651 unsigned EMUSHORT *x, *y;
5652 enum machine_mode mode;
5653{
5654 unsigned EMUSHORT i;
5655 unsigned EMUSHORT *p;
5656 int r;
5657
5658 p = x;
5659 *y = 0;
5660 if (*p++)
5661 *y = 0x8000;
5662 i = *p++;
5663 if (i == 0)
5664 {
5665 *y++ = 0;
5666 *y++ = 0;
5667 if (mode != SFmode)
5668 {
5669 *y++ = 0;
5670 *y++ = 0;
5671 }
5672 return;
5673 }
5674 r = i & 0x3;
5675 i >>= 2;
5676 if (i > 0x7f)
5677 {
5678 *y++ |= 0x7fff;
5679 *y++ = 0xffff;
5680 if (mode != SFmode)
5681 {
5682 *y++ = 0xffff;
5683 *y++ = 0xffff;
5684 }
5685#ifdef ERANGE
5686 errno = ERANGE;
5687#endif
5688 return;
5689 }
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)
5696 {
5697 *y++ = x[M + 2];
5698 *y++ = x[M + 3];
5699 }
5700}
5701#endif /* IBM */
5702
5703/* Output a binary NaN bit pattern in the target machine's format. */
5704
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
5709TFMODE_NAN;
5710#else
5711#ifdef IEEE
5712unsigned EMUSHORT TFbignan[8] =
5713 {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
5714unsigned EMUSHORT TFlittlenan[8] = {0, 0, 0, 0, 0, 0, 0x8000, 0xffff};
5715#endif
5716#endif
5717
5718#ifdef XFMODE_NAN
5719XFMODE_NAN;
5720#else
5721#ifdef IEEE
5722unsigned EMUSHORT XFbignan[6] =
5723 {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
5724unsigned EMUSHORT XFlittlenan[6] = {0, 0, 0, 0xc000, 0xffff, 0};
5725#endif
5726#endif
5727
5728#ifdef DFMODE_NAN
5729DFMODE_NAN;
5730#else
5731#ifdef IEEE
5732unsigned EMUSHORT DFbignan[4] = {0x7fff, 0xffff, 0xffff, 0xffff};
5733unsigned EMUSHORT DFlittlenan[4] = {0, 0, 0, 0xfff8};
5734#endif
5735#endif
5736
5737#ifdef SFMODE_NAN
5738SFMODE_NAN;
5739#else
5740#ifdef IEEE
5741unsigned EMUSHORT SFbignan[2] = {0x7fff, 0xffff};
5742unsigned EMUSHORT SFlittlenan[2] = {0, 0xffc0};
5743#endif
5744#endif
5745
5746
5747static void
5748make_nan (nan, sign, mode)
5749 unsigned EMUSHORT *nan;
5750 int sign;
5751 enum machine_mode mode;
5752{
5753 int n;
5754 unsigned EMUSHORT *p;
5755
5756 switch (mode)
5757 {
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 ();
5793 }
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;
5800}
5801
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. */
5805
5806REAL_VALUE_TYPE
5807ereal_from_float (f)
5808 HOST_WIDE_INT f;
5809{
5810 REAL_VALUE_TYPE r;
5811 unsigned EMUSHORT s[2];
5812 unsigned EMUSHORT e[NE];
5813
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)
5817 {
5818 s[0] = (unsigned EMUSHORT) (f >> 16);
5819 s[1] = (unsigned EMUSHORT) f;
5820 }
5821 else
5822 {
5823 s[0] = (unsigned EMUSHORT) f;
5824 s[1] = (unsigned EMUSHORT) (f >> 16);
5825 }
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;
5831}
5832
5833
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.
5837
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. */
5842
5843REAL_VALUE_TYPE
5844ereal_from_double (d)
5845 HOST_WIDE_INT d[];
5846{
5847 REAL_VALUE_TYPE r;
5848 unsigned EMUSHORT s[4];
5849 unsigned EMUSHORT e[NE];
5850
5851 /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces. */
5852 if (REAL_WORDS_BIG_ENDIAN)
5853 {
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
5866 }
5867 else
5868 {
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
5879 }
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;
5885}
5886
5887
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. */
5891
5892static void
5893uditoe (di, e)
5894 unsigned EMUSHORT *di; /* Address of the 64-bit int. */
5895 unsigned EMUSHORT *e;
5896{
5897 unsigned EMUSHORT yi[NI];
5898 int k;
5899
5900 ecleaz (yi);
5901 if (WORDS_BIG_ENDIAN)
5902 {
5903 for (k = M; k < M + 4; k++)
5904 yi[k] = *di++;
5905 }
5906 else
5907 {
5908 for (k = M + 3; k >= M; k--)
5909 yi[k] = *di++;
5910 }
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);
5917}
5918
5919/* Convert target computer signed 64-bit integer to e-type. */
5920
5921static void
5922ditoe (di, e)
5923 unsigned EMUSHORT *di; /* Address of the 64-bit int. */
5924 unsigned EMUSHORT *e;
5925{
5926 unsigned EMULONG acc;
5927 unsigned EMUSHORT yi[NI];
5928 unsigned EMUSHORT carry;
5929 int k, sign;
5930
5931 ecleaz (yi);
5932 if (WORDS_BIG_ENDIAN)
5933 {
5934 for (k = M; k < M + 4; k++)
5935 yi[k] = *di++;
5936 }
5937 else
5938 {
5939 for (k = M + 3; k >= M; k--)
5940 yi[k] = *di++;
5941 }
5942 /* Take absolute value */
5943 sign = 0;
5944 if (yi[M] & 0x8000)
5945 {
5946 sign = 1;
5947 carry = 0;
5948 for (k = M + 3; k >= M; k--)
5949 {
5950 acc = (unsigned EMULONG) (~yi[k] & 0xffff) + carry;
5951 yi[k] = acc;
5952 carry = 0;
5953 if (acc & 0x10000)
5954 carry = 1;
5955 }
5956 }
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);
5965}
5966
5967
5968/* Convert e-type to unsigned 64-bit int. */
5969
5970static void
5971etoudi (x, i)
5972 unsigned EMUSHORT *x;
5973 unsigned EMUSHORT *i;
5974{
5975 unsigned EMUSHORT xi[NI];
5976 int j, k;
5977
5978 emovi (x, xi);
5979 if (xi[0])
5980 {
5981 xi[M] = 0;
5982 goto noshift;
5983 }
5984 k = (int) xi[E] - (EXONE - 1);
5985 if (k <= 0)
5986 {
5987 for (j = 0; j < 4; j++)
5988 *i++ = 0;
5989 return;
5990 }
5991 if (k > 64)
5992 {
5993 for (j = 0; j < 4; j++)
5994 *i++ = 0xffff;
5995 if (extra_warnings)
5996 warning ("overflow on truncation to integer");
5997 return;
5998 }
5999 if (k > 16)
6000 {
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
6010 {
6011 i += 3;
6012 *i-- = xi[M];
6013 }
6014 k -= j;
6015 do
6016 {
6017 eshup6 (xi);
6018 if (WORDS_BIG_ENDIAN)
6019 *i++ = xi[M];
6020 else
6021 *i-- = xi[M];
6022 }
6023 while ((k -= 16) > 0);
6024 }
6025 else
6026 {
6027 /* shift not more than 16 bits */
6028 eshift (xi, k);
6029
6030noshift:
6031
6032 if (WORDS_BIG_ENDIAN)
6033 {
6034 i += 3;
6035 *i-- = xi[M];
6036 *i-- = 0;
6037 *i-- = 0;
6038 *i = 0;
6039 }
6040 else
6041 {
6042 *i++ = xi[M];
6043 *i++ = 0;
6044 *i++ = 0;
6045 *i = 0;
6046 }
6047 }
6048}
6049
6050
6051/* Convert e-type to signed 64-bit int. */
6052
6053static void
6054etodi (x, i)
6055 unsigned EMUSHORT *x;
6056 unsigned EMUSHORT *i;
6057{
6058 unsigned EMULONG acc;
6059 unsigned EMUSHORT xi[NI];
6060 unsigned EMUSHORT carry;
6061 unsigned EMUSHORT *isave;
6062 int j, k;
6063
6064 emovi (x, xi);
6065 k = (int) xi[E] - (EXONE - 1);
6066 if (k <= 0)
6067 {
6068 for (j = 0; j < 4; j++)
6069 *i++ = 0;
6070 return;
6071 }
6072 if (k > 64)
6073 {
6074 for (j = 0; j < 4; j++)
6075 *i++ = 0xffff;
6076 if (extra_warnings)
6077 warning ("overflow on truncation to integer");
6078 return;
6079 }
6080 isave = i;
6081 if (k > 16)
6082 {
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
6092 {
6093 i += 3;
6094 *i-- = xi[M];
6095 }
6096 k -= j;
6097 do
6098 {
6099 eshup6 (xi);
6100 if (WORDS_BIG_ENDIAN)
6101 *i++ = xi[M];
6102 else
6103 *i-- = xi[M];
6104 }
6105 while ((k -= 16) > 0);
6106 }
6107 else
6108 {
6109 /* shift not more than 16 bits */
6110 eshift (xi, k);
6111
6112 if (WORDS_BIG_ENDIAN)
6113 {
6114 i += 3;
6115 *i = xi[M];
6116 *i-- = 0;
6117 *i-- = 0;
6118 *i = 0;
6119 }
6120 else
6121 {
6122 *i++ = xi[M];
6123 *i++ = 0;
6124 *i++ = 0;
6125 *i = 0;
6126 }
6127 }
6128 /* Negate if negative */
6129 if (xi[0])
6130 {
6131 carry = 0;
6132 if (WORDS_BIG_ENDIAN)
6133 isave += 3;
6134 for (k = 0; k < 4; k++)
6135 {
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;
6144 }
6145 }
6146}
6147
6148
6149/* Longhand square root routine. */
6150
6151
6152static int esqinited = 0;
6153static unsigned short sqrndbit[NI];
6154
6155static void
6156esqrt (x, y)
6157 unsigned EMUSHORT *x, *y;
6158{
6159 unsigned EMUSHORT temp[NI], num[NI], sq[NI], xx[NI];
6160 EMULONG m, exp;
6161 int i, j, k, n, nlups;
6162
6163 if (esqinited == 0)
6164 {
6165 ecleaz (sqrndbit);
6166 sqrndbit[NI - 2] = 1;
6167 esqinited = 1;
6168 }
6169 /* Check for arg <= 0 */
6170 i = ecmp (x, ezero);
6171 if (i <= 0)
6172 {
6173 if (i == -1)
6174 {
6175 mtherr ("esqrt", DOMAIN);
6176 eclear (y);
6177 }
6178 else
6179 emov (x, y);
6180 return;
6181 }
6182
6183#ifdef INFINITY
6184 if (eisinf (x))
6185 {
6186 eclear (y);
6187 einfin (y);
6188 return;
6189 }
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);
6196
6197 /* Divide exponent by 2 */
6198 m -= 0x3ffe;
6199 exp = (unsigned short) ((m / 2) + 0x3ffe);
6200
6201 /* Adjust if exponent odd */
6202 if ((m & 1) != 0)
6203 {
6204 if (m > 0)
6205 exp += 1;
6206 eshdn1 (xx);
6207 }
6208
6209 ecleaz (sq);
6210 ecleaz (num);
6211 n = 8; /* get 8 bits of result per inner loop */
6212 nlups = rndprc;
6213 j = 0;
6214
6215 while (nlups > 0)
6216 {
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++)
6224 {
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)
6237 {
6238 esubm (temp, num);
6239 sq[NI - 2] |= 1;
6240 }
6241 }
6242 nlups -= n;
6243 j += 1;
6244 }
6245
6246 /* Adjust for extra, roundoff loop done. */
6247 exp += (NBITS - 1) - rndprc;
6248
6249 /* Sticky bit = 1 if the remainder is nonzero. */
6250 k = 0;
6251 for (i = 3; i < NI; i++)
6252 k |= (int) num[i];
6253
6254 /* Renormalize and round off. */
6255 emdnorm (sq, k, 0, exp, 64);
6256 emovo (sq, y);
6257}
6258#endif /* EMU_NON_COMPILE not defined */
6259\f
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. */
6263
6264int
6265significand_size (mode)
6266 enum machine_mode mode;
6267{
6268
6269/* Don't test the modes, but their sizes, lest this
6270 code won't work for BITS_PER_UNIT != 8 . */
6271
6272switch (GET_MODE_BITSIZE (mode))
6273 {
6274 case 32:
6275 return 24;
6276
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
6291
6292 case 96:
6293 return 64;
6294 case 128:
6295 return 113;
6296
6297 default:
6298 abort ();
6299 }
6300}
GNU Compiler Collection
This page took 0.304842 seconds and 5 git commands to generate.