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