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