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