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