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