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