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