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