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