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