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