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