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