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