]> gcc.gnu.org Git - gcc.git/blame - gcc/real.c
git://gcc.gnu.org / gcc.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
pa.h (ASM_DECLARE_FUNCTION_NAME): Handle invisible reference arguments passed in...
[gcc.git] / gcc / real.c
CommitLineData
985b6196 1/* real.c - implementation of REAL_ARITHMETIC, REAL_VALUE_ATOF,
29e11dab 2 and support for XFmode IEEE extended real floating point arithmetic.
82e974d4 3 Copyright (C) 1993, 1994, 1995 Free Software Foundation, Inc.
c764eafd 4 Contributed by Stephen L. Moshier (moshier@world.std.com).
985b6196
RS		 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, 675 Mass Ave, Cambridge, MA 02139, USA. */
21
22#include <stdio.h>
64685ffa 23#include <errno.h>
985b6196 24#include "config.h"
985b6196
RS		 25#include "tree.h"
26
64685ffa
RS		 27#ifndef errno
28extern int errno;
29#endif
30
985b6196
RS		 31/* To enable support of XFmode extended real floating point, define
32LONG_DOUBLE_TYPE_SIZE 96 in the tm.h file (m68k.h or i386.h).
33
842fbaaa 34To support cross compilation between IEEE, VAX and IBM floating
985b6196
RS		 35point formats, define REAL_ARITHMETIC in the tm.h file.
36
37In either case the machine files (tm.h) must not contain any code
38that tries to use host floating point arithmetic to convert
39REAL_VALUE_TYPEs from `double' to `float', pass them to fprintf,
40etc. In cross-compile situations a REAL_VALUE_TYPE may not
41be intelligible to the host computer's native arithmetic.
42
43The emulator defaults to the host's floating point format so that
44its decimal conversion functions can be used if desired (see
45real.h).
46
47The first part of this file interfaces gcc to ieee.c, which is a
48floating point arithmetic suite that was not written with gcc in
49mind. The interface is followed by ieee.c itself and related
50items. Avoid changing ieee.c unless you have suitable test
51programs available. A special version of the PARANOIA floating
52point arithmetic tester, modified for this purpose, can be found
53on usc.edu : /pub/C-numanal/ieeetest.zoo. Some tutorial
54information on ieee.c is given in my book: S. L. Moshier,
55_Methods and Programs for Mathematical Functions_, Prentice-Hall
56or Simon & Schuster Int'l, 1989. A library of XFmode elementary
57transcendental functions can be obtained by ftp from
58research.att.com: netlib/cephes/ldouble.shar.Z */
775ba35d 59\f
985b6196 60/* Type of computer arithmetic.
f76b9db2
ILT		 61 Only one of DEC, IBM, IEEE, or UNK should get defined.
62
63 `IEEE', when FLOAT_WORDS_BIG_ENDIAN is non-zero, refers generically
64 to big-endian IEEE floating-point data structure. This definition
65 should work in SFmode `float' type and DFmode `double' type on
66 virtually all big-endian IEEE machines. If LONG_DOUBLE_TYPE_SIZE
67 has been defined to be 96, then IEEE also invokes the particular
68 XFmode (`long double' type) data structure used by the Motorola
69 680x0 series processors.
70
71 `IEEE', when FLOAT_WORDS_BIG_ENDIAN is zero, refers generally to
72 little-endian IEEE machines. In this case, if LONG_DOUBLE_TYPE_SIZE
73 has been defined to be 96, then IEEE also invokes the particular
74 XFmode `long double' data structure used by the Intel 80x86 series
75 processors.
66b6d60b
RS		 76
77 `DEC' refers specifically to the Digital Equipment Corp PDP-11
78 and VAX floating point data structure. This model currently
79 supports no type wider than DFmode.
80
842fbaaa
JW		 81 `IBM' refers specifically to the IBM System/370 and compatible
82 floating point data structure. This model currently supports
83 no type wider than DFmode. The IBM conversions were contributed by
84 frank@atom.ansto.gov.au (Frank Crawford).
85
66b6d60b
RS		 86 If LONG_DOUBLE_TYPE_SIZE = 64 (the default, unless tm.h defines it)
87 then `long double' and `double' are both implemented, but they
88 both mean DFmode. In this case, the software floating-point
89 support available here is activated by writing
90 #define REAL_ARITHMETIC
91 in tm.h.
92
93 The case LONG_DOUBLE_TYPE_SIZE = 128 activates TFmode support
842fbaaa 94 and may deactivate XFmode since `long double' is used to refer
b51ab098
RK		 95 to both modes.
96
97 The macros FLOAT_WORDS_BIG_ENDIAN, HOST_FLOAT_WORDS_BIG_ENDIAN,
98 contributed by Richard Earnshaw <Richard.Earnshaw@cl.cam.ac.uk>,
99 separate the floating point unit's endian-ness from that of
100 the integer addressing. This permits one to define a big-endian
101 FPU on a little-endian machine (e.g., ARM). An extension to
102 BYTES_BIG_ENDIAN may be required for some machines in the future.
103 These optional macros may be defined in tm.h. In real.h, they
104 default to WORDS_BIG_ENDIAN, etc., so there is no need to define
105 them for any normal host or target machine on which the floats
106 and the integers have the same endian-ness. */
107
66b6d60b
RS		 108
109/* The following converts gcc macros into the ones used by this file. */
110
985b6196
RS		 111/* REAL_ARITHMETIC defined means that macros in real.h are
112 defined to call emulator functions. */
113#ifdef REAL_ARITHMETIC
114
115#if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
116/* PDP-11, Pro350, VAX: */
117#define DEC 1
118#else /* it's not VAX */
842fbaaa
JW		 119#if TARGET_FLOAT_FORMAT == IBM_FLOAT_FORMAT
120/* IBM System/370 style */
121#define IBM 1
122#else /* it's also not an IBM */
985b6196 123#if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
f76b9db2 124#define IEEE
985b6196
RS		 125#else /* it's not IEEE either */
126/* UNKnown arithmetic. We don't support this and can't go on. */
127unknown arithmetic type
128#define UNK 1
129#endif /* not IEEE */
842fbaaa 130#endif /* not IBM */
985b6196
RS		 131#endif /* not VAX */
132
133#else
134/* REAL_ARITHMETIC not defined means that the *host's* data
135 structure will be used. It may differ by endian-ness from the
136 target machine's structure and will get its ends swapped
137 accordingly (but not here). Probably only the decimal <-> binary
138 functions in this file will actually be used in this case. */
defb5dab 139
985b6196
RS		 140#if HOST_FLOAT_FORMAT == VAX_FLOAT_FORMAT
141#define DEC 1
142#else /* it's not VAX */
842fbaaa
JW		 143#if HOST_FLOAT_FORMAT == IBM_FLOAT_FORMAT
144/* IBM System/370 style */
145#define IBM 1
146#else /* it's also not an IBM */
985b6196 147#if HOST_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
f76b9db2 148#define IEEE
985b6196
RS		 149#else /* it's not IEEE either */
150unknown arithmetic type
151#define UNK 1
152#endif /* not IEEE */
842fbaaa 153#endif /* not IBM */
985b6196
RS		 154#endif /* not VAX */
155
156#endif /* REAL_ARITHMETIC not defined */
157
66b6d60b
RS		 158/* Define INFINITY for support of infinity.
159 Define NANS for support of Not-a-Number's (NaN's). */
842fbaaa 160#if !defined(DEC) && !defined(IBM)
985b6196 161#define INFINITY
66b6d60b 162#define NANS
985b6196
RS		 163#endif
164
66b6d60b
RS		 165/* Support of NaNs requires support of infinity. */
166#ifdef NANS
167#ifndef INFINITY
168#define INFINITY
169#endif
170#endif
775ba35d 171\f
985b6196
RS		 172/* Find a host integer type that is at least 16 bits wide,
173 and another type at least twice whatever that size is. */
174
175#if HOST_BITS_PER_CHAR >= 16
176#define EMUSHORT char
177#define EMUSHORT_SIZE HOST_BITS_PER_CHAR
178#define EMULONG_SIZE (2 * HOST_BITS_PER_CHAR)
179#else
180#if HOST_BITS_PER_SHORT >= 16
181#define EMUSHORT short
182#define EMUSHORT_SIZE HOST_BITS_PER_SHORT
183#define EMULONG_SIZE (2 * HOST_BITS_PER_SHORT)
184#else
185#if HOST_BITS_PER_INT >= 16
186#define EMUSHORT int
187#define EMUSHORT_SIZE HOST_BITS_PER_INT
188#define EMULONG_SIZE (2 * HOST_BITS_PER_INT)
189#else
190#if HOST_BITS_PER_LONG >= 16
191#define EMUSHORT long
192#define EMUSHORT_SIZE HOST_BITS_PER_LONG
193#define EMULONG_SIZE (2 * HOST_BITS_PER_LONG)
194#else
195/* You will have to modify this program to have a smaller unit size. */
196#define EMU_NON_COMPILE
197#endif
198#endif
199#endif
200#endif
201
202#if HOST_BITS_PER_SHORT >= EMULONG_SIZE
203#define EMULONG short
204#else
205#if HOST_BITS_PER_INT >= EMULONG_SIZE
206#define EMULONG int
207#else
208#if HOST_BITS_PER_LONG >= EMULONG_SIZE
209#define EMULONG long
210#else
211#if HOST_BITS_PER_LONG_LONG >= EMULONG_SIZE
212#define EMULONG long long int
213#else
214/* You will have to modify this program to have a smaller unit size. */
215#define EMU_NON_COMPILE
216#endif
217#endif
218#endif
219#endif
220
221
222/* The host interface doesn't work if no 16-bit size exists. */
223#if EMUSHORT_SIZE != 16
224#define EMU_NON_COMPILE
225#endif
226
227/* OK to continue compilation. */
228#ifndef EMU_NON_COMPILE
229
230/* Construct macros to translate between REAL_VALUE_TYPE and e type.
231 In GET_REAL and PUT_REAL, r and e are pointers.
232 A REAL_VALUE_TYPE is guaranteed to occupy contiguous locations
233 in memory, with no holes. */
234
235#if LONG_DOUBLE_TYPE_SIZE == 96
842fbaaa
JW		 236/* Number of 16 bit words in external e type format */
237#define NE 6
238#define MAXDECEXP 4932
239#define MINDECEXP -4956
4ba46f43
RK		 240#define GET_REAL(r,e) bcopy ((char *) r, (char *) e, 2*NE)
241#define PUT_REAL(e,r) bcopy ((char *) e, (char *) r, 2*NE)
985b6196 242#else /* no XFmode */
842fbaaa
JW		 243#if LONG_DOUBLE_TYPE_SIZE == 128
244#define NE 10
245#define MAXDECEXP 4932
246#define MINDECEXP -4977
4ba46f43
RK		 247#define GET_REAL(r,e) bcopy ((char *) r, (char *) e, 2*NE)
248#define PUT_REAL(e,r) bcopy ((char *) e, (char *) r, 2*NE)
842fbaaa
JW		 249#else
250#define NE 6
251#define MAXDECEXP 4932
252#define MINDECEXP -4956
985b6196
RS		 253#ifdef REAL_ARITHMETIC
254/* Emulator uses target format internally
255 but host stores it in host endian-ness. */
256
f76b9db2
ILT		 257#define GET_REAL(r,e) \
258do { \
259 if (HOST_FLOAT_WORDS_BIG_ENDIAN == FLOAT_WORDS_BIG_ENDIAN) \
260 e53toe ((unsigned EMUSHORT*) (r), (e)); \
261 else \
262 { \
263 unsigned EMUSHORT w[4]; \
264 w[3] = ((EMUSHORT *) r)[0]; \
265 w[2] = ((EMUSHORT *) r)[1]; \
266 w[1] = ((EMUSHORT *) r)[2]; \
267 w[0] = ((EMUSHORT *) r)[3]; \
268 e53toe (w, (e)); \
269 } \
270 } while (0)
271
272#define PUT_REAL(e,r) \
273do { \
274 if (HOST_FLOAT_WORDS_BIG_ENDIAN == FLOAT_WORDS_BIG_ENDIAN) \
275 etoe53 ((e), (unsigned EMUSHORT *) (r)); \
276 else \
277 { \
278 unsigned EMUSHORT w[4]; \
279 etoe53 ((e), w); \
280 *((EMUSHORT *) r) = w[3]; \
281 *((EMUSHORT *) r + 1) = w[2]; \
282 *((EMUSHORT *) r + 2) = w[1]; \
283 *((EMUSHORT *) r + 3) = w[0]; \
284 } \
285 } while (0)
985b6196
RS		 286
287#else /* not REAL_ARITHMETIC */
288
289/* emulator uses host format */
a0353055
RK		 290#define GET_REAL(r,e) e53toe ((unsigned EMUSHORT *) (r), (e))
291#define PUT_REAL(e,r) etoe53 ((e), (unsigned EMUSHORT *) (r))
985b6196
RS		 292
293#endif /* not REAL_ARITHMETIC */
842fbaaa 294#endif /* not TFmode */
985b6196
RS		 295#endif /* no XFmode */
296
842fbaaa
JW		 297
298/* Number of 16 bit words in internal format */
299#define NI (NE+3)
300
301/* Array offset to exponent */
302#define E 1
303
304/* Array offset to high guard word */
305#define M 2
306
307/* Number of bits of precision */
308#define NBITS ((NI-4)*16)
309
310/* Maximum number of decimal digits in ASCII conversion
311 * = NBITS*log10(2)
312 */
313#define NDEC (NBITS*8/27)
314
315/* The exponent of 1.0 */
316#define EXONE (0x3fff)
317
64685ffa 318extern int extra_warnings;
9d1bd99c
MM		 319extern unsigned EMUSHORT ezero[], ehalf[], eone[], etwo[];
320extern unsigned EMUSHORT elog2[], esqrt2[];
a0353055
RK		 321
322static void endian PROTO((unsigned EMUSHORT *, long *,
323 enum machine_mode));
324static void eclear PROTO((unsigned EMUSHORT *));
325static void emov PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
326static void eabs PROTO((unsigned EMUSHORT *));
327static void eneg PROTO((unsigned EMUSHORT *));
328static int eisneg PROTO((unsigned EMUSHORT *));
329static int eisinf PROTO((unsigned EMUSHORT *));
330static int eisnan PROTO((unsigned EMUSHORT *));
331static void einfin PROTO((unsigned EMUSHORT *));
332static void enan PROTO((unsigned EMUSHORT *, int));
333static void emovi PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
334static void emovo PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
335static void ecleaz PROTO((unsigned EMUSHORT *));
336static void ecleazs PROTO((unsigned EMUSHORT *));
337static void emovz PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
338static void einan PROTO((unsigned EMUSHORT *));
339static int eiisnan PROTO((unsigned EMUSHORT *));
340static int eiisneg PROTO((unsigned EMUSHORT *));
341static void eiinfin PROTO((unsigned EMUSHORT *));
342static int eiisinf PROTO((unsigned EMUSHORT *));
343static int ecmpm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
344static void eshdn1 PROTO((unsigned EMUSHORT *));
345static void eshup1 PROTO((unsigned EMUSHORT *));
346static void eshdn8 PROTO((unsigned EMUSHORT *));
347static void eshup8 PROTO((unsigned EMUSHORT *));
348static void eshup6 PROTO((unsigned EMUSHORT *));
349static void eshdn6 PROTO((unsigned EMUSHORT *));
350static void eaddm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));\f
351static void esubm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
c92d992a 352static void m16m PROTO((unsigned int, unsigned short *,
a0353055
RK		 353 unsigned short *));
354static int edivm PROTO((unsigned short *, unsigned short *));
355static int emulm PROTO((unsigned short *, unsigned short *));
356static void emdnorm PROTO((unsigned EMUSHORT *, int, int, EMULONG, int));
357static void esub PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
358 unsigned EMUSHORT *));
359static void eadd PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
360 unsigned EMUSHORT *));
361static void eadd1 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
362 unsigned EMUSHORT *));
363static void ediv PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
364 unsigned EMUSHORT *));
365static void emul PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
366 unsigned EMUSHORT *));
367static void e53toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
368static void e64toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
369static void e113toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
370static void e24toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
371static void etoe113 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
372static void toe113 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
373static void etoe64 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
374static void toe64 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
375static void etoe53 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
376static void toe53 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
377static void etoe24 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
378static void toe24 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
379static int ecmp PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
380static void eround PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
381static void ltoe PROTO((HOST_WIDE_INT *, unsigned EMUSHORT *));
382static void ultoe PROTO((unsigned HOST_WIDE_INT *, unsigned EMUSHORT *));
383static void eifrac PROTO((unsigned EMUSHORT *, HOST_WIDE_INT *,
384 unsigned EMUSHORT *));
385static void euifrac PROTO((unsigned EMUSHORT *, unsigned HOST_WIDE_INT *,
386 unsigned EMUSHORT *));
387static int eshift PROTO((unsigned EMUSHORT *, int));
388static int enormlz PROTO((unsigned EMUSHORT *));
389static void e24toasc PROTO((unsigned EMUSHORT *, char *, int));
390static void e53toasc PROTO((unsigned EMUSHORT *, char *, int));
391static void e64toasc PROTO((unsigned EMUSHORT *, char *, int));
392static void e113toasc PROTO((unsigned EMUSHORT *, char *, int));
393static void etoasc PROTO((unsigned EMUSHORT *, char *, int));
394static void asctoe24 PROTO((char *, unsigned EMUSHORT *));
395static void asctoe53 PROTO((char *, unsigned EMUSHORT *));
396static void asctoe64 PROTO((char *, unsigned EMUSHORT *));
397static void asctoe113 PROTO((char *, unsigned EMUSHORT *));
398static void asctoe PROTO((char *, unsigned EMUSHORT *));
399static void asctoeg PROTO((char *, unsigned EMUSHORT *, int));
400static void efloor PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
401static void efrexp PROTO((unsigned EMUSHORT *, int *,
402 unsigned EMUSHORT *));
403static void eldexp PROTO((unsigned EMUSHORT *, int, unsigned EMUSHORT *));
404static void eremain PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
405 unsigned EMUSHORT *));
406static void eiremain PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
407static void mtherr PROTO((char *, int));
408static void dectoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
409static void etodec PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
410static void todec PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
411static void ibmtoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
412 enum machine_mode));
413static void etoibm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
414 enum machine_mode));
415static void toibm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
416 enum machine_mode));
417static void make_nan PROTO((unsigned EMUSHORT *, int, enum machine_mode));
418static void uditoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
419static void ditoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
420static void etoudi PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
421static void etodi PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
422static void esqrt PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
775ba35d 423\f
b51ab098
RK		 424/* Copy 32-bit numbers obtained from array containing 16-bit numbers,
425 swapping ends if required, into output array of longs. The
426 result is normally passed to fprintf by the ASM_OUTPUT_ macros. */
a0353055
RK		 427
428static void
985b6196
RS		 429endian (e, x, mode)
430 unsigned EMUSHORT e[];
431 long x[];
432 enum machine_mode mode;
433{
434 unsigned long th, t;
435
f76b9db2 436 if (FLOAT_WORDS_BIG_ENDIAN)
985b6196 437 {
f76b9db2
ILT		 438 switch (mode)
439 {
985b6196 440
f76b9db2
ILT		 441 case TFmode:
442 /* Swap halfwords in the fourth long. */
443 th = (unsigned long) e[6] & 0xffff;
444 t = (unsigned long) e[7] & 0xffff;
445 t |= th << 16;
446 x[3] = (long) t;
447
448 case XFmode:
449
450 /* Swap halfwords in the third long. */
451 th = (unsigned long) e[4] & 0xffff;
452 t = (unsigned long) e[5] & 0xffff;
453 t |= th << 16;
454 x[2] = (long) t;
455 /* fall into the double case */
456
457 case DFmode:
458
459 /* swap halfwords in the second word */
460 th = (unsigned long) e[2] & 0xffff;
461 t = (unsigned long) e[3] & 0xffff;
462 t |= th << 16;
463 x[1] = (long) t;
464 /* fall into the float case */
465
466 case HFmode:
467 case SFmode:
468
469 /* swap halfwords in the first word */
470 th = (unsigned long) e[0] & 0xffff;
471 t = (unsigned long) e[1] & 0xffff;
472 t |= th << 16;
473 x[0] = t;
474 break;
985b6196 475
f76b9db2
ILT		 476 default:
477 abort ();
478 }
985b6196 479 }
f76b9db2 480 else
985b6196 481 {
f76b9db2 482 /* Pack the output array without swapping. */
985b6196 483
f76b9db2
ILT		 484 switch (mode)
485 {
985b6196 486
f76b9db2
ILT		 487 case TFmode:
488
489 /* Pack the fourth long. */
490 th = (unsigned long) e[7] & 0xffff;
491 t = (unsigned long) e[6] & 0xffff;
492 t |= th << 16;
493 x[3] = (long) t;
494
495 case XFmode:
496
497 /* Pack the third long.
498 Each element of the input REAL_VALUE_TYPE array has 16 useful bits
499 in it. */
500 th = (unsigned long) e[5] & 0xffff;
501 t = (unsigned long) e[4] & 0xffff;
502 t |= th << 16;
503 x[2] = (long) t;
504 /* fall into the double case */
505
506 case DFmode:
507
508 /* pack the second long */
509 th = (unsigned long) e[3] & 0xffff;
510 t = (unsigned long) e[2] & 0xffff;
511 t |= th << 16;
512 x[1] = (long) t;
513 /* fall into the float case */
514
515 case HFmode:
516 case SFmode:
517
518 /* pack the first long */
519 th = (unsigned long) e[1] & 0xffff;
520 t = (unsigned long) e[0] & 0xffff;
521 t |= th << 16;
522 x[0] = t;
523 break;
985b6196 524
f76b9db2
ILT		 525 default:
526 abort ();
527 }
985b6196 528 }
985b6196
RS		 529}
530
531
defb5dab 532/* This is the implementation of the REAL_ARITHMETIC macro. */
a0353055 533
985b6196
RS		 534void
535earith (value, icode, r1, r2)
536 REAL_VALUE_TYPE *value;
537 int icode;
538 REAL_VALUE_TYPE *r1;
539 REAL_VALUE_TYPE *r2;
540{
541 unsigned EMUSHORT d1[NE], d2[NE], v[NE];
542 enum tree_code code;
543
544 GET_REAL (r1, d1);
545 GET_REAL (r2, d2);
66b6d60b
RS		 546#ifdef NANS
547/* Return NaN input back to the caller. */
548 if (eisnan (d1))
549 {
550 PUT_REAL (d1, value);
551 return;
552 }
553 if (eisnan (d2))
554 {
555 PUT_REAL (d2, value);
556 return;
557 }
558#endif
985b6196
RS		 559 code = (enum tree_code) icode;
560 switch (code)
561 {
562 case PLUS_EXPR:
563 eadd (d2, d1, v);
564 break;
565
566 case MINUS_EXPR:
567 esub (d2, d1, v); /* d1 - d2 */
568 break;
569
570 case MULT_EXPR:
571 emul (d2, d1, v);
572 break;
573
574 case RDIV_EXPR:
575#ifndef REAL_INFINITY
576 if (ecmp (d2, ezero) == 0)
66b6d60b
RS		 577 {
578#ifdef NANS
29e11dab 579 enan (v, eisneg (d1) ^ eisneg (d2));
66b6d60b
RS		 580 break;
581#else
985b6196 582 abort ();
66b6d60b
RS		 583#endif
584 }
985b6196
RS		 585#endif
586 ediv (d2, d1, v); /* d1/d2 */
587 break;
588
589 case MIN_EXPR: /* min (d1,d2) */
590 if (ecmp (d1, d2) < 0)
591 emov (d1, v);
592 else
593 emov (d2, v);
594 break;
595
596 case MAX_EXPR: /* max (d1,d2) */
597 if (ecmp (d1, d2) > 0)
598 emov (d1, v);
599 else
600 emov (d2, v);
601 break;
602 default:
603 emov (ezero, v);
604 break;
605 }
606PUT_REAL (v, value);
607}
608
609
defb5dab
RK		 610/* Truncate REAL_VALUE_TYPE toward zero to signed HOST_WIDE_INT.
611 implements REAL_VALUE_RNDZINT (x) (etrunci (x)). */
612
985b6196
RS		 613REAL_VALUE_TYPE
614etrunci (x)
615 REAL_VALUE_TYPE x;
616{
617 unsigned EMUSHORT f[NE], g[NE];
618 REAL_VALUE_TYPE r;
b51ab098 619 HOST_WIDE_INT l;
985b6196
RS		 620
621 GET_REAL (&x, g);
66b6d60b
RS		 622#ifdef NANS
623 if (eisnan (g))
624 return (x);
625#endif
985b6196
RS		 626 eifrac (g, &l, f);
627 ltoe (&l, g);
628 PUT_REAL (g, &r);
629 return (r);
630}
631
632
defb5dab
RK		 633/* Truncate REAL_VALUE_TYPE toward zero to unsigned HOST_WIDE_INT;
634 implements REAL_VALUE_UNSIGNED_RNDZINT (x) (etruncui (x)). */
635
985b6196
RS		 636REAL_VALUE_TYPE
637etruncui (x)
638 REAL_VALUE_TYPE x;
639{
640 unsigned EMUSHORT f[NE], g[NE];
641 REAL_VALUE_TYPE r;
b51ab098 642 unsigned HOST_WIDE_INT l;
985b6196
RS		 643
644 GET_REAL (&x, g);
66b6d60b
RS		 645#ifdef NANS
646 if (eisnan (g))
647 return (x);
648#endif
985b6196
RS		 649 euifrac (g, &l, f);
650 ultoe (&l, g);
651 PUT_REAL (g, &r);
652 return (r);
653}
654
655
defb5dab
RK		 656/* This is the REAL_VALUE_ATOF function. It converts a decimal string to
657 binary, rounding off as indicated by the machine_mode argument. Then it
658 promotes the rounded value to REAL_VALUE_TYPE. */
659
985b6196
RS		 660REAL_VALUE_TYPE
661ereal_atof (s, t)
662 char *s;
663 enum machine_mode t;
664{
665 unsigned EMUSHORT tem[NE], e[NE];
666 REAL_VALUE_TYPE r;
667
668 switch (t)
669 {
bfbc6416 670 case HFmode:
985b6196
RS		 671 case SFmode:
672 asctoe24 (s, tem);
673 e24toe (tem, e);
674 break;
675 case DFmode:
676 asctoe53 (s, tem);
677 e53toe (tem, e);
678 break;
679 case XFmode:
680 asctoe64 (s, tem);
681 e64toe (tem, e);
682 break;
842fbaaa
JW		 683 case TFmode:
684 asctoe113 (s, tem);
685 e113toe (tem, e);
686 break;
985b6196
RS		 687 default:
688 asctoe (s, e);
689 }
690 PUT_REAL (e, &r);
691 return (r);
692}
693
694
defb5dab
RK		 695/* Expansion of REAL_NEGATE. */
696
985b6196
RS		 697REAL_VALUE_TYPE
698ereal_negate (x)
699 REAL_VALUE_TYPE x;
700{
701 unsigned EMUSHORT e[NE];
702 REAL_VALUE_TYPE r;
703
704 GET_REAL (&x, e);
705 eneg (e);
706 PUT_REAL (e, &r);
707 return (r);
708}
709
710
defb5dab
RK		 711/* Round real toward zero to HOST_WIDE_INT;
712 implements REAL_VALUE_FIX (x). */
713
b51ab098 714HOST_WIDE_INT
842fbaaa 715efixi (x)
985b6196
RS		 716 REAL_VALUE_TYPE x;
717{
718 unsigned EMUSHORT f[NE], g[NE];
b51ab098 719 HOST_WIDE_INT l;
985b6196
RS		 720
721 GET_REAL (&x, f);
66b6d60b
RS		 722#ifdef NANS
723 if (eisnan (f))
724 {
725 warning ("conversion from NaN to int");
726 return (-1);
727 }
728#endif
842fbaaa
JW		 729 eifrac (f, &l, g);
730 return l;
985b6196
RS		 731}
732
842fbaaa 733/* Round real toward zero to unsigned HOST_WIDE_INT
defb5dab
RK		 734 implements REAL_VALUE_UNSIGNED_FIX (x).
735 Negative input returns zero. */
736
b51ab098 737unsigned HOST_WIDE_INT
842fbaaa 738efixui (x)
985b6196
RS		 739 REAL_VALUE_TYPE x;
740{
741 unsigned EMUSHORT f[NE], g[NE];
b51ab098 742 unsigned HOST_WIDE_INT l;
985b6196
RS		 743
744 GET_REAL (&x, f);
66b6d60b
RS		 745#ifdef NANS
746 if (eisnan (f))
747 {
748 warning ("conversion from NaN to unsigned int");
749 return (-1);
750 }
751#endif
842fbaaa
JW		 752 euifrac (f, &l, g);
753 return l;
985b6196
RS		 754}
755
756
defb5dab
RK		 757/* REAL_VALUE_FROM_INT macro. */
758
985b6196
RS		 759void
760ereal_from_int (d, i, j)
761 REAL_VALUE_TYPE *d;
b51ab098 762 HOST_WIDE_INT i, j;
985b6196
RS		 763{
764 unsigned EMUSHORT df[NE], dg[NE];
b51ab098 765 HOST_WIDE_INT low, high;
985b6196
RS		 766 int sign;
767
768 sign = 0;
769 low = i;
770 if ((high = j) < 0)
771 {
772 sign = 1;
773 /* complement and add 1 */
774 high = ~high;
775 if (low)
776 low = -low;
777 else
778 high += 1;
779 }
b51ab098 780 eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
60e61165 781 ultoe ((unsigned HOST_WIDE_INT *) &high, dg);
985b6196 782 emul (dg, df, dg);
60e61165 783 ultoe ((unsigned HOST_WIDE_INT *) &low, df);
985b6196
RS		 784 eadd (df, dg, dg);
785 if (sign)
786 eneg (dg);
787 PUT_REAL (dg, d);
788}
789
790
defb5dab 791/* REAL_VALUE_FROM_UNSIGNED_INT macro. */
a0353055 792
985b6196
RS		 793void
794ereal_from_uint (d, i, j)
795 REAL_VALUE_TYPE *d;
b51ab098 796 unsigned HOST_WIDE_INT i, j;
985b6196
RS		 797{
798 unsigned EMUSHORT df[NE], dg[NE];
b51ab098 799 unsigned HOST_WIDE_INT low, high;
985b6196
RS		 800
801 low = i;
802 high = j;
b51ab098 803 eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
985b6196
RS		 804 ultoe (&high, dg);
805 emul (dg, df, dg);
806 ultoe (&low, df);
807 eadd (df, dg, dg);
808 PUT_REAL (dg, d);
809}
810
811
defb5dab
RK		 812/* REAL_VALUE_TO_INT macro. */
813
985b6196
RS		 814void
815ereal_to_int (low, high, rr)
b51ab098 816 HOST_WIDE_INT *low, *high;
985b6196
RS		 817 REAL_VALUE_TYPE rr;
818{
819 unsigned EMUSHORT d[NE], df[NE], dg[NE], dh[NE];
820 int s;
821
822 GET_REAL (&rr, d);
66b6d60b 823#ifdef NANS
970491df 824 if (eisnan (d))
66b6d60b
RS		 825 {
826 warning ("conversion from NaN to int");
827 *low = -1;
828 *high = -1;
829 return;
830 }
831#endif
985b6196
RS		 832 /* convert positive value */
833 s = 0;
834 if (eisneg (d))
835 {
836 eneg (d);
837 s = 1;
838 }
b51ab098 839 eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
985b6196 840 ediv (df, d, dg); /* dg = d / 2^32 is the high word */
60e61165 841 euifrac (dg, (unsigned HOST_WIDE_INT *) high, dh);
985b6196 842 emul (df, dh, dg); /* fractional part is the low word */
60e61165 843 euifrac (dg, (unsigned HOST_WIDE_INT *)low, dh);
985b6196
RS		 844 if (s)
845 {
846 /* complement and add 1 */
847 *high = ~(*high);
848 if (*low)
849 *low = -(*low);
850 else
851 *high += 1;
852 }
853}
854
855
defb5dab
RK		 856/* REAL_VALUE_LDEXP macro. */
857
985b6196
RS		 858REAL_VALUE_TYPE
859ereal_ldexp (x, n)
860 REAL_VALUE_TYPE x;
861 int n;
862{
863 unsigned EMUSHORT e[NE], y[NE];
864 REAL_VALUE_TYPE r;
865
866 GET_REAL (&x, e);
66b6d60b
RS		 867#ifdef NANS
868 if (eisnan (e))
869 return (x);
870#endif
985b6196
RS		 871 eldexp (e, n, y);
872 PUT_REAL (y, &r);
873 return (r);
874}
875
876/* These routines are conditionally compiled because functions
defb5dab
RK		 877 of the same names may be defined in fold-const.c. */
878
985b6196
RS		 879#ifdef REAL_ARITHMETIC
880
881/* Check for infinity in a REAL_VALUE_TYPE. */
defb5dab 882
985b6196
RS		 883int
884target_isinf (x)
885 REAL_VALUE_TYPE x;
886{
887 unsigned EMUSHORT e[NE];
888
889#ifdef INFINITY
890 GET_REAL (&x, e);
891 return (eisinf (e));
892#else
893 return 0;
894#endif
895}
896
897
66b6d60b 898/* Check whether a REAL_VALUE_TYPE item is a NaN. */
985b6196
RS		 899
900int
901target_isnan (x)
902 REAL_VALUE_TYPE x;
903{
9d72da33
RS		 904 unsigned EMUSHORT e[NE];
905
66b6d60b 906#ifdef NANS
9d72da33
RS		 907 GET_REAL (&x, e);
908 return (eisnan (e));
66b6d60b 909#else
985b6196 910 return (0);
66b6d60b 911#endif
985b6196
RS		 912}
913
914
66b6d60b 915/* Check for a negative REAL_VALUE_TYPE number.
defb5dab 916 This just checks the sign bit, so that -0 counts as negative. */
985b6196
RS		 917
918int
919target_negative (x)
920 REAL_VALUE_TYPE x;
921{
281bb5e4 922 return ereal_isneg (x);
985b6196
RS		 923}
924
925/* Expansion of REAL_VALUE_TRUNCATE.
defb5dab
RK		 926 The result is in floating point, rounded to nearest or even. */
927
985b6196
RS		 928REAL_VALUE_TYPE
929real_value_truncate (mode, arg)
930 enum machine_mode mode;
931 REAL_VALUE_TYPE arg;
932{
933 unsigned EMUSHORT e[NE], t[NE];
934 REAL_VALUE_TYPE r;
935
936 GET_REAL (&arg, e);
66b6d60b
RS		 937#ifdef NANS
938 if (eisnan (e))
939 return (arg);
940#endif
985b6196
RS		 941 eclear (t);
942 switch (mode)
943 {
842fbaaa
JW		 944 case TFmode:
945 etoe113 (e, t);
946 e113toe (t, t);
947 break;
948
985b6196
RS		 949 case XFmode:
950 etoe64 (e, t);
951 e64toe (t, t);
952 break;
953
954 case DFmode:
955 etoe53 (e, t);
956 e53toe (t, t);
957 break;
958
bfbc6416 959 case HFmode:
985b6196
RS		 960 case SFmode:
961 etoe24 (e, t);
962 e24toe (t, t);
963 break;
964
965 case SImode:
f8ece317 966 r = etrunci (arg);
985b6196
RS		 967 return (r);
968
0de689b7
RK		 969 /* If an unsupported type was requested, presume that
970 the machine files know something useful to do with
971 the unmodified value. */
defb5dab 972
985b6196 973 default:
0de689b7 974 return (arg);
985b6196
RS		 975 }
976 PUT_REAL (t, &r);
977 return (r);
978}
979
980#endif /* REAL_ARITHMETIC defined */
981
775ba35d
RS		 982/* Used for debugging--print the value of R in human-readable format
983 on stderr. */
984
985void
986debug_real (r)
987 REAL_VALUE_TYPE r;
988{
989 char dstr[30];
990
991 REAL_VALUE_TO_DECIMAL (r, "%.20g", dstr);
992 fprintf (stderr, "%s", dstr);
993}
994
995\f
985b6196
RS		 996/* Target values are arrays of host longs. A long is guaranteed
997 to be at least 32 bits wide. */
842fbaaa
JW		 998
999/* 128-bit long double */
defb5dab 1000
842fbaaa
JW		 1001void
1002etartdouble (r, l)
1003 REAL_VALUE_TYPE r;
1004 long l[];
1005{
1006 unsigned EMUSHORT e[NE];
1007
1008 GET_REAL (&r, e);
1009 etoe113 (e, e);
1010 endian (e, l, TFmode);
1011}
1012
1013/* 80-bit long double */
defb5dab 1014
985b6196
RS		 1015void
1016etarldouble (r, l)
1017 REAL_VALUE_TYPE r;
1018 long l[];
1019{
1020 unsigned EMUSHORT e[NE];
1021
1022 GET_REAL (&r, e);
1023 etoe64 (e, e);
1024 endian (e, l, XFmode);
1025}
1026
1027void
1028etardouble (r, l)
1029 REAL_VALUE_TYPE r;
1030 long l[];
1031{
1032 unsigned EMUSHORT e[NE];
1033
1034 GET_REAL (&r, e);
1035 etoe53 (e, e);
1036 endian (e, l, DFmode);
1037}
1038
1039long
1040etarsingle (r)
1041 REAL_VALUE_TYPE r;
1042{
1043 unsigned EMUSHORT e[NE];
60e61165 1044 long l;
985b6196
RS		 1045
1046 GET_REAL (&r, e);
1047 etoe24 (e, e);
1048 endian (e, &l, SFmode);
1049 return ((long) l);
1050}
1051
1052void
1053ereal_to_decimal (x, s)
1054 REAL_VALUE_TYPE x;
1055 char *s;
1056{
1057 unsigned EMUSHORT e[NE];
1058
1059 GET_REAL (&x, e);
1060 etoasc (e, s, 20);
1061}
1062
1063int
1064ereal_cmp (x, y)
1065 REAL_VALUE_TYPE x, y;
1066{
1067 unsigned EMUSHORT ex[NE], ey[NE];
1068
1069 GET_REAL (&x, ex);
1070 GET_REAL (&y, ey);
1071 return (ecmp (ex, ey));
1072}
1073
1074int
1075ereal_isneg (x)
1076 REAL_VALUE_TYPE x;
1077{
1078 unsigned EMUSHORT ex[NE];
1079
1080 GET_REAL (&x, ex);
1081 return (eisneg (ex));
1082}
1083
1084/* End of REAL_ARITHMETIC interface */
775ba35d 1085\f
defb5dab
RK		 1086/*
1087 Extended precision IEEE binary floating point arithmetic routines
1088
1089 Numbers are stored in C language as arrays of 16-bit unsigned
1090 short integers. The arguments of the routines are pointers to
1091 the arrays.
1092
1093 External e type data structure, simulates Intel 8087 chip
1094 temporary real format but possibly with a larger significand:
1095
1096 NE-1 significand words (least significant word first,
1097 most significant bit is normally set)
1098 exponent (value = EXONE for 1.0,
1099 top bit is the sign)
1100
1101
1102 Internal data structure of a number (a "word" is 16 bits):
1103
1104 ei[0] sign word (0 for positive, 0xffff for negative)
1105 ei[1] biased exponent (value = EXONE for the number 1.0)
1106 ei[2] high guard word (always zero after normalization)
1107 ei[3]
1108 to ei[NI-2] significand (NI-4 significand words,
1109 most significant word first,
1110 most significant bit is set)
1111 ei[NI-1] low guard word (0x8000 bit is rounding place)
1112
1113
1114
1115 Routines for external format numbers
1116
1117 asctoe (string, e) ASCII string to extended double e type
1118 asctoe64 (string, &d) ASCII string to long double
1119 asctoe53 (string, &d) ASCII string to double
1120 asctoe24 (string, &f) ASCII string to single
1121 asctoeg (string, e, prec) ASCII string to specified precision
1122 e24toe (&f, e) IEEE single precision to e type
1123 e53toe (&d, e) IEEE double precision to e type
1124 e64toe (&d, e) IEEE long double precision to e type
1125 e113toe (&d, e) 128-bit long double precision to e type
1126 eabs (e) absolute value
1127 eadd (a, b, c) c = b + a
1128 eclear (e) e = 0
1129 ecmp (a, b) Returns 1 if a > b, 0 if a == b,
1130 -1 if a < b, -2 if either a or b is a NaN.
1131 ediv (a, b, c) c = b / a
1132 efloor (a, b) truncate to integer, toward -infinity
1133 efrexp (a, exp, s) extract exponent and significand
1134 eifrac (e, &l, frac) e to HOST_WIDE_INT and e type fraction
1135 euifrac (e, &l, frac) e to unsigned HOST_WIDE_INT and e type fraction
1136 einfin (e) set e to infinity, leaving its sign alone
1137 eldexp (a, n, b) multiply by 2**n
1138 emov (a, b) b = a
1139 emul (a, b, c) c = b * a
1140 eneg (e) e = -e
1141 eround (a, b) b = nearest integer value to a
1142 esub (a, b, c) c = b - a
1143 e24toasc (&f, str, n) single to ASCII string, n digits after decimal
1144 e53toasc (&d, str, n) double to ASCII string, n digits after decimal
1145 e64toasc (&d, str, n) 80-bit long double to ASCII string
1146 e113toasc (&d, str, n) 128-bit long double to ASCII string
1147 etoasc (e, str, n) e to ASCII string, n digits after decimal
1148 etoe24 (e, &f) convert e type to IEEE single precision
1149 etoe53 (e, &d) convert e type to IEEE double precision
1150 etoe64 (e, &d) convert e type to IEEE long double precision
1151 ltoe (&l, e) HOST_WIDE_INT to e type
1152 ultoe (&l, e) unsigned HOST_WIDE_INT to e type
1153 eisneg (e) 1 if sign bit of e != 0, else 0
1154 eisinf (e) 1 if e has maximum exponent (non-IEEE)
1155 or is infinite (IEEE)
1156 eisnan (e) 1 if e is a NaN
1157
1158
1159 Routines for internal format numbers
1160
1161 eaddm (ai, bi) add significands, bi = bi + ai
1162 ecleaz (ei) ei = 0
1163 ecleazs (ei) set ei = 0 but leave its sign alone
1164 ecmpm (ai, bi) compare significands, return 1, 0, or -1
1165 edivm (ai, bi) divide significands, bi = bi / ai
1166 emdnorm (ai,l,s,exp) normalize and round off
1167 emovi (a, ai) convert external a to internal ai
1168 emovo (ai, a) convert internal ai to external a
1169 emovz (ai, bi) bi = ai, low guard word of bi = 0
1170 emulm (ai, bi) multiply significands, bi = bi * ai
1171 enormlz (ei) left-justify the significand
1172 eshdn1 (ai) shift significand and guards down 1 bit
1173 eshdn8 (ai) shift down 8 bits
1174 eshdn6 (ai) shift down 16 bits
1175 eshift (ai, n) shift ai n bits up (or down if n < 0)
1176 eshup1 (ai) shift significand and guards up 1 bit
1177 eshup8 (ai) shift up 8 bits
1178 eshup6 (ai) shift up 16 bits
1179 esubm (ai, bi) subtract significands, bi = bi - ai
1180 eiisinf (ai) 1 if infinite
1181 eiisnan (ai) 1 if a NaN
1182 eiisneg (ai) 1 if sign bit of ai != 0, else 0
1183 einan (ai) set ai = NaN
1184 eiinfin (ai) set ai = infinity
1185
1186 The result is always normalized and rounded to NI-4 word precision
1187 after each arithmetic operation.
1188
1189 Exception flags are NOT fully supported.
1190
1191 Signaling NaN's are NOT supported; they are treated the same
1192 as quiet NaN's.
1193
1194 Define INFINITY for support of infinity; otherwise a
1195 saturation arithmetic is implemented.
1196
1197 Define NANS for support of Not-a-Number items; otherwise the
1198 arithmetic will never produce a NaN output, and might be confused
1199 by a NaN input.
1200 If NaN's are supported, the output of `ecmp (a,b)' is -2 if
1201 either a or b is a NaN. This means asking `if (ecmp (a,b) < 0)'
1202 may not be legitimate. Use `if (ecmp (a,b) == -1)' for `less than'
1203 if in doubt.
1204
1205 Denormals are always supported here where appropriate (e.g., not
1206 for conversion to DEC numbers). */
1207
1208/* Definitions for error codes that are passed to the common error handling
1209 routine mtherr.
1210
1211 For Digital Equipment PDP-11 and VAX computers, certain
1212 IBM systems, and others that use numbers with a 56-bit
1213 significand, the symbol DEC should be defined. In this
1214 mode, most floating point constants are given as arrays
1215 of octal integers to eliminate decimal to binary conversion
1216 errors that might be introduced by the compiler.
1217
1218 For computers, such as IBM PC, that follow the IEEE
1219 Standard for Binary Floating Point Arithmetic (ANSI/IEEE
1220 Std 754-1985), the symbol IBMPC or MIEEE should be defined.
1221 These numbers have 53-bit significands. In this mode, constants
1222 are provided as arrays of hexadecimal 16 bit integers.
f76b9db2
ILT		 1223 [This has been changed to instead check the preprocessor macros IEEE
1224 and FLOAT_WORDS_BIG_ENDIAN].
defb5dab
RK		 1225
1226 To accommodate other types of computer arithmetic, all
1227 constants are also provided in a normal decimal radix
1228 which one can hope are correctly converted to a suitable
1229 format by the available C language compiler. To invoke
1230 this mode, the symbol UNK is defined.
1231
1232 An important difference among these modes is a predefined
1233 set of machine arithmetic constants for each. The numbers
1234 MACHEP (the machine roundoff error), MAXNUM (largest number
1235 represented), and several other parameters are preset by
1236 the configuration symbol. Check the file const.c to
1237 ensure that these values are correct for your computer.
1238
1239 For ANSI C compatibility, define ANSIC equal to 1. Currently
1240 this affects only the atan2 function and others that use it. */
985b6196 1241
e8650b8f 1242/* Constant definitions for math error conditions. */
985b6196
RS		 1243
1244#define DOMAIN 1 /* argument domain error */
1245#define SING 2 /* argument singularity */
1246#define OVERFLOW 3 /* overflow range error */
1247#define UNDERFLOW 4 /* underflow range error */
1248#define TLOSS 5 /* total loss of precision */
1249#define PLOSS 6 /* partial loss of precision */
66b6d60b 1250#define INVALID 7 /* NaN-producing operation */
985b6196 1251
985b6196
RS		 1252/* e type constants used by high precision check routines */
1253
842fbaaa 1254#if LONG_DOUBLE_TYPE_SIZE == 128
985b6196
RS		 1255/* 0.0 */
1256unsigned EMUSHORT ezero[NE] =
842fbaaa
JW		 1257 {0x0000, 0x0000, 0x0000, 0x0000,
1258 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,};
985b6196
RS		 1259extern unsigned EMUSHORT ezero[];
1260
1261/* 5.0E-1 */
1262unsigned EMUSHORT ehalf[NE] =
842fbaaa
JW		 1263 {0x0000, 0x0000, 0x0000, 0x0000,
1264 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3ffe,};
985b6196
RS		 1265extern unsigned EMUSHORT ehalf[];
1266
1267/* 1.0E0 */
1268unsigned EMUSHORT eone[NE] =
842fbaaa
JW		 1269 {0x0000, 0x0000, 0x0000, 0x0000,
1270 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,};
985b6196
RS		 1271extern unsigned EMUSHORT eone[];
1272
1273/* 2.0E0 */
1274unsigned EMUSHORT etwo[NE] =
842fbaaa
JW		 1275 {0x0000, 0x0000, 0x0000, 0x0000,
1276 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4000,};
985b6196
RS		 1277extern unsigned EMUSHORT etwo[];
1278
1279/* 3.2E1 */
1280unsigned EMUSHORT e32[NE] =
842fbaaa
JW		 1281 {0x0000, 0x0000, 0x0000, 0x0000,
1282 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4004,};
985b6196
RS		 1283extern unsigned EMUSHORT e32[];
1284
1285/* 6.93147180559945309417232121458176568075500134360255E-1 */
1286unsigned EMUSHORT elog2[NE] =
842fbaaa
JW		 1287 {0x40f3, 0xf6af, 0x03f2, 0xb398,
1288 0xc9e3, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
985b6196
RS		 1289extern unsigned EMUSHORT elog2[];
1290
1291/* 1.41421356237309504880168872420969807856967187537695E0 */
1292unsigned EMUSHORT esqrt2[NE] =
842fbaaa
JW		 1293 {0x1d6f, 0xbe9f, 0x754a, 0x89b3,
1294 0x597d, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
985b6196
RS		 1295extern unsigned EMUSHORT esqrt2[];
1296
985b6196
RS		 1297/* 3.14159265358979323846264338327950288419716939937511E0 */
1298unsigned EMUSHORT epi[NE] =
842fbaaa 1299 {0x2902, 0x1cd1, 0x80dc, 0x628b,
985b6196
RS		 1300 0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1301extern unsigned EMUSHORT epi[];
1302
842fbaaa
JW		 1303#else
1304/* LONG_DOUBLE_TYPE_SIZE is other than 128 */
1305unsigned EMUSHORT ezero[NE] =
1306 {0, 0000000, 0000000, 0000000, 0000000, 0000000,};
1307unsigned EMUSHORT ehalf[NE] =
1308 {0, 0000000, 0000000, 0000000, 0100000, 0x3ffe,};
1309unsigned EMUSHORT eone[NE] =
1310 {0, 0000000, 0000000, 0000000, 0100000, 0x3fff,};
1311unsigned EMUSHORT etwo[NE] =
1312 {0, 0000000, 0000000, 0000000, 0100000, 0040000,};
1313unsigned EMUSHORT e32[NE] =
1314 {0, 0000000, 0000000, 0000000, 0100000, 0040004,};
1315unsigned EMUSHORT elog2[NE] =
1316 {0xc9e4, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1317unsigned EMUSHORT esqrt2[NE] =
1318 {0x597e, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1319unsigned EMUSHORT epi[NE] =
1320 {0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1321#endif
985b6196
RS		 1322
1323
1324
1325/* Control register for rounding precision.
defb5dab
RK		 1326 This can be set to 113 (if NE=10), 80 (if NE=6), 64, 56, 53, or 24 bits. */
1327
985b6196
RS		 1328int rndprc = NBITS;
1329extern int rndprc;
1330
defb5dab 1331/* Clear out entire external format number. */
985b6196 1332
a0353055 1333static void
985b6196
RS		 1334eclear (x)
1335 register unsigned EMUSHORT *x;
1336{
1337 register int i;
1338
1339 for (i = 0; i < NE; i++)
1340 *x++ = 0;
1341}
1342
1343
1344
defb5dab 1345/* Move external format number from a to b. */
985b6196 1346
a0353055 1347static void
985b6196
RS		 1348emov (a, b)
1349 register unsigned EMUSHORT *a, *b;
1350{
1351 register int i;
1352
1353 for (i = 0; i < NE; i++)
1354 *b++ = *a++;
1355}
1356
1357
defb5dab 1358/* Absolute value of external format number. */
985b6196 1359
a0353055 1360static void
985b6196 1361eabs (x)
a0353055 1362 unsigned EMUSHORT x[];
985b6196 1363{
a0353055
RK		 1364 /* sign is top bit of last word of external format */
1365 x[NE - 1] &= 0x7fff;
985b6196
RS		 1366}
1367
defb5dab 1368/* Negate external format number. */
985b6196 1369
a0353055 1370static void
985b6196
RS		 1371eneg (x)
1372 unsigned EMUSHORT x[];
1373{
1374
1375 x[NE - 1] ^= 0x8000; /* Toggle the sign bit */
1376}
1377
1378
1379
defb5dab
RK		 1380/* Return 1 if sign bit of external format number is nonzero, else zero. */
1381
a0353055 1382static int
985b6196
RS		 1383eisneg (x)
1384 unsigned EMUSHORT x[];
1385{
1386
1387 if (x[NE - 1] & 0x8000)
1388 return (1);
1389 else
1390 return (0);
1391}
1392
1393
defb5dab 1394/* Return 1 if external format number is infinity, else return zero. */
a0353055
RK		 1395
1396static int
985b6196
RS		 1397eisinf (x)
1398 unsigned EMUSHORT x[];
1399{
1400
66b6d60b
RS		 1401#ifdef NANS
1402 if (eisnan (x))
1403 return (0);
1404#endif
985b6196
RS		 1405 if ((x[NE - 1] & 0x7fff) == 0x7fff)
1406 return (1);
1407 else
1408 return (0);
1409}
1410
1411
defb5dab
RK		 1412/* Check if e-type number is not a number. The bit pattern is one that we
1413 defined, so we know for sure how to detect it. */
66b6d60b 1414
a0353055 1415static int
66b6d60b
RS		 1416eisnan (x)
1417 unsigned EMUSHORT x[];
1418{
66b6d60b
RS		 1419#ifdef NANS
1420 int i;
defb5dab
RK		 1421
1422 /* NaN has maximum exponent */
66b6d60b
RS		 1423 if ((x[NE - 1] & 0x7fff) != 0x7fff)
1424 return (0);
defb5dab 1425 /* ... and non-zero significand field. */
66b6d60b
RS		 1426 for (i = 0; i < NE - 1; i++)
1427 {
1428 if (*x++ != 0)
1429 return (1);
1430 }
1431#endif
defb5dab 1432
66b6d60b
RS		 1433 return (0);
1434}
1435
1436/* Fill external format number with infinity pattern (IEEE)
842fbaaa 1437 or largest possible number (non-IEEE). */
985b6196 1438
a0353055 1439static void
985b6196
RS		 1440einfin (x)
1441 register unsigned EMUSHORT *x;
1442{
1443 register int i;
1444
1445#ifdef INFINITY
1446 for (i = 0; i < NE - 1; i++)
1447 *x++ = 0;
1448 *x |= 32767;
1449#else
1450 for (i = 0; i < NE - 1; i++)
1451 *x++ = 0xffff;
1452 *x |= 32766;
1453 if (rndprc < NBITS)
1454 {
842fbaaa
JW		 1455 if (rndprc == 113)
1456 {
1457 *(x - 9) = 0;
1458 *(x - 8) = 0;
1459 }
985b6196
RS		 1460 if (rndprc == 64)
1461 {
1462 *(x - 5) = 0;
1463 }
1464 if (rndprc == 53)
1465 {
1466 *(x - 4) = 0xf800;
1467 }
1468 else
1469 {
1470 *(x - 4) = 0;
1471 *(x - 3) = 0;
1472 *(x - 2) = 0xff00;
1473 }
1474 }
1475#endif
1476}
1477
1478
66b6d60b
RS		 1479/* Output an e-type NaN.
1480 This generates Intel's quiet NaN pattern for extended real.
1481 The exponent is 7fff, the leading mantissa word is c000. */
1482
a0353055 1483static void
29e11dab 1484enan (x, sign)
66b6d60b 1485 register unsigned EMUSHORT *x;
29e11dab 1486 int sign;
66b6d60b
RS		 1487{
1488 register int i;
1489
1490 for (i = 0; i < NE - 2; i++)
1491 *x++ = 0;
1492 *x++ = 0xc000;
29e11dab 1493 *x = (sign << 15) | 0x7fff;
66b6d60b
RS		 1494}
1495
985b6196 1496
defb5dab
RK		 1497/* Move in external format number, converting it to internal format. */
1498
a0353055 1499static void
985b6196
RS		 1500emovi (a, b)
1501 unsigned EMUSHORT *a, *b;
1502{
1503 register unsigned EMUSHORT *p, *q;
1504 int i;
1505
1506 q = b;
1507 p = a + (NE - 1); /* point to last word of external number */
1508 /* get the sign bit */
1509 if (*p & 0x8000)
1510 *q++ = 0xffff;
1511 else
1512 *q++ = 0;
1513 /* get the exponent */
1514 *q = *p--;
1515 *q++ &= 0x7fff; /* delete the sign bit */
1516#ifdef INFINITY
1517 if ((*(q - 1) & 0x7fff) == 0x7fff)
1518 {
66b6d60b
RS		 1519#ifdef NANS
1520 if (eisnan (a))
1521 {
1522 *q++ = 0;
1523 for (i = 3; i < NI; i++)
1524 *q++ = *p--;
1525 return;
1526 }
1527#endif
defb5dab 1528
985b6196
RS		 1529 for (i = 2; i < NI; i++)
1530 *q++ = 0;
1531 return;
1532 }
1533#endif
defb5dab 1534
985b6196
RS		 1535 /* clear high guard word */
1536 *q++ = 0;
1537 /* move in the significand */
1538 for (i = 0; i < NE - 1; i++)
1539 *q++ = *p--;
1540 /* clear low guard word */
1541 *q = 0;
1542}
1543
1544
defb5dab
RK		 1545/* Move internal format number out, converting it to external format. */
1546
a0353055 1547static void
985b6196
RS		 1548emovo (a, b)
1549 unsigned EMUSHORT *a, *b;
1550{
1551 register unsigned EMUSHORT *p, *q;
1552 unsigned EMUSHORT i;
239b043b 1553 int j;
985b6196
RS		 1554
1555 p = a;
1556 q = b + (NE - 1); /* point to output exponent */
1557 /* combine sign and exponent */
1558 i = *p++;
1559 if (i)
1560 *q-- = *p++ | 0x8000;
1561 else
1562 *q-- = *p++;
1563#ifdef INFINITY
1564 if (*(p - 1) == 0x7fff)
1565 {
66b6d60b
RS		 1566#ifdef NANS
1567 if (eiisnan (a))
1568 {
29e11dab 1569 enan (b, eiisneg (a));
66b6d60b
RS		 1570 return;
1571 }
1572#endif
985b6196 1573 einfin (b);
842fbaaa 1574 return;
985b6196
RS		 1575 }
1576#endif
1577 /* skip over guard word */
1578 ++p;
1579 /* move the significand */
239b043b 1580 for (j = 0; j < NE - 1; j++)
985b6196
RS		 1581 *q-- = *p++;
1582}
1583
defb5dab 1584/* Clear out internal format number. */
985b6196 1585
a0353055 1586static void
985b6196
RS		 1587ecleaz (xi)
1588 register unsigned EMUSHORT *xi;
1589{
1590 register int i;
1591
1592 for (i = 0; i < NI; i++)
1593 *xi++ = 0;
1594}
1595
1596
defb5dab 1597/* Same, but don't touch the sign. */
985b6196 1598
a0353055 1599static void
985b6196
RS		 1600ecleazs (xi)
1601 register unsigned EMUSHORT *xi;
1602{
1603 register int i;
1604
1605 ++xi;
1606 for (i = 0; i < NI - 1; i++)
1607 *xi++ = 0;
1608}
1609
1610
1611
defb5dab 1612/* Move internal format number from a to b. */
a0353055
RK		 1613
1614static void
985b6196
RS		 1615emovz (a, b)
1616 register unsigned EMUSHORT *a, *b;
1617{
1618 register int i;
1619
1620 for (i = 0; i < NI - 1; i++)
1621 *b++ = *a++;
1622 /* clear low guard word */
1623 *b = 0;
1624}
1625
66b6d60b
RS		 1626/* Generate internal format NaN.
1627 The explicit pattern for this is maximum exponent and
defb5dab 1628 top two significant bits set. */
66b6d60b 1629
a0353055 1630static void
66b6d60b
RS		 1631einan (x)
1632 unsigned EMUSHORT x[];
1633{
1634
1635 ecleaz (x);
1636 x[E] = 0x7fff;
1637 x[M + 1] = 0xc000;
1638}
1639
1640/* Return nonzero if internal format number is a NaN. */
1641
a0353055 1642static int
66b6d60b
RS		 1643eiisnan (x)
1644 unsigned EMUSHORT x[];
1645{
1646 int i;
1647
1648 if ((x[E] & 0x7fff) == 0x7fff)
1649 {
1650 for (i = M + 1; i < NI; i++)
1651 {
1652 if (x[i] != 0)
1653 return (1);
1654 }
1655 }
1656 return (0);
1657}
1658
29e11dab
RK		 1659/* Return nonzero if sign of internal format number is nonzero. */
1660
a0353055 1661static int
29e11dab
RK		 1662eiisneg (x)
1663 unsigned EMUSHORT x[];
1664{
1665
1666 return x[0] != 0;
1667}
1668
66b6d60b
RS		 1669/* Fill internal format number with infinity pattern.
1670 This has maximum exponent and significand all zeros. */
1671
a0353055 1672static void
66b6d60b
RS		 1673eiinfin (x)
1674 unsigned EMUSHORT x[];
1675{
1676
1677 ecleaz (x);
1678 x[E] = 0x7fff;
1679}
1680
1681/* Return nonzero if internal format number is infinite. */
1682
a0353055 1683static int
66b6d60b
RS		 1684eiisinf (x)
1685 unsigned EMUSHORT x[];
1686{
1687
1688#ifdef NANS
1689 if (eiisnan (x))
1690 return (0);
1691#endif
1692 if ((x[E] & 0x7fff) == 0x7fff)
1693 return (1);
1694 return (0);
1695}
1696
985b6196 1697
defb5dab
RK		 1698/* Compare significands of numbers in internal format.
1699 Guard words are included in the comparison.
1700
1701 Returns +1 if a > b
1702 0 if a == b
1703 -1 if a < b */
a0353055
RK		 1704
1705static int
985b6196
RS		 1706ecmpm (a, b)
1707 register unsigned EMUSHORT *a, *b;
1708{
1709 int i;
1710
1711 a += M; /* skip up to significand area */
1712 b += M;
1713 for (i = M; i < NI; i++)
1714 {
1715 if (*a++ != *b++)
1716 goto difrnt;
1717 }
1718 return (0);
1719
1720 difrnt:
1721 if (*(--a) > *(--b))
1722 return (1);
1723 else
1724 return (-1);
1725}
1726
1727
defb5dab 1728/* Shift significand down by 1 bit. */
985b6196 1729
a0353055 1730static void
985b6196
RS		 1731eshdn1 (x)
1732 register unsigned EMUSHORT *x;
1733{
1734 register unsigned EMUSHORT bits;
1735 int i;
1736
1737 x += M; /* point to significand area */
1738
1739 bits = 0;
1740 for (i = M; i < NI; i++)
1741 {
1742 if (*x & 1)
1743 bits |= 1;
1744 *x >>= 1;
1745 if (bits & 2)
1746 *x |= 0x8000;
1747 bits <<= 1;
1748 ++x;
1749 }
1750}
1751
1752
1753
defb5dab 1754/* Shift significand up by 1 bit. */
985b6196 1755
a0353055 1756static void
985b6196
RS		 1757eshup1 (x)
1758 register unsigned EMUSHORT *x;
1759{
1760 register unsigned EMUSHORT bits;
1761 int i;
1762
1763 x += NI - 1;
1764 bits = 0;
1765
1766 for (i = M; i < NI; i++)
1767 {
1768 if (*x & 0x8000)
1769 bits |= 1;
1770 *x <<= 1;
1771 if (bits & 2)
1772 *x |= 1;
1773 bits <<= 1;
1774 --x;
1775 }
1776}
1777
1778
defb5dab 1779/* Shift significand down by 8 bits. */
985b6196 1780
a0353055 1781static void
985b6196
RS		 1782eshdn8 (x)
1783 register unsigned EMUSHORT *x;
1784{
1785 register unsigned EMUSHORT newbyt, oldbyt;
1786 int i;
1787
1788 x += M;
1789 oldbyt = 0;
1790 for (i = M; i < NI; i++)
1791 {
1792 newbyt = *x << 8;
1793 *x >>= 8;
1794 *x |= oldbyt;
1795 oldbyt = newbyt;
1796 ++x;
1797 }
1798}
1799
defb5dab 1800/* Shift significand up by 8 bits. */
985b6196 1801
a0353055 1802static void
985b6196
RS		 1803eshup8 (x)
1804 register unsigned EMUSHORT *x;
1805{
1806 int i;
1807 register unsigned EMUSHORT newbyt, oldbyt;
1808
1809 x += NI - 1;
1810 oldbyt = 0;
1811
1812 for (i = M; i < NI; i++)
1813 {
1814 newbyt = *x >> 8;
1815 *x <<= 8;
1816 *x |= oldbyt;
1817 oldbyt = newbyt;
1818 --x;
1819 }
1820}
1821
defb5dab 1822/* Shift significand up by 16 bits. */
985b6196 1823
a0353055 1824static void
985b6196
RS		 1825eshup6 (x)
1826 register unsigned EMUSHORT *x;
1827{
1828 int i;
1829 register unsigned EMUSHORT *p;
1830
1831 p = x + M;
1832 x += M + 1;
1833
1834 for (i = M; i < NI - 1; i++)
1835 *p++ = *x++;
1836
1837 *p = 0;
1838}
1839
defb5dab 1840/* Shift significand down by 16 bits. */
985b6196 1841
a0353055 1842static void
985b6196
RS		 1843eshdn6 (x)
1844 register unsigned EMUSHORT *x;
1845{
1846 int i;
1847 register unsigned EMUSHORT *p;
1848
1849 x += NI - 1;
1850 p = x + 1;
1851
1852 for (i = M; i < NI - 1; i++)
1853 *(--p) = *(--x);
1854
1855 *(--p) = 0;
1856}
1857\f
defb5dab 1858/* Add significands. x + y replaces y. */
985b6196 1859
a0353055 1860static void
985b6196
RS		 1861eaddm (x, y)
1862 unsigned EMUSHORT *x, *y;
1863{
1864 register unsigned EMULONG a;
1865 int i;
1866 unsigned int carry;
1867
1868 x += NI - 1;
1869 y += NI - 1;
1870 carry = 0;
1871 for (i = M; i < NI; i++)
1872 {
1873 a = (unsigned EMULONG) (*x) + (unsigned EMULONG) (*y) + carry;
1874 if (a & 0x10000)
1875 carry = 1;
1876 else
1877 carry = 0;
1878 *y = (unsigned EMUSHORT) a;
1879 --x;
1880 --y;
1881 }
1882}
1883
defb5dab 1884/* Subtract significands. y - x replaces y. */
985b6196 1885
a0353055 1886static void
985b6196
RS		 1887esubm (x, y)
1888 unsigned EMUSHORT *x, *y;
1889{
1890 unsigned EMULONG a;
1891 int i;
1892 unsigned int carry;
1893
1894 x += NI - 1;
1895 y += NI - 1;
1896 carry = 0;
1897 for (i = M; i < NI; i++)
1898 {
1899 a = (unsigned EMULONG) (*y) - (unsigned EMULONG) (*x) - carry;
1900 if (a & 0x10000)
1901 carry = 1;
1902 else
1903 carry = 0;
1904 *y = (unsigned EMUSHORT) a;
1905 --x;
1906 --y;
1907 }
1908}
1909
1910
985b6196
RS		 1911static unsigned EMUSHORT equot[NI];
1912
842fbaaa
JW		 1913
1914#if 0
1915/* Radix 2 shift-and-add versions of multiply and divide */
1916
1917
1918/* Divide significands */
1919
985b6196
RS		 1920int
1921edivm (den, num)
1922 unsigned EMUSHORT den[], num[];
1923{
1924 int i;
1925 register unsigned EMUSHORT *p, *q;
1926 unsigned EMUSHORT j;
1927
1928 p = &equot[0];
1929 *p++ = num[0];
1930 *p++ = num[1];
1931
1932 for (i = M; i < NI; i++)
1933 {
1934 *p++ = 0;
1935 }
1936
defb5dab
RK		 1937 /* Use faster compare and subtraction if denominator has only 15 bits of
1938 significance. */
1939
985b6196
RS		 1940 p = &den[M + 2];
1941 if (*p++ == 0)
1942 {
1943 for (i = M + 3; i < NI; i++)
1944 {
1945 if (*p++ != 0)
1946 goto fulldiv;
1947 }
1948 if ((den[M + 1] & 1) != 0)
1949 goto fulldiv;
1950 eshdn1 (num);
1951 eshdn1 (den);
1952
1953 p = &den[M + 1];
1954 q = &num[M + 1];
1955
1956 for (i = 0; i < NBITS + 2; i++)
1957 {
1958 if (*p <= *q)
1959 {
1960 *q -= *p;
1961 j = 1;
1962 }
1963 else
1964 {
1965 j = 0;
1966 }
1967 eshup1 (equot);
1968 equot[NI - 2] |= j;
1969 eshup1 (num);
1970 }
1971 goto divdon;
1972 }
1973
defb5dab
RK		 1974 /* The number of quotient bits to calculate is NBITS + 1 scaling guard
1975 bit + 1 roundoff bit. */
1976
985b6196
RS		 1977 fulldiv:
1978
1979 p = &equot[NI - 2];
1980 for (i = 0; i < NBITS + 2; i++)
1981 {
1982 if (ecmpm (den, num) <= 0)
1983 {
1984 esubm (den, num);
1985 j = 1; /* quotient bit = 1 */
1986 }
1987 else
1988 j = 0;
1989 eshup1 (equot);
1990 *p |= j;
1991 eshup1 (num);
1992 }
1993
1994 divdon:
1995
1996 eshdn1 (equot);
1997 eshdn1 (equot);
1998
1999 /* test for nonzero remainder after roundoff bit */
2000 p = &num[M];
2001 j = 0;
2002 for (i = M; i < NI; i++)
2003 {
2004 j |= *p++;
2005 }
2006 if (j)
2007 j = 1;
2008
2009
2010 for (i = 0; i < NI; i++)
2011 num[i] = equot[i];
2012 return ((int) j);
2013}
2014
2015
2016/* Multiply significands */
2017int
2018emulm (a, b)
2019 unsigned EMUSHORT a[], b[];
2020{
2021 unsigned EMUSHORT *p, *q;
2022 int i, j, k;
2023
2024 equot[0] = b[0];
2025 equot[1] = b[1];
2026 for (i = M; i < NI; i++)
2027 equot[i] = 0;
2028
2029 p = &a[NI - 2];
2030 k = NBITS;
defb5dab 2031 while (*p == 0) /* significand is not supposed to be zero */
985b6196
RS		 2032 {
2033 eshdn6 (a);
2034 k -= 16;
2035 }
2036 if ((*p & 0xff) == 0)
2037 {
2038 eshdn8 (a);
2039 k -= 8;
2040 }
2041
2042 q = &equot[NI - 1];
2043 j = 0;
2044 for (i = 0; i < k; i++)
2045 {
2046 if (*p & 1)
2047 eaddm (b, equot);
2048 /* remember if there were any nonzero bits shifted out */
2049 if (*q & 1)
2050 j |= 1;
2051 eshdn1 (a);
2052 eshdn1 (equot);
2053 }
2054
2055 for (i = 0; i < NI; i++)
2056 b[i] = equot[i];
2057
2058 /* return flag for lost nonzero bits */
2059 return (j);
2060}
2061
842fbaaa
JW		 2062#else
2063
2064/* Radix 65536 versions of multiply and divide */
2065
2066
2067/* Multiply significand of e-type number b
a0353055 2068 by 16-bit quantity a, e-type result to c. */
842fbaaa 2069
a0353055 2070static void
242cef1e 2071m16m (a, b, c)
c92d992a 2072 unsigned int a;
242cef1e 2073 unsigned short b[], c[];
842fbaaa 2074{
242cef1e
RS		 2075 register unsigned short *pp;
2076 register unsigned long carry;
2077 unsigned short *ps;
2078 unsigned short p[NI];
2079 unsigned long aa, m;
2080 int i;
2081
2082 aa = a;
2083 pp = &p[NI-2];
2084 *pp++ = 0;
2085 *pp = 0;
2086 ps = &b[NI-1];
2087
2088 for (i=M+1; i<NI; i++)
2089 {
2090 if (*ps == 0)
842fbaaa 2091 {
242cef1e
RS		 2092 --ps;
2093 --pp;
2094 *(pp-1) = 0;
842fbaaa 2095 }
242cef1e
RS		 2096 else
2097 {
2098 m = (unsigned long) aa * *ps--;
2099 carry = (m & 0xffff) + *pp;
2100 *pp-- = (unsigned short)carry;
2101 carry = (carry >> 16) + (m >> 16) + *pp;
2102 *pp = (unsigned short)carry;
2103 *(pp-1) = carry >> 16;
2104 }
2105 }
2106 for (i=M; i<NI; i++)
2107 c[i] = p[i];
842fbaaa
JW		 2108}
2109
2110
2111/* Divide significands. Neither the numerator nor the denominator
242cef1e 2112 is permitted to have its high guard word nonzero. */
842fbaaa 2113
a0353055 2114static int
242cef1e
RS		 2115edivm (den, num)
2116 unsigned short den[], num[];
842fbaaa 2117{
242cef1e
RS		 2118 int i;
2119 register unsigned short *p;
2120 unsigned long tnum;
2121 unsigned short j, tdenm, tquot;
2122 unsigned short tprod[NI+1];
842fbaaa 2123
242cef1e
RS		 2124 p = &equot[0];
2125 *p++ = num[0];
2126 *p++ = num[1];
842fbaaa 2127
242cef1e
RS		 2128 for (i=M; i<NI; i++)
2129 {
2130 *p++ = 0;
2131 }
2132 eshdn1 (num);
2133 tdenm = den[M+1];
2134 for (i=M; i<NI; i++)
2135 {
2136 /* Find trial quotient digit (the radix is 65536). */
2137 tnum = (((unsigned long) num[M]) << 16) + num[M+1];
2138
2139 /* Do not execute the divide instruction if it will overflow. */
2140 if ((tdenm * 0xffffL) < tnum)
2141 tquot = 0xffff;
2142 else
2143 tquot = tnum / tdenm;
2144 /* Multiply denominator by trial quotient digit. */
c92d992a 2145 m16m ((unsigned int)tquot, den, tprod);
242cef1e
RS		 2146 /* The quotient digit may have been overestimated. */
2147 if (ecmpm (tprod, num) > 0)
842fbaaa 2148 {
242cef1e
RS		 2149 tquot -= 1;
2150 esubm (den, tprod);
2151 if (ecmpm (tprod, num) > 0)
2152 {
2153 tquot -= 1;
2154 esubm (den, tprod);
2155 }
842fbaaa 2156 }
242cef1e
RS		 2157 esubm (tprod, num);
2158 equot[i] = tquot;
2159 eshup6(num);
2160 }
2161 /* test for nonzero remainder after roundoff bit */
2162 p = &num[M];
2163 j = 0;
2164 for (i=M; i<NI; i++)
2165 {
2166 j |= *p++;
2167 }
2168 if (j)
2169 j = 1;
842fbaaa 2170
242cef1e
RS		 2171 for (i=0; i<NI; i++)
2172 num[i] = equot[i];
842fbaaa 2173
242cef1e 2174 return ((int)j);
842fbaaa
JW		 2175}
2176
2177
2178
2179/* Multiply significands */
a0353055 2180static int
242cef1e
RS		 2181emulm (a, b)
2182 unsigned short a[], b[];
842fbaaa 2183{
242cef1e
RS		 2184 unsigned short *p, *q;
2185 unsigned short pprod[NI];
2186 unsigned short j;
2187 int i;
2188
2189 equot[0] = b[0];
2190 equot[1] = b[1];
2191 for (i=M; i<NI; i++)
2192 equot[i] = 0;
2193
2194 j = 0;
2195 p = &a[NI-1];
2196 q = &equot[NI-1];
2197 for (i=M+1; i<NI; i++)
2198 {
2199 if (*p == 0)
842fbaaa 2200 {
242cef1e
RS		 2201 --p;
2202 }
2203 else
2204 {
c92d992a 2205 m16m ((unsigned int) *p--, b, pprod);
242cef1e 2206 eaddm(pprod, equot);
842fbaaa 2207 }
242cef1e
RS		 2208 j |= *q;
2209 eshdn6(equot);
2210 }
842fbaaa 2211
242cef1e
RS		 2212 for (i=0; i<NI; i++)
2213 b[i] = equot[i];
842fbaaa 2214
242cef1e
RS		 2215 /* return flag for lost nonzero bits */
2216 return ((int)j);
842fbaaa
JW		 2217}
2218#endif
985b6196
RS		 2219
2220
defb5dab 2221/* Normalize and round off.
985b6196 2222
defb5dab
RK		 2223 The internal format number to be rounded is "s".
2224 Input "lost" indicates whether or not the number is exact.
2225 This is the so-called sticky bit.
2226
2227 Input "subflg" indicates whether the number was obtained
2228 by a subtraction operation. In that case if lost is nonzero
2229 then the number is slightly smaller than indicated.
2230
2231 Input "exp" is the biased exponent, which may be negative.
2232 the exponent field of "s" is ignored but is replaced by
2233 "exp" as adjusted by normalization and rounding.
2234
2235 Input "rcntrl" is the rounding control.
2236
2237 For future reference: In order for emdnorm to round off denormal
842fbaaa
JW		 2238 significands at the right point, the input exponent must be
2239 adjusted to be the actual value it would have after conversion to
2240 the final floating point type. This adjustment has been
2241 implemented for all type conversions (etoe53, etc.) and decimal
2242 conversions, but not for the arithmetic functions (eadd, etc.).
2243 Data types having standard 15-bit exponents are not affected by
2244 this, but SFmode and DFmode are affected. For example, ediv with
2245 rndprc = 24 will not round correctly to 24-bit precision if the
2246 result is denormal. */
2247
985b6196
RS		 2248static int rlast = -1;
2249static int rw = 0;
2250static unsigned EMUSHORT rmsk = 0;
2251static unsigned EMUSHORT rmbit = 0;
2252static unsigned EMUSHORT rebit = 0;
2253static int re = 0;
2254static unsigned EMUSHORT rbit[NI];
2255
a0353055 2256static void
985b6196
RS		 2257emdnorm (s, lost, subflg, exp, rcntrl)
2258 unsigned EMUSHORT s[];
2259 int lost;
2260 int subflg;
2261 EMULONG exp;
2262 int rcntrl;
2263{
2264 int i, j;
2265 unsigned EMUSHORT r;
2266
2267 /* Normalize */
2268 j = enormlz (s);
2269
2270 /* a blank significand could mean either zero or infinity. */
2271#ifndef INFINITY
2272 if (j > NBITS)
2273 {
2274 ecleazs (s);
2275 return;
2276 }
2277#endif
2278 exp -= j;
2279#ifndef INFINITY
2280 if (exp >= 32767L)
2281 goto overf;
2282#else
2283 if ((j > NBITS) && (exp < 32767))
2284 {
2285 ecleazs (s);
2286 return;
2287 }
2288#endif
2289 if (exp < 0L)
2290 {
2291 if (exp > (EMULONG) (-NBITS - 1))
2292 {
2293 j = (int) exp;
2294 i = eshift (s, j);
2295 if (i)
2296 lost = 1;
2297 }
2298 else
2299 {
2300 ecleazs (s);
2301 return;
2302 }
2303 }
2304 /* Round off, unless told not to by rcntrl. */
2305 if (rcntrl == 0)
2306 goto mdfin;
2307 /* Set up rounding parameters if the control register changed. */
2308 if (rndprc != rlast)
2309 {
2310 ecleaz (rbit);
2311 switch (rndprc)
2312 {
2313 default:
2314 case NBITS:
2315 rw = NI - 1; /* low guard word */
2316 rmsk = 0xffff;
2317 rmbit = 0x8000;
842fbaaa 2318 re = rw - 1;
985b6196
RS		 2319 rebit = 1;
2320 break;
842fbaaa
JW		 2321 case 113:
2322 rw = 10;
2323 rmsk = 0x7fff;
2324 rmbit = 0x4000;
2325 rebit = 0x8000;
2326 re = rw;
2327 break;
985b6196
RS		 2328 case 64:
2329 rw = 7;
2330 rmsk = 0xffff;
2331 rmbit = 0x8000;
985b6196
RS		 2332 re = rw - 1;
2333 rebit = 1;
2334 break;
842fbaaa 2335 /* For DEC or IBM arithmetic */
985b6196
RS		 2336 case 56:
2337 rw = 6;
2338 rmsk = 0xff;
2339 rmbit = 0x80;
985b6196 2340 rebit = 0x100;
842fbaaa 2341 re = rw;
985b6196
RS		 2342 break;
2343 case 53:
2344 rw = 6;
2345 rmsk = 0x7ff;
2346 rmbit = 0x0400;
985b6196 2347 rebit = 0x800;
842fbaaa 2348 re = rw;
985b6196
RS		 2349 break;
2350 case 24:
2351 rw = 4;
2352 rmsk = 0xff;
2353 rmbit = 0x80;
985b6196 2354 rebit = 0x100;
842fbaaa 2355 re = rw;
985b6196
RS		 2356 break;
2357 }
842fbaaa 2358 rbit[re] = rebit;
985b6196
RS		 2359 rlast = rndprc;
2360 }
2361
842fbaaa
JW		 2362 /* Shift down 1 temporarily if the data structure has an implied
2363 most significant bit and the number is denormal. */
2364 if ((exp <= 0) && (rndprc != 64) && (rndprc != NBITS))
985b6196 2365 {
842fbaaa
JW		 2366 lost |= s[NI - 1] & 1;
2367 eshdn1 (s);
985b6196 2368 }
842fbaaa
JW		 2369 /* Clear out all bits below the rounding bit,
2370 remembering in r if any were nonzero. */
2371 r = s[rw] & rmsk;
2372 if (rndprc < NBITS)
985b6196 2373 {
985b6196
RS		 2374 i = rw + 1;
2375 while (i < NI)
2376 {
2377 if (s[i])
2378 r |= 1;
2379 s[i] = 0;
2380 ++i;
2381 }
985b6196 2382 }
afb817fd 2383 s[rw] &= ~rmsk;
985b6196
RS		 2384 if ((r & rmbit) != 0)
2385 {
2386 if (r == rmbit)
2387 {
2388 if (lost == 0)
2389 { /* round to even */
2390 if ((s[re] & rebit) == 0)
2391 goto mddone;
2392 }
2393 else
2394 {
2395 if (subflg != 0)
2396 goto mddone;
2397 }
2398 }
2399 eaddm (rbit, s);
2400 }
2401 mddone:
842fbaaa 2402 if ((exp <= 0) && (rndprc != 64) && (rndprc != NBITS))
985b6196
RS		 2403 {
2404 eshup1 (s);
2405 }
2406 if (s[2] != 0)
2407 { /* overflow on roundoff */
2408 eshdn1 (s);
2409 exp += 1;
2410 }
2411 mdfin:
2412 s[NI - 1] = 0;
2413 if (exp >= 32767L)
2414 {
2415#ifndef INFINITY
2416 overf:
2417#endif
2418#ifdef INFINITY
2419 s[1] = 32767;
2420 for (i = 2; i < NI - 1; i++)
2421 s[i] = 0;
64685ffa
RS		 2422 if (extra_warnings)
2423 warning ("floating point overflow");
985b6196
RS		 2424#else
2425 s[1] = 32766;
2426 s[2] = 0;
2427 for (i = M + 1; i < NI - 1; i++)
2428 s[i] = 0xffff;
2429 s[NI - 1] = 0;
842fbaaa 2430 if ((rndprc < 64) || (rndprc == 113))
985b6196
RS		 2431 {
2432 s[rw] &= ~rmsk;
2433 if (rndprc == 24)
2434 {
2435 s[5] = 0;
2436 s[6] = 0;
2437 }
2438 }
2439#endif
2440 return;
2441 }
2442 if (exp < 0)
2443 s[1] = 0;
2444 else
2445 s[1] = (unsigned EMUSHORT) exp;
2446}
2447
2448
2449
defb5dab 2450/* Subtract external format numbers. */
985b6196
RS		 2451
2452static int subflg = 0;
2453
a0353055 2454static void
985b6196
RS		 2455esub (a, b, c)
2456 unsigned EMUSHORT *a, *b, *c;
2457{
2458
66b6d60b
RS		 2459#ifdef NANS
2460 if (eisnan (a))
2461 {
2462 emov (a, c);
2463 return;
2464 }
2465 if (eisnan (b))
2466 {
2467 emov (b, c);
2468 return;
2469 }
2470/* Infinity minus infinity is a NaN.
2471 Test for subtracting infinities of the same sign. */
2472 if (eisinf (a) && eisinf (b)
2473 && ((eisneg (a) ^ eisneg (b)) == 0))
2474 {
2475 mtherr ("esub", INVALID);
29e11dab 2476 enan (c, 0);
66b6d60b
RS		 2477 return;
2478 }
2479#endif
985b6196
RS		 2480 subflg = 1;
2481 eadd1 (a, b, c);
2482}
2483
2484
defb5dab 2485/* Add. */
a0353055
RK		 2486
2487static void
985b6196
RS		 2488eadd (a, b, c)
2489 unsigned EMUSHORT *a, *b, *c;
2490{
2491
66b6d60b
RS		 2492#ifdef NANS
2493/* NaN plus anything is a NaN. */
2494 if (eisnan (a))
2495 {
2496 emov (a, c);
2497 return;
2498 }
2499 if (eisnan (b))
2500 {
2501 emov (b, c);
2502 return;
2503 }
2504/* Infinity minus infinity is a NaN.
2505 Test for adding infinities of opposite signs. */
2506 if (eisinf (a) && eisinf (b)
2507 && ((eisneg (a) ^ eisneg (b)) != 0))
2508 {
2509 mtherr ("esub", INVALID);
29e11dab 2510 enan (c, 0);
66b6d60b
RS		 2511 return;
2512 }
2513#endif
985b6196
RS		 2514 subflg = 0;
2515 eadd1 (a, b, c);
2516}
2517
a0353055 2518static void
985b6196
RS		 2519eadd1 (a, b, c)
2520 unsigned EMUSHORT *a, *b, *c;
2521{
2522 unsigned EMUSHORT ai[NI], bi[NI], ci[NI];
2523 int i, lost, j, k;
2524 EMULONG lt, lta, ltb;
2525
2526#ifdef INFINITY
2527 if (eisinf (a))
2528 {
2529 emov (a, c);
2530 if (subflg)
2531 eneg (c);
2532 return;
2533 }
2534 if (eisinf (b))
2535 {
2536 emov (b, c);
2537 return;
2538 }
2539#endif
2540 emovi (a, ai);
2541 emovi (b, bi);
2542 if (subflg)
2543 ai[0] = ~ai[0];
2544
2545 /* compare exponents */
2546 lta = ai[E];
2547 ltb = bi[E];
2548 lt = lta - ltb;
2549 if (lt > 0L)
2550 { /* put the larger number in bi */
2551 emovz (bi, ci);
2552 emovz (ai, bi);
2553 emovz (ci, ai);
2554 ltb = bi[E];
2555 lt = -lt;
2556 }
2557 lost = 0;
2558 if (lt != 0L)
2559 {
2560 if (lt < (EMULONG) (-NBITS - 1))
2561 goto done; /* answer same as larger addend */
2562 k = (int) lt;
2563 lost = eshift (ai, k); /* shift the smaller number down */
2564 }
2565 else
2566 {
2567 /* exponents were the same, so must compare significands */
2568 i = ecmpm (ai, bi);
2569 if (i == 0)
2570 { /* the numbers are identical in magnitude */
2571 /* if different signs, result is zero */
2572 if (ai[0] != bi[0])
2573 {
2574 eclear (c);
2575 return;
2576 }
2577 /* if same sign, result is double */
2578 /* double denomalized tiny number */
2579 if ((bi[E] == 0) && ((bi[3] & 0x8000) == 0))
2580 {
2581 eshup1 (bi);
2582 goto done;
2583 }
2584 /* add 1 to exponent unless both are zero! */
2585 for (j = 1; j < NI - 1; j++)
2586 {
2587 if (bi[j] != 0)
2588 {
2589 /* This could overflow, but let emovo take care of that. */
2590 ltb += 1;
2591 break;
2592 }
2593 }
2594 bi[E] = (unsigned EMUSHORT) ltb;
2595 goto done;
2596 }
2597 if (i > 0)
2598 { /* put the larger number in bi */
2599 emovz (bi, ci);
2600 emovz (ai, bi);
2601 emovz (ci, ai);
2602 }
2603 }
2604 if (ai[0] == bi[0])
2605 {
2606 eaddm (ai, bi);
2607 subflg = 0;
2608 }
2609 else
2610 {
2611 esubm (ai, bi);
2612 subflg = 1;
2613 }
2614 emdnorm (bi, lost, subflg, ltb, 64);
2615
2616 done:
2617 emovo (bi, c);
2618}
2619
2620
2621
defb5dab 2622/* Divide. */
a0353055
RK		 2623
2624static void
985b6196
RS		 2625ediv (a, b, c)
2626 unsigned EMUSHORT *a, *b, *c;
2627{
2628 unsigned EMUSHORT ai[NI], bi[NI];
2629 int i;
2630 EMULONG lt, lta, ltb;
2631
66b6d60b
RS		 2632#ifdef NANS
2633/* Return any NaN input. */
2634 if (eisnan (a))
2635 {
2636 emov (a, c);
2637 return;
2638 }
2639 if (eisnan (b))
2640 {
2641 emov (b, c);
2642 return;
2643 }
2644/* Zero over zero, or infinity over infinity, is a NaN. */
2645 if (((ecmp (a, ezero) == 0) && (ecmp (b, ezero) == 0))
2646 || (eisinf (a) && eisinf (b)))
2647 {
2648 mtherr ("ediv", INVALID);
29e11dab 2649 enan (c, eisneg (a) ^ eisneg (b));
66b6d60b
RS		 2650 return;
2651 }
2652#endif
2653/* Infinity over anything else is infinity. */
985b6196
RS		 2654#ifdef INFINITY
2655 if (eisinf (b))
2656 {
2657 if (eisneg (a) ^ eisneg (b))
2658 *(c + (NE - 1)) = 0x8000;
2659 else
2660 *(c + (NE - 1)) = 0;
2661 einfin (c);
2662 return;
2663 }
66b6d60b 2664/* Anything else over infinity is zero. */
985b6196
RS		 2665 if (eisinf (a))
2666 {
2667 eclear (c);
2668 return;
2669 }
2670#endif
2671 emovi (a, ai);
2672 emovi (b, bi);
2673 lta = ai[E];
2674 ltb = bi[E];
2675 if (bi[E] == 0)
2676 { /* See if numerator is zero. */
2677 for (i = 1; i < NI - 1; i++)
2678 {
2679 if (bi[i] != 0)
2680 {
2681 ltb -= enormlz (bi);
2682 goto dnzro1;
2683 }
2684 }
2685 eclear (c);
2686 return;
2687 }
2688 dnzro1:
2689
2690 if (ai[E] == 0)
2691 { /* possible divide by zero */
2692 for (i = 1; i < NI - 1; i++)
2693 {
2694 if (ai[i] != 0)
2695 {
2696 lta -= enormlz (ai);
2697 goto dnzro2;
2698 }
2699 }
2700 if (ai[0] == bi[0])
2701 *(c + (NE - 1)) = 0;
2702 else
2703 *(c + (NE - 1)) = 0x8000;
66b6d60b
RS		 2704/* Divide by zero is not an invalid operation.
2705 It is a divide-by-zero operation! */
985b6196
RS		 2706 einfin (c);
2707 mtherr ("ediv", SING);
2708 return;
2709 }
2710 dnzro2:
2711
2712 i = edivm (ai, bi);
2713 /* calculate exponent */
2714 lt = ltb - lta + EXONE;
2715 emdnorm (bi, i, 0, lt, 64);
2716 /* set the sign */
2717 if (ai[0] == bi[0])
2718 bi[0] = 0;
2719 else
2720 bi[0] = 0Xffff;
2721 emovo (bi, c);
2722}
2723
2724
2725
defb5dab 2726/* Multiply. */
a0353055
RK		 2727
2728static void
985b6196
RS		 2729emul (a, b, c)
2730 unsigned EMUSHORT *a, *b, *c;
2731{
2732 unsigned EMUSHORT ai[NI], bi[NI];
2733 int i, j;
2734 EMULONG lt, lta, ltb;
2735
66b6d60b
RS		 2736#ifdef NANS
2737/* NaN times anything is the same NaN. */
2738 if (eisnan (a))
2739 {
2740 emov (a, c);
2741 return;
2742 }
2743 if (eisnan (b))
2744 {
2745 emov (b, c);
2746 return;
2747 }
2748/* Zero times infinity is a NaN. */
2749 if ((eisinf (a) && (ecmp (b, ezero) == 0))
2750 || (eisinf (b) && (ecmp (a, ezero) == 0)))
2751 {
2752 mtherr ("emul", INVALID);
29e11dab 2753 enan (c, eisneg (a) ^ eisneg (b));
66b6d60b
RS		 2754 return;
2755 }
2756#endif
2757/* Infinity times anything else is infinity. */
985b6196
RS		 2758#ifdef INFINITY
2759 if (eisinf (a) || eisinf (b))
2760 {
2761 if (eisneg (a) ^ eisneg (b))
2762 *(c + (NE - 1)) = 0x8000;
2763 else
2764 *(c + (NE - 1)) = 0;
2765 einfin (c);
2766 return;
2767 }
2768#endif
2769 emovi (a, ai);
2770 emovi (b, bi);
2771 lta = ai[E];
2772 ltb = bi[E];
2773 if (ai[E] == 0)
2774 {
2775 for (i = 1; i < NI - 1; i++)
2776 {
2777 if (ai[i] != 0)
2778 {
2779 lta -= enormlz (ai);
2780 goto mnzer1;
2781 }
2782 }
2783 eclear (c);
2784 return;
2785 }
2786 mnzer1:
2787
2788 if (bi[E] == 0)
2789 {
2790 for (i = 1; i < NI - 1; i++)
2791 {
2792 if (bi[i] != 0)
2793 {
2794 ltb -= enormlz (bi);
2795 goto mnzer2;
2796 }
2797 }
2798 eclear (c);
2799 return;
2800 }
2801 mnzer2:
2802
2803 /* Multiply significands */
2804 j = emulm (ai, bi);
2805 /* calculate exponent */
2806 lt = lta + ltb - (EXONE - 1);
2807 emdnorm (bi, j, 0, lt, 64);
2808 /* calculate sign of product */
2809 if (ai[0] == bi[0])
2810 bi[0] = 0;
2811 else
2812 bi[0] = 0xffff;
2813 emovo (bi, c);
2814}
2815
2816
2817
2818
defb5dab 2819/* Convert IEEE double precision to e type. */
a0353055
RK		 2820
2821static void
66b6d60b
RS		 2822e53toe (pe, y)
2823 unsigned EMUSHORT *pe, *y;
985b6196
RS		 2824{
2825#ifdef DEC
2826
66b6d60b 2827 dectoe (pe, y); /* see etodec.c */
985b6196
RS		 2828
2829#else
842fbaaa
JW		 2830#ifdef IBM
2831
2832 ibmtoe (pe, y, DFmode);
985b6196 2833
842fbaaa 2834#else
985b6196 2835 register unsigned EMUSHORT r;
66b6d60b 2836 register unsigned EMUSHORT *e, *p;
985b6196
RS		 2837 unsigned EMUSHORT yy[NI];
2838 int denorm, k;
2839
66b6d60b 2840 e = pe;
985b6196
RS		 2841 denorm = 0; /* flag if denormalized number */
2842 ecleaz (yy);
f76b9db2
ILT		 2843 if (! FLOAT_WORDS_BIG_ENDIAN)
2844 e += 3;
985b6196
RS		 2845 r = *e;
2846 yy[0] = 0;
2847 if (r & 0x8000)
2848 yy[0] = 0xffff;
2849 yy[M] = (r & 0x0f) | 0x10;
2850 r &= ~0x800f; /* strip sign and 4 significand bits */
2851#ifdef INFINITY
2852 if (r == 0x7ff0)
2853 {
66b6d60b 2854#ifdef NANS
f76b9db2 2855 if (! FLOAT_WORDS_BIG_ENDIAN)
66b6d60b 2856 {
f76b9db2
ILT		 2857 if (((pe[3] & 0xf) != 0) || (pe[2] != 0)
2858 || (pe[1] != 0) || (pe[0] != 0))
2859 {
2860 enan (y, yy[0] != 0);
2861 return;
2862 }
66b6d60b 2863 }
f76b9db2 2864 else
66b6d60b 2865 {
f76b9db2
ILT		 2866 if (((pe[0] & 0xf) != 0) || (pe[1] != 0)
2867 || (pe[2] != 0) || (pe[3] != 0))
2868 {
2869 enan (y, yy[0] != 0);
2870 return;
2871 }
66b6d60b 2872 }
66b6d60b 2873#endif /* NANS */
dca821e1 2874 eclear (y);
985b6196 2875 einfin (y);
dca821e1 2876 if (yy[0])
985b6196
RS		 2877 eneg (y);
2878 return;
2879 }
66b6d60b 2880#endif /* INFINITY */
985b6196
RS		 2881 r >>= 4;
2882 /* If zero exponent, then the significand is denormalized.
defb5dab
RK		 2883 So take back the understood high significand bit. */
2884
985b6196
RS		 2885 if (r == 0)
2886 {
2887 denorm = 1;
2888 yy[M] &= ~0x10;
2889 }
2890 r += EXONE - 01777;
2891 yy[E] = r;
2892 p = &yy[M + 1];
f76b9db2
ILT		 2893#ifdef IEEE
2894 if (! FLOAT_WORDS_BIG_ENDIAN)
2895 {
2896 *p++ = *(--e);
2897 *p++ = *(--e);
2898 *p++ = *(--e);
2899 }
2900 else
2901 {
2902 ++e;
2903 *p++ = *e++;
2904 *p++ = *e++;
2905 *p++ = *e++;
2906 }
985b6196 2907#endif
64685ffa 2908 eshift (yy, -5);
985b6196
RS		 2909 if (denorm)
2910 { /* if zero exponent, then normalize the significand */
2911 if ((k = enormlz (yy)) > NBITS)
2912 ecleazs (yy);
2913 else
2914 yy[E] -= (unsigned EMUSHORT) (k - 1);
2915 }
2916 emovo (yy, y);
842fbaaa 2917#endif /* not IBM */
985b6196
RS		 2918#endif /* not DEC */
2919}
2920
a0353055 2921static void
66b6d60b
RS		 2922e64toe (pe, y)
2923 unsigned EMUSHORT *pe, *y;
985b6196
RS		 2924{
2925 unsigned EMUSHORT yy[NI];
66b6d60b 2926 unsigned EMUSHORT *e, *p, *q;
985b6196
RS		 2927 int i;
2928
66b6d60b 2929 e = pe;
985b6196
RS		 2930 p = yy;
2931 for (i = 0; i < NE - 5; i++)
2932 *p++ = 0;
842fbaaa 2933/* This precision is not ordinarily supported on DEC or IBM. */
985b6196
RS		 2934#ifdef DEC
2935 for (i = 0; i < 5; i++)
2936 *p++ = *e++;
2937#endif
842fbaaa
JW		 2938#ifdef IBM
2939 p = &yy[0] + (NE - 1);
2940 *p-- = *e++;
2941 ++e;
2942 for (i = 0; i < 5; i++)
2943 *p-- = *e++;
2944#endif
f76b9db2
ILT		 2945#ifdef IEEE
2946 if (! FLOAT_WORDS_BIG_ENDIAN)
2947 {
2948 for (i = 0; i < 5; i++)
2949 *p++ = *e++;
2950 }
2951 else
2952 {
2953 p = &yy[0] + (NE - 1);
2954 *p-- = *e++;
2955 ++e;
2956 for (i = 0; i < 4; i++)
2957 *p-- = *e++;
2958 }
985b6196 2959#endif
985b6196 2960#ifdef INFINITY
82e974d4
RK		 2961 /* Point to the exponent field and check max exponent cases. */
2962 p = &yy[NE - 1];
985b6196
RS		 2963 if (*p == 0x7fff)
2964 {
66b6d60b 2965#ifdef NANS
f76b9db2 2966 if (! FLOAT_WORDS_BIG_ENDIAN)
66b6d60b 2967 {
f76b9db2 2968 for (i = 0; i < 4; i++)
66b6d60b 2969 {
82e974d4
RK		 2970 if ((i != 3 && pe[i] != 0)
2971 /* Anything but 0x8000 here, including 0, is a NaN. */
2972 || (i == 3 && pe[i] != 0x8000))
f76b9db2
ILT		 2973 {
2974 enan (y, (*p & 0x8000) != 0);
2975 return;
2976 }
66b6d60b
RS		 2977 }
2978 }
f76b9db2 2979 else
66b6d60b 2980 {
f76b9db2 2981 for (i = 1; i <= 4; i++)
66b6d60b 2982 {
f76b9db2
ILT		 2983 if (pe[i] != 0)
2984 {
2985 enan (y, (*p & 0x8000) != 0);
2986 return;
2987 }
66b6d60b
RS		 2988 }
2989 }
66b6d60b 2990#endif /* NANS */
dca821e1 2991 eclear (y);
985b6196
RS		 2992 einfin (y);
2993 if (*p & 0x8000)
2994 eneg (y);
2995 return;
2996 }
66b6d60b 2997#endif /* INFINITY */
82e974d4
RK		 2998 p = yy;
2999 q = y;
985b6196
RS		 3000 for (i = 0; i < NE; i++)
3001 *q++ = *p++;
3002}
3003
3004
a0353055 3005static void
842fbaaa 3006e113toe (pe, y)
66b6d60b 3007 unsigned EMUSHORT *pe, *y;
985b6196
RS		 3008{
3009 register unsigned EMUSHORT r;
842fbaaa 3010 unsigned EMUSHORT *e, *p;
985b6196 3011 unsigned EMUSHORT yy[NI];
842fbaaa 3012 int denorm, i;
985b6196 3013
66b6d60b 3014 e = pe;
842fbaaa 3015 denorm = 0;
985b6196 3016 ecleaz (yy);
f76b9db2
ILT		 3017#ifdef IEEE
3018 if (! FLOAT_WORDS_BIG_ENDIAN)
3019 e += 7;
985b6196
RS		 3020#endif
3021 r = *e;
3022 yy[0] = 0;
3023 if (r & 0x8000)
3024 yy[0] = 0xffff;
842fbaaa 3025 r &= 0x7fff;
985b6196 3026#ifdef INFINITY
842fbaaa 3027 if (r == 0x7fff)
985b6196 3028 {
66b6d60b 3029#ifdef NANS
f76b9db2 3030 if (! FLOAT_WORDS_BIG_ENDIAN)
66b6d60b 3031 {
f76b9db2 3032 for (i = 0; i < 7; i++)
842fbaaa 3033 {
f76b9db2
ILT		 3034 if (pe[i] != 0)
3035 {
3036 enan (y, yy[0] != 0);
3037 return;
3038 }
842fbaaa 3039 }
66b6d60b 3040 }
f76b9db2 3041 else
66b6d60b 3042 {
f76b9db2 3043 for (i = 1; i < 8; i++)
842fbaaa 3044 {
f76b9db2
ILT		 3045 if (pe[i] != 0)
3046 {
3047 enan (y, yy[0] != 0);
3048 return;
3049 }
842fbaaa 3050 }
66b6d60b 3051 }
842fbaaa 3052#endif /* NANS */
dca821e1 3053 eclear (y);
985b6196 3054 einfin (y);
dca821e1 3055 if (yy[0])
985b6196
RS		 3056 eneg (y);
3057 return;
3058 }
66b6d60b 3059#endif /* INFINITY */
985b6196
RS		 3060 yy[E] = r;
3061 p = &yy[M + 1];
f76b9db2
ILT		 3062#ifdef IEEE
3063 if (! FLOAT_WORDS_BIG_ENDIAN)
3064 {
3065 for (i = 0; i < 7; i++)
3066 *p++ = *(--e);
3067 }
3068 else
3069 {
3070 ++e;
3071 for (i = 0; i < 7; i++)
3072 *p++ = *e++;
3073 }
985b6196 3074#endif
842fbaaa
JW		 3075/* If denormal, remove the implied bit; else shift down 1. */
3076 if (r == 0)
3077 {
3078 yy[M] = 0;
3079 }
3080 else
3081 {
3082 yy[M] = 1;
3083 eshift (yy, -1);
3084 }
3085 emovo (yy, y);
3086}
3087
3088
defb5dab 3089/* Convert IEEE single precision to e type. */
a0353055
RK		 3090
3091static void
842fbaaa
JW		 3092e24toe (pe, y)
3093 unsigned EMUSHORT *pe, *y;
3094{
3095#ifdef IBM
3096
3097 ibmtoe (pe, y, SFmode);
3098
3099#else
3100 register unsigned EMUSHORT r;
3101 register unsigned EMUSHORT *e, *p;
3102 unsigned EMUSHORT yy[NI];
3103 int denorm, k;
3104
3105 e = pe;
3106 denorm = 0; /* flag if denormalized number */
3107 ecleaz (yy);
f76b9db2
ILT		 3108#ifdef IEEE
3109 if (! FLOAT_WORDS_BIG_ENDIAN)
3110 e += 1;
842fbaaa
JW		 3111#endif
3112#ifdef DEC
3113 e += 1;
3114#endif
3115 r = *e;
3116 yy[0] = 0;
3117 if (r & 0x8000)
3118 yy[0] = 0xffff;
3119 yy[M] = (r & 0x7f) | 0200;
3120 r &= ~0x807f; /* strip sign and 7 significand bits */
3121#ifdef INFINITY
3122 if (r == 0x7f80)
3123 {
3124#ifdef NANS
f76b9db2 3125 if (FLOAT_WORDS_BIG_ENDIAN)
842fbaaa 3126 {
f76b9db2
ILT		 3127 if (((pe[0] & 0x7f) != 0) || (pe[1] != 0))
3128 {
3129 enan (y, yy[0] != 0);
3130 return;
3131 }
842fbaaa 3132 }
f76b9db2 3133 else
842fbaaa 3134 {
f76b9db2
ILT		 3135 if (((pe[1] & 0x7f) != 0) || (pe[0] != 0))
3136 {
3137 enan (y, yy[0] != 0);
3138 return;
3139 }
842fbaaa 3140 }
842fbaaa
JW		 3141#endif /* NANS */
3142 eclear (y);
3143 einfin (y);
3144 if (yy[0])
3145 eneg (y);
3146 return;
3147 }
3148#endif /* INFINITY */
3149 r >>= 7;
3150 /* If zero exponent, then the significand is denormalized.
defb5dab 3151 So take back the understood high significand bit. */
842fbaaa
JW		 3152 if (r == 0)
3153 {
3154 denorm = 1;
3155 yy[M] &= ~0200;
3156 }
3157 r += EXONE - 0177;
3158 yy[E] = r;
3159 p = &yy[M + 1];
842fbaaa
JW		 3160#ifdef DEC
3161 *p++ = *(--e);
3162#endif
f76b9db2
ILT		 3163#ifdef IEEE
3164 if (! FLOAT_WORDS_BIG_ENDIAN)
3165 *p++ = *(--e);
3166 else
3167 {
3168 ++e;
3169 *p++ = *e++;
3170 }
842fbaaa
JW		 3171#endif
3172 eshift (yy, -8);
3173 if (denorm)
3174 { /* if zero exponent, then normalize the significand */
3175 if ((k = enormlz (yy)) > NBITS)
3176 ecleazs (yy);
3177 else
3178 yy[E] -= (unsigned EMUSHORT) (k - 1);
985b6196
RS		 3179 }
3180 emovo (yy, y);
842fbaaa
JW		 3181#endif /* not IBM */
3182}
3183
3184
a0353055 3185static void
842fbaaa
JW		 3186etoe113 (x, e)
3187 unsigned EMUSHORT *x, *e;
3188{
3189 unsigned EMUSHORT xi[NI];
3190 EMULONG exp;
3191 int rndsav;
3192
3193#ifdef NANS
3194 if (eisnan (x))
3195 {
29e11dab 3196 make_nan (e, eisneg (x), TFmode);
842fbaaa
JW		 3197 return;
3198 }
3199#endif
3200 emovi (x, xi);
3201 exp = (EMULONG) xi[E];
3202#ifdef INFINITY
3203 if (eisinf (x))
3204 goto nonorm;
3205#endif
3206 /* round off to nearest or even */
3207 rndsav = rndprc;
3208 rndprc = 113;
3209 emdnorm (xi, 0, 0, exp, 64);
3210 rndprc = rndsav;
3211 nonorm:
3212 toe113 (xi, e);
985b6196
RS		 3213}
3214
defb5dab 3215/* Move out internal format to ieee long double */
a0353055 3216
842fbaaa
JW		 3217static void
3218toe113 (a, b)
3219 unsigned EMUSHORT *a, *b;
3220{
3221 register unsigned EMUSHORT *p, *q;
3222 unsigned EMUSHORT i;
3223
3224#ifdef NANS
3225 if (eiisnan (a))
3226 {
29e11dab 3227 make_nan (b, eiisneg (a), TFmode);
842fbaaa
JW		 3228 return;
3229 }
3230#endif
3231 p = a;
f76b9db2
ILT		 3232 if (FLOAT_WORDS_BIG_ENDIAN)
3233 q = b;
3234 else
3235 q = b + 7; /* point to output exponent */
842fbaaa
JW		 3236
3237 /* If not denormal, delete the implied bit. */
3238 if (a[E] != 0)
3239 {
3240 eshup1 (a);
3241 }
3242 /* combine sign and exponent */
3243 i = *p++;
f76b9db2
ILT		 3244 if (FLOAT_WORDS_BIG_ENDIAN)
3245 {
3246 if (i)
3247 *q++ = *p++ | 0x8000;
3248 else
3249 *q++ = *p++;
3250 }
842fbaaa 3251 else
f76b9db2
ILT		 3252 {
3253 if (i)
3254 *q-- = *p++ | 0x8000;
3255 else
3256 *q-- = *p++;
3257 }
842fbaaa
JW		 3258 /* skip over guard word */
3259 ++p;
3260 /* move the significand */
f76b9db2
ILT		 3261 if (FLOAT_WORDS_BIG_ENDIAN)
3262 {
3263 for (i = 0; i < 7; i++)
3264 *q++ = *p++;
3265 }
3266 else
3267 {
3268 for (i = 0; i < 7; i++)
3269 *q-- = *p++;
3270 }
842fbaaa 3271}
985b6196 3272
a0353055 3273static void
985b6196
RS		 3274etoe64 (x, e)
3275 unsigned EMUSHORT *x, *e;
3276{
3277 unsigned EMUSHORT xi[NI];
3278 EMULONG exp;
3279 int rndsav;
3280
66b6d60b
RS		 3281#ifdef NANS
3282 if (eisnan (x))
3283 {
29e11dab 3284 make_nan (e, eisneg (x), XFmode);
66b6d60b
RS		 3285 return;
3286 }
3287#endif
985b6196
RS		 3288 emovi (x, xi);
3289 /* adjust exponent for offset */
3290 exp = (EMULONG) xi[E];
3291#ifdef INFINITY
3292 if (eisinf (x))
3293 goto nonorm;
3294#endif
3295 /* round off to nearest or even */
3296 rndsav = rndprc;
3297 rndprc = 64;
3298 emdnorm (xi, 0, 0, exp, 64);
3299 rndprc = rndsav;
3300 nonorm:
3301 toe64 (xi, e);
3302}
3303
a0353055 3304
defb5dab
RK		 3305/* Move out internal format to ieee long double. */
3306
985b6196
RS		 3307static void
3308toe64 (a, b)
3309 unsigned EMUSHORT *a, *b;
3310{
3311 register unsigned EMUSHORT *p, *q;
3312 unsigned EMUSHORT i;
3313
66b6d60b
RS		 3314#ifdef NANS
3315 if (eiisnan (a))
3316 {
29e11dab 3317 make_nan (b, eiisneg (a), XFmode);
66b6d60b
RS		 3318 return;
3319 }
3320#endif
985b6196 3321 p = a;
f76b9db2 3322#ifdef IBM
985b6196 3323 q = b;
f76b9db2
ILT		 3324#endif
3325#ifdef DEC
3326 q = b + 4;
3327#endif
3328#ifdef IEEE
3329 if (FLOAT_WORDS_BIG_ENDIAN)
3330 q = b;
3331 else
3332 {
3333 q = b + 4; /* point to output exponent */
985b6196 3334#if LONG_DOUBLE_TYPE_SIZE == 96
f76b9db2
ILT		 3335 /* Clear the last two bytes of 12-byte Intel format */
3336 *(q+1) = 0;
985b6196 3337#endif
f76b9db2 3338 }
985b6196
RS		 3339#endif
3340
3341 /* combine sign and exponent */
3342 i = *p++;
f76b9db2 3343#ifdef IBM
985b6196
RS		 3344 if (i)
3345 *q++ = *p++ | 0x8000;
3346 else
3347 *q++ = *p++;
3348 *q++ = 0;
f76b9db2
ILT		 3349#endif
3350#ifdef DEC
985b6196
RS		 3351 if (i)
3352 *q-- = *p++ | 0x8000;
3353 else
3354 *q-- = *p++;
f76b9db2
ILT		 3355#endif
3356#ifdef IEEE
3357 if (FLOAT_WORDS_BIG_ENDIAN)
3358 {
3359 if (i)
3360 *q++ = *p++ | 0x8000;
3361 else
3362 *q++ = *p++;
3363 *q++ = 0;
3364 }
3365 else
3366 {
3367 if (i)
3368 *q-- = *p++ | 0x8000;
3369 else
3370 *q-- = *p++;
3371 }
985b6196
RS		 3372#endif
3373 /* skip over guard word */
3374 ++p;
3375 /* move the significand */
f76b9db2 3376#ifdef IBM
985b6196
RS		 3377 for (i = 0; i < 4; i++)
3378 *q++ = *p++;
f76b9db2
ILT		 3379#endif
3380#ifdef DEC
985b6196
RS		 3381 for (i = 0; i < 4; i++)
3382 *q-- = *p++;
3383#endif
f76b9db2
ILT		 3384#ifdef IEEE
3385 if (FLOAT_WORDS_BIG_ENDIAN)
3386 {
3387 for (i = 0; i < 4; i++)
3388 *q++ = *p++;
3389 }
3390 else
3391 {
82e974d4
RK		 3392#ifdef INFINITY
3393 if (eiisinf (a))
3394 {
3395 /* Intel long double infinity significand. */
3396 *q-- = 0x8000;
3397 *q-- = 0;
3398 *q-- = 0;
3399 *q = 0;
3400 return;
3401 }
3402#endif
f76b9db2
ILT		 3403 for (i = 0; i < 4; i++)
3404 *q-- = *p++;
3405 }
3406#endif
985b6196
RS		 3407}
3408
3409
defb5dab 3410/* e type to IEEE double precision. */
985b6196
RS		 3411
3412#ifdef DEC
3413
a0353055 3414static void
985b6196
RS		 3415etoe53 (x, e)
3416 unsigned EMUSHORT *x, *e;
3417{
3418 etodec (x, e); /* see etodec.c */
3419}
3420
3421static void
3422toe53 (x, y)
3423 unsigned EMUSHORT *x, *y;
3424{
3425 todec (x, y);
3426}
3427
3428#else
842fbaaa
JW		 3429#ifdef IBM
3430
008f0d36 3431static void
842fbaaa
JW		 3432etoe53 (x, e)
3433 unsigned EMUSHORT *x, *e;
3434{
3435 etoibm (x, e, DFmode);
3436}
3437
3438static void
3439toe53 (x, y)
3440 unsigned EMUSHORT *x, *y;
3441{
3442 toibm (x, y, DFmode);
3443}
3444
3445#else /* it's neither DEC nor IBM */
985b6196 3446
008f0d36 3447static void
985b6196
RS		 3448etoe53 (x, e)
3449 unsigned EMUSHORT *x, *e;
3450{
3451 unsigned EMUSHORT xi[NI];
3452 EMULONG exp;
3453 int rndsav;
3454
66b6d60b
RS		 3455#ifdef NANS
3456 if (eisnan (x))
3457 {
29e11dab 3458 make_nan (e, eisneg (x), DFmode);
66b6d60b
RS		 3459 return;
3460 }
3461#endif
985b6196
RS		 3462 emovi (x, xi);
3463 /* adjust exponent for offsets */
3464 exp = (EMULONG) xi[E] - (EXONE - 0x3ff);
3465#ifdef INFINITY
3466 if (eisinf (x))
3467 goto nonorm;
3468#endif
3469 /* round off to nearest or even */
3470 rndsav = rndprc;
3471 rndprc = 53;
3472 emdnorm (xi, 0, 0, exp, 64);
3473 rndprc = rndsav;
3474 nonorm:
3475 toe53 (xi, e);
3476}
3477
3478
3479static void
3480toe53 (x, y)
3481 unsigned EMUSHORT *x, *y;
3482{
3483 unsigned EMUSHORT i;
3484 unsigned EMUSHORT *p;
3485
66b6d60b
RS		 3486#ifdef NANS
3487 if (eiisnan (x))
3488 {
29e11dab 3489 make_nan (y, eiisneg (x), DFmode);
66b6d60b
RS		 3490 return;
3491 }
3492#endif
985b6196 3493 p = &x[0];
f76b9db2
ILT		 3494#ifdef IEEE
3495 if (! FLOAT_WORDS_BIG_ENDIAN)
3496 y += 3;
985b6196
RS		 3497#endif
3498 *y = 0; /* output high order */
3499 if (*p++)
3500 *y = 0x8000; /* output sign bit */
3501
3502 i = *p++;
3503 if (i >= (unsigned int) 2047)
3504 { /* Saturate at largest number less than infinity. */
3505#ifdef INFINITY
3506 *y |= 0x7ff0;
f76b9db2
ILT		 3507 if (! FLOAT_WORDS_BIG_ENDIAN)
3508 {
3509 *(--y) = 0;
3510 *(--y) = 0;
3511 *(--y) = 0;
3512 }
3513 else
3514 {
3515 ++y;
3516 *y++ = 0;
3517 *y++ = 0;
3518 *y++ = 0;
3519 }
985b6196
RS		 3520#else
3521 *y |= (unsigned EMUSHORT) 0x7fef;
f76b9db2
ILT		 3522 if (! FLOAT_WORDS_BIG_ENDIAN)
3523 {
3524 *(--y) = 0xffff;
3525 *(--y) = 0xffff;
3526 *(--y) = 0xffff;
3527 }
3528 else
3529 {
3530 ++y;
3531 *y++ = 0xffff;
3532 *y++ = 0xffff;
3533 *y++ = 0xffff;
3534 }
985b6196
RS		 3535#endif
3536 return;
3537 }
3538 if (i == 0)
3539 {
64685ffa 3540 eshift (x, 4);
985b6196
RS		 3541 }
3542 else
3543 {
3544 i <<= 4;
64685ffa 3545 eshift (x, 5);
985b6196
RS		 3546 }
3547 i |= *p++ & (unsigned EMUSHORT) 0x0f; /* *p = xi[M] */
3548 *y |= (unsigned EMUSHORT) i; /* high order output already has sign bit set */
f76b9db2
ILT		 3549 if (! FLOAT_WORDS_BIG_ENDIAN)
3550 {
3551 *(--y) = *p++;
3552 *(--y) = *p++;
3553 *(--y) = *p;
3554 }
3555 else
3556 {
3557 ++y;
3558 *y++ = *p++;
3559 *y++ = *p++;
3560 *y++ = *p++;
3561 }
985b6196
RS		 3562}
3563
842fbaaa 3564#endif /* not IBM */
985b6196
RS		 3565#endif /* not DEC */
3566
3567
3568
defb5dab
RK		 3569/* e type to IEEE single precision. */
3570
842fbaaa
JW		 3571#ifdef IBM
3572
008f0d36 3573static void
842fbaaa
JW		 3574etoe24 (x, e)
3575 unsigned EMUSHORT *x, *e;
3576{
3577 etoibm (x, e, SFmode);
3578}
3579
3580static void
3581toe24 (x, y)
3582 unsigned EMUSHORT *x, *y;
3583{
3584 toibm (x, y, SFmode);
3585}
3586
3587#else
3588
008f0d36 3589static void
985b6196
RS		 3590etoe24 (x, e)
3591 unsigned EMUSHORT *x, *e;
3592{
3593 EMULONG exp;
3594 unsigned EMUSHORT xi[NI];
3595 int rndsav;
3596
66b6d60b
RS		 3597#ifdef NANS
3598 if (eisnan (x))
3599 {
29e11dab 3600 make_nan (e, eisneg (x), SFmode);
66b6d60b
RS		 3601 return;
3602 }
3603#endif
985b6196
RS		 3604 emovi (x, xi);
3605 /* adjust exponent for offsets */
3606 exp = (EMULONG) xi[E] - (EXONE - 0177);
3607#ifdef INFINITY
3608 if (eisinf (x))
3609 goto nonorm;
3610#endif
3611 /* round off to nearest or even */
3612 rndsav = rndprc;
3613 rndprc = 24;
3614 emdnorm (xi, 0, 0, exp, 64);
3615 rndprc = rndsav;
3616 nonorm:
3617 toe24 (xi, e);
3618}
3619
3620static void
3621toe24 (x, y)
3622 unsigned EMUSHORT *x, *y;
3623{
3624 unsigned EMUSHORT i;
3625 unsigned EMUSHORT *p;
3626
66b6d60b
RS		 3627#ifdef NANS
3628 if (eiisnan (x))
3629 {
29e11dab 3630 make_nan (y, eiisneg (x), SFmode);
66b6d60b
RS		 3631 return;
3632 }
3633#endif
985b6196 3634 p = &x[0];
f76b9db2
ILT		 3635#ifdef IEEE
3636 if (! FLOAT_WORDS_BIG_ENDIAN)
3637 y += 1;
985b6196
RS		 3638#endif
3639#ifdef DEC
3640 y += 1;
3641#endif
3642 *y = 0; /* output high order */
3643 if (*p++)
3644 *y = 0x8000; /* output sign bit */
3645
3646 i = *p++;
64685ffa 3647/* Handle overflow cases. */
985b6196 3648 if (i >= 255)
64685ffa 3649 {
985b6196
RS		 3650#ifdef INFINITY
3651 *y |= (unsigned EMUSHORT) 0x7f80;
985b6196
RS		 3652#ifdef DEC
3653 *(--y) = 0;
3654#endif
f76b9db2
ILT		 3655#ifdef IEEE
3656 if (! FLOAT_WORDS_BIG_ENDIAN)
3657 *(--y) = 0;
3658 else
3659 {
3660 ++y;
3661 *y = 0;
3662 }
985b6196 3663#endif
64685ffa 3664#else /* no INFINITY */
985b6196 3665 *y |= (unsigned EMUSHORT) 0x7f7f;
985b6196
RS		 3666#ifdef DEC
3667 *(--y) = 0xffff;
3668#endif
f76b9db2
ILT		 3669#ifdef IEEE
3670 if (! FLOAT_WORDS_BIG_ENDIAN)
3671 *(--y) = 0xffff;
3672 else
3673 {
3674 ++y;
3675 *y = 0xffff;
3676 }
985b6196 3677#endif
64685ffa
RS		 3678#ifdef ERANGE
3679 errno = ERANGE;
985b6196 3680#endif
64685ffa 3681#endif /* no INFINITY */
985b6196
RS		 3682 return;
3683 }
3684 if (i == 0)
3685 {
64685ffa 3686 eshift (x, 7);
985b6196
RS		 3687 }
3688 else
3689 {
3690 i <<= 7;
64685ffa 3691 eshift (x, 8);
985b6196
RS		 3692 }
3693 i |= *p++ & (unsigned EMUSHORT) 0x7f; /* *p = xi[M] */
3694 *y |= i; /* high order output already has sign bit set */
985b6196
RS		 3695#ifdef DEC
3696 *(--y) = *p;
3697#endif
f76b9db2
ILT		 3698#ifdef IEEE
3699 if (! FLOAT_WORDS_BIG_ENDIAN)
3700 *(--y) = *p;
3701 else
3702 {
3703 ++y;
3704 *y = *p;
3705 }
985b6196
RS		 3706#endif
3707}
842fbaaa 3708#endif /* not IBM */
985b6196 3709
defb5dab
RK		 3710/* Compare two e type numbers.
3711 Return +1 if a > b
3712 0 if a == b
3713 -1 if a < b
3714 -2 if either a or b is a NaN. */
a0353055
RK		 3715
3716static int
985b6196
RS		 3717ecmp (a, b)
3718 unsigned EMUSHORT *a, *b;
3719{
3720 unsigned EMUSHORT ai[NI], bi[NI];
3721 register unsigned EMUSHORT *p, *q;
3722 register int i;
3723 int msign;
3724
66b6d60b
RS		 3725#ifdef NANS
3726 if (eisnan (a) || eisnan (b))
3727 return (-2);
3728#endif
985b6196
RS		 3729 emovi (a, ai);
3730 p = ai;
3731 emovi (b, bi);
3732 q = bi;
3733
3734 if (*p != *q)
3735 { /* the signs are different */
3736 /* -0 equals + 0 */
3737 for (i = 1; i < NI - 1; i++)
3738 {
3739 if (ai[i] != 0)
3740 goto nzro;
3741 if (bi[i] != 0)
3742 goto nzro;
3743 }
3744 return (0);
3745 nzro:
3746 if (*p == 0)
3747 return (1);
3748 else
3749 return (-1);
3750 }
3751 /* both are the same sign */
3752 if (*p == 0)
3753 msign = 1;
3754 else
3755 msign = -1;
3756 i = NI - 1;
3757 do
3758 {
3759 if (*p++ != *q++)
3760 {
3761 goto diff;
3762 }
3763 }
3764 while (--i > 0);
3765
3766 return (0); /* equality */
3767
3768
3769
3770 diff:
3771
3772 if (*(--p) > *(--q))
3773 return (msign); /* p is bigger */
3774 else
3775 return (-msign); /* p is littler */
3776}
3777
3778
3779
3780
defb5dab 3781/* Find nearest integer to x = floor (x + 0.5). */
a0353055
RK		 3782
3783static void
985b6196
RS		 3784eround (x, y)
3785 unsigned EMUSHORT *x, *y;
3786{
3787 eadd (ehalf, x, y);
3788 efloor (y, y);
3789}
3790
3791
3792
3793
defb5dab 3794/* Convert HOST_WIDE_INT to e type. */
a0353055
RK		 3795
3796static void
985b6196 3797ltoe (lp, y)
b51ab098
RK		 3798 HOST_WIDE_INT *lp;
3799 unsigned EMUSHORT *y;
985b6196
RS		 3800{
3801 unsigned EMUSHORT yi[NI];
b51ab098 3802 unsigned HOST_WIDE_INT ll;
985b6196
RS		 3803 int k;
3804
3805 ecleaz (yi);
3806 if (*lp < 0)
3807 {
3808 /* make it positive */
b51ab098 3809 ll = (unsigned HOST_WIDE_INT) (-(*lp));
985b6196
RS		 3810 yi[0] = 0xffff; /* put correct sign in the e type number */
3811 }
3812 else
3813 {
b51ab098 3814 ll = (unsigned HOST_WIDE_INT) (*lp);
985b6196
RS		 3815 }
3816 /* move the long integer to yi significand area */
b51ab098 3817#if HOST_BITS_PER_WIDE_INT == 64
7729f1ca
RS		 3818 yi[M] = (unsigned EMUSHORT) (ll >> 48);
3819 yi[M + 1] = (unsigned EMUSHORT) (ll >> 32);
3820 yi[M + 2] = (unsigned EMUSHORT) (ll >> 16);
3821 yi[M + 3] = (unsigned EMUSHORT) ll;
3822 yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
3823#else
985b6196
RS		 3824 yi[M] = (unsigned EMUSHORT) (ll >> 16);
3825 yi[M + 1] = (unsigned EMUSHORT) ll;
985b6196 3826 yi[E] = EXONE + 15; /* exponent if normalize shift count were 0 */
7729f1ca
RS		 3827#endif
3828
985b6196
RS		 3829 if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
3830 ecleaz (yi); /* it was zero */
3831 else
3832 yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
3833 emovo (yi, y); /* output the answer */
3834}
3835
defb5dab 3836/* Convert unsigned HOST_WIDE_INT to e type. */
a0353055
RK		 3837
3838static void
985b6196 3839ultoe (lp, y)
b51ab098
RK		 3840 unsigned HOST_WIDE_INT *lp;
3841 unsigned EMUSHORT *y;
985b6196
RS		 3842{
3843 unsigned EMUSHORT yi[NI];
b51ab098 3844 unsigned HOST_WIDE_INT ll;
985b6196
RS		 3845 int k;
3846
3847 ecleaz (yi);
3848 ll = *lp;
3849
3850 /* move the long integer to ayi significand area */
b51ab098 3851#if HOST_BITS_PER_WIDE_INT == 64
7729f1ca
RS		 3852 yi[M] = (unsigned EMUSHORT) (ll >> 48);
3853 yi[M + 1] = (unsigned EMUSHORT) (ll >> 32);
3854 yi[M + 2] = (unsigned EMUSHORT) (ll >> 16);
3855 yi[M + 3] = (unsigned EMUSHORT) ll;
3856 yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
3857#else
985b6196
RS		 3858 yi[M] = (unsigned EMUSHORT) (ll >> 16);
3859 yi[M + 1] = (unsigned EMUSHORT) ll;
985b6196 3860 yi[E] = EXONE + 15; /* exponent if normalize shift count were 0 */
7729f1ca
RS		 3861#endif
3862
985b6196
RS		 3863 if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
3864 ecleaz (yi); /* it was zero */
3865 else
3866 yi[E] -= (unsigned EMUSHORT) k; /* subtract shift count from exponent */
3867 emovo (yi, y); /* output the answer */
3868}
3869
3870
c764eafd
RK		 3871/* Find signed HOST_WIDE_INT integer and floating point fractional
3872 parts of e-type (packed internal format) floating point input X.
3873 The integer output I has the sign of the input, except that
3874 positive overflow is permitted if FIXUNS_TRUNC_LIKE_FIX_TRUNC.
3875 The output e-type fraction FRAC is the positive fractional
3876 part of abs (X). */
985b6196 3877
a0353055 3878static void
985b6196
RS		 3879eifrac (x, i, frac)
3880 unsigned EMUSHORT *x;
b51ab098 3881 HOST_WIDE_INT *i;
985b6196
RS		 3882 unsigned EMUSHORT *frac;
3883{
3884 unsigned EMUSHORT xi[NI];
7729f1ca 3885 int j, k;
b51ab098 3886 unsigned HOST_WIDE_INT ll;
985b6196
RS		 3887
3888 emovi (x, xi);
3889 k = (int) xi[E] - (EXONE - 1);
3890 if (k <= 0)
3891 {
3892 /* if exponent <= 0, integer = 0 and real output is fraction */
3893 *i = 0L;
3894 emovo (xi, frac);
3895 return;
3896 }
b51ab098 3897 if (k > (HOST_BITS_PER_WIDE_INT - 1))
985b6196 3898 {
7729f1ca
RS		 3899 /* long integer overflow: output large integer
3900 and correct fraction */
985b6196 3901 if (xi[0])
b51ab098 3902 *i = ((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1);
985b6196 3903 else
c764eafd
RK		 3904 {
3905#ifdef FIXUNS_TRUNC_LIKE_FIX_TRUNC
3906 /* In this case, let it overflow and convert as if unsigned. */
3907 euifrac (x, &ll, frac);
3908 *i = (HOST_WIDE_INT) ll;
3909 return;
3910#else
3911 /* In other cases, return the largest positive integer. */
3912 *i = (((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1)) - 1;
3913#endif
3914 }
64685ffa
RS		 3915 eshift (xi, k);
3916 if (extra_warnings)
3917 warning ("overflow on truncation to integer");
985b6196 3918 }
7729f1ca 3919 else if (k > 16)
985b6196 3920 {
7729f1ca
RS		 3921 /* Shift more than 16 bits: first shift up k-16 mod 16,
3922 then shift up by 16's. */
3923 j = k - ((k >> 4) << 4);
3924 eshift (xi, j);
3925 ll = xi[M];
3926 k -= j;
3927 do
3928 {
3929 eshup6 (xi);
3930 ll = (ll << 16) | xi[M];
3931 }
3932 while ((k -= 16) > 0);
3933 *i = ll;
3934 if (xi[0])
3935 *i = -(*i);
3936 }
3937 else
842fbaaa
JW		 3938 {
3939 /* shift not more than 16 bits */
3940 eshift (xi, k);
b51ab098 3941 *i = (HOST_WIDE_INT) xi[M] & 0xffff;
842fbaaa
JW		 3942 if (xi[0])
3943 *i = -(*i);
3944 }
985b6196
RS		 3945 xi[0] = 0;
3946 xi[E] = EXONE - 1;
3947 xi[M] = 0;
3948 if ((k = enormlz (xi)) > NBITS)
3949 ecleaz (xi);
3950 else
3951 xi[E] -= (unsigned EMUSHORT) k;
3952
3953 emovo (xi, frac);
3954}
3955
3956
b51ab098 3957/* Find unsigned HOST_WIDE_INT integer and floating point fractional parts.
7729f1ca
RS		 3958 A negative e type input yields integer output = 0
3959 but correct fraction. */
985b6196 3960
a0353055 3961static void
985b6196
RS		 3962euifrac (x, i, frac)
3963 unsigned EMUSHORT *x;
b51ab098 3964 unsigned HOST_WIDE_INT *i;
985b6196
RS		 3965 unsigned EMUSHORT *frac;
3966{
b51ab098 3967 unsigned HOST_WIDE_INT ll;
985b6196 3968 unsigned EMUSHORT xi[NI];
7729f1ca 3969 int j, k;
985b6196
RS		 3970
3971 emovi (x, xi);
3972 k = (int) xi[E] - (EXONE - 1);
3973 if (k <= 0)
3974 {
3975 /* if exponent <= 0, integer = 0 and argument is fraction */
3976 *i = 0L;
3977 emovo (xi, frac);
3978 return;
3979 }
b51ab098 3980 if (k > HOST_BITS_PER_WIDE_INT)
985b6196 3981 {
7729f1ca
RS		 3982 /* Long integer overflow: output large integer
3983 and correct fraction.
3984 Note, the BSD microvax compiler says that ~(0UL)
3985 is a syntax error. */
985b6196 3986 *i = ~(0L);
64685ffa
RS		 3987 eshift (xi, k);
3988 if (extra_warnings)
3989 warning ("overflow on truncation to unsigned integer");
985b6196 3990 }
7729f1ca 3991 else if (k > 16)
985b6196 3992 {
7729f1ca
RS		 3993 /* Shift more than 16 bits: first shift up k-16 mod 16,
3994 then shift up by 16's. */
3995 j = k - ((k >> 4) << 4);
3996 eshift (xi, j);
3997 ll = xi[M];
3998 k -= j;
3999 do
4000 {
4001 eshup6 (xi);
4002 ll = (ll << 16) | xi[M];
4003 }
4004 while ((k -= 16) > 0);
4005 *i = ll;
4006 }
4007 else
4008 {
4009 /* shift not more than 16 bits */
64685ffa 4010 eshift (xi, k);
b51ab098 4011 *i = (HOST_WIDE_INT) xi[M] & 0xffff;
985b6196
RS		 4012 }
4013
7729f1ca 4014 if (xi[0]) /* A negative value yields unsigned integer 0. */
985b6196 4015 *i = 0L;
842fbaaa 4016
985b6196
RS		 4017 xi[0] = 0;
4018 xi[E] = EXONE - 1;
4019 xi[M] = 0;
4020 if ((k = enormlz (xi)) > NBITS)
4021 ecleaz (xi);
4022 else
4023 xi[E] -= (unsigned EMUSHORT) k;
4024
4025 emovo (xi, frac);
4026}
4027
4028
4029
defb5dab 4030/* Shift significand area up or down by the number of bits given by SC. */
a0353055
RK		 4031
4032static int
985b6196
RS		 4033eshift (x, sc)
4034 unsigned EMUSHORT *x;
4035 int sc;
4036{
4037 unsigned EMUSHORT lost;
4038 unsigned EMUSHORT *p;
4039
4040 if (sc == 0)
4041 return (0);
4042
4043 lost = 0;
4044 p = x + NI - 1;
4045
4046 if (sc < 0)
4047 {
4048 sc = -sc;
4049 while (sc >= 16)
4050 {
4051 lost |= *p; /* remember lost bits */
4052 eshdn6 (x);
4053 sc -= 16;
4054 }
4055
4056 while (sc >= 8)
4057 {
4058 lost |= *p & 0xff;
4059 eshdn8 (x);
4060 sc -= 8;
4061 }
4062
4063 while (sc > 0)
4064 {
4065 lost |= *p & 1;
4066 eshdn1 (x);
4067 sc -= 1;
4068 }
4069 }
4070 else
4071 {
4072 while (sc >= 16)
4073 {
4074 eshup6 (x);
4075 sc -= 16;
4076 }
4077
4078 while (sc >= 8)
4079 {
4080 eshup8 (x);
4081 sc -= 8;
4082 }
4083
4084 while (sc > 0)
4085 {
4086 eshup1 (x);
4087 sc -= 1;
4088 }
4089 }
4090 if (lost)
4091 lost = 1;
4092 return ((int) lost);
4093}
4094
4095
4096
defb5dab
RK		 4097/* Shift normalize the significand area pointed to by argument.
4098 Shift count (up = positive) is returned. */
a0353055
RK		 4099
4100static int
985b6196
RS		 4101enormlz (x)
4102 unsigned EMUSHORT x[];
4103{
4104 register unsigned EMUSHORT *p;
4105 int sc;
4106
4107 sc = 0;
4108 p = &x[M];
4109 if (*p != 0)
4110 goto normdn;
4111 ++p;
4112 if (*p & 0x8000)
4113 return (0); /* already normalized */
4114 while (*p == 0)
4115 {
4116 eshup6 (x);
4117 sc += 16;
defb5dab 4118
985b6196 4119 /* With guard word, there are NBITS+16 bits available.
defb5dab 4120 Return true if all are zero. */
985b6196
RS		 4121 if (sc > NBITS)
4122 return (sc);
4123 }
4124 /* see if high byte is zero */
4125 while ((*p & 0xff00) == 0)
4126 {
4127 eshup8 (x);
4128 sc += 8;
4129 }
4130 /* now shift 1 bit at a time */
4131 while ((*p & 0x8000) == 0)
4132 {
4133 eshup1 (x);
4134 sc += 1;
4135 if (sc > NBITS)
4136 {
4137 mtherr ("enormlz", UNDERFLOW);
4138 return (sc);
4139 }
4140 }
4141 return (sc);
4142
4143 /* Normalize by shifting down out of the high guard word
4144 of the significand */
4145 normdn:
4146
4147 if (*p & 0xff00)
4148 {
4149 eshdn8 (x);
4150 sc -= 8;
4151 }
4152 while (*p != 0)
4153 {
4154 eshdn1 (x);
4155 sc -= 1;
4156
4157 if (sc < -NBITS)
4158 {
4159 mtherr ("enormlz", OVERFLOW);
4160 return (sc);
4161 }
4162 }
4163 return (sc);
4164}
4165
4166
4167
4168
4169/* Convert e type number to decimal format ASCII string.
defb5dab 4170 The constants are for 64 bit precision. */
985b6196
RS		 4171
4172#define NTEN 12
4173#define MAXP 4096
4174
842fbaaa
JW		 4175#if LONG_DOUBLE_TYPE_SIZE == 128
4176static unsigned EMUSHORT etens[NTEN + 1][NE] =
4177{
4178 {0x6576, 0x4a92, 0x804a, 0x153f,
4179 0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,}, /* 10**4096 */
4180 {0x6a32, 0xce52, 0x329a, 0x28ce,
4181 0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,}, /* 10**2048 */
4182 {0x526c, 0x50ce, 0xf18b, 0x3d28,
4183 0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
4184 {0x9c66, 0x58f8, 0xbc50, 0x5c54,
4185 0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
4186 {0x851e, 0xeab7, 0x98fe, 0x901b,
4187 0xddbb, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
4188 {0x0235, 0x0137, 0x36b1, 0x336c,
4189 0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
4190 {0x50f8, 0x25fb, 0xc76b, 0x6b71,
4191 0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
4192 {0x0000, 0x0000, 0x0000, 0x0000,
4193 0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
4194 {0x0000, 0x0000, 0x0000, 0x0000,
4195 0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
4196 {0x0000, 0x0000, 0x0000, 0x0000,
4197 0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
4198 {0x0000, 0x0000, 0x0000, 0x0000,
4199 0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
4200 {0x0000, 0x0000, 0x0000, 0x0000,
4201 0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
4202 {0x0000, 0x0000, 0x0000, 0x0000,
4203 0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,}, /* 10**1 */
4204};
4205
4206static unsigned EMUSHORT emtens[NTEN + 1][NE] =
4207{
4208 {0x2030, 0xcffc, 0xa1c3, 0x8123,
4209 0x2de3, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,}, /* 10**-4096 */
4210 {0x8264, 0xd2cb, 0xf2ea, 0x12d4,
4211 0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,}, /* 10**-2048 */
4212 {0xf53f, 0xf698, 0x6bd3, 0x0158,
4213 0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
4214 {0xe731, 0x04d4, 0xe3f2, 0xd332,
4215 0x7132, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
4216 {0xa23e, 0x5308, 0xfefb, 0x1155,
4217 0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
4218 {0xe26d, 0xdbde, 0xd05d, 0xb3f6,
4219 0xac7c, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
4220 {0x2a20, 0x6224, 0x47b3, 0x98d7,
4221 0x3f23, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
4222 {0x0b5b, 0x4af2, 0xa581, 0x18ed,
4223 0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
4224 {0xbf71, 0xa9b3, 0x7989, 0xbe68,
4225 0x4c2e, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
4226 {0x3d4d, 0x7c3d, 0x36ba, 0x0d2b,
4227 0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
4228 {0xc155, 0xa4a8, 0x404e, 0x6113,
4229 0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
4230 {0xd70a, 0x70a3, 0x0a3d, 0xa3d7,
4231 0x3d70, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
4232 {0xcccd, 0xcccc, 0xcccc, 0xcccc,
4233 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,}, /* 10**-1 */
4234};
4235#else
4236/* LONG_DOUBLE_TYPE_SIZE is other than 128 */
985b6196
RS		 4237static unsigned EMUSHORT etens[NTEN + 1][NE] =
4238{
4239 {0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,}, /* 10**4096 */
4240 {0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,}, /* 10**2048 */
4241 {0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
4242 {0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
4243 {0xddbc, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
4244 {0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
4245 {0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
4246 {0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
4247 {0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
4248 {0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
4249 {0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
4250 {0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
4251 {0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,}, /* 10**1 */
4252};
4253
4254static unsigned EMUSHORT emtens[NTEN + 1][NE] =
4255{
4256 {0x2de4, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,}, /* 10**-4096 */
4257 {0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,}, /* 10**-2048 */
4258 {0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
4259 {0x7133, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
4260 {0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
4261 {0xac7d, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
4262 {0x3f24, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
4263 {0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
4264 {0x4c2f, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
4265 {0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
4266 {0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
4267 {0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
4268 {0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,}, /* 10**-1 */
4269};
842fbaaa 4270#endif
985b6196 4271
a0353055 4272static void
985b6196
RS		 4273e24toasc (x, string, ndigs)
4274 unsigned EMUSHORT x[];
4275 char *string;
4276 int ndigs;
4277{
4278 unsigned EMUSHORT w[NI];
4279
985b6196
RS		 4280 e24toe (x, w);
4281 etoasc (w, string, ndigs);
4282}
4283
4284
a0353055 4285static void
985b6196
RS		 4286e53toasc (x, string, ndigs)
4287 unsigned EMUSHORT x[];
4288 char *string;
4289 int ndigs;
4290{
4291 unsigned EMUSHORT w[NI];
4292
985b6196
RS		 4293 e53toe (x, w);
4294 etoasc (w, string, ndigs);
4295}
4296
4297
a0353055 4298static void
985b6196
RS		 4299e64toasc (x, string, ndigs)
4300 unsigned EMUSHORT x[];
4301 char *string;
4302 int ndigs;
4303{
4304 unsigned EMUSHORT w[NI];
4305
985b6196
RS		 4306 e64toe (x, w);
4307 etoasc (w, string, ndigs);
4308}
4309
a0353055 4310static void
842fbaaa
JW		 4311e113toasc (x, string, ndigs)
4312 unsigned EMUSHORT x[];
4313 char *string;
4314 int ndigs;
4315{
4316 unsigned EMUSHORT w[NI];
4317
4318 e113toe (x, w);
4319 etoasc (w, string, ndigs);
4320}
4321
985b6196
RS		 4322
4323static char wstring[80]; /* working storage for ASCII output */
4324
a0353055 4325static void
985b6196
RS		 4326etoasc (x, string, ndigs)
4327 unsigned EMUSHORT x[];
4328 char *string;
4329 int ndigs;
4330{
4331 EMUSHORT digit;
4332 unsigned EMUSHORT y[NI], t[NI], u[NI], w[NI];
4333 unsigned EMUSHORT *p, *r, *ten;
4334 unsigned EMUSHORT sign;
4335 int i, j, k, expon, rndsav;
4336 char *s, *ss;
4337 unsigned EMUSHORT m;
4338
66b6d60b
RS		 4339
4340 rndsav = rndprc;
985b6196
RS		 4341 ss = string;
4342 s = wstring;
66b6d60b
RS		 4343 *ss = '\0';
4344 *s = '\0';
4345#ifdef NANS
4346 if (eisnan (x))
4347 {
4348 sprintf (wstring, " NaN ");
4349 goto bxit;
4350 }
4351#endif
985b6196
RS		 4352 rndprc = NBITS; /* set to full precision */
4353 emov (x, y); /* retain external format */
4354 if (y[NE - 1] & 0x8000)
4355 {
4356 sign = 0xffff;
4357 y[NE - 1] &= 0x7fff;
4358 }
4359 else
4360 {
4361 sign = 0;
4362 }
4363 expon = 0;
4364 ten = &etens[NTEN][0];
4365 emov (eone, t);
4366 /* Test for zero exponent */
4367 if (y[NE - 1] == 0)
4368 {
4369 for (k = 0; k < NE - 1; k++)
4370 {
4371 if (y[k] != 0)
4372 goto tnzro; /* denormalized number */
4373 }
43b55a67 4374 goto isone; /* valid all zeros */
985b6196
RS		 4375 }
4376 tnzro:
4377
66b6d60b 4378 /* Test for infinity. */
985b6196
RS		 4379 if (y[NE - 1] == 0x7fff)
4380 {
4381 if (sign)
4382 sprintf (wstring, " -Infinity ");
4383 else
4384 sprintf (wstring, " Infinity ");
4385 goto bxit;
4386 }
4387
4388 /* Test for exponent nonzero but significand denormalized.
4389 * This is an error condition.
4390 */
4391 if ((y[NE - 1] != 0) && ((y[NE - 2] & 0x8000) == 0))
4392 {
4393 mtherr ("etoasc", DOMAIN);
4394 sprintf (wstring, "NaN");
4395 goto bxit;
4396 }
4397
4398 /* Compare to 1.0 */
4399 i = ecmp (eone, y);
4400 if (i == 0)
4401 goto isone;
4402
66b6d60b
RS		 4403 if (i == -2)
4404 abort ();
4405
985b6196
RS		 4406 if (i < 0)
4407 { /* Number is greater than 1 */
4408 /* Convert significand to an integer and strip trailing decimal zeros. */
4409 emov (y, u);
4410 u[NE - 1] = EXONE + NBITS - 1;
4411
4412 p = &etens[NTEN - 4][0];
4413 m = 16;
4414 do
4415 {
4416 ediv (p, u, t);
4417 efloor (t, w);
4418 for (j = 0; j < NE - 1; j++)
4419 {
4420 if (t[j] != w[j])
4421 goto noint;
4422 }
4423 emov (t, u);
4424 expon += (int) m;
4425 noint:
4426 p += NE;
4427 m >>= 1;
4428 }
4429 while (m != 0);
4430
4431 /* Rescale from integer significand */
4432 u[NE - 1] += y[NE - 1] - (unsigned int) (EXONE + NBITS - 1);
4433 emov (u, y);
4434 /* Find power of 10 */
4435 emov (eone, t);
4436 m = MAXP;
4437 p = &etens[0][0];
66b6d60b 4438 /* An unordered compare result shouldn't happen here. */
985b6196
RS		 4439 while (ecmp (ten, u) <= 0)
4440 {
4441 if (ecmp (p, u) <= 0)
4442 {
4443 ediv (p, u, u);
4444 emul (p, t, t);
4445 expon += (int) m;
4446 }
4447 m >>= 1;
4448 if (m == 0)
4449 break;
4450 p += NE;
4451 }
4452 }
4453 else
4454 { /* Number is less than 1.0 */
4455 /* Pad significand with trailing decimal zeros. */
4456 if (y[NE - 1] == 0)
4457 {
4458 while ((y[NE - 2] & 0x8000) == 0)
4459 {
4460 emul (ten, y, y);
4461 expon -= 1;
4462 }
4463 }
4464 else
4465 {
4466 emovi (y, w);
4467 for (i = 0; i < NDEC + 1; i++)
4468 {
4469 if ((w[NI - 1] & 0x7) != 0)
4470 break;
4471 /* multiply by 10 */
4472 emovz (w, u);
4473 eshdn1 (u);
4474 eshdn1 (u);
4475 eaddm (w, u);
4476 u[1] += 3;
4477 while (u[2] != 0)
4478 {
4479 eshdn1 (u);
4480 u[1] += 1;
4481 }
4482 if (u[NI - 1] != 0)
4483 break;
4484 if (eone[NE - 1] <= u[1])
4485 break;
4486 emovz (u, w);
4487 expon -= 1;
4488 }
4489 emovo (w, y);
4490 }
4491 k = -MAXP;
4492 p = &emtens[0][0];
4493 r = &etens[0][0];
4494 emov (y, w);
4495 emov (eone, t);
4496 while (ecmp (eone, w) > 0)
4497 {
4498 if (ecmp (p, w) >= 0)
4499 {
4500 emul (r, w, w);
4501 emul (r, t, t);
4502 expon += k;
4503 }
4504 k /= 2;
4505 if (k == 0)
4506 break;
4507 p += NE;
4508 r += NE;
4509 }
4510 ediv (t, eone, t);
4511 }
4512 isone:
4513 /* Find the first (leading) digit. */
4514 emovi (t, w);
4515 emovz (w, t);
4516 emovi (y, w);
4517 emovz (w, y);
4518 eiremain (t, y);
4519 digit = equot[NI - 1];
4520 while ((digit == 0) && (ecmp (y, ezero) != 0))
4521 {
4522 eshup1 (y);
4523 emovz (y, u);
4524 eshup1 (u);
4525 eshup1 (u);
4526 eaddm (u, y);
4527 eiremain (t, y);
4528 digit = equot[NI - 1];
4529 expon -= 1;
4530 }
4531 s = wstring;
4532 if (sign)
4533 *s++ = '-';
4534 else
4535 *s++ = ' ';
985b6196
RS		 4536 /* Examine number of digits requested by caller. */
4537 if (ndigs < 0)
4538 ndigs = 0;
4539 if (ndigs > NDEC)
4540 ndigs = NDEC;
64685ffa
RS		 4541 if (digit == 10)
4542 {
4543 *s++ = '1';
4544 *s++ = '.';
4545 if (ndigs > 0)
4546 {
4547 *s++ = '0';
4548 ndigs -= 1;
4549 }
4550 expon += 1;
4551 }
4552 else
4553 {
242cef1e 4554 *s++ = (char)digit + '0';
64685ffa
RS		 4555 *s++ = '.';
4556 }
985b6196
RS		 4557 /* Generate digits after the decimal point. */
4558 for (k = 0; k <= ndigs; k++)
4559 {
4560 /* multiply current number by 10, without normalizing */
4561 eshup1 (y);
4562 emovz (y, u);
4563 eshup1 (u);
4564 eshup1 (u);
4565 eaddm (u, y);
4566 eiremain (t, y);
4567 *s++ = (char) equot[NI - 1] + '0';
4568 }
4569 digit = equot[NI - 1];
4570 --s;
4571 ss = s;
4572 /* round off the ASCII string */
4573 if (digit > 4)
4574 {
4575 /* Test for critical rounding case in ASCII output. */
4576 if (digit == 5)
4577 {
4578 emovo (y, t);
4579 if (ecmp (t, ezero) != 0)
4580 goto roun; /* round to nearest */
4581 if ((*(s - 1) & 1) == 0)
4582 goto doexp; /* round to even */
4583 }
4584 /* Round up and propagate carry-outs */
4585 roun:
4586 --s;
4587 k = *s & 0x7f;
4588 /* Carry out to most significant digit? */
4589 if (k == '.')
4590 {
4591 --s;
4592 k = *s;
4593 k += 1;
4594 *s = (char) k;
4595 /* Most significant digit carries to 10? */
4596 if (k > '9')
4597 {
4598 expon += 1;
4599 *s = '1';
4600 }
4601 goto doexp;
4602 }
4603 /* Round up and carry out from less significant digits */
4604 k += 1;
4605 *s = (char) k;
4606 if (k > '9')
4607 {
4608 *s = '0';
4609 goto roun;
4610 }
4611 }
4612 doexp:
4613 /*
4614 if (expon >= 0)
4615 sprintf (ss, "e+%d", expon);
4616 else
4617 sprintf (ss, "e%d", expon);
4618 */
4619 sprintf (ss, "e%d", expon);
4620 bxit:
4621 rndprc = rndsav;
4622 /* copy out the working string */
4623 s = string;
4624 ss = wstring;
4625 while (*ss == ' ') /* strip possible leading space */
4626 ++ss;
4627 while ((*s++ = *ss++) != '\0')
4628 ;
4629}
4630
4631
defb5dab 4632/* Convert ASCII string to quadruple precision floating point
985b6196 4633
defb5dab
RK		 4634 Numeric input is free field decimal number with max of 15 digits with or
4635 without decimal point entered as ASCII from teletype. Entering E after
4636 the number followed by a second number causes the second number to be
4637 interpreted as a power of 10 to be multiplied by the first number
4638 (i.e., "scientific" notation). */
985b6196
RS		 4639
4640/* ASCII to single */
a0353055
RK		 4641
4642static void
985b6196
RS		 4643asctoe24 (s, y)
4644 char *s;
4645 unsigned EMUSHORT *y;
4646{
4647 asctoeg (s, y, 24);
4648}
4649
4650
4651/* ASCII to double */
a0353055
RK		 4652
4653static void
985b6196
RS		 4654asctoe53 (s, y)
4655 char *s;
4656 unsigned EMUSHORT *y;
4657{
842fbaaa 4658#if defined(DEC) || defined(IBM)
985b6196
RS		 4659 asctoeg (s, y, 56);
4660#else
4661 asctoeg (s, y, 53);
4662#endif
4663}
4664
4665
4666/* ASCII to long double */
a0353055
RK		 4667
4668static void
985b6196
RS		 4669asctoe64 (s, y)
4670 char *s;
4671 unsigned EMUSHORT *y;
4672{
4673 asctoeg (s, y, 64);
4674}
4675
842fbaaa 4676/* ASCII to 128-bit long double */
a0353055
RK		 4677
4678static void
842fbaaa
JW		 4679asctoe113 (s, y)
4680 char *s;
4681 unsigned EMUSHORT *y;
4682{
4683 asctoeg (s, y, 113);
4684}
4685
985b6196 4686/* ASCII to super double */
defb5dab 4687
a0353055 4688static void
985b6196
RS		 4689asctoe (s, y)
4690 char *s;
4691 unsigned EMUSHORT *y;
4692{
4693 asctoeg (s, y, NBITS);
4694}
4695
985b6196 4696
d73e9b8d 4697/* ASCII to e type, with specified rounding precision = oprec. */
defb5dab 4698
a0353055 4699static void
985b6196
RS		 4700asctoeg (ss, y, oprec)
4701 char *ss;
4702 unsigned EMUSHORT *y;
4703 int oprec;
4704{
4705 unsigned EMUSHORT yy[NI], xt[NI], tt[NI];
4706 int esign, decflg, sgnflg, nexp, exp, prec, lost;
4707 int k, trail, c, rndsav;
4708 EMULONG lexp;
4709 unsigned EMUSHORT nsign, *p;
d73e9b8d 4710 char *sp, *s, *lstr;
985b6196
RS		 4711
4712 /* Copy the input string. */
d73e9b8d 4713 lstr = (char *) alloca (strlen (ss) + 1);
985b6196
RS		 4714 s = ss;
4715 while (*s == ' ') /* skip leading spaces */
4716 ++s;
4717 sp = lstr;
a9456cd3
RS		 4718 while ((*sp++ = *s++) != '\0')
4719 ;
985b6196
RS		 4720 s = lstr;
4721
4722 rndsav = rndprc;
4723 rndprc = NBITS; /* Set to full precision */
4724 lost = 0;
4725 nsign = 0;
4726 decflg = 0;
4727 sgnflg = 0;
4728 nexp = 0;
4729 exp = 0;
4730 prec = 0;
4731 ecleaz (yy);
4732 trail = 0;
4733
4734 nxtcom:
4735 k = *s - '0';
4736 if ((k >= 0) && (k <= 9))
4737 {
4738 /* Ignore leading zeros */
4739 if ((prec == 0) && (decflg == 0) && (k == 0))
4740 goto donchr;
4741 /* Identify and strip trailing zeros after the decimal point. */
4742 if ((trail == 0) && (decflg != 0))
4743 {
4744 sp = s;
4745 while ((*sp >= '0') && (*sp <= '9'))
4746 ++sp;
4747 /* Check for syntax error */
4748 c = *sp & 0x7f;
4749 if ((c != 'e') && (c != 'E') && (c != '\0')
4750 && (c != '\n') && (c != '\r') && (c != ' ')
4751 && (c != ','))
4752 goto error;
4753 --sp;
4754 while (*sp == '0')
4755 *sp-- = 'z';
4756 trail = 1;
4757 if (*s == 'z')
4758 goto donchr;
4759 }
defb5dab 4760
985b6196 4761 /* If enough digits were given to more than fill up the yy register,
defb5dab
RK		 4762 continuing until overflow into the high guard word yy[2]
4763 guarantees that there will be a roundoff bit at the top
4764 of the low guard word after normalization. */
4765
985b6196
RS		 4766 if (yy[2] == 0)
4767 {
4768 if (decflg)
4769 nexp += 1; /* count digits after decimal point */
4770 eshup1 (yy); /* multiply current number by 10 */
4771 emovz (yy, xt);
4772 eshup1 (xt);
4773 eshup1 (xt);
4774 eaddm (xt, yy);
4775 ecleaz (xt);
4776 xt[NI - 2] = (unsigned EMUSHORT) k;
4777 eaddm (xt, yy);
4778 }
4779 else
4780 {
d73e9b8d 4781 /* Mark any lost non-zero digit. */
985b6196 4782 lost |= k;
d73e9b8d
RS		 4783 /* Count lost digits before the decimal point. */
4784 if (decflg == 0)
4785 nexp -= 1;
985b6196
RS		 4786 }
4787 prec += 1;
4788 goto donchr;
4789 }
4790
4791 switch (*s)
4792 {
4793 case 'z':
4794 break;
4795 case 'E':
4796 case 'e':
4797 goto expnt;
4798 case '.': /* decimal point */
4799 if (decflg)
4800 goto error;
4801 ++decflg;
4802 break;
4803 case '-':
4804 nsign = 0xffff;
4805 if (sgnflg)
4806 goto error;
4807 ++sgnflg;
4808 break;
4809 case '+':
4810 if (sgnflg)
4811 goto error;
4812 ++sgnflg;
4813 break;
4814 case ',':
4815 case ' ':
4816 case '\0':
4817 case '\n':
4818 case '\r':
4819 goto daldone;
4820 case 'i':
4821 case 'I':
64685ffa 4822 goto infinite;
985b6196
RS		 4823 default:
4824 error:
66b6d60b
RS		 4825#ifdef NANS
4826 einan (yy);
4827#else
985b6196 4828 mtherr ("asctoe", DOMAIN);
66b6d60b
RS		 4829 eclear (yy);
4830#endif
985b6196
RS		 4831 goto aexit;
4832 }
4833 donchr:
4834 ++s;
4835 goto nxtcom;
4836
4837 /* Exponent interpretation */
4838 expnt:
4839
4840 esign = 1;
4841 exp = 0;
4842 ++s;
4843 /* check for + or - */
4844 if (*s == '-')
4845 {
4846 esign = -1;
4847 ++s;
4848 }
4849 if (*s == '+')
4850 ++s;
4851 while ((*s >= '0') && (*s <= '9'))
4852 {
4853 exp *= 10;
4854 exp += *s++ - '0';
842fbaaa 4855 if (exp > -(MINDECEXP))
64685ffa
RS		 4856 {
4857 if (esign < 0)
4858 goto zero;
4859 else
4860 goto infinite;
4861 }
985b6196
RS		 4862 }
4863 if (esign < 0)
4864 exp = -exp;
842fbaaa 4865 if (exp > MAXDECEXP)
64685ffa
RS		 4866 {
4867 infinite:
4868 ecleaz (yy);
4869 yy[E] = 0x7fff; /* infinity */
4870 goto aexit;
4871 }
842fbaaa 4872 if (exp < MINDECEXP)
64685ffa
RS		 4873 {
4874 zero:
4875 ecleaz (yy);
4876 goto aexit;
4877 }
985b6196
RS		 4878
4879 daldone:
4880 nexp = exp - nexp;
4881 /* Pad trailing zeros to minimize power of 10, per IEEE spec. */
4882 while ((nexp > 0) && (yy[2] == 0))
4883 {
4884 emovz (yy, xt);
4885 eshup1 (xt);
4886 eshup1 (xt);
4887 eaddm (yy, xt);
4888 eshup1 (xt);
4889 if (xt[2] != 0)
4890 break;
4891 nexp -= 1;
4892 emovz (xt, yy);
4893 }
4894 if ((k = enormlz (yy)) > NBITS)
4895 {
4896 ecleaz (yy);
4897 goto aexit;
4898 }
4899 lexp = (EXONE - 1 + NBITS) - k;
4900 emdnorm (yy, lost, 0, lexp, 64);
985b6196 4901
defb5dab
RK		 4902 /* Convert to external format:
4903
4904 Multiply by 10**nexp. If precision is 64 bits,
4905 the maximum relative error incurred in forming 10**n
4906 for 0 <= n <= 324 is 8.2e-20, at 10**180.
4907 For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
4908 For 0 >= n >= -999, it is -1.55e-19 at 10**-435. */
985b6196 4909
985b6196
RS		 4910 lexp = yy[E];
4911 if (nexp == 0)
4912 {
4913 k = 0;
4914 goto expdon;
4915 }
4916 esign = 1;
4917 if (nexp < 0)
4918 {
4919 nexp = -nexp;
4920 esign = -1;
4921 if (nexp > 4096)
defb5dab
RK		 4922 {
4923 /* Punt. Can't handle this without 2 divides. */
985b6196
RS		 4924 emovi (etens[0], tt);
4925 lexp -= tt[E];
4926 k = edivm (tt, yy);
4927 lexp += EXONE;
4928 nexp -= 4096;
4929 }
4930 }
4931 p = &etens[NTEN][0];
4932 emov (eone, xt);
4933 exp = 1;
4934 do
4935 {
4936 if (exp & nexp)
4937 emul (p, xt, xt);
4938 p -= NE;
4939 exp = exp + exp;
4940 }
4941 while (exp <= MAXP);
4942
4943 emovi (xt, tt);
4944 if (esign < 0)
4945 {
4946 lexp -= tt[E];
4947 k = edivm (tt, yy);
4948 lexp += EXONE;
4949 }
4950 else
4951 {
4952 lexp += tt[E];
4953 k = emulm (tt, yy);
4954 lexp -= EXONE - 1;
4955 }
4956
4957 expdon:
4958
4959 /* Round and convert directly to the destination type */
4960 if (oprec == 53)
4961 lexp -= EXONE - 0x3ff;
842fbaaa
JW		 4962#ifdef IBM
4963 else if (oprec == 24 || oprec == 56)
4964 lexp -= EXONE - (0x41 << 2);
4965#else
985b6196
RS		 4966 else if (oprec == 24)
4967 lexp -= EXONE - 0177;
842fbaaa 4968#endif
985b6196
RS		 4969#ifdef DEC
4970 else if (oprec == 56)
4971 lexp -= EXONE - 0201;
4972#endif
4973 rndprc = oprec;
4974 emdnorm (yy, k, 0, lexp, 64);
4975
4976 aexit:
4977
4978 rndprc = rndsav;
4979 yy[0] = nsign;
4980 switch (oprec)
4981 {
4982#ifdef DEC
4983 case 56:
4984 todec (yy, y); /* see etodec.c */
4985 break;
842fbaaa
JW		 4986#endif
4987#ifdef IBM
4988 case 56:
4989 toibm (yy, y, DFmode);
4990 break;
985b6196
RS		 4991#endif
4992 case 53:
4993 toe53 (yy, y);
4994 break;
4995 case 24:
4996 toe24 (yy, y);
4997 break;
4998 case 64:
4999 toe64 (yy, y);
5000 break;
842fbaaa
JW		 5001 case 113:
5002 toe113 (yy, y);
5003 break;
985b6196
RS		 5004 case NBITS:
5005 emovo (yy, y);
5006 break;
5007 }
5008}
5009
5010
5011
defb5dab
RK		 5012/* y = largest integer not greater than x (truncated toward minus infinity) */
5013
985b6196
RS		 5014static unsigned EMUSHORT bmask[] =
5015{
5016 0xffff,
5017 0xfffe,
5018 0xfffc,
5019 0xfff8,
5020 0xfff0,
5021 0xffe0,
5022 0xffc0,
5023 0xff80,
5024 0xff00,
5025 0xfe00,
5026 0xfc00,
5027 0xf800,
5028 0xf000,
5029 0xe000,
5030 0xc000,
5031 0x8000,
5032 0x0000,
5033};
5034
a0353055 5035static void
985b6196
RS		 5036efloor (x, y)
5037 unsigned EMUSHORT x[], y[];
5038{
5039 register unsigned EMUSHORT *p;
5040 int e, expon, i;
5041 unsigned EMUSHORT f[NE];
5042
5043 emov (x, f); /* leave in external format */
5044 expon = (int) f[NE - 1];
5045 e = (expon & 0x7fff) - (EXONE - 1);
5046 if (e <= 0)
5047 {
5048 eclear (y);
5049 goto isitneg;
5050 }
5051 /* number of bits to clear out */
5052 e = NBITS - e;
5053 emov (f, y);
5054 if (e <= 0)
5055 return;
5056
5057 p = &y[0];
5058 while (e >= 16)
5059 {
5060 *p++ = 0;
5061 e -= 16;
5062 }
5063 /* clear the remaining bits */
5064 *p &= bmask[e];
5065 /* truncate negatives toward minus infinity */
5066 isitneg:
5067
5068 if ((unsigned EMUSHORT) expon & (unsigned EMUSHORT) 0x8000)
5069 {
5070 for (i = 0; i < NE - 1; i++)
5071 {
5072 if (f[i] != y[i])
5073 {
5074 esub (eone, y, y);
5075 break;
5076 }
5077 }
5078 }
5079}
5080
5081
defb5dab
RK		 5082/* Returns s and exp such that s * 2**exp = x and .5 <= s < 1.
5083 For example, 1.1 = 0.55 * 2**1
5084 Handles denormalized numbers properly using long integer exp. */
a0353055
RK		 5085
5086static void
985b6196
RS		 5087efrexp (x, exp, s)
5088 unsigned EMUSHORT x[];
5089 int *exp;
5090 unsigned EMUSHORT s[];
5091{
5092 unsigned EMUSHORT xi[NI];
5093 EMULONG li;
5094
5095 emovi (x, xi);
5096 li = (EMULONG) ((EMUSHORT) xi[1]);
5097
5098 if (li == 0)
5099 {
5100 li -= enormlz (xi);
5101 }
5102 xi[1] = 0x3ffe;
5103 emovo (xi, s);
5104 *exp = (int) (li - 0x3ffe);
5105}
5106
5107
5108
defb5dab 5109/* Return y = x * 2**pwr2. */
a0353055
RK		 5110
5111static void
985b6196
RS		 5112eldexp (x, pwr2, y)
5113 unsigned EMUSHORT x[];
5114 int pwr2;
5115 unsigned EMUSHORT y[];
5116{
5117 unsigned EMUSHORT xi[NI];
5118 EMULONG li;
5119 int i;
5120
5121 emovi (x, xi);
5122 li = xi[1];
5123 li += pwr2;
5124 i = 0;
5125 emdnorm (xi, i, i, li, 64);
5126 emovo (xi, y);
5127}
5128
5129
5130/* c = remainder after dividing b by a
defb5dab 5131 Least significant integer quotient bits left in equot[]. */
a0353055
RK		 5132
5133static void
985b6196
RS		 5134eremain (a, b, c)
5135 unsigned EMUSHORT a[], b[], c[];
5136{
5137 unsigned EMUSHORT den[NI], num[NI];
5138
66b6d60b 5139#ifdef NANS
242cef1e
RS		 5140 if (eisinf (b)
5141 || (ecmp (a, ezero) == 0)
5142 || eisnan (a)
5143 || eisnan (b))
66b6d60b 5144 {
29e11dab 5145 enan (c, 0);
66b6d60b
RS		 5146 return;
5147 }
5148#endif
985b6196
RS		 5149 if (ecmp (a, ezero) == 0)
5150 {
5151 mtherr ("eremain", SING);
5152 eclear (c);
5153 return;
5154 }
5155 emovi (a, den);
5156 emovi (b, num);
5157 eiremain (den, num);
5158 /* Sign of remainder = sign of quotient */
5159 if (a[0] == b[0])
5160 num[0] = 0;
5161 else
5162 num[0] = 0xffff;
5163 emovo (num, c);
5164}
5165
a0353055 5166static void
985b6196
RS		 5167eiremain (den, num)
5168 unsigned EMUSHORT den[], num[];
5169{
5170 EMULONG ld, ln;
5171 unsigned EMUSHORT j;
5172
5173 ld = den[E];
5174 ld -= enormlz (den);
5175 ln = num[E];
5176 ln -= enormlz (num);
5177 ecleaz (equot);
5178 while (ln >= ld)
5179 {
5180 if (ecmpm (den, num) <= 0)
5181 {
5182 esubm (den, num);
5183 j = 1;
5184 }
5185 else
5186 {
5187 j = 0;
5188 }
5189 eshup1 (equot);
5190 equot[NI - 1] |= j;
5191 eshup1 (num);
5192 ln -= 1;
5193 }
5194 emdnorm (num, 0, 0, ln, 0);
5195}
5196
defb5dab
RK		 5197/* This routine may be called to report one of the following
5198 error conditions (in the include file mconf.h).
5199
5200 Mnemonic Value Significance
5201
5202 DOMAIN 1 argument domain error
5203 SING 2 function singularity
5204 OVERFLOW 3 overflow range error
5205 UNDERFLOW 4 underflow range error
5206 TLOSS 5 total loss of precision
5207 PLOSS 6 partial loss of precision
5208 INVALID 7 NaN - producing operation
5209 EDOM 33 Unix domain error code
5210 ERANGE 34 Unix range error code
5211
5212 The default version of the file prints the function name,
5213 passed to it by the pointer fctnam, followed by the
5214 error condition. The display is directed to the standard
5215 output device. The routine then returns to the calling
5216 program. Users may wish to modify the program to abort by
5217 calling exit under severe error conditions such as domain
5218 errors.
5219
5220 Since all error conditions pass control to this function,
5221 the display may be easily changed, eliminated, or directed
5222 to an error logging device. */
5223
5224/* Note: the order of appearance of the following messages is bound to the
5225 error codes defined above. */
985b6196 5226
66b6d60b
RS		 5227#define NMSGS 8
5228static char *ermsg[NMSGS] =
985b6196
RS		 5229{
5230 "unknown", /* error code 0 */
5231 "domain", /* error code 1 */
5232 "singularity", /* et seq. */
5233 "overflow",
5234 "underflow",
5235 "total loss of precision",
66b6d60b
RS		 5236 "partial loss of precision",
5237 "invalid operation"
985b6196
RS		 5238};
5239
5240int merror = 0;
5241extern int merror;
5242
a0353055 5243static void
985b6196
RS		 5244mtherr (name, code)
5245 char *name;
5246 int code;
5247{
5248 char errstr[80];
5249
defb5dab
RK		 5250 /* Display string passed by calling program, which is supposed to be the
5251 name of the function in which the error occurred.
5252
5253 Display error message defined by the code argument. */
985b6196 5254
66b6d60b 5255 if ((code <= 0) || (code >= NMSGS))
985b6196 5256 code = 0;
a8d78514 5257 sprintf (errstr, " %s %s error", name, ermsg[code]);
64685ffa
RS		 5258 if (extra_warnings)
5259 warning (errstr);
985b6196
RS		 5260 /* Set global error message word */
5261 merror = code + 1;
985b6196
RS		 5262}
5263
842fbaaa 5264#ifdef DEC
defb5dab 5265/* Convert DEC double precision to e type. */
a0353055
RK		 5266
5267static void
985b6196
RS		 5268dectoe (d, e)
5269 unsigned EMUSHORT *d;
5270 unsigned EMUSHORT *e;
5271{
5272 unsigned EMUSHORT y[NI];
5273 register unsigned EMUSHORT r, *p;
5274
5275 ecleaz (y); /* start with a zero */
5276 p = y; /* point to our number */
5277 r = *d; /* get DEC exponent word */
5278 if (*d & (unsigned int) 0x8000)
5279 *p = 0xffff; /* fill in our sign */
5280 ++p; /* bump pointer to our exponent word */
5281 r &= 0x7fff; /* strip the sign bit */
5282 if (r == 0) /* answer = 0 if high order DEC word = 0 */
5283 goto done;
5284
5285
5286 r >>= 7; /* shift exponent word down 7 bits */
5287 r += EXONE - 0201; /* subtract DEC exponent offset */
5288 /* add our e type exponent offset */
5289 *p++ = r; /* to form our exponent */
5290
5291 r = *d++; /* now do the high order mantissa */
5292 r &= 0177; /* strip off the DEC exponent and sign bits */
5293 r |= 0200; /* the DEC understood high order mantissa bit */
5294 *p++ = r; /* put result in our high guard word */
5295
5296 *p++ = *d++; /* fill in the rest of our mantissa */
5297 *p++ = *d++;
5298 *p = *d;
5299
5300 eshdn8 (y); /* shift our mantissa down 8 bits */
5301 done:
5302 emovo (y, e);
5303}
5304
5305
5306
5307/*
5308; convert e type to DEC double precision
5309; double d;
5310; EMUSHORT e[NE];
5311; etodec (e, &d);
5312*/
985b6196 5313
a0353055 5314static void
985b6196
RS		 5315etodec (x, d)
5316 unsigned EMUSHORT *x, *d;
5317{
5318 unsigned EMUSHORT xi[NI];
842fbaaa
JW		 5319 EMULONG exp;
5320 int rndsav;
985b6196
RS		 5321
5322 emovi (x, xi);
5323 exp = (EMULONG) xi[E] - (EXONE - 0201); /* adjust exponent for offsets */
5324/* round off to nearest or even */
5325 rndsav = rndprc;
5326 rndprc = 56;
5327 emdnorm (xi, 0, 0, exp, 64);
5328 rndprc = rndsav;
5329 todec (xi, d);
5330}
5331
a0353055 5332static void
985b6196
RS		 5333todec (x, y)
5334 unsigned EMUSHORT *x, *y;
5335{
5336 unsigned EMUSHORT i;
5337 unsigned EMUSHORT *p;
5338
5339 p = x;
5340 *y = 0;
5341 if (*p++)
5342 *y = 0100000;
5343 i = *p++;
5344 if (i == 0)
5345 {
5346 *y++ = 0;
5347 *y++ = 0;
5348 *y++ = 0;
5349 *y++ = 0;
5350 return;
5351 }
5352 if (i > 0377)
5353 {
5354 *y++ |= 077777;
5355 *y++ = 0xffff;
5356 *y++ = 0xffff;
5357 *y++ = 0xffff;
64685ffa
RS		 5358#ifdef ERANGE
5359 errno = ERANGE;
5360#endif
985b6196
RS		 5361 return;
5362 }
5363 i &= 0377;
5364 i <<= 7;
5365 eshup8 (x);
5366 x[M] &= 0177;
5367 i |= x[M];
5368 *y++ |= i;
5369 *y++ = x[M + 1];
5370 *y++ = x[M + 2];
5371 *y++ = x[M + 3];
5372}
842fbaaa
JW		 5373#endif /* DEC */
5374
5375#ifdef IBM
defb5dab 5376/* Convert IBM single/double precision to e type. */
a0353055
RK		 5377
5378static void
842fbaaa
JW		 5379ibmtoe (d, e, mode)
5380 unsigned EMUSHORT *d;
5381 unsigned EMUSHORT *e;
5382 enum machine_mode mode;
5383{
5384 unsigned EMUSHORT y[NI];
5385 register unsigned EMUSHORT r, *p;
5386 int rndsav;
5387
5388 ecleaz (y); /* start with a zero */
5389 p = y; /* point to our number */
5390 r = *d; /* get IBM exponent word */
5391 if (*d & (unsigned int) 0x8000)
5392 *p = 0xffff; /* fill in our sign */
5393 ++p; /* bump pointer to our exponent word */
5394 r &= 0x7f00; /* strip the sign bit */
5395 r >>= 6; /* shift exponent word down 6 bits */
5396 /* in fact shift by 8 right and 2 left */
5397 r += EXONE - (0x41 << 2); /* subtract IBM exponent offset */
5398 /* add our e type exponent offset */
5399 *p++ = r; /* to form our exponent */
5400
5401 *p++ = *d++ & 0xff; /* now do the high order mantissa */
5402 /* strip off the IBM exponent and sign bits */
5403 if (mode != SFmode) /* there are only 2 words in SFmode */
5404 {
5405 *p++ = *d++; /* fill in the rest of our mantissa */
5406 *p++ = *d++;
5407 }
5408 *p = *d;
5409
5410 if (y[M] == 0 && y[M+1] == 0 && y[M+2] == 0 && y[M+3] == 0)
5411 y[0] = y[E] = 0;
5412 else
5413 y[E] -= 5 + enormlz (y); /* now normalise the mantissa */
5414 /* handle change in RADIX */
5415 emovo (y, e);
5416}
5417
985b6196 5418
985b6196 5419
defb5dab 5420/* Convert e type to IBM single/double precision. */
842fbaaa 5421
a0353055 5422static void
842fbaaa
JW		 5423etoibm (x, d, mode)
5424 unsigned EMUSHORT *x, *d;
5425 enum machine_mode mode;
5426{
5427 unsigned EMUSHORT xi[NI];
5428 EMULONG exp;
5429 int rndsav;
5430
5431 emovi (x, xi);
5432 exp = (EMULONG) xi[E] - (EXONE - (0x41 << 2)); /* adjust exponent for offsets */
5433 /* round off to nearest or even */
5434 rndsav = rndprc;
5435 rndprc = 56;
5436 emdnorm (xi, 0, 0, exp, 64);
5437 rndprc = rndsav;
5438 toibm (xi, d, mode);
5439}
5440
a0353055 5441static void
842fbaaa
JW		 5442toibm (x, y, mode)
5443 unsigned EMUSHORT *x, *y;
5444 enum machine_mode mode;
5445{
5446 unsigned EMUSHORT i;
5447 unsigned EMUSHORT *p;
5448 int r;
5449
5450 p = x;
5451 *y = 0;
5452 if (*p++)
5453 *y = 0x8000;
5454 i = *p++;
5455 if (i == 0)
5456 {
5457 *y++ = 0;
5458 *y++ = 0;
5459 if (mode != SFmode)
5460 {
5461 *y++ = 0;
5462 *y++ = 0;
5463 }
5464 return;
5465 }
5466 r = i & 0x3;
5467 i >>= 2;
5468 if (i > 0x7f)
5469 {
5470 *y++ |= 0x7fff;
5471 *y++ = 0xffff;
5472 if (mode != SFmode)
5473 {
5474 *y++ = 0xffff;
5475 *y++ = 0xffff;
5476 }
5477#ifdef ERANGE
5478 errno = ERANGE;
5479#endif
5480 return;
5481 }
5482 i &= 0x7f;
5483 *y |= (i << 8);
5484 eshift (x, r + 5);
5485 *y++ |= x[M];
5486 *y++ = x[M + 1];
5487 if (mode != SFmode)
5488 {
5489 *y++ = x[M + 2];
5490 *y++ = x[M + 3];
5491 }
5492}
5493#endif /* IBM */
66b6d60b
RS		 5494
5495/* Output a binary NaN bit pattern in the target machine's format. */
5496
5497/* If special NaN bit patterns are required, define them in tm.h
5498 as arrays of unsigned 16-bit shorts. Otherwise, use the default
5499 patterns here. */
7729f1ca
RS		 5500#ifdef TFMODE_NAN
5501TFMODE_NAN;
5502#else
f76b9db2
ILT		 5503#ifdef IEEE
5504unsigned EMUSHORT TFbignan[8] =
66b6d60b 5505 {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
f76b9db2 5506unsigned EMUSHORT TFlittlenan[8] = {0, 0, 0, 0, 0, 0, 0x8000, 0xffff};
66b6d60b
RS		 5507#endif
5508#endif
5509
7729f1ca
RS		 5510#ifdef XFMODE_NAN
5511XFMODE_NAN;
5512#else
f76b9db2
ILT		 5513#ifdef IEEE
5514unsigned EMUSHORT XFbignan[6] =
5515 {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
5516unsigned EMUSHORT XFlittlenan[6] = {0, 0, 0, 0xc000, 0xffff, 0};
66b6d60b
RS		 5517#endif
5518#endif
5519
7729f1ca
RS		 5520#ifdef DFMODE_NAN
5521DFMODE_NAN;
5522#else
f76b9db2
ILT		 5523#ifdef IEEE
5524unsigned EMUSHORT DFbignan[4] = {0x7fff, 0xffff, 0xffff, 0xffff};
5525unsigned EMUSHORT DFlittlenan[4] = {0, 0, 0, 0xfff8};
66b6d60b
RS		 5526#endif
5527#endif
5528
7729f1ca
RS		 5529#ifdef SFMODE_NAN
5530SFMODE_NAN;
5531#else
f76b9db2
ILT		 5532#ifdef IEEE
5533unsigned EMUSHORT SFbignan[2] = {0x7fff, 0xffff};
5534unsigned EMUSHORT SFlittlenan[2] = {0, 0xffc0};
66b6d60b
RS		 5535#endif
5536#endif
5537
5538
a0353055 5539static void
29e11dab 5540make_nan (nan, sign, mode)
a0353055
RK		 5541 unsigned EMUSHORT *nan;
5542 int sign;
5543 enum machine_mode mode;
66b6d60b 5544{
29e11dab 5545 int n;
66b6d60b
RS		 5546 unsigned EMUSHORT *p;
5547
5548 switch (mode)
5549 {
5550/* Possibly the `reserved operand' patterns on a VAX can be
5551 used like NaN's, but probably not in the same way as IEEE. */
842fbaaa 5552#if !defined(DEC) && !defined(IBM)
66b6d60b
RS		 5553 case TFmode:
5554 n = 8;
f76b9db2
ILT		 5555 if (FLOAT_WORDS_BIG_ENDIAN)
5556 p = TFbignan;
5557 else
5558 p = TFlittlenan;
66b6d60b
RS		 5559 break;
5560 case XFmode:
5561 n = 6;
f76b9db2
ILT		 5562 if (FLOAT_WORDS_BIG_ENDIAN)
5563 p = XFbignan;
5564 else
5565 p = XFlittlenan;
66b6d60b
RS		 5566 break;
5567 case DFmode:
5568 n = 4;
f76b9db2
ILT		 5569 if (FLOAT_WORDS_BIG_ENDIAN)
5570 p = DFbignan;
5571 else
5572 p = DFlittlenan;
66b6d60b 5573 break;
bfbc6416 5574 case HFmode:
66b6d60b
RS		 5575 case SFmode:
5576 n = 2;
f76b9db2
ILT		 5577 if (FLOAT_WORDS_BIG_ENDIAN)
5578 p = SFbignan;
5579 else
5580 p = SFlittlenan;
66b6d60b
RS		 5581 break;
5582#endif
5583 default:
5584 abort ();
5585 }
f76b9db2
ILT		 5586 if (FLOAT_WORDS_BIG_ENDIAN)
5587 *nan++ = (sign << 15) | *p++;
29e11dab 5588 while (--n != 0)
66b6d60b 5589 *nan++ = *p++;
f76b9db2
ILT		 5590 if (! FLOAT_WORDS_BIG_ENDIAN)
5591 *nan = (sign << 15) | *p;
66b6d60b
RS		 5592}
5593
b31c244f
RS		 5594/* Convert an SFmode target `float' value to a REAL_VALUE_TYPE.
5595 This is the inverse of the function `etarsingle' invoked by
5596 REAL_VALUE_TO_TARGET_SINGLE. */
5597
5598REAL_VALUE_TYPE
5599ereal_from_float (f)
04ae9e4c 5600 HOST_WIDE_INT f;
b31c244f
RS		 5601{
5602 REAL_VALUE_TYPE r;
5603 unsigned EMUSHORT s[2];
5604 unsigned EMUSHORT e[NE];
5605
5606 /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
5607 This is the inverse operation to what the function `endian' does. */
f76b9db2
ILT		 5608 if (FLOAT_WORDS_BIG_ENDIAN)
5609 {
5610 s[0] = (unsigned EMUSHORT) (f >> 16);
5611 s[1] = (unsigned EMUSHORT) f;
5612 }
5613 else
5614 {
5615 s[0] = (unsigned EMUSHORT) f;
5616 s[1] = (unsigned EMUSHORT) (f >> 16);
5617 }
b31c244f
RS		 5618 /* Convert and promote the target float to E-type. */
5619 e24toe (s, e);
5620 /* Output E-type to REAL_VALUE_TYPE. */
5621 PUT_REAL (e, &r);
5622 return r;
5623}
5624
842fbaaa 5625
b31c244f
RS		 5626/* Convert a DFmode target `double' value to a REAL_VALUE_TYPE.
5627 This is the inverse of the function `etardouble' invoked by
5628 REAL_VALUE_TO_TARGET_DOUBLE.
5629
04ae9e4c
RK		 5630 The DFmode is stored as an array of HOST_WIDE_INT in the target's
5631 data format, with no holes in the bit packing. The first element
b31c244f
RS		 5632 of the input array holds the bits that would come first in the
5633 target computer's memory. */
5634
5635REAL_VALUE_TYPE
5636ereal_from_double (d)
04ae9e4c 5637 HOST_WIDE_INT d[];
b31c244f
RS		 5638{
5639 REAL_VALUE_TYPE r;
5640 unsigned EMUSHORT s[4];
5641 unsigned EMUSHORT e[NE];
5642
04ae9e4c 5643 /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces. */
f76b9db2
ILT		 5644 if (FLOAT_WORDS_BIG_ENDIAN)
5645 {
5646 s[0] = (unsigned EMUSHORT) (d[0] >> 16);
5647 s[1] = (unsigned EMUSHORT) d[0];
60e61165 5648#if HOST_BITS_PER_WIDE_INT == 32
f76b9db2
ILT		 5649 s[2] = (unsigned EMUSHORT) (d[1] >> 16);
5650 s[3] = (unsigned EMUSHORT) d[1];
60e61165 5651#else
f76b9db2
ILT		 5652 /* In this case the entire target double is contained in the
5653 first array element. The second element of the input is
5654 ignored. */
5655 s[2] = (unsigned EMUSHORT) (d[0] >> 48);
5656 s[3] = (unsigned EMUSHORT) (d[0] >> 32);
60e61165 5657#endif
f76b9db2
ILT		 5658 }
5659 else
5660 {
5661 /* Target float words are little-endian. */
5662 s[0] = (unsigned EMUSHORT) d[0];
5663 s[1] = (unsigned EMUSHORT) (d[0] >> 16);
60e61165 5664#if HOST_BITS_PER_WIDE_INT == 32
f76b9db2
ILT		 5665 s[2] = (unsigned EMUSHORT) d[1];
5666 s[3] = (unsigned EMUSHORT) (d[1] >> 16);
60e61165 5667#else
f76b9db2
ILT		 5668 s[2] = (unsigned EMUSHORT) (d[0] >> 32);
5669 s[3] = (unsigned EMUSHORT) (d[0] >> 48);
b31c244f 5670#endif
f76b9db2 5671 }
b31c244f
RS		 5672 /* Convert target double to E-type. */
5673 e53toe (s, e);
5674 /* Output E-type to REAL_VALUE_TYPE. */
5675 PUT_REAL (e, &r);
5676 return r;
5677}
842fbaaa
JW		 5678
5679
b51ab098
RK		 5680/* Convert target computer unsigned 64-bit integer to e-type.
5681 The endian-ness of DImode follows the convention for integers,
5682 so we use WORDS_BIG_ENDIAN here, not FLOAT_WORDS_BIG_ENDIAN. */
842fbaaa 5683
a0353055 5684static void
842fbaaa
JW		 5685uditoe (di, e)
5686 unsigned EMUSHORT *di; /* Address of the 64-bit int. */
5687 unsigned EMUSHORT *e;
5688{
5689 unsigned EMUSHORT yi[NI];
5690 int k;
5691
5692 ecleaz (yi);
f76b9db2
ILT		 5693 if (WORDS_BIG_ENDIAN)
5694 {
5695 for (k = M; k < M + 4; k++)
5696 yi[k] = *di++;
5697 }
5698 else
5699 {
5700 for (k = M + 3; k >= M; k--)
5701 yi[k] = *di++;
5702 }
842fbaaa
JW		 5703 yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
5704 if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
5705 ecleaz (yi); /* it was zero */
5706 else
5707 yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
5708 emovo (yi, e);
5709}
5710
5711/* Convert target computer signed 64-bit integer to e-type. */
5712
a0353055 5713static void
842fbaaa
JW		 5714ditoe (di, e)
5715 unsigned EMUSHORT *di; /* Address of the 64-bit int. */
5716 unsigned EMUSHORT *e;
5717{
5718 unsigned EMULONG acc;
5719 unsigned EMUSHORT yi[NI];
5720 unsigned EMUSHORT carry;
5721 int k, sign;
5722
5723 ecleaz (yi);
f76b9db2
ILT		 5724 if (WORDS_BIG_ENDIAN)
5725 {
5726 for (k = M; k < M + 4; k++)
5727 yi[k] = *di++;
5728 }
5729 else
5730 {
5731 for (k = M + 3; k >= M; k--)
5732 yi[k] = *di++;
5733 }
842fbaaa
JW		 5734 /* Take absolute value */
5735 sign = 0;
5736 if (yi[M] & 0x8000)
5737 {
5738 sign = 1;
5739 carry = 0;
5740 for (k = M + 3; k >= M; k--)
5741 {
5742 acc = (unsigned EMULONG) (~yi[k] & 0xffff) + carry;
5743 yi[k] = acc;
5744 carry = 0;
5745 if (acc & 0x10000)
5746 carry = 1;
5747 }
5748 }
5749 yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
5750 if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
5751 ecleaz (yi); /* it was zero */
5752 else
5753 yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
5754 emovo (yi, e);
5755 if (sign)
5756 eneg (e);
5757}
5758
5759
5760/* Convert e-type to unsigned 64-bit int. */
5761
008f0d36
RK		 5762static void
5763etoudi (x, i)
842fbaaa
JW		 5764 unsigned EMUSHORT *x;
5765 unsigned EMUSHORT *i;
5766{
5767 unsigned EMUSHORT xi[NI];
5768 int j, k;
5769
5770 emovi (x, xi);
5771 if (xi[0])
5772 {
5773 xi[M] = 0;
5774 goto noshift;
5775 }
5776 k = (int) xi[E] - (EXONE - 1);
5777 if (k <= 0)
5778 {
5779 for (j = 0; j < 4; j++)
5780 *i++ = 0;
5781 return;
5782 }
5783 if (k > 64)
5784 {
5785 for (j = 0; j < 4; j++)
5786 *i++ = 0xffff;
5787 if (extra_warnings)
5788 warning ("overflow on truncation to integer");
5789 return;
5790 }
5791 if (k > 16)
5792 {
5793 /* Shift more than 16 bits: first shift up k-16 mod 16,
5794 then shift up by 16's. */
5795 j = k - ((k >> 4) << 4);
5796 if (j == 0)
5797 j = 16;
5798 eshift (xi, j);
f76b9db2
ILT		 5799 if (WORDS_BIG_ENDIAN)
5800 *i++ = xi[M];
5801 else
5802 {
5803 i += 3;
5804 *i-- = xi[M];
5805 }
842fbaaa
JW		 5806 k -= j;
5807 do
5808 {
5809 eshup6 (xi);
f76b9db2
ILT		 5810 if (WORDS_BIG_ENDIAN)
5811 *i++ = xi[M];
5812 else
5813 *i-- = xi[M];
842fbaaa
JW		 5814 }
5815 while ((k -= 16) > 0);
5816 }
5817 else
5818 {
5819 /* shift not more than 16 bits */
5820 eshift (xi, k);
5821
5822noshift:
5823
f76b9db2
ILT		 5824 if (WORDS_BIG_ENDIAN)
5825 {
5826 i += 3;
5827 *i-- = xi[M];
5828 *i-- = 0;
5829 *i-- = 0;
5830 *i = 0;
5831 }
5832 else
5833 {
5834 *i++ = xi[M];
5835 *i++ = 0;
5836 *i++ = 0;
5837 *i = 0;
5838 }
842fbaaa
JW		 5839 }
5840}
5841
5842
5843/* Convert e-type to signed 64-bit int. */
5844
a0353055 5845static void
842fbaaa
JW		 5846etodi (x, i)
5847 unsigned EMUSHORT *x;
5848 unsigned EMUSHORT *i;
5849{
5850 unsigned EMULONG acc;
5851 unsigned EMUSHORT xi[NI];
5852 unsigned EMUSHORT carry;
5853 unsigned EMUSHORT *isave;
5854 int j, k;
5855
5856 emovi (x, xi);
5857 k = (int) xi[E] - (EXONE - 1);
5858 if (k <= 0)
5859 {
5860 for (j = 0; j < 4; j++)
5861 *i++ = 0;
5862 return;
5863 }
5864 if (k > 64)
5865 {
5866 for (j = 0; j < 4; j++)
5867 *i++ = 0xffff;
5868 if (extra_warnings)
5869 warning ("overflow on truncation to integer");
5870 return;
5871 }
5872 isave = i;
5873 if (k > 16)
5874 {
5875 /* Shift more than 16 bits: first shift up k-16 mod 16,
5876 then shift up by 16's. */
5877 j = k - ((k >> 4) << 4);
5878 if (j == 0)
5879 j = 16;
5880 eshift (xi, j);
f76b9db2
ILT		 5881 if (WORDS_BIG_ENDIAN)
5882 *i++ = xi[M];
5883 else
5884 {
5885 i += 3;
5886 *i-- = xi[M];
5887 }
842fbaaa
JW		 5888 k -= j;
5889 do
5890 {
5891 eshup6 (xi);
f76b9db2
ILT		 5892 if (WORDS_BIG_ENDIAN)
5893 *i++ = xi[M];
5894 else
5895 *i-- = xi[M];
842fbaaa
JW		 5896 }
5897 while ((k -= 16) > 0);
5898 }
5899 else
5900 {
5901 /* shift not more than 16 bits */
5902 eshift (xi, k);
5903
f76b9db2
ILT		 5904 if (WORDS_BIG_ENDIAN)
5905 {
5906 i += 3;
5907 *i = xi[M];
5908 *i-- = 0;
5909 *i-- = 0;
5910 *i = 0;
5911 }
5912 else
5913 {
5914 *i++ = xi[M];
5915 *i++ = 0;
5916 *i++ = 0;
5917 *i = 0;
5918 }
842fbaaa
JW		 5919 }
5920 /* Negate if negative */
5921 if (xi[0])
5922 {
5923 carry = 0;
f76b9db2
ILT		 5924 if (WORDS_BIG_ENDIAN)
5925 isave += 3;
842fbaaa
JW		 5926 for (k = 0; k < 4; k++)
5927 {
5928 acc = (unsigned EMULONG) (~(*isave) & 0xffff) + carry;
f76b9db2
ILT		 5929 if (WORDS_BIG_ENDIAN)
5930 *isave-- = acc;
5931 else
5932 *isave++ = acc;
842fbaaa
JW		 5933 carry = 0;
5934 if (acc & 0x10000)
5935 carry = 1;
5936 }
5937 }
5938}
5939
5940
5941/* Longhand square root routine. */
5942
5943
5944static int esqinited = 0;
5945static unsigned short sqrndbit[NI];
5946
a0353055 5947static void
842fbaaa
JW		 5948esqrt (x, y)
5949 unsigned EMUSHORT *x, *y;
5950{
5951 unsigned EMUSHORT temp[NI], num[NI], sq[NI], xx[NI];
5952 EMULONG m, exp;
5953 int i, j, k, n, nlups;
5954
5955 if (esqinited == 0)
5956 {
5957 ecleaz (sqrndbit);
5958 sqrndbit[NI - 2] = 1;
5959 esqinited = 1;
5960 }
5961 /* Check for arg <= 0 */
5962 i = ecmp (x, ezero);
5963 if (i <= 0)
5964 {
29e11dab 5965 if (i == -1)
842fbaaa 5966 {
29e11dab
RK		 5967 mtherr ("esqrt", DOMAIN);
5968 eclear (y);
842fbaaa 5969 }
29e11dab
RK		 5970 else
5971 emov (x, y);
842fbaaa
JW		 5972 return;
5973 }
5974
5975#ifdef INFINITY
5976 if (eisinf (x))
5977 {
5978 eclear (y);
5979 einfin (y);
5980 return;
5981 }
5982#endif
5983 /* Bring in the arg and renormalize if it is denormal. */
5984 emovi (x, xx);
5985 m = (EMULONG) xx[1]; /* local long word exponent */
5986 if (m == 0)
5987 m -= enormlz (xx);
5988
5989 /* Divide exponent by 2 */
5990 m -= 0x3ffe;
5991 exp = (unsigned short) ((m / 2) + 0x3ffe);
5992
5993 /* Adjust if exponent odd */
5994 if ((m & 1) != 0)
5995 {
5996 if (m > 0)
5997 exp += 1;
5998 eshdn1 (xx);
5999 }
6000
6001 ecleaz (sq);
6002 ecleaz (num);
6003 n = 8; /* get 8 bits of result per inner loop */
6004 nlups = rndprc;
6005 j = 0;
6006
6007 while (nlups > 0)
6008 {
6009 /* bring in next word of arg */
6010 if (j < NE)
6011 num[NI - 1] = xx[j + 3];
6012 /* Do additional bit on last outer loop, for roundoff. */
6013 if (nlups <= 8)
6014 n = nlups + 1;
6015 for (i = 0; i < n; i++)
6016 {
6017 /* Next 2 bits of arg */
6018 eshup1 (num);
6019 eshup1 (num);
6020 /* Shift up answer */
6021 eshup1 (sq);
6022 /* Make trial divisor */
6023 for (k = 0; k < NI; k++)
6024 temp[k] = sq[k];
6025 eshup1 (temp);
6026 eaddm (sqrndbit, temp);
6027 /* Subtract and insert answer bit if it goes in */
6028 if (ecmpm (temp, num) <= 0)
6029 {
6030 esubm (temp, num);
6031 sq[NI - 2] |= 1;
6032 }
6033 }
6034 nlups -= n;
6035 j += 1;
6036 }
6037
6038 /* Adjust for extra, roundoff loop done. */
6039 exp += (NBITS - 1) - rndprc;
6040
6041 /* Sticky bit = 1 if the remainder is nonzero. */
6042 k = 0;
6043 for (i = 3; i < NI; i++)
6044 k |= (int) num[i];
6045
6046 /* Renormalize and round off. */
6047 emdnorm (sq, k, 0, exp, 64);
6048 emovo (sq, y);
6049}
985b6196 6050#endif /* EMU_NON_COMPILE not defined */
8ddae348
RK		 6051\f
6052/* Return the binary precision of the significand for a given
6053 floating point mode. The mode can hold an integer value
6054 that many bits wide, without losing any bits. */
6055
6056int
6057significand_size (mode)
6058 enum machine_mode mode;
6059{
6060
6061switch (mode)
6062 {
6063 case SFmode:
6064 return 24;
6065
6066 case DFmode:
6067#if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
6068 return 53;
6069#else
6070#if TARGET_FLOAT_FORMAT == IBM_FLOAT_FORMAT
6071 return 56;
6072#else
6073#if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
6074 return 56;
6075#else
6076 abort ();
6077#endif
6078#endif
6079#endif
6080
6081 case XFmode:
6082 return 64;
6083 case TFmode:
6084 return 113;
6085
6086 default:
6087 abort ();
6088 }
6089}
GNU Compiler Collection
This page took 1.187608 seconds and 5 git commands to generate.