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