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