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