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