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