1 /* real.c - software floating point emulation.
2 Copyright (C) 1993, 1994, 1995, 1996, 1997, 1998, 1999,
3 2000, 2002, 2003, 2004, 2005 Free Software Foundation, Inc.
4 Contributed by Stephen L. Moshier (moshier@world.std.com).
5 Re-written by Richard Henderson <rth@redhat.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, 51 Franklin Street, Fifth Floor, Boston, MA
22 02110-1301, USA. */
23
24 #include "config.h"
25 #include "system.h"
26 #include "coretypes.h"
27 #include "tm.h"
28 #include "tree.h"
29 #include "toplev.h"
30 #include "real.h"
31 #include "tm_p.h"
32
33 /* The floating point model used internally is not exactly IEEE 754
34 compliant, and close to the description in the ISO C99 standard,
35 section 5.2.4.2.2 Characteristics of floating types.
36
37 Specifically
38
39 x = s * b^e * \sum_{k=1}^p f_k * b^{-k}
40
41 where
42 s = sign (+- 1)
43 b = base or radix, here always 2
44 e = exponent
45 p = precision (the number of base-b digits in the significand)
46 f_k = the digits of the significand.
47
48 We differ from typical IEEE 754 encodings in that the entire
49 significand is fractional. Normalized significands are in the
50 range [0.5, 1.0).
51
52 A requirement of the model is that P be larger than the largest
53 supported target floating-point type by at least 2 bits. This gives
54 us proper rounding when we truncate to the target type. In addition,
55 E must be large enough to hold the smallest supported denormal number
56 in a normalized form.
57
58 Both of these requirements are easily satisfied. The largest target
59 significand is 113 bits; we store at least 160. The smallest
60 denormal number fits in 17 exponent bits; we store 27.
61
62 Note that the decimal string conversion routines are sensitive to
63 rounding errors. Since the raw arithmetic routines do not themselves
64 have guard digits or rounding, the computation of 10**exp can
65 accumulate more than a few digits of error. The previous incarnation
66 of real.c successfully used a 144-bit fraction; given the current
67 layout of REAL_VALUE_TYPE we're forced to expand to at least 160 bits.
68
69 Target floating point models that use base 16 instead of base 2
70 (i.e. IBM 370), are handled during round_for_format, in which we
71 canonicalize the exponent to be a multiple of 4 (log2(16)), and
72 adjust the significand to match. */
73
74
75 #if 0
76 /* Used to classify two numbers simultaneously. */
77 #define CLASS2(A, B) ((A) << 2 | (B))
78
79 #if HOST_BITS_PER_LONG != 64 && HOST_BITS_PER_LONG != 32
80 #error "Some constant folding done by hand to avoid shift count warnings"
81 #endif
82
83 static void get_zero (REAL_VALUE_TYPE *, int);
84 static void get_canonical_qnan (REAL_VALUE_TYPE *, int);
85 static void get_canonical_snan (REAL_VALUE_TYPE *, int);
86 static void get_inf (REAL_VALUE_TYPE *, int);
87 static bool sticky_rshift_significand (REAL_VALUE_TYPE *,
88 const REAL_VALUE_TYPE *, unsigned int);
89 static void rshift_significand (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
90 unsigned int);
91 static void lshift_significand (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
92 unsigned int);
93 static void lshift_significand_1 (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *);
94 static bool add_significands (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *,
95 const REAL_VALUE_TYPE *);
96 static bool sub_significands (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
97 const REAL_VALUE_TYPE *, int);
98 static void neg_significand (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *);
99 static int cmp_significands (const REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *);
100 static int cmp_significand_0 (const REAL_VALUE_TYPE *);
101 static void set_significand_bit (REAL_VALUE_TYPE *, unsigned int);
102 static void clear_significand_bit (REAL_VALUE_TYPE *, unsigned int);
103 static bool test_significand_bit (REAL_VALUE_TYPE *, unsigned int);
104 static void clear_significand_below (REAL_VALUE_TYPE *, unsigned int);
105 static bool div_significands (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
106 const REAL_VALUE_TYPE *);
107 static void normalize (REAL_VALUE_TYPE *);
108
109 static bool do_add (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
110 const REAL_VALUE_TYPE *, int);
111 static bool do_multiply (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
112 const REAL_VALUE_TYPE *);
113 static bool do_divide (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
114 const REAL_VALUE_TYPE *);
115 static int do_compare (const REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *, int);
116 static void do_fix_trunc (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *);
117
118 static unsigned long rtd_divmod (REAL_VALUE_TYPE *, REAL_VALUE_TYPE *);
119
120 static const REAL_VALUE_TYPE * ten_to_ptwo (int);
121 static const REAL_VALUE_TYPE * ten_to_mptwo (int);
122 static const REAL_VALUE_TYPE * real_digit (int);
123 static void times_pten (REAL_VALUE_TYPE *, int);
124
125 static void round_for_format (const struct real_format *, REAL_VALUE_TYPE *);
126 #endif /* 0 */
127
128 #if 0
129 /* Initialize R with a positive zero. */
130
131 static inline void
132 get_zero (REAL_VALUE_TYPE *r, int sign)
133 {
134 memset (r, 0, sizeof (*r));
135 r->sign = sign;
136 }
137
138 /* Initialize R with the canonical quiet NaN. */
139
140 static inline void
141 get_canonical_qnan (REAL_VALUE_TYPE *r, int sign)
142 {
143 memset (r, 0, sizeof (*r));
144 r->cl = rvc_nan;
145 r->sign = sign;
146 r->canonical = 1;
147 }
148
149 static inline void
150 get_canonical_snan (REAL_VALUE_TYPE *r, int sign)
151 {
152 memset (r, 0, sizeof (*r));
153 r->cl = rvc_nan;
154 r->sign = sign;
155 r->signalling = 1;
156 r->canonical = 1;
157 }
158
159 static inline void
160 get_inf (REAL_VALUE_TYPE *r, int sign)
161 {
162 memset (r, 0, sizeof (*r));
163 r->cl = rvc_inf;
164 r->sign = sign;
165 }
166
167
168 /* Right-shift the significand of A by N bits; put the result in the
169 significand of R. If any one bits are shifted out, return true. */
170
171 static bool
172 sticky_rshift_significand (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
173 unsigned int n)
174 {
175 unsigned long sticky = 0;
176 unsigned int i, ofs = 0;
177
178 if (n >= HOST_BITS_PER_LONG)
179 {
180 for (i = 0, ofs = n / HOST_BITS_PER_LONG; i < ofs; ++i)
181 sticky |= a->sig[i];
182 n &= HOST_BITS_PER_LONG - 1;
183 }
184
185 if (n != 0)
186 {
187 sticky |= a->sig[ofs] & (((unsigned long)1 << n) - 1);
188 for (i = 0; i < SIGSZ; ++i)
189 {
190 r->sig[i]
191 = (((ofs + i >= SIGSZ ? 0 : a->sig[ofs + i]) >> n)
192 | ((ofs + i + 1 >= SIGSZ ? 0 : a->sig[ofs + i + 1])
193 << (HOST_BITS_PER_LONG - n)));
194 }
195 }
196 else
197 {
198 for (i = 0; ofs + i < SIGSZ; ++i)
199 r->sig[i] = a->sig[ofs + i];
200 for (; i < SIGSZ; ++i)
201 r->sig[i] = 0;
202 }
203
204 return sticky != 0;
205 }
206
207 /* Right-shift the significand of A by N bits; put the result in the
208 significand of R. */
209
210 static void
211 rshift_significand (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
212 unsigned int n)
213 {
214 unsigned int i, ofs = n / HOST_BITS_PER_LONG;
215
216 n &= HOST_BITS_PER_LONG - 1;
217 if (n != 0)
218 {
219 for (i = 0; i < SIGSZ; ++i)
220 {
221 r->sig[i]
222 = (((ofs + i >= SIGSZ ? 0 : a->sig[ofs + i]) >> n)
223 | ((ofs + i + 1 >= SIGSZ ? 0 : a->sig[ofs + i + 1])
224 << (HOST_BITS_PER_LONG - n)));
225 }
226 }
227 else
228 {
229 for (i = 0; ofs + i < SIGSZ; ++i)
230 r->sig[i] = a->sig[ofs + i];
231 for (; i < SIGSZ; ++i)
232 r->sig[i] = 0;
233 }
234 }
235
236 /* Left-shift the significand of A by N bits; put the result in the
237 significand of R. */
238
239 static void
240 lshift_significand (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
241 unsigned int n)
242 {
243 unsigned int i, ofs = n / HOST_BITS_PER_LONG;
244
245 n &= HOST_BITS_PER_LONG - 1;
246 if (n == 0)
247 {
248 for (i = 0; ofs + i < SIGSZ; ++i)
249 r->sig[SIGSZ-1-i] = a->sig[SIGSZ-1-i-ofs];
250 for (; i < SIGSZ; ++i)
251 r->sig[SIGSZ-1-i] = 0;
252 }
253 else
254 for (i = 0; i < SIGSZ; ++i)
255 {
256 r->sig[SIGSZ-1-i]
257 = (((ofs + i >= SIGSZ ? 0 : a->sig[SIGSZ-1-i-ofs]) << n)
258 | ((ofs + i + 1 >= SIGSZ ? 0 : a->sig[SIGSZ-1-i-ofs-1])
259 >> (HOST_BITS_PER_LONG - n)));
260 }
261 }
262
263 /* Likewise, but N is specialized to 1. */
264
265 static inline void
266 lshift_significand_1 (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a)
267 {
268 unsigned int i;
269
270 for (i = SIGSZ - 1; i > 0; --i)
271 r->sig[i] = (a->sig[i] << 1) | (a->sig[i-1] >> (HOST_BITS_PER_LONG - 1));
272 r->sig[0] = a->sig[0] << 1;
273 }
274
275 /* Add the significands of A and B, placing the result in R. Return
276 true if there was carry out of the most significant word. */
277
278 static inline bool
279 add_significands (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
280 const REAL_VALUE_TYPE *b)
281 {
282 bool carry = false;
283 int i;
284
285 for (i = 0; i < SIGSZ; ++i)
286 {
287 unsigned long ai = a->sig[i];
288 unsigned long ri = ai + b->sig[i];
289
290 if (carry)
291 {
292 carry = ri < ai;
293 carry |= ++ri == 0;
294 }
295 else
296 carry = ri < ai;
297
298 r->sig[i] = ri;
299 }
300
301 return carry;
302 }
303
304 /* Subtract the significands of A and B, placing the result in R. CARRY is
305 true if there's a borrow incoming to the least significant word.
306 Return true if there was borrow out of the most significant word. */
307
308 static inline bool
309 sub_significands (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
310 const REAL_VALUE_TYPE *b, int carry)
311 {
312 int i;
313
314 for (i = 0; i < SIGSZ; ++i)
315 {
316 unsigned long ai = a->sig[i];
317 unsigned long ri = ai - b->sig[i];
318
319 if (carry)
320 {
321 carry = ri > ai;
322 carry |= ~--ri == 0;
323 }
324 else
325 carry = ri > ai;
326
327 r->sig[i] = ri;
328 }
329
330 return carry;
331 }
332
333 /* Negate the significand A, placing the result in R. */
334
335 static inline void
336 neg_significand (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a)
337 {
338 bool carry = true;
339 int i;
340
341 for (i = 0; i < SIGSZ; ++i)
342 {
343 unsigned long ri, ai = a->sig[i];
344
345 if (carry)
346 {
347 if (ai)
348 {
349 ri = -ai;
350 carry = false;
351 }
352 else
353 ri = ai;
354 }
355 else
356 ri = ~ai;
357
358 r->sig[i] = ri;
359 }
360 }
361
362 /* Compare significands. Return tri-state vs zero. */
363
364 static inline int
365 cmp_significands (const REAL_VALUE_TYPE *a, const REAL_VALUE_TYPE *b)
366 {
367 int i;
368
369 for (i = SIGSZ - 1; i >= 0; --i)
370 {
371 unsigned long ai = a->sig[i];
372 unsigned long bi = b->sig[i];
373
374 if (ai > bi)
375 return 1;
376 if (ai < bi)
377 return -1;
378 }
379
380 return 0;
381 }
382
383 /* Return true if A is nonzero. */
384
385 static inline int
386 cmp_significand_0 (const REAL_VALUE_TYPE *a)
387 {
388 int i;
389
390 for (i = SIGSZ - 1; i >= 0; --i)
391 if (a->sig[i])
392 return 1;
393
394 return 0;
395 }
396
397 /* Set bit N of the significand of R. */
398
399 static inline void
400 set_significand_bit (REAL_VALUE_TYPE *r, unsigned int n)
401 {
402 r->sig[n / HOST_BITS_PER_LONG]
403 |= (unsigned long)1 << (n % HOST_BITS_PER_LONG);
404 }
405
406 /* Clear bit N of the significand of R. */
407
408 static inline void
409 clear_significand_bit (REAL_VALUE_TYPE *r, unsigned int n)
410 {
411 r->sig[n / HOST_BITS_PER_LONG]
412 &= ~((unsigned long)1 << (n % HOST_BITS_PER_LONG));
413 }
414
415 /* Test bit N of the significand of R. */
416
417 static inline bool
418 test_significand_bit (REAL_VALUE_TYPE *r, unsigned int n)
419 {
420 /* ??? Compiler bug here if we return this expression directly.
421 The conversion to bool strips the "&1" and we wind up testing
422 e.g. 2 != 0 -> true. Seen in gcc version 3.2 20020520. */
423 int t = (r->sig[n / HOST_BITS_PER_LONG] >> (n % HOST_BITS_PER_LONG)) & 1;
424 return t;
425 }
426
427 /* Clear bits 0..N-1 of the significand of R. */
428
429 static void
430 clear_significand_below (REAL_VALUE_TYPE *r, unsigned int n)
431 {
432 int i, w = n / HOST_BITS_PER_LONG;
433
434 for (i = 0; i < w; ++i)
435 r->sig[i] = 0;
436
437 r->sig[w] &= ~(((unsigned long)1 << (n % HOST_BITS_PER_LONG)) - 1);
438 }
439
440 /* Divide the significands of A and B, placing the result in R. Return
441 true if the division was inexact. */
442
443 static inline bool
444 div_significands (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
445 const REAL_VALUE_TYPE *b)
446 {
447 REAL_VALUE_TYPE u;
448 int i, bit = SIGNIFICAND_BITS - 1;
449 unsigned long msb, inexact;
450
451 u = *a;
452 memset (r->sig, 0, sizeof (r->sig));
453
454 msb = 0;
455 goto start;
456 do
457 {
458 msb = u.sig[SIGSZ-1] & SIG_MSB;
459 lshift_significand_1 (&u, &u);
460 start:
461 if (msb || cmp_significands (&u, b) >= 0)
462 {
463 sub_significands (&u, &u, b, 0);
464 set_significand_bit (r, bit);
465 }
466 }
467 while (--bit >= 0);
468
469 for (i = 0, inexact = 0; i < SIGSZ; i++)
470 inexact |= u.sig[i];
471
472 return inexact != 0;
473 }
474
475 /* Adjust the exponent and significand of R such that the most
476 significant bit is set. We underflow to zero and overflow to
477 infinity here, without denormals. (The intermediate representation
478 exponent is large enough to handle target denormals normalized.) */
479
480 static void
481 normalize (REAL_VALUE_TYPE *r)
482 {
483 int shift = 0, exp;
484 int i, j;
485
486 /* Find the first word that is nonzero. */
487 for (i = SIGSZ - 1; i >= 0; i--)
488 if (r->sig[i] == 0)
489 shift += HOST_BITS_PER_LONG;
490 else
491 break;
492
493 /* Zero significand flushes to zero. */
494 if (i < 0)
495 {
496 r->cl = rvc_zero;
497 SET_REAL_EXP (r, 0);
498 return;
499 }
500
501 /* Find the first bit that is nonzero. */
502 for (j = 0; ; j++)
503 if (r->sig[i] & ((unsigned long)1 << (HOST_BITS_PER_LONG - 1 - j)))
504 break;
505 shift += j;
506
507 if (shift > 0)
508 {
509 exp = REAL_EXP (r) - shift;
510 if (exp > MAX_EXP)
511 get_inf (r, r->sign);
512 else if (exp < -MAX_EXP)
513 get_zero (r, r->sign);
514 else
515 {
516 SET_REAL_EXP (r, exp);
517 lshift_significand (r, r, shift);
518 }
519 }
520 }
521
522 /* Calculate R = A + (SUBTRACT_P ? -B : B). Return true if the
523 result may be inexact due to a loss of precision. */
524
525 static bool
526 do_add (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
527 const REAL_VALUE_TYPE *b, int subtract_p)
528 {
529 int dexp, sign, exp;
530 REAL_VALUE_TYPE t;
531 bool inexact = false;
532
533 /* Determine if we need to add or subtract. */
534 sign = a->sign;
535 subtract_p = (sign ^ b->sign) ^ subtract_p;
536
537 switch (CLASS2 (a->cl, b->cl))
538 {
539 case CLASS2 (rvc_zero, rvc_zero):
540 /* -0 + -0 = -0, -0 - +0 = -0; all other cases yield +0. */
541 get_zero (r, sign & !subtract_p);
542 return false;
543
544 case CLASS2 (rvc_zero, rvc_normal):
545 case CLASS2 (rvc_zero, rvc_inf):
546 case CLASS2 (rvc_zero, rvc_nan):
547 /* 0 + ANY = ANY. */
548 case CLASS2 (rvc_normal, rvc_nan):
549 case CLASS2 (rvc_inf, rvc_nan):
550 case CLASS2 (rvc_nan, rvc_nan):
551 /* ANY + NaN = NaN. */
552 case CLASS2 (rvc_normal, rvc_inf):
553 /* R + Inf = Inf. */
554 *r = *b;
555 r->sign = sign ^ subtract_p;
556 return false;
557
558 case CLASS2 (rvc_normal, rvc_zero):
559 case CLASS2 (rvc_inf, rvc_zero):
560 case CLASS2 (rvc_nan, rvc_zero):
561 /* ANY + 0 = ANY. */
562 case CLASS2 (rvc_nan, rvc_normal):
563 case CLASS2 (rvc_nan, rvc_inf):
564 /* NaN + ANY = NaN. */
565 case CLASS2 (rvc_inf, rvc_normal):
566 /* Inf + R = Inf. */
567 *r = *a;
568 return false;
569
570 case CLASS2 (rvc_inf, rvc_inf):
571 if (subtract_p)
572 /* Inf - Inf = NaN. */
573 get_canonical_qnan (r, 0);
574 else
575 /* Inf + Inf = Inf. */
576 *r = *a;
577 return false;
578
579 case CLASS2 (rvc_normal, rvc_normal):
580 break;
581
582 default:
583 gcc_unreachable ();
584 }
585
586 /* Swap the arguments such that A has the larger exponent. */
587 dexp = REAL_EXP (a) - REAL_EXP (b);
588 if (dexp < 0)
589 {
590 const REAL_VALUE_TYPE *t;
591 t = a, a = b, b = t;
592 dexp = -dexp;
593 sign ^= subtract_p;
594 }
595 exp = REAL_EXP (a);
596
597 /* If the exponents are not identical, we need to shift the
598 significand of B down. */
599 if (dexp > 0)
600 {
601 /* If the exponents are too far apart, the significands
602 do not overlap, which makes the subtraction a noop. */
603 if (dexp >= SIGNIFICAND_BITS)
604 {
605 *r = *a;
606 r->sign = sign;
607 return true;
608 }
609
610 inexact |= sticky_rshift_significand (&t, b, dexp);
611 b = &t;
612 }
613
614 if (subtract_p)
615 {
616 if (sub_significands (r, a, b, inexact))
617 {
618 /* We got a borrow out of the subtraction. That means that
619 A and B had the same exponent, and B had the larger
620 significand. We need to swap the sign and negate the
621 significand. */
622 sign ^= 1;
623 neg_significand (r, r);
624 }
625 }
626 else
627 {
628 if (add_significands (r, a, b))
629 {
630 /* We got carry out of the addition. This means we need to
631 shift the significand back down one bit and increase the
632 exponent. */
633 inexact |= sticky_rshift_significand (r, r, 1);
634 r->sig[SIGSZ-1] |= SIG_MSB;
635 if (++exp > MAX_EXP)
636 {
637 get_inf (r, sign);
638 return true;
639 }
640 }
641 }
642
643 r->cl = rvc_normal;
644 r->sign = sign;
645 SET_REAL_EXP (r, exp);
646 /* Zero out the remaining fields. */
647 r->signalling = 0;
648 r->canonical = 0;
649
650 /* Re-normalize the result. */
651 normalize (r);
652
653 /* Special case: if the subtraction results in zero, the result
654 is positive. */
655 if (r->cl == rvc_zero)
656 r->sign = 0;
657 else
658 r->sig[0] |= inexact;
659
660 return inexact;
661 }
662
663 /* Calculate R = A * B. Return true if the result may be inexact. */
664
665 static bool
666 do_multiply (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
667 const REAL_VALUE_TYPE *b)
668 {
669 REAL_VALUE_TYPE u, t, *rr;
670 unsigned int i, j, k;
671 int sign = a->sign ^ b->sign;
672 bool inexact = false;
673
674 switch (CLASS2 (a->cl, b->cl))
675 {
676 case CLASS2 (rvc_zero, rvc_zero):
677 case CLASS2 (rvc_zero, rvc_normal):
678 case CLASS2 (rvc_normal, rvc_zero):
679 /* +-0 * ANY = 0 with appropriate sign. */
680 get_zero (r, sign);
681 return false;
682
683 case CLASS2 (rvc_zero, rvc_nan):
684 case CLASS2 (rvc_normal, rvc_nan):
685 case CLASS2 (rvc_inf, rvc_nan):
686 case CLASS2 (rvc_nan, rvc_nan):
687 /* ANY * NaN = NaN. */
688 *r = *b;
689 r->sign = sign;
690 return false;
691
692 case CLASS2 (rvc_nan, rvc_zero):
693 case CLASS2 (rvc_nan, rvc_normal):
694 case CLASS2 (rvc_nan, rvc_inf):
695 /* NaN * ANY = NaN. */
696 *r = *a;
697 r->sign = sign;
698 return false;
699
700 case CLASS2 (rvc_zero, rvc_inf):
701 case CLASS2 (rvc_inf, rvc_zero):
702 /* 0 * Inf = NaN */
703 get_canonical_qnan (r, sign);
704 return false;
705
706 case CLASS2 (rvc_inf, rvc_inf):
707 case CLASS2 (rvc_normal, rvc_inf):
708 case CLASS2 (rvc_inf, rvc_normal):
709 /* Inf * Inf = Inf, R * Inf = Inf */
710 get_inf (r, sign);
711 return false;
712
713 case CLASS2 (rvc_normal, rvc_normal):
714 break;
715
716 default:
717 gcc_unreachable ();
718 }
719
720 if (r == a || r == b)
721 rr = &t;
722 else
723 rr = r;
724 get_zero (rr, 0);
725
726 /* Collect all the partial products. Since we don't have sure access
727 to a widening multiply, we split each long into two half-words.
728
729 Consider the long-hand form of a four half-word multiplication:
730
731 A B C D
732 * E F G H
733 --------------
734 DE DF DG DH
735 CE CF CG CH
736 BE BF BG BH
737 AE AF AG AH
738
739 We construct partial products of the widened half-word products
740 that are known to not overlap, e.g. DF+DH. Each such partial
741 product is given its proper exponent, which allows us to sum them
742 and obtain the finished product. */
743
744 for (i = 0; i < SIGSZ * 2; ++i)
745 {
746 unsigned long ai = a->sig[i / 2];
747 if (i & 1)
748 ai >>= HOST_BITS_PER_LONG / 2;
749 else
750 ai &= ((unsigned long)1 << (HOST_BITS_PER_LONG / 2)) - 1;
751
752 if (ai == 0)
753 continue;
754
755 for (j = 0; j < 2; ++j)
756 {
757 int exp = (REAL_EXP (a) - (2*SIGSZ-1-i)*(HOST_BITS_PER_LONG/2)
758 + (REAL_EXP (b) - (1-j)*(HOST_BITS_PER_LONG/2)));
759
760 if (exp > MAX_EXP)
761 {
762 get_inf (r, sign);
763 return true;
764 }
765 if (exp < -MAX_EXP)
766 {
767 /* Would underflow to zero, which we shouldn't bother adding. */
768 inexact = true;
769 continue;
770 }
771
772 memset (&u, 0, sizeof (u));
773 u.cl = rvc_normal;
774 SET_REAL_EXP (&u, exp);
775
776 for (k = j; k < SIGSZ * 2; k += 2)
777 {
778 unsigned long bi = b->sig[k / 2];
779 if (k & 1)
780 bi >>= HOST_BITS_PER_LONG / 2;
781 else
782 bi &= ((unsigned long)1 << (HOST_BITS_PER_LONG / 2)) - 1;
783
784 u.sig[k / 2] = ai * bi;
785 }
786
787 normalize (&u);
788 inexact |= do_add (rr, rr, &u, 0);
789 }
790 }
791
792 rr->sign = sign;
793 if (rr != r)
794 *r = t;
795
796 return inexact;
797 }
798
799 /* Calculate R = A / B. Return true if the result may be inexact. */
800
801 static bool
802 do_divide (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
803 const REAL_VALUE_TYPE *b)
804 {
805 int exp, sign = a->sign ^ b->sign;
806 REAL_VALUE_TYPE t, *rr;
807 bool inexact;
808
809 switch (CLASS2 (a->cl, b->cl))
810 {
811 case CLASS2 (rvc_zero, rvc_zero):
812 /* 0 / 0 = NaN. */
813 case CLASS2 (rvc_inf, rvc_inf):
814 /* Inf / Inf = NaN. */
815 get_canonical_qnan (r, sign);
816 return false;
817
818 case CLASS2 (rvc_zero, rvc_normal):
819 case CLASS2 (rvc_zero, rvc_inf):
820 /* 0 / ANY = 0. */
821 case CLASS2 (rvc_normal, rvc_inf):
822 /* R / Inf = 0. */
823 get_zero (r, sign);
824 return false;
825
826 case CLASS2 (rvc_normal, rvc_zero):
827 /* R / 0 = Inf. */
828 case CLASS2 (rvc_inf, rvc_zero):
829 /* Inf / 0 = Inf. */
830 get_inf (r, sign);
831 return false;
832
833 case CLASS2 (rvc_zero, rvc_nan):
834 case CLASS2 (rvc_normal, rvc_nan):
835 case CLASS2 (rvc_inf, rvc_nan):
836 case CLASS2 (rvc_nan, rvc_nan):
837 /* ANY / NaN = NaN. */
838 *r = *b;
839 r->sign = sign;
840 return false;
841
842 case CLASS2 (rvc_nan, rvc_zero):
843 case CLASS2 (rvc_nan, rvc_normal):
844 case CLASS2 (rvc_nan, rvc_inf):
845 /* NaN / ANY = NaN. */
846 *r = *a;
847 r->sign = sign;
848 return false;
849
850 case CLASS2 (rvc_inf, rvc_normal):
851 /* Inf / R = Inf. */
852 get_inf (r, sign);
853 return false;
854
855 case CLASS2 (rvc_normal, rvc_normal):
856 break;
857
858 default:
859 gcc_unreachable ();
860 }
861
862 if (r == a || r == b)
863 rr = &t;
864 else
865 rr = r;
866
867 /* Make sure all fields in the result are initialized. */
868 get_zero (rr, 0);
869 rr->cl = rvc_normal;
870 rr->sign = sign;
871
872 exp = REAL_EXP (a) - REAL_EXP (b) + 1;
873 if (exp > MAX_EXP)
874 {
875 get_inf (r, sign);
876 return true;
877 }
878 if (exp < -MAX_EXP)
879 {
880 get_zero (r, sign);
881 return true;
882 }
883 SET_REAL_EXP (rr, exp);
884
885 inexact = div_significands (rr, a, b);
886
887 /* Re-normalize the result. */
888 normalize (rr);
889 rr->sig[0] |= inexact;
890
891 if (rr != r)
892 *r = t;
893
894 return inexact;
895 }
896
897 /* Return a tri-state comparison of A vs B. Return NAN_RESULT if
898 one of the two operands is a NaN. */
899
900 static int
901 do_compare (const REAL_VALUE_TYPE *a, const REAL_VALUE_TYPE *b,
902 int nan_result)
903 {
904 int ret;
905
906 switch (CLASS2 (a->cl, b->cl))
907 {
908 case CLASS2 (rvc_zero, rvc_zero):
909 /* Sign of zero doesn't matter for compares. */
910 return 0;
911
912 case CLASS2 (rvc_inf, rvc_zero):
913 case CLASS2 (rvc_inf, rvc_normal):
914 case CLASS2 (rvc_normal, rvc_zero):
915 return (a->sign ? -1 : 1);
916
917 case CLASS2 (rvc_inf, rvc_inf):
918 return -a->sign - -b->sign;
919
920 case CLASS2 (rvc_zero, rvc_normal):
921 case CLASS2 (rvc_zero, rvc_inf):
922 case CLASS2 (rvc_normal, rvc_inf):
923 return (b->sign ? 1 : -1);
924
925 case CLASS2 (rvc_zero, rvc_nan):
926 case CLASS2 (rvc_normal, rvc_nan):
927 case CLASS2 (rvc_inf, rvc_nan):
928 case CLASS2 (rvc_nan, rvc_nan):
929 case CLASS2 (rvc_nan, rvc_zero):
930 case CLASS2 (rvc_nan, rvc_normal):
931 case CLASS2 (rvc_nan, rvc_inf):
932 return nan_result;
933
934 case CLASS2 (rvc_normal, rvc_normal):
935 break;
936
937 default:
938 gcc_unreachable ();
939 }
940
941 if (a->sign != b->sign)
942 return -a->sign - -b->sign;
943
944 if (REAL_EXP (a) > REAL_EXP (b))
945 ret = 1;
946 else if (REAL_EXP (a) < REAL_EXP (b))
947 ret = -1;
948 else
949 ret = cmp_significands (a, b);
950
951 return (a->sign ? -ret : ret);
952 }
953
954 /* Return A truncated to an integral value toward zero. */
955
956 static void
957 do_fix_trunc (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a)
958 {
959 *r = *a;
960
961 switch (r->cl)
962 {
963 case rvc_zero:
964 case rvc_inf:
965 case rvc_nan:
966 break;
967
968 case rvc_normal:
969 if (REAL_EXP (r) <= 0)
970 get_zero (r, r->sign);
971 else if (REAL_EXP (r) < SIGNIFICAND_BITS)
972 clear_significand_below (r, SIGNIFICAND_BITS - REAL_EXP (r));
973 break;
974
975 default:
976 gcc_unreachable ();
977 }
978 }
979 #endif /* 0 */
980
981 /* Add 1 unit in the last place to a positive SMAP II BCD float. */
bcdpadd1ulp(smap_bcd_float * op)982 static void bcdpadd1ulp(smap_bcd_float *op)
983 {
984 int i;
985 op->mantissa++; /* This will give a digit of 10 in some cases. */
986 for (i=0;i<60;i+=4) { /* for each digit except the first one */
987 int d=(op->mantissa>>i)&15;
988 if (d<10) /* if the digit is <10, we're done */
989 break;
990 op->mantissa += 6ull<<i; /* subtract 10 from the digit, add 1 to the next digit */
991 }
992 /* Now, we should have a carry only in one case: the mantissa was
993 9999999999999999. So now, we have A000000000000000, where A is 10 in a
994 single digit. We can safely truncate the last digit because it is 0 (so
995 there is no risk of double-rounding here). */
996 if (op->mantissa >= 0xA000000000000000ull) { /* if we have a carry */
997 gcc_assert (op->mantissa == 0xA000000000000000ull); /* sanity check */
998 op->mantissa = 0x1000000000000000ull; /* A -> 10, drop the last digit */
999 op->exponent++; /* adjust the exponent */
1000 if (op->exponent > 0x4000+16383) /* exponent too large, overflow to +infinity */
1001 *op = POSITIVE_INF;
1002 }
1003 }
1004
1005 /* Add 2 positive SMAP II BCD floats. */
bcdppadd(smap_bcd_float op0,smap_bcd_float op1)1006 static smap_bcd_float bcdppadd(smap_bcd_float op0, smap_bcd_float op1)
1007 {
1008 int lastdigit=0;
1009
1010 if (op0.exponent < op1.exponent) {
1011 /* If op0 does not contribute to the result, avoid shift count overflow by
1012 returning op1 immediately. */
1013 if (op1.exponent-op0.exponent>16) return op1;
1014 /* Adjust op0 exponent to fit op1. */
1015 op0.mantissa >>= ((op1.exponent-op0.exponent-1)<<2);
1016 lastdigit = (op0.mantissa&15); /* save the last digit of op0 */
1017 op0.mantissa >>= 4; /* drop the last digit of op0 */
1018 op0.exponent = op1.exponent; /* adjust the exponent of op0 */
1019 } else if (op1.exponent < op0.exponent) {
1020 /* If op1 does not contribute to the result, avoid shift count overflow by
1021 returning op0 immediately. */
1022 if (op0.exponent-op1.exponent>16) return op0;
1023 /* Adjust op1 exponent to fit op0. */
1024 op1.mantissa >>= ((op0.exponent-op1.exponent-1)<<2);
1025 lastdigit = (op1.mantissa&15); /* save the last digit of op1 */
1026 op1.mantissa >>= 4; /* drop the last digit of op1 */
1027 op1.exponent = op0.exponent; /* adjust the exponent of op1 */
1028 }
1029
1030 /* Now do the addition in the BCD-coded mantissa. */
1031 {
1032 int i,d,carry=0;
1033 smap_bcd_float result={0,0};result.exponent=op0.exponent;
1034 for (i=0;i<64;i+=4) { /* for each digit */
1035 d = (op0.mantissa&15)+(op1.mantissa&15)+carry;
1036 carry = (d>=10); /* handle carry */
1037 if (carry) d -= 10;
1038 /* We are done with this digit of the mantissa. */
1039 op0.mantissa >>= 4;
1040 op1.mantissa >>= 4;
1041 /* Store it into the resulting mantissa. */
1042 result.mantissa += ((unsigned long long)d)<<i;
1043 }
1044 if (carry) { /* The mantissa overflowed, so we need to adjust the exponent. */
1045 lastdigit = (result.mantissa&15); /* save the last digit of result,
1046 dropping the previously saved last
1047 digit */
1048 result.mantissa >>= 4; /* drop the last digit of result */
1049 result.exponent++; /* adjust the exponent of result */
1050 if (result.exponent > 0x4000+16383) /* exponent too large, overflow to */
1051 return POSITIVE_INF; /* +infinity */
1052 /* prepend the first digit */
1053 result.mantissa += ((unsigned long long)carry)<<60;
1054 }
1055 if (lastdigit>=5) /* now do the correct rounding */
1056 bcdpadd1ulp(&result);
1057 return result;
1058 }
1059 }
1060
1061 /* Subtract 2 positive SMAP II BCD floats. */
bcdppsub(smap_bcd_float op0,smap_bcd_float op1)1062 static smap_bcd_float bcdppsub(smap_bcd_float op0, smap_bcd_float op1)
1063 {
1064 unsigned long long lastdigits=0;
1065
1066 if (REAL_VALUES_LESS(op0,op1)) /* if op0<op1, subtract the other way round */
1067 return REAL_VALUE_NEGATE(bcdppsub(op1,op0));
1068 else if (REAL_VALUES_IDENTICAL(op0,op1)) /* if op0=op1, return unsigned 0 */
1069 return UNSIGNED_ZERO;
1070
1071 /* Now, we can assume op0>op1. */
1072 if (op1.exponent < op0.exponent) {
1073 /* If op1 does not contribute to the result, avoid shift count overflow by
1074 returning op0 immediately. */
1075 if (op0.exponent-op1.exponent>16) return op0;
1076 /* Adjust op1 exponent to fit op0.
1077 Save all dropped digits, so we can round correctly. */
1078 lastdigits = op1.mantissa << (64-((op0.exponent-op1.exponent)<<2));
1079 if (op0.exponent-op1.exponent==16)
1080 op1.mantissa = 0; /* special-cased because shifts by 64 aren't allowed */
1081 else
1082 op1.mantissa >>= ((op0.exponent-op1.exponent)<<2);
1083 op1.exponent = op0.exponent; /* adjust the exponent of op1 */
1084 }
1085
1086 /* Now do the subtraction in the BCD-coded mantissa. */
1087 {
1088 int i,d,carry=0;
1089 smap_bcd_float result={0,0};result.exponent=op0.exponent;
1090 for (i=0;i<64;i+=4) { /* for each digit */
1091 d = (op0.mantissa&15)-((op1.mantissa&15)+carry);
1092 carry = (d<0); /* handle carry */
1093 if (carry) d += 10;
1094 /* We are done with this digit of the mantissa. */
1095 op0.mantissa >>= 4;
1096 op1.mantissa >>= 4;
1097 /* Store it into the resulting mantissa. */
1098 result.mantissa += ((unsigned long long)d)<<i;
1099 }
1100 gcc_assert (!carry); /* Carry should be 0 here! */
1101 /* Normalize and return the result (which is always positive here, since we
1102 assumed op1>op0). Handle rounding and extra digits during normalization.
1103 */
1104 while ((result.mantissa<0x1000000000000000ull)
1105 || ((result.mantissa==0x1000000000000000ull)
1106 && ((lastdigits>0x5000000000000000ull)
1107 || ((lastdigits>0x500000000000000ull)
1108 && (result.exponent>0x4000-16383))))) {
1109 /* while mantissa<1 */
1110 int c = (result.mantissa==0x1000000000000000ull); /* used by the sanity
1111 check below */
1112 result.exponent--; /* decrease exponent by 1 */
1113 if (result.exponent<0x4000-16383) /* exponent too small, underflow to +0 */
1114 return POSITIVE_ZERO;
1115 result.mantissa <<= 4; /* left-shift mantissa */
1116 carry = !!(lastdigits>>60);
1117 result.mantissa -= (lastdigits>>60)+(carry*6); /* Subtract extra digit. If
1118 there is a carry, add 10 to the digit, subtract 1 from the next digit.
1119 This can give digits of 15 in some cases. */
1120 for (i=4;i<64;i+=4) { /* for each digit except the last one */
1121 int d=(result.mantissa>>i)&15;
1122 if (d<10) /* if the digit is <10, we're done */
1123 break;
1124 result.mantissa -= 6ull<<i; /* add 10 to the digit, subtract 1 from the
1125 next digit */
1126 /* There should be no carry in the first digit except if the mantissa
1127 was exactly 1. */
1128 gcc_assert (c || result.mantissa<0xA000000000000000ull);
1129 }
1130 lastdigits <<= 4; /* We handled an extra digit. */
1131 }
1132 /* Now do the rounding. */
1133 if (lastdigits>0x5000000000000000ull) {
1134 int c = !result.mantissa; /* used by the sanity check below */
1135 /* This corner case should have been handled above. */
1136 gcc_assert (result.mantissa!=0x1000000000000000ull);
1137 result.mantissa--; /* Subtract 1. This can give digits of 15 in some cases. */
1138 for (i=0;i<64;i+=4) { /* for each digit */
1139 int d=(result.mantissa>>i)&15;
1140 if (d<10) /* if the digit is <10, we're done */
1141 break;
1142 result.mantissa -= 6ull<<i; /* add 10 to the digit, subtract 1 from the
1143 next digit */
1144 /* There should be no carry in the first digit except if the mantissa was
1145 0 (which actually means 10^16 here). */
1146 gcc_assert (c || result.mantissa<0xA000000000000000ull);
1147 }
1148 }
1149 return result;
1150 }
1151 }
1152
1153 /* Add 2 SMAP II BCD floats. */
bcdadd(smap_bcd_float op0,smap_bcd_float op1)1154 static smap_bcd_float bcdadd(smap_bcd_float op0, smap_bcd_float op1)
1155 {
1156 if (REAL_VALUE_ISNAN(op0) || REAL_VALUE_ISNAN(op1)) /* keep NAN */
1157 return NAN;
1158 else if (REAL_VALUE_ISINF(op0) && REAL_VALUE_ISINF(op1)) { /* both operands
1159 are infinity */
1160 if (!REAL_VALUES_IDENTICAL(op0,op1)
1161 || REAL_VALUES_IDENTICAL(op0,UNSIGNED_INF)) /* differing signs or both
1162 unsigned yield NAN */
1163 return NAN;
1164 else /* both positive or both negative */
1165 return op0;
1166 } else if (REAL_VALUE_ISINF(op0)) /* op0=inf, so op0+op1=inf+op1=inf=op0 */
1167 return op0;
1168 else if (REAL_VALUE_ISINF(op1)) /* op1=inf, so op0+op1=op0+inf=inf=op1 */
1169 return op1;
1170 else if (REAL_VALUE_ISZERO(op0) && REAL_VALUE_ISZERO(op1)) { /* both operands
1171 are 0 */
1172 if (REAL_VALUES_IDENTICAL(op0,op1)) /* both positive, both negative or both
1173 unsigned */
1174 return op0;
1175 else /* differing signs yield unsigned zero */
1176 return UNSIGNED_ZERO;
1177 } else if (REAL_VALUE_ISZERO(op0)) /* op0=0, so op0+op1=0+op1=op1 */
1178 return op1;
1179 else if (REAL_VALUE_ISZERO(op1)) /* op1=0, so op0+op1=op0+0=op0 */
1180 return op0;
1181 else if (REAL_VALUE_ISPOSITIVE(op0) && REAL_VALUE_ISPOSITIVE(op1))
1182 return bcdppadd(op0,op1);
1183 else if (REAL_VALUE_ISPOSITIVE(op0)) /* and op1 negative */
1184 return bcdppsub(op0,REAL_VALUE_NEGATE(op1));
1185 else if (REAL_VALUE_ISPOSITIVE(op1)) /* and op0 negative */
1186 return bcdppsub(op1,REAL_VALUE_NEGATE(op0));
1187 else /* both negative */
1188 return REAL_VALUE_NEGATE(bcdppadd(REAL_VALUE_NEGATE(op0),
1189 REAL_VALUE_NEGATE(op1)));
1190 }
1191
1192 /* Subtract 2 SMAP II BCD floats. */
bcdsub(smap_bcd_float op0,smap_bcd_float op1)1193 static smap_bcd_float bcdsub(smap_bcd_float op0, smap_bcd_float op1)
1194 {
1195 return bcdadd(op0,REAL_VALUE_NEGATE(op1));
1196 }
1197
1198 /* Multiply 2 positive SMAP II BCD floats. */
bcdppmul(smap_bcd_float op0,smap_bcd_float op1)1199 static smap_bcd_float bcdppmul(smap_bcd_float op0, smap_bcd_float op1)
1200 {
1201 /* Compute the result in 2 binary parts. The upper 16 decimal digits and the
1202 lower ones. */
1203 unsigned long long resulth=0, resultl=0;
1204 int i,j,k,d0,d1,d32,exponent;
1205 unsigned long long factor=1ull;
1206 for (i=0;i<64;i+=4,factor*=10ull) { /* for each result digit <16 */
1207 for (j=0;j<64;j+=4) { /* for each digit of op0 */
1208 k = i-j; /* corresponding op1 digit */
1209 if (k<0 || k>=64) continue; /* digit out of range */
1210 d0 = (op0.mantissa>>j)&15; /* jth digit of op0 */
1211 d1 = (op1.mantissa>>k)&15; /* kth digit of op0 */
1212 resultl += (unsigned long long)(d0*d1)*factor;
1213 while (resultl>=10000000000000000ull/*10^16*/) { /* carry into resulth */
1214 resultl -= 10000000000000000ull/*10^16*/;
1215 resulth++;
1216 }
1217 }
1218 }
1219 for (factor=1ull;i<128;i+=4,factor*=10ull) { /* for each result digit >=16 */
1220 for (j=0;j<64;j+=4) { /* for each digit of op0 */
1221 k = i-j; /* corresponding op1 digit */
1222 if (k<0 || k>=64) continue; /* digit out of range */
1223 d0 = (op0.mantissa>>j)&15; /* jth digit of op0 */
1224 d1 = (op1.mantissa>>k)&15; /* kth digit of op0 */
1225 resulth += (unsigned long long)(d0*d1)*factor;
1226 }
1227 }
1228
1229 /* resultl should always be <10^16 */
1230 gcc_assert (resultl<10000000000000000ull);
1231
1232 /* Because of normalization, the result has either 31 or 32 digits. */
1233 d32 = (resulth>=1000000000000000ull/*10^15*/
1234 || (resulth==999999999999999ull/*10^15-1*/
1235 && resultl>=9500000000000000ull/*9.5*10^15*/));
1236 if (!d32) { /* if we have only 15 digits in resulth, take one from resultl */
1237 resulth = resulth*10ull+resultl/1000000000000000ull/*10^15*/;
1238 resultl = (resultl%1000000000000000ull/*10^15*/)*10ull;
1239 }
1240 if (resultl>=5000000000000000ull/*5*10^15*/) /* round resultl into resulth */
1241 resulth++;
1242
1243 /* Now compute the exponent. */
1244 exponent=op0.exponent+op1.exponent-0x4000+d32;
1245 if (exponent>0x4000+16383) /* exponent too large, overflow to +infinity */
1246 return POSITIVE_INF;
1247 if (exponent<0x4000-16383) /* exponent too small, underflow to +0 */
1248 return POSITIVE_ZERO;
1249
1250 /* Now convert resulth into a BCD mantissa. */
1251 {
1252 unsigned long long d;
1253 smap_bcd_float result={0,0};result.exponent=exponent;
1254 for (i=0;i<64;i+=4) { /* for each digit */
1255 d = resulth%10ull; /* Extract the digit. */
1256 resulth /= 10ull; /* We are done with this digit of the mantissa. */
1257 result.mantissa += d<<i; /* Store it into the resulting mantissa. */
1258 }
1259 return result;
1260 }
1261 }
1262
1263 /* Multiply 2 SMAP II BCD floats. */
bcdmul(smap_bcd_float op0,smap_bcd_float op1)1264 static smap_bcd_float bcdmul(smap_bcd_float op0, smap_bcd_float op1)
1265 {
1266 if (REAL_VALUE_ISNAN(op0) || REAL_VALUE_ISNAN(op1)
1267 || (REAL_VALUE_ISINF(op0) && REAL_VALUE_ISZERO(op1))
1268 || (REAL_VALUE_ISZERO(op0) && REAL_VALUE_ISINF(op1))) /* keep NAN,
1269 0*inf=NAN */
1270 return NAN;
1271 else if (REAL_VALUES_IDENTICAL(op0,UNSIGNED_INF)
1272 || REAL_VALUES_IDENTICAL(op1,UNSIGNED_INF)) /* unsigned inf * non-0 =
1273 unsigned inf */
1274 return UNSIGNED_INF;
1275 else if ((REAL_VALUES_IDENTICAL(op0,POSITIVE_INF)
1276 && REAL_VALUES_IDENTICAL(op1,POSITIVE_INF)) /* +inf * +inf = +inf */
1277 || (REAL_VALUES_IDENTICAL(op0,NEGATIVE_INF)
1278 && REAL_VALUES_IDENTICAL(op1,NEGATIVE_INF))) /* -inf * -inf = +inf */
1279 return POSITIVE_INF;
1280 else if ((REAL_VALUES_IDENTICAL(op0,POSITIVE_INF)
1281 && REAL_VALUES_IDENTICAL(op1,NEGATIVE_INF)) /* +inf * -inf = -inf */
1282 || (REAL_VALUES_IDENTICAL(op0,NEGATIVE_INF)
1283 && REAL_VALUES_IDENTICAL(op1,POSITIVE_INF))) /* -inf * +inf = -inf */
1284 return NEGATIVE_INF;
1285 /* Now we can assume that at least 1 value is finite, and that we don't have
1286 an unsigned infinity, a NAN or an inf*0 indeterminate form. */
1287 else if ((REAL_VALUES_IDENTICAL(op0,POSITIVE_INF)
1288 && REAL_VALUE_ISPOSITIVE(op1)) /* +inf * +finite = +inf */
1289 || (REAL_VALUE_ISPOSITIVE(op0)
1290 && REAL_VALUES_IDENTICAL(op1,POSITIVE_INF)) /* +finite * +inf = +inf */
1291 || (REAL_VALUES_IDENTICAL(op0,NEGATIVE_INF)
1292 && REAL_VALUE_ISNEGATIVE(op1)) /* -inf * -finite = +inf */
1293 || (REAL_VALUE_ISNEGATIVE(op0)
1294 && REAL_VALUES_IDENTICAL(op1,NEGATIVE_INF))) /* -finite * -inf = +inf */
1295 return POSITIVE_INF;
1296 else if ((REAL_VALUES_IDENTICAL(op0,POSITIVE_INF)
1297 && REAL_VALUE_ISNEGATIVE(op1)) /* +inf * -finite = -inf */
1298 || (REAL_VALUE_ISPOSITIVE(op0)
1299 && REAL_VALUES_IDENTICAL(op1,NEGATIVE_INF)) /* +finite * -inf = -inf */
1300 || (REAL_VALUES_IDENTICAL(op0,NEGATIVE_INF)
1301 && REAL_VALUE_ISPOSITIVE(op1)) /* -inf * +finite = -inf */
1302 || (REAL_VALUE_ISNEGATIVE(op0)
1303 && REAL_VALUES_IDENTICAL(op1,POSITIVE_INF))) /* -finite * +inf = -inf */
1304 return NEGATIVE_INF;
1305 /* Now we can assume that both values are finite. */
1306 else if (REAL_VALUES_IDENTICAL(op0,UNSIGNED_ZERO)
1307 || REAL_VALUES_IDENTICAL(op1,UNSIGNED_ZERO)) /* unsigned 0 * finite =
1308 unsigned 0 */
1309 return UNSIGNED_ZERO;
1310 else if ((REAL_VALUES_IDENTICAL(op0,POSITIVE_ZERO)
1311 && REAL_VALUES_IDENTICAL(op1,POSITIVE_ZERO)) /* +0 * +0 = +0 */
1312 || (REAL_VALUES_IDENTICAL(op0,NEGATIVE_ZERO)
1313 && REAL_VALUES_IDENTICAL(op1,NEGATIVE_ZERO))) /* -0 * -0 = +0 */
1314 return POSITIVE_ZERO;
1315 else if ((REAL_VALUES_IDENTICAL(op0,POSITIVE_ZERO)
1316 && REAL_VALUES_IDENTICAL(op1,NEGATIVE_ZERO)) /* +0 * -0 = -0 */
1317 || (REAL_VALUES_IDENTICAL(op0,NEGATIVE_ZERO)
1318 && REAL_VALUES_IDENTICAL(op1,POSITIVE_ZERO))) /* -0 * +0 = -0 */
1319 return NEGATIVE_ZERO;
1320 /* Now we can assume that both values are finite, at least 1 value is non-0,
1321 and that we don't have an unsigned 0. */
1322 else if ((REAL_VALUES_IDENTICAL(op0,POSITIVE_ZERO)
1323 && REAL_VALUE_ISPOSITIVE(op1)) /* +0 * +finite = +0 */
1324 || (REAL_VALUE_ISPOSITIVE(op0)
1325 && REAL_VALUES_IDENTICAL(op1,POSITIVE_ZERO)) /* +finite * +0 = +0 */
1326 || (REAL_VALUES_IDENTICAL(op0,NEGATIVE_ZERO)
1327 && REAL_VALUE_ISNEGATIVE(op1)) /* -0 * -finite = +0 */
1328 || (REAL_VALUE_ISNEGATIVE(op0)
1329 && REAL_VALUES_IDENTICAL(op1,NEGATIVE_ZERO))) /* -finite * -0 = +0 */
1330 return POSITIVE_ZERO;
1331 else if ((REAL_VALUES_IDENTICAL(op0,POSITIVE_ZERO)
1332 && REAL_VALUE_ISNEGATIVE(op1)) /* +0 * -finite = -0 */
1333 || (REAL_VALUE_ISPOSITIVE(op0)
1334 && REAL_VALUES_IDENTICAL(op1,NEGATIVE_ZERO)) /* +finite * -0 = -0 */
1335 || (REAL_VALUES_IDENTICAL(op0,NEGATIVE_ZERO)
1336 && REAL_VALUE_ISPOSITIVE(op1)) /* -0 * +finite = -0 */
1337 || (REAL_VALUE_ISNEGATIVE(op0)
1338 && REAL_VALUES_IDENTICAL(op1,POSITIVE_ZERO))) /* -finite * +0 = -0 */
1339 return NEGATIVE_ZERO;
1340 /* Now we can assume that both values are finite and non-0. */
1341 else if (REAL_VALUE_ISPOSITIVE(op0) && REAL_VALUE_ISPOSITIVE(op1))
1342 return bcdppmul(op0,op1);
1343 else if (REAL_VALUE_ISPOSITIVE(op0)) /* and op1 negative */
1344 return REAL_VALUE_NEGATE(bcdppmul(op0,REAL_VALUE_NEGATE(op1)));
1345 else if (REAL_VALUE_ISPOSITIVE(op1)) /* and op0 negative */
1346 return REAL_VALUE_NEGATE(bcdppmul(REAL_VALUE_NEGATE(op0),op1));
1347 else /* both negative */
1348 return bcdppmul(REAL_VALUE_NEGATE(op0),REAL_VALUE_NEGATE(op1));
1349 }
1350
1351 /* Divide 2 positive SMAP II BCD floats. Indicate if the computed result was exact. */
bcdppdiv(smap_bcd_float op0,smap_bcd_float op1,int * exactresult)1352 static smap_bcd_float bcdppdiv(smap_bcd_float op0, smap_bcd_float op1, int *exactresult)
1353 {
1354 /* The dividend is represented multiplied by 10^16, and in 2 binary parts. The
1355 upper 64 bits and the lower ones. The divisor is represented in binary, in
1356 2 parts to allow shifting. Compute the result mantissa in binary. */
1357 unsigned long long dividendh, dividendl=0, divisorh, divisorl=0, resultm=0;
1358 int i,d17,exponent;
1359 unsigned long long factor=1ull;
1360 /* Convert the mantissa of op0 to binary. */
1361 for (i=0;i<64;i+=4,factor*=10ull) /* for each digit of op0 */
1362 dividendl += (unsigned long long)((op0.mantissa>>i)&15)*factor;
1363 /* Multiply the dividend with 10^16. */
1364 dividendh = dividendl>>48; /* multiply first 16 bits with 2^16 */
1365 dividendh *= 152587890625ull; /* and with 5^16 */
1366 dividendl &= 0xFFFFFFFFFFFFull; /* remove them from dividendl */
1367 dividendl *= 15625ull; /* multiply the remaining 48 bits with 5^6
1368 we still have to multiply them by 2^16*5^10 */
1369 dividendh += (dividendl>>48)*9765625ull; /* multiply first 16 bits with 2^16*5^10 */
1370 dividendl &= 0xFFFFFFFFFFFFull; /* remove them from dividendl */
1371 dividendl *= 15625ull; /* multiply the remaining 48 bits with 5^6
1372 we still have to multiply them by 2^16*5^4 */
1373 dividendh += (dividendl>>48)*625ull; /* multiply first 16 bits with 2^16*5^4 */
1374 dividendl &= 0xFFFFFFFFFFFFull; /* remove them from dividendl */
1375 dividendl *= 625ull; /* multiply the remaining 48 bits with 5^4
1376 we still have to multiply them by 2^16 */
1377 dividendh += dividendl>>48; /* multiply first 16 bits with 2^16 */
1378 dividendl <<= 16; /* multiply the remaining 48 bits with 2^16 */
1379 /* Convert the mantissa of op1 to binary. */
1380 for (i=0,factor=1ull;i<64;i+=4,factor*=10ull) /* for each digit of op1 */
1381 divisorl += (unsigned long long)((op1.mantissa>>i)&15)*factor;
1382 /* Multiply the divisor with 2^56. We know that, due to normalization, the
1383 result is always <10^17, which is smaller than 2^57, so we don't have to go
1384 through the full 128-bit division. */
1385 divisorh = divisorl>>(64-56);
1386 divisorl <<= 56;
1387 /* Now do the 128-bit division. */
1388 for (i=56;i>=0;i--) {
1389 /* Shift the result to the left. */
1390 resultm <<= 1;
1391 /* Check if the divisor fits into the dividend. */
1392 if (divisorh<dividendh || (divisorh==dividendh && divisorl<=dividendl)) {
1393 /* Add 1 to the result. */
1394 resultm++;
1395 /* Subtract the divisor from the dividend. */
1396 if (dividendl<divisorl) /* handle carry, use unsigned wraparound overflow */
1397 dividendh--;
1398 dividendl -= divisorl; /* now do the subtraction */
1399 dividendh -= divisorh;
1400 }
1401 if (i) {
1402 /* Shift the divisor to the right. */
1403 divisorl = ((divisorh&1)<<63)+(divisorl>>1);
1404 divisorh >>= 1;
1405 }
1406 }
1407 /* dividendl now contains the remainder. It is always smaller than the
1408 divisor, so it always fits into 64 bits. divisorl now contains the original
1409 divisor. */
1410
1411 if (exactresult)
1412 *exactresult = !dividendl; /* if there is a remainder, the result sure is
1413 not exact, otherwise, let's assume it is for
1414 a moment */
1415
1416 /* Because of normalization, the result has either 16 or 17 digits. */
1417 d17 = (resultm>=10000000000000000ull/*10^16*/
1418 || (resultm==9999999999999999ull/*10^16-1*/
1419 && (dividendl<<1)>=divisorl) /* 2r>=d <=> r>=d/2 */);
1420 if (d17) { /* if we have 17 digits in the result, drop one */
1421 if (exactresult && (resultm%10ull)) /* if we are about to drop a non-0
1422 digit, the result is not exact
1423 anymore */
1424 *exactresult = 0;
1425 resultm = (resultm+5ull)/10ull; /* add 5 for correct rounding */
1426 } else {
1427 if ((dividendl<<1)>=divisorl /* r>=d/2 */) /* round remainder into result */
1428 resultm++;
1429 }
1430
1431 /* Now compute the exponent. */
1432 exponent=op0.exponent-op1.exponent+0x4000-(!d17);
1433 if (exactresult && (exponent>0x4000+16383 || exponent<0x4000-16383))
1434 /* if we overflowed, the result is not exact anymore */
1435 *exactresult = 0;
1436 if (exponent>0x4000+16383) /* exponent too large, overflow to +infinity */
1437 return POSITIVE_INF;
1438 if (exponent<0x4000-16383) /* exponent too small, underflow to +0 */
1439 return POSITIVE_ZERO;
1440
1441 /* Now convert resulth into a BCD mantissa. */
1442 {
1443 unsigned long long d;
1444 smap_bcd_float result={0,0};result.exponent=exponent;
1445 for (i=0;i<64;i+=4) { /* for each digit */
1446 d = resultm%10ull; /* Extract the digit. */
1447 resultm /= 10ull; /* We are done with this digit of the mantissa. */
1448 result.mantissa += d<<i; /* Store it into the resulting mantissa. */
1449 }
1450 return result;
1451 }
1452 }
1453
1454 /* Divide 2 SMAP II BCD floats. Indicate if the computed result was exact. */
bcddiv(smap_bcd_float op0,smap_bcd_float op1,int * exactresult)1455 static smap_bcd_float bcddiv(smap_bcd_float op0, smap_bcd_float op1, int *exactresult)
1456 {
1457 if (exactresult) *exactresult=1; /* special cases are all exact */
1458 if (REAL_VALUE_ISNAN(op0) || REAL_VALUE_ISNAN(op1)
1459 || (REAL_VALUE_ISINF(op0) && REAL_VALUE_ISINF(op1))
1460 || (REAL_VALUE_ISZERO(op0) && REAL_VALUE_ISZERO(op1))) /* keep NAN,
1461 0/0=inf/inf=NAN */
1462 return NAN;
1463 else if (REAL_VALUES_IDENTICAL(op0,UNSIGNED_INF)
1464 || REAL_VALUES_IDENTICAL(op1,UNSIGNED_ZERO)) /* unsigned inf / finite
1465 = non-0 / unsigned 0 = unsigned inf */
1466 return UNSIGNED_INF;
1467 else if ((REAL_VALUES_IDENTICAL(op0,POSITIVE_INF)
1468 && REAL_VALUES_IDENTICAL(op1,POSITIVE_ZERO)) /* +inf / +0 = +inf */
1469 || (REAL_VALUES_IDENTICAL(op0,NEGATIVE_INF)
1470 && REAL_VALUES_IDENTICAL(op1,NEGATIVE_ZERO))) /* -inf / -0 = +inf */
1471 return POSITIVE_INF;
1472 else if ((REAL_VALUES_IDENTICAL(op0,POSITIVE_INF)
1473 && REAL_VALUES_IDENTICAL(op1,NEGATIVE_ZERO)) /* +inf / -0 = -inf */
1474 || (REAL_VALUES_IDENTICAL(op0,NEGATIVE_INF)
1475 && REAL_VALUES_IDENTICAL(op1,POSITIVE_ZERO))) /* -inf / +0 = -inf */
1476 return NEGATIVE_INF;
1477 /* Now we can assume that at least 1 of op0 and 1/op1 is finite, and that we
1478 don't have op0 = unsigned inf, op1 = unsigned 0, a NAN or an inf/inf or 0/0
1479 indeterminate form. */
1480 else if ((REAL_VALUES_IDENTICAL(op0,POSITIVE_INF)
1481 && REAL_VALUE_ISPOSITIVE(op1)) /* +inf / +non-0 = +inf */
1482 || (REAL_VALUE_ISPOSITIVE(op0)
1483 && REAL_VALUES_IDENTICAL(op1,POSITIVE_ZERO)) /* +finite / +0 = +inf */
1484 || (REAL_VALUES_IDENTICAL(op0,NEGATIVE_INF)
1485 && REAL_VALUE_ISNEGATIVE(op1)) /* -inf / -non-0 = +inf */
1486 || (REAL_VALUE_ISNEGATIVE(op0)
1487 && REAL_VALUES_IDENTICAL(op1,NEGATIVE_ZERO))) /* -finite / -0 = +inf */
1488 return POSITIVE_INF;
1489 else if ((REAL_VALUES_IDENTICAL(op0,POSITIVE_INF)
1490 && REAL_VALUE_ISNEGATIVE(op1)) /* +inf / -non-0 = -inf */
1491 || (REAL_VALUE_ISPOSITIVE(op0)
1492 && REAL_VALUES_IDENTICAL(op1,NEGATIVE_ZERO)) /* +finite / -0 = -inf */
1493 || (REAL_VALUES_IDENTICAL(op0,NEGATIVE_INF)
1494 && REAL_VALUE_ISPOSITIVE(op1)) /* -inf / +non-0 = -inf */
1495 || (REAL_VALUE_ISNEGATIVE(op0)
1496 && REAL_VALUES_IDENTICAL(op1,POSITIVE_ZERO))) /* -finite / +0 = -inf */
1497 return NEGATIVE_INF;
1498 /* Now we can assume that both of op0 and 1/op1 are finite. */
1499 else if (REAL_VALUES_IDENTICAL(op0,UNSIGNED_ZERO)
1500 || REAL_VALUES_IDENTICAL(op1,UNSIGNED_INF)) /* unsigned 0 / non-0 =
1501 finite / unsigned inf = unsigned 0 */
1502 return UNSIGNED_ZERO;
1503 else if ((REAL_VALUES_IDENTICAL(op0,POSITIVE_ZERO)
1504 && REAL_VALUES_IDENTICAL(op1,POSITIVE_INF)) /* +0 / +inf = +0 */
1505 || (REAL_VALUES_IDENTICAL(op0,NEGATIVE_ZERO)
1506 && REAL_VALUES_IDENTICAL(op1,NEGATIVE_INF))) /* -0 / -inf = +0 */
1507 return POSITIVE_ZERO;
1508 else if ((REAL_VALUES_IDENTICAL(op0,POSITIVE_ZERO)
1509 && REAL_VALUES_IDENTICAL(op1,NEGATIVE_INF)) /* +0 / -inf = -0 */
1510 || (REAL_VALUES_IDENTICAL(op0,NEGATIVE_ZERO)
1511 && REAL_VALUES_IDENTICAL(op1,POSITIVE_INF))) /* -0 / +inf = -0 */
1512 return NEGATIVE_ZERO;
1513 /* Now we can assume that both of op0 and 1/op1 are finite, at least 1 of them
1514 is non-0, and that neither of them is an unsigned 0. */
1515 else if ((REAL_VALUES_IDENTICAL(op0,POSITIVE_ZERO)
1516 && REAL_VALUE_ISPOSITIVE(op1)) /* +0 / +non-0 = +0 */
1517 || (REAL_VALUE_ISPOSITIVE(op0)
1518 && REAL_VALUES_IDENTICAL(op1,POSITIVE_INF)) /* +finite / +inf = +0 */
1519 || (REAL_VALUES_IDENTICAL(op0,NEGATIVE_ZERO)
1520 && REAL_VALUE_ISNEGATIVE(op1)) /* -0 / -non-0 = +0 */
1521 || (REAL_VALUE_ISNEGATIVE(op0)
1522 && REAL_VALUES_IDENTICAL(op1,NEGATIVE_INF))) /* -finite / -inf = +0 */
1523 return POSITIVE_ZERO;
1524 else if ((REAL_VALUES_IDENTICAL(op0,POSITIVE_ZERO)
1525 && REAL_VALUE_ISNEGATIVE(op1)) /* +0 / -non-0 = -0 */
1526 || (REAL_VALUE_ISPOSITIVE(op0)
1527 && REAL_VALUES_IDENTICAL(op1,NEGATIVE_INF)) /* +finite / -inf = -0 */
1528 || (REAL_VALUES_IDENTICAL(op0,NEGATIVE_ZERO)
1529 && REAL_VALUE_ISPOSITIVE(op1)) /* -0 / +non-0 = -0 */
1530 || (REAL_VALUE_ISNEGATIVE(op0)
1531 && REAL_VALUES_IDENTICAL(op1,POSITIVE_INF))) /* -finite / +inf = -0 */
1532 return NEGATIVE_ZERO;
1533 /* Now we can assume that both values are finite and non-0. */
1534 else if (REAL_VALUE_ISPOSITIVE(op0) && REAL_VALUE_ISPOSITIVE(op1))
1535 return bcdppdiv(op0,op1,exactresult);
1536 else if (REAL_VALUE_ISPOSITIVE(op0)) /* and op1 negative */
1537 return REAL_VALUE_NEGATE(bcdppdiv(op0,REAL_VALUE_NEGATE(op1),exactresult));
1538 else if (REAL_VALUE_ISPOSITIVE(op1)) /* and op0 negative */
1539 return REAL_VALUE_NEGATE(bcdppdiv(REAL_VALUE_NEGATE(op0),op1,exactresult));
1540 else /* both negative */
1541 return bcdppdiv(REAL_VALUE_NEGATE(op0),REAL_VALUE_NEGATE(op1),exactresult);
1542 }
1543
1544 /* Compute min of 2 SMAP II BCD floats. */
bcdmin(smap_bcd_float op0,smap_bcd_float op1)1545 static smap_bcd_float bcdmin(smap_bcd_float op0, smap_bcd_float op1)
1546 {
1547 if (REAL_VALUE_ISNANUINF(op0) || REAL_VALUE_ISNANUINF(op1)) /* keep NAN,
1548 UNSIGNED_INF is neither smaller nor larger than the other */
1549 return NAN;
1550 else if (REAL_VALUES_IDENTICAL(op0,NEGATIVE_INF)
1551 || REAL_VALUES_IDENTICAL(op1,NEGATIVE_INF)) /* keep -infinity */
1552 return NEGATIVE_INF;
1553 else if (REAL_VALUES_IDENTICAL(op0,POSITIVE_INF)) /* all numbers are smaller
1554 than +infinity */
1555 return op1;
1556 else if (REAL_VALUES_IDENTICAL(op1,POSITIVE_INF)) /* all numbers are smaller
1557 than +infinity */
1558 return op0;
1559 /* Now, we can assume that all values are finite. */
1560 else if (REAL_VALUES_LESS(op0,op1)) /* if op0<op1, return op0 */
1561 return op0;
1562 else /* if op0>=op1, return op1 */
1563 return op1;
1564 }
1565
1566 /* Compute max of 2 SMAP II BCD floats. */
bcdmax(smap_bcd_float op0,smap_bcd_float op1)1567 static smap_bcd_float bcdmax(smap_bcd_float op0, smap_bcd_float op1)
1568 {
1569 if (REAL_VALUE_ISNANUINF(op0) || REAL_VALUE_ISNANUINF(op1)) /* keep NAN,
1570 UNSIGNED_INF is neither smaller nor larger than the other */
1571 return NAN;
1572 else if (REAL_VALUES_IDENTICAL(op0,POSITIVE_INF)
1573 || REAL_VALUES_IDENTICAL(op1,POSITIVE_INF)) /* keep +infinity */
1574 return POSITIVE_INF;
1575 else if (REAL_VALUES_IDENTICAL(op0,NEGATIVE_INF)) /* all numbers are smaller
1576 than -infinity */
1577 return op1;
1578 else if (REAL_VALUES_IDENTICAL(op1,NEGATIVE_INF)) /* all numbers are smaller
1579 than -infinity */
1580 return op0;
1581 /* Now, we can assume that all values are finite. */
1582 else if (REAL_VALUES_LESS(op0,op1)) /* if op0<op1, return op1 */
1583 return op1;
1584 else /* if op0>=op1, return op0 */
1585 return op0;
1586 }
1587
1588 /* Perform the binary or unary operation described by CODE.
1589 For a unary operation, leave OP1 NULL. This function returns
1590 true if the result may be inexact due to loss of precision. */
1591
1592 bool
real_arithmetic(REAL_VALUE_TYPE * r,int icode,const REAL_VALUE_TYPE * op0,const REAL_VALUE_TYPE * op1)1593 real_arithmetic (REAL_VALUE_TYPE *r, int icode, const REAL_VALUE_TYPE *op0,
1594 const REAL_VALUE_TYPE *op1)
1595 {
1596 switch (icode)
1597 {
1598 case PLUS_EXPR:
1599 *r = bcdadd (*op0, *op1);
1600 return true;
1601
1602 case MINUS_EXPR:
1603 *r = bcdsub (*op0, *op1);
1604 return true;
1605
1606 case MULT_EXPR:
1607 *r = bcdmul (*op0, *op1);
1608 return true;
1609
1610 case RDIV_EXPR:
1611 *r = bcddiv (*op0, *op1, NULL);
1612 return true;
1613
1614 case MIN_EXPR:
1615 *r = bcdmin (*op0, *op1);
1616 return false;
1617
1618 case MAX_EXPR:
1619 *r = bcdmax (*op0, *op1);
1620 return false;
1621
1622 default:
1623 gcc_unreachable ();
1624 }
1625
1626 #if 0
1627 enum tree_code code = icode;
1628
1629 switch (code)
1630 {
1631 case PLUS_EXPR:
1632 return do_add (r, op0, op1, 0);
1633
1634 case MINUS_EXPR:
1635 return do_add (r, op0, op1, 1);
1636
1637 case MULT_EXPR:
1638 return do_multiply (r, op0, op1);
1639
1640 case RDIV_EXPR:
1641 return do_divide (r, op0, op1);
1642
1643 case MIN_EXPR:
1644 if (op1->cl == rvc_nan)
1645 *r = *op1;
1646 else if (do_compare (op0, op1, -1) < 0)
1647 *r = *op0;
1648 else
1649 *r = *op1;
1650 break;
1651
1652 case MAX_EXPR:
1653 if (op1->cl == rvc_nan)
1654 *r = *op1;
1655 else if (do_compare (op0, op1, 1) < 0)
1656 *r = *op1;
1657 else
1658 *r = *op0;
1659 break;
1660
1661 case NEGATE_EXPR:
1662 *r = *op0;
1663 r->sign ^= 1;
1664 break;
1665
1666 case ABS_EXPR:
1667 *r = *op0;
1668 r->sign = 0;
1669 break;
1670
1671 case FIX_TRUNC_EXPR:
1672 do_fix_trunc (r, op0);
1673 break;
1674
1675 default:
1676 gcc_unreachable ();
1677 }
1678 #endif /* 0 */
1679 return false;
1680 }
1681
1682 #if 0
1683 /* Legacy. Similar, but return the result directly. */
1684
1685 REAL_VALUE_TYPE
1686 real_arithmetic2 (int icode, const REAL_VALUE_TYPE *op0,
1687 const REAL_VALUE_TYPE *op1)
1688 {
1689 REAL_VALUE_TYPE r;
1690 real_arithmetic (&r, icode, op0, op1);
1691 return r;
1692 }
1693 #endif /* 0 */
1694
1695 bool
real_compare(int icode,const REAL_VALUE_TYPE * op0,const REAL_VALUE_TYPE * op1)1696 real_compare (int icode, const REAL_VALUE_TYPE *op0,
1697 const REAL_VALUE_TYPE *op1)
1698 {
1699 enum tree_code code = icode;
1700
1701 switch (code)
1702 {
1703 case LT_EXPR:
1704 return REAL_VALUES_LESS (*op0, *op1);
1705 case LE_EXPR:
1706 return !real_compare (UNGT_EXPR, op0, op1);
1707 case GT_EXPR:
1708 return real_compare (LT_EXPR, op1, op0);
1709 case GE_EXPR:
1710 return !real_compare (UNLT_EXPR, op0, op1);
1711 case EQ_EXPR:
1712 return REAL_VALUES_EQUAL (*op0, *op1);
1713 case NE_EXPR:
1714 return !real_compare (EQ_EXPR, op0, op1);
1715 case UNORDERED_EXPR:
1716 return REAL_VALUE_ISNANUINF (*op0) || REAL_VALUE_ISNANUINF (*op1);
1717 case ORDERED_EXPR:
1718 return !real_compare (UNORDERED_EXPR, op0, op1);
1719 case UNLT_EXPR:
1720 return real_compare (UNORDERED_EXPR, op0, op1) || real_compare (LT_EXPR, op0, op1);
1721 case UNLE_EXPR:
1722 return !real_compare (GT_EXPR, op0, op1);
1723 case UNGT_EXPR:
1724 return real_compare (UNORDERED_EXPR, op0, op1) || real_compare (GT_EXPR, op0, op1);
1725 case UNGE_EXPR:
1726 return !real_compare (LT_EXPR, op0, op1);
1727 case UNEQ_EXPR:
1728 return real_compare (UNORDERED_EXPR, op0, op1) || real_compare (EQ_EXPR, op0, op1);
1729 case LTGT_EXPR:
1730 return !real_compare (UNEQ_EXPR, op0, op1);
1731
1732 default:
1733 gcc_unreachable ();
1734 }
1735 }
1736
1737 #if 0
1738 /* Return floor log2(R). */
1739
1740 int
1741 real_exponent (const REAL_VALUE_TYPE *r)
1742 {
1743 switch (r->cl)
1744 {
1745 case rvc_zero:
1746 return 0;
1747 case rvc_inf:
1748 case rvc_nan:
1749 return (unsigned int)-1 >> 1;
1750 case rvc_normal:
1751 return REAL_EXP (r);
1752 default:
1753 gcc_unreachable ();
1754 }
1755 }
1756
1757 /* R = OP0 * 2**EXP. */
1758
1759 void
1760 real_ldexp (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *op0, int exp)
1761 {
1762 *r = *op0;
1763 switch (r->cl)
1764 {
1765 case rvc_zero:
1766 case rvc_inf:
1767 case rvc_nan:
1768 break;
1769
1770 case rvc_normal:
1771 exp += REAL_EXP (op0);
1772 if (exp > MAX_EXP)
1773 get_inf (r, r->sign);
1774 else if (exp < -MAX_EXP)
1775 get_zero (r, r->sign);
1776 else
1777 SET_REAL_EXP (r, exp);
1778 break;
1779
1780 default:
1781 gcc_unreachable ();
1782 }
1783 }
1784
1785 /* Determine whether a floating-point value X is infinite. */
1786
1787 bool
1788 real_isinf (const REAL_VALUE_TYPE *r)
1789 {
1790 return (r->cl == rvc_inf);
1791 }
1792
1793 /* Determine whether a floating-point value X is a NaN. */
1794
1795 bool
1796 real_isnan (const REAL_VALUE_TYPE *r)
1797 {
1798 return (r->cl == rvc_nan);
1799 }
1800
1801 /* Determine whether a floating-point value X is negative. */
1802
1803 bool
1804 real_isneg (const REAL_VALUE_TYPE *r)
1805 {
1806 return r->sign;
1807 }
1808
1809 /* Determine whether a floating-point value X is minus zero. */
1810
1811 bool
1812 real_isnegzero (const REAL_VALUE_TYPE *r)
1813 {
1814 return r->sign && r->cl == rvc_zero;
1815 }
1816
1817 /* Compare two floating-point objects for bitwise identity. */
1818
1819 bool
1820 real_identical (const REAL_VALUE_TYPE *a, const REAL_VALUE_TYPE *b)
1821 {
1822 int i;
1823
1824 if (a->cl != b->cl)
1825 return false;
1826 if (a->sign != b->sign)
1827 return false;
1828
1829 switch (a->cl)
1830 {
1831 case rvc_zero:
1832 case rvc_inf:
1833 return true;
1834
1835 case rvc_normal:
1836 if (REAL_EXP (a) != REAL_EXP (b))
1837 return false;
1838 break;
1839
1840 case rvc_nan:
1841 if (a->signalling != b->signalling)
1842 return false;
1843 /* The significand is ignored for canonical NaNs. */
1844 if (a->canonical || b->canonical)
1845 return a->canonical == b->canonical;
1846 break;
1847
1848 default:
1849 gcc_unreachable ();
1850 }
1851
1852 for (i = 0; i < SIGSZ; ++i)
1853 if (a->sig[i] != b->sig[i])
1854 return false;
1855
1856 return true;
1857 }
1858 #endif /* 0 */
1859
1860 /* Try to change R into its exact multiplicative inverse in machine
1861 mode MODE. Return true if successful. */
1862
1863 bool
exact_real_inverse(enum machine_mode mode,REAL_VALUE_TYPE * r)1864 exact_real_inverse (enum machine_mode mode, REAL_VALUE_TYPE *r)
1865 {
1866 int exactresult=0;
1867 smap_bcd_float one={0x4000,0x1000000000000000ull}, invr;
1868
1869 invr = bcddiv(one,*r,&exactresult);
1870 if (exactresult) *r=invr;
1871 return exactresult;
1872 #if 0
1873 const REAL_VALUE_TYPE *one = real_digit (1);
1874 REAL_VALUE_TYPE u;
1875 int i;
1876
1877 if (r->cl != rvc_normal)
1878 return false;
1879
1880 /* Check for a power of two: all significand bits zero except the MSB. */
1881 for (i = 0; i < SIGSZ-1; ++i)
1882 if (r->sig[i] != 0)
1883 return false;
1884 if (r->sig[SIGSZ-1] != SIG_MSB)
1885 return false;
1886
1887 /* Find the inverse and truncate to the required mode. */
1888 do_divide (&u, one, r);
1889 real_convert (&u, mode, &u);
1890
1891 /* The rounding may have overflowed. */
1892 if (u.cl != rvc_normal)
1893 return false;
1894 for (i = 0; i < SIGSZ-1; ++i)
1895 if (u.sig[i] != 0)
1896 return false;
1897 if (u.sig[SIGSZ-1] != SIG_MSB)
1898 return false;
1899
1900 *r = u;
1901 return true;
1902 #endif /* 0 */
1903 }
1904
1905 #if 0
1906 /* Render R as an integer. */
1907
1908 HOST_WIDE_INT
1909 real_to_integer (const REAL_VALUE_TYPE *r)
1910 {
1911 unsigned HOST_WIDE_INT i;
1912
1913 switch (r->cl)
1914 {
1915 case rvc_zero:
1916 underflow:
1917 return 0;
1918
1919 case rvc_inf:
1920 case rvc_nan:
1921 overflow:
1922 i = (unsigned HOST_WIDE_INT) 1 << (HOST_BITS_PER_WIDE_INT - 1);
1923 if (!r->sign)
1924 i--;
1925 return i;
1926
1927 case rvc_normal:
1928 if (REAL_EXP (r) <= 0)
1929 goto underflow;
1930 /* Only force overflow for unsigned overflow. Signed overflow is
1931 undefined, so it doesn't matter what we return, and some callers
1932 expect to be able to use this routine for both signed and
1933 unsigned conversions. */
1934 if (REAL_EXP (r) > HOST_BITS_PER_WIDE_INT)
1935 goto overflow;
1936
1937 if (HOST_BITS_PER_WIDE_INT == HOST_BITS_PER_LONG)
1938 i = r->sig[SIGSZ-1];
1939 else
1940 {
1941 gcc_assert (HOST_BITS_PER_WIDE_INT == 2 * HOST_BITS_PER_LONG);
1942 i = r->sig[SIGSZ-1];
1943 i = i << (HOST_BITS_PER_LONG - 1) << 1;
1944 i |= r->sig[SIGSZ-2];
1945 }
1946
1947 i >>= HOST_BITS_PER_WIDE_INT - REAL_EXP (r);
1948
1949 if (r->sign)
1950 i = -i;
1951 return i;
1952
1953 default:
1954 gcc_unreachable ();
1955 }
1956 }
1957 #endif /* 0 */
1958
1959 /* Likewise, but to an integer pair, HI+LOW. */
1960
1961 void
real_to_integer2(HOST_WIDE_INT * plow,HOST_WIDE_INT * phigh,const REAL_VALUE_TYPE * r)1962 real_to_integer2 (HOST_WIDE_INT *plow, HOST_WIDE_INT *phigh,
1963 const REAL_VALUE_TYPE *r)
1964 {
1965 unsigned short exponent = (r->exponent & 0x7fff);
1966 unsigned HOST_WIDE_INT low = 0, high = 0;
1967 int e, i;
1968 unsigned long long mantissa = r->mantissa;
1969
1970 /* Return 0 for transfinite, zero or |*r|<1 */
1971 if (!REAL_VALUE_ISFINITE(*r) || REAL_VALUE_ISZERO(*r) || exponent < 0x4000) {
1972 *plow = *phigh = 0;
1973 return;
1974 }
1975
1976 /* Compute the effective exponent */
1977 e = (int)exponent - (0x4000+15);
1978
1979 /* Delete any digits following the decimal point right now */
1980 if (e < 0) {
1981 mantissa >>= ((-e)<<2);
1982 e = 0;
1983 } else if (e > 2*HOST_BITS_PER_WIDE_INT) {
1984 /* We'll multiply by 10^(2*HBPWI)=5^(2*HBPWI)*2^(2*HBPWI). This will zero
1985 out all digits. So avoid the long loop and return 0 right now. */
1986 *plow = *phigh = 0;
1987 return;
1988 }
1989
1990 /* Convert the mantissa from BCD to binary */
1991 for (i=60;i>=0;i-=4) {
1992 /* Multiply by 10 using x*10=(x<<3)+(x<<1) */
1993 high = (high<<3) + (low>>(HOST_BITS_PER_WIDE_INT-3)) /* (x<<3) */
1994 + (high<<1) + (low>>(HOST_BITS_PER_WIDE_INT-1)) /* (x<<1) */
1995 + ((low<<3) + (low<<1) < (low<<1)); /* (carry) */
1996 low = (low<<3) + (low<<1);
1997 low += (mantissa>>i)&15;
1998 }
1999
2000 /* Multiply by 10^e */
2001 for (i=0;i<e;i++) {
2002 /* Multiply by 10 using x*10=(x<<3)+(x<<1) */
2003 high = (high<<3) + (low>>(HOST_BITS_PER_WIDE_INT-3)) /* (x<<3) */
2004 + (high<<1) + (low>>(HOST_BITS_PER_WIDE_INT-1)) /* (x<<1) */
2005 + ((low<<3) + (low<<1) < (low<<1)); /* (carry) */
2006 low = (low<<3) + (low<<1);
2007 }
2008
2009 /* Negate the result if *r was negative */
2010 if (r->exponent >= 0x8000) {
2011 if (low == 0)
2012 high = -high;
2013 else
2014 low = -low, high = ~high;
2015 }
2016
2017 #if 0
2018 REAL_VALUE_TYPE t;
2019 HOST_WIDE_INT low, high;
2020 int exp;
2021
2022 switch (r->cl)
2023 {
2024 case rvc_zero:
2025 underflow:
2026 low = high = 0;
2027 break;
2028
2029 case rvc_inf:
2030 case rvc_nan:
2031 overflow:
2032 high = (unsigned HOST_WIDE_INT) 1 << (HOST_BITS_PER_WIDE_INT - 1);
2033 if (r->sign)
2034 low = 0;
2035 else
2036 {
2037 high--;
2038 low = -1;
2039 }
2040 break;
2041
2042 case rvc_normal:
2043 exp = REAL_EXP (r);
2044 if (exp <= 0)
2045 goto underflow;
2046 /* Only force overflow for unsigned overflow. Signed overflow is
2047 undefined, so it doesn't matter what we return, and some callers
2048 expect to be able to use this routine for both signed and
2049 unsigned conversions. */
2050 if (exp > 2*HOST_BITS_PER_WIDE_INT)
2051 goto overflow;
2052
2053 rshift_significand (&t, r, 2*HOST_BITS_PER_WIDE_INT - exp);
2054 if (HOST_BITS_PER_WIDE_INT == HOST_BITS_PER_LONG)
2055 {
2056 high = t.sig[SIGSZ-1];
2057 low = t.sig[SIGSZ-2];
2058 }
2059 else
2060 {
2061 gcc_assert (HOST_BITS_PER_WIDE_INT == 2*HOST_BITS_PER_LONG);
2062 high = t.sig[SIGSZ-1];
2063 high = high << (HOST_BITS_PER_LONG - 1) << 1;
2064 high |= t.sig[SIGSZ-2];
2065
2066 low = t.sig[SIGSZ-3];
2067 low = low << (HOST_BITS_PER_LONG - 1) << 1;
2068 low |= t.sig[SIGSZ-4];
2069 }
2070
2071 if (r->sign)
2072 {
2073 if (low == 0)
2074 high = -high;
2075 else
2076 low = -low, high = ~high;
2077 }
2078 break;
2079
2080 default:
2081 gcc_unreachable ();
2082 }
2083 #endif /* 0 */
2084
2085 *plow = low;
2086 *phigh = high;
2087 }
2088
2089 #if 0
2090 /* A subroutine of real_to_decimal. Compute the quotient and remainder
2091 of NUM / DEN. Return the quotient and place the remainder in NUM.
2092 It is expected that NUM / DEN are close enough that the quotient is
2093 small. */
2094
2095 static unsigned long
2096 rtd_divmod (REAL_VALUE_TYPE *num, REAL_VALUE_TYPE *den)
2097 {
2098 unsigned long q, msb;
2099 int expn = REAL_EXP (num), expd = REAL_EXP (den);
2100
2101 if (expn < expd)
2102 return 0;
2103
2104 q = msb = 0;
2105 goto start;
2106 do
2107 {
2108 msb = num->sig[SIGSZ-1] & SIG_MSB;
2109 q <<= 1;
2110 lshift_significand_1 (num, num);
2111 start:
2112 if (msb || cmp_significands (num, den) >= 0)
2113 {
2114 sub_significands (num, num, den, 0);
2115 q |= 1;
2116 }
2117 }
2118 while (--expn >= expd);
2119
2120 SET_REAL_EXP (num, expd);
2121 normalize (num);
2122
2123 return q;
2124 }
2125
2126 /* Render R as a decimal floating point constant. Emit DIGITS significant
2127 digits in the result, bounded by BUF_SIZE. If DIGITS is 0, choose the
2128 maximum for the representation. If CROP_TRAILING_ZEROS, strip trailing
2129 zeros. */
2130
2131 #define M_LOG10_2 0.30102999566398119521
2132
2133 void
2134 real_to_decimal (char *str, const REAL_VALUE_TYPE *r_orig, size_t buf_size,
2135 size_t digits, int crop_trailing_zeros)
2136 {
2137 const REAL_VALUE_TYPE *one, *ten;
2138 REAL_VALUE_TYPE r, pten, u, v;
2139 int dec_exp, cmp_one, digit;
2140 size_t max_digits;
2141 char *p, *first, *last;
2142 bool sign;
2143
2144 r = *r_orig;
2145 switch (r.cl)
2146 {
2147 case rvc_zero:
2148 strcpy (str, (r.sign ? "-0.0" : "0.0"));
2149 return;
2150 case rvc_normal:
2151 break;
2152 case rvc_inf:
2153 strcpy (str, (r.sign ? "-Inf" : "+Inf"));
2154 return;
2155 case rvc_nan:
2156 /* ??? Print the significand as well, if not canonical? */
2157 strcpy (str, (r.sign ? "-NaN" : "+NaN"));
2158 return;
2159 default:
2160 gcc_unreachable ();
2161 }
2162
2163 /* Bound the number of digits printed by the size of the representation. */
2164 max_digits = SIGNIFICAND_BITS * M_LOG10_2;
2165 if (digits == 0 || digits > max_digits)
2166 digits = max_digits;
2167
2168 /* Estimate the decimal exponent, and compute the length of the string it
2169 will print as. Be conservative and add one to account for possible
2170 overflow or rounding error. */
2171 dec_exp = REAL_EXP (&r) * M_LOG10_2;
2172 for (max_digits = 1; dec_exp ; max_digits++)
2173 dec_exp /= 10;
2174
2175 /* Bound the number of digits printed by the size of the output buffer. */
2176 max_digits = buf_size - 1 - 1 - 2 - max_digits - 1;
2177 gcc_assert (max_digits <= buf_size);
2178 if (digits > max_digits)
2179 digits = max_digits;
2180
2181 one = real_digit (1);
2182 ten = ten_to_ptwo (0);
2183
2184 sign = r.sign;
2185 r.sign = 0;
2186
2187 dec_exp = 0;
2188 pten = *one;
2189
2190 cmp_one = do_compare (&r, one, 0);
2191 if (cmp_one > 0)
2192 {
2193 int m;
2194
2195 /* Number is greater than one. Convert significand to an integer
2196 and strip trailing decimal zeros. */
2197
2198 u = r;
2199 SET_REAL_EXP (&u, SIGNIFICAND_BITS - 1);
2200
2201 /* Largest M, such that 10**2**M fits within SIGNIFICAND_BITS. */
2202 m = floor_log2 (max_digits);
2203
2204 /* Iterate over the bits of the possible powers of 10 that might
2205 be present in U and eliminate them. That is, if we find that
2206 10**2**M divides U evenly, keep the division and increase
2207 DEC_EXP by 2**M. */
2208 do
2209 {
2210 REAL_VALUE_TYPE t;
2211
2212 do_divide (&t, &u, ten_to_ptwo (m));
2213 do_fix_trunc (&v, &t);
2214 if (cmp_significands (&v, &t) == 0)
2215 {
2216 u = t;
2217 dec_exp += 1 << m;
2218 }
2219 }
2220 while (--m >= 0);
2221
2222 /* Revert the scaling to integer that we performed earlier. */
2223 SET_REAL_EXP (&u, REAL_EXP (&u) + REAL_EXP (&r)
2224 - (SIGNIFICAND_BITS - 1));
2225 r = u;
2226
2227 /* Find power of 10. Do this by dividing out 10**2**M when
2228 this is larger than the current remainder. Fill PTEN with
2229 the power of 10 that we compute. */
2230 if (REAL_EXP (&r) > 0)
2231 {
2232 m = floor_log2 ((int)(REAL_EXP (&r) * M_LOG10_2)) + 1;
2233 do
2234 {
2235 const REAL_VALUE_TYPE *ptentwo = ten_to_ptwo (m);
2236 if (do_compare (&u, ptentwo, 0) >= 0)
2237 {
2238 do_divide (&u, &u, ptentwo);
2239 do_multiply (&pten, &pten, ptentwo);
2240 dec_exp += 1 << m;
2241 }
2242 }
2243 while (--m >= 0);
2244 }
2245 else
2246 /* We managed to divide off enough tens in the above reduction
2247 loop that we've now got a negative exponent. Fall into the
2248 less-than-one code to compute the proper value for PTEN. */
2249 cmp_one = -1;
2250 }
2251 if (cmp_one < 0)
2252 {
2253 int m;
2254
2255 /* Number is less than one. Pad significand with leading
2256 decimal zeros. */
2257
2258 v = r;
2259 while (1)
2260 {
2261 /* Stop if we'd shift bits off the bottom. */
2262 if (v.sig[0] & 7)
2263 break;
2264
2265 do_multiply (&u, &v, ten);
2266
2267 /* Stop if we're now >= 1. */
2268 if (REAL_EXP (&u) > 0)
2269 break;
2270
2271 v = u;
2272 dec_exp -= 1;
2273 }
2274 r = v;
2275
2276 /* Find power of 10. Do this by multiplying in P=10**2**M when
2277 the current remainder is smaller than 1/P. Fill PTEN with the
2278 power of 10 that we compute. */
2279 m = floor_log2 ((int)(-REAL_EXP (&r) * M_LOG10_2)) + 1;
2280 do
2281 {
2282 const REAL_VALUE_TYPE *ptentwo = ten_to_ptwo (m);
2283 const REAL_VALUE_TYPE *ptenmtwo = ten_to_mptwo (m);
2284
2285 if (do_compare (&v, ptenmtwo, 0) <= 0)
2286 {
2287 do_multiply (&v, &v, ptentwo);
2288 do_multiply (&pten, &pten, ptentwo);
2289 dec_exp -= 1 << m;
2290 }
2291 }
2292 while (--m >= 0);
2293
2294 /* Invert the positive power of 10 that we've collected so far. */
2295 do_divide (&pten, one, &pten);
2296 }
2297
2298 p = str;
2299 if (sign)
2300 *p++ = '-';
2301 first = p++;
2302
2303 /* At this point, PTEN should contain the nearest power of 10 smaller
2304 than R, such that this division produces the first digit.
2305
2306 Using a divide-step primitive that returns the complete integral
2307 remainder avoids the rounding error that would be produced if
2308 we were to use do_divide here and then simply multiply by 10 for
2309 each subsequent digit. */
2310
2311 digit = rtd_divmod (&r, &pten);
2312
2313 /* Be prepared for error in that division via underflow ... */
2314 if (digit == 0 && cmp_significand_0 (&r))
2315 {
2316 /* Multiply by 10 and try again. */
2317 do_multiply (&r, &r, ten);
2318 digit = rtd_divmod (&r, &pten);
2319 dec_exp -= 1;
2320 gcc_assert (digit != 0);
2321 }
2322
2323 /* ... or overflow. */
2324 if (digit == 10)
2325 {
2326 *p++ = '1';
2327 if (--digits > 0)
2328 *p++ = '0';
2329 dec_exp += 1;
2330 }
2331 else
2332 {
2333 gcc_assert (digit <= 10);
2334 *p++ = digit + '0';
2335 }
2336
2337 /* Generate subsequent digits. */
2338 while (--digits > 0)
2339 {
2340 do_multiply (&r, &r, ten);
2341 digit = rtd_divmod (&r, &pten);
2342 *p++ = digit + '0';
2343 }
2344 last = p;
2345
2346 /* Generate one more digit with which to do rounding. */
2347 do_multiply (&r, &r, ten);
2348 digit = rtd_divmod (&r, &pten);
2349
2350 /* Round the result. */
2351 if (digit == 5)
2352 {
2353 /* Round to nearest. If R is nonzero there are additional
2354 nonzero digits to be extracted. */
2355 if (cmp_significand_0 (&r))
2356 digit++;
2357 /* Round to even. */
2358 else if ((p[-1] - '0') & 1)
2359 digit++;
2360 }
2361 if (digit > 5)
2362 {
2363 while (p > first)
2364 {
2365 digit = *--p;
2366 if (digit == '9')
2367 *p = '0';
2368 else
2369 {
2370 *p = digit + 1;
2371 break;
2372 }
2373 }
2374
2375 /* Carry out of the first digit. This means we had all 9's and
2376 now have all 0's. "Prepend" a 1 by overwriting the first 0. */
2377 if (p == first)
2378 {
2379 first[1] = '1';
2380 dec_exp++;
2381 }
2382 }
2383
2384 /* Insert the decimal point. */
2385 first[0] = first[1];
2386 first[1] = '.';
2387
2388 /* If requested, drop trailing zeros. Never crop past "1.0". */
2389 if (crop_trailing_zeros)
2390 while (last > first + 3 && last[-1] == '0')
2391 last--;
2392
2393 /* Append the exponent. */
2394 sprintf (last, "e%+d", dec_exp);
2395 }
2396
2397 /* Render R as a hexadecimal floating point constant. Emit DIGITS
2398 significant digits in the result, bounded by BUF_SIZE. If DIGITS is 0,
2399 choose the maximum for the representation. If CROP_TRAILING_ZEROS,
2400 strip trailing zeros. */
2401
2402 void
2403 real_to_hexadecimal (char *str, const REAL_VALUE_TYPE *r, size_t buf_size,
2404 size_t digits, int crop_trailing_zeros)
2405 {
2406 int i, j, exp = REAL_EXP (r);
2407 char *p, *first;
2408 char exp_buf[16];
2409 size_t max_digits;
2410
2411 switch (r->cl)
2412 {
2413 case rvc_zero:
2414 exp = 0;
2415 break;
2416 case rvc_normal:
2417 break;
2418 case rvc_inf:
2419 strcpy (str, (r->sign ? "-Inf" : "+Inf"));
2420 return;
2421 case rvc_nan:
2422 /* ??? Print the significand as well, if not canonical? */
2423 strcpy (str, (r->sign ? "-NaN" : "+NaN"));
2424 return;
2425 default:
2426 gcc_unreachable ();
2427 }
2428
2429 if (digits == 0)
2430 digits = SIGNIFICAND_BITS / 4;
2431
2432 /* Bound the number of digits printed by the size of the output buffer. */
2433
2434 sprintf (exp_buf, "p%+d", exp);
2435 max_digits = buf_size - strlen (exp_buf) - r->sign - 4 - 1;
2436 gcc_assert (max_digits <= buf_size);
2437 if (digits > max_digits)
2438 digits = max_digits;
2439
2440 p = str;
2441 if (r->sign)
2442 *p++ = '-';
2443 *p++ = '0';
2444 *p++ = 'x';
2445 *p++ = '0';
2446 *p++ = '.';
2447 first = p;
2448
2449 for (i = SIGSZ - 1; i >= 0; --i)
2450 for (j = HOST_BITS_PER_LONG - 4; j >= 0; j -= 4)
2451 {
2452 *p++ = "0123456789abcdef"[(r->sig[i] >> j) & 15];
2453 if (--digits == 0)
2454 goto out;
2455 }
2456
2457 out:
2458 if (crop_trailing_zeros)
2459 while (p > first + 1 && p[-1] == '0')
2460 p--;
2461
2462 sprintf (p, "p%+d", exp);
2463 }
2464 #endif /* 0 */
2465
real_value_dtof(REAL_VALUE_TYPE * r,const char * string)2466 static void real_value_dtof (REAL_VALUE_TYPE *r, const char *string)
2467 {
2468 const char *strpart = string;
2469 unsigned char state = 0;
2470 unsigned long long mpmul = 0x1000000000000000ull;
2471 signed short expshift = -1;
2472 unsigned short exp = 0;
2473 unsigned short expmul = 1;
2474 unsigned short expadd = 0;
2475 int endrounding = -1;
2476 r->mantissa = 0;
2477 while (strpart && *strpart)
2478 {
2479 switch (state)
2480 {
2481 case 0:
2482 if (*strpart == '.')
2483 {
2484 state = 1;
2485 break;
2486 }
2487 else if (*strpart == '+')
2488 expadd = 0;
2489 else if (*strpart == '-')
2490 expadd = 0x8000;
2491 else if (((*strpart) >= '0' && (*strpart) <= '9') && ((r->mantissa) || (*strpart != '0')))
2492 expshift++;
2493 else if ((*strpart == 'i') || (*strpart == 'I'))
2494 {
2495 if (expadd == 0x8000)
2496 *r = NEGATIVE_INF;
2497 else
2498 *r = POSITIVE_INF;
2499 expadd = 0x3FFF;
2500 strpart = 0;
2501 }
2502 else if ((*strpart == 'n') || (*strpart == 'N'))
2503 {
2504 *r = NAN;
2505 expadd = 0x3FFF;
2506 strpart = 0;
2507 }
2508 case 1:
2509 if ((*strpart == 'e') || (*strpart == 'E'))
2510 state = 2;
2511 else if ((*strpart) >= '0' && (*strpart) <= '9')
2512 {
2513 if (mpmul)
2514 r->mantissa |= (*strpart - '0') * mpmul;
2515 else if (endrounding < 0)
2516 endrounding = (*strpart >= '5');
2517 if (r->mantissa)
2518 mpmul /= 0x10;
2519 else if (state)
2520 expshift--;
2521 }
2522 break;
2523 case 2:
2524 if (*strpart == '+')
2525 expmul = 1;
2526 else if (*strpart == '-')
2527 expmul = -1;
2528 if ((*strpart) >= '0' && (*strpart) <= '9')
2529 {
2530 exp *= 10;
2531 exp += (*strpart - '0');
2532 }
2533 break;
2534 }
2535 if (strpart)
2536 strpart++;
2537 else
2538 break;
2539 }
2540 if (!(r->mantissa))
2541 *r = ZERO;
2542 else if (expadd != 0x3FFF)
2543 {
2544 r->exponent = exp * expmul + 0x4000 + expshift;
2545 if (r->exponent>0x4000+16383) /* exponent too large, overflow to +infinity */
2546 *r = POSITIVE_INF;
2547 else if (r->exponent<0x4000-16383) /* exponent too small, underflow to +0 */
2548 *r = POSITIVE_ZERO;
2549 else if (endrounding > 0)
2550 bcdpadd1ulp (r);
2551 r->exponent += expadd;
2552 }
2553 }
2554
2555 typedef struct {
2556 int ndigits;
2557 unsigned char *digits;
2558 long effexp;
2559 } arbprec_decimal;
2560
arbprec_pack(arbprec_decimal * r)2561 static void arbprec_pack (arbprec_decimal *r)
2562 {
2563 unsigned char *p = r->digits;
2564 if (r->ndigits) {
2565 while (p < r->digits + r->ndigits && !*p) p++;
2566 r->ndigits -= p - r->digits;
2567 memmove (r->digits, p, r->ndigits);
2568 if (r->ndigits) {
2569 p = r->digits + r->ndigits - 1;
2570 while (p > r->digits && !*p) {
2571 p--;
2572 r->effexp++;
2573 }
2574 r->ndigits = (p + 1) - r->digits;
2575 }
2576 r->digits = xrealloc (r->digits, r->ndigits);
2577 }
2578 }
2579
arbprec_mul2(arbprec_decimal * r)2580 static void arbprec_mul2 (arbprec_decimal *r)
2581 {
2582 unsigned char *p;
2583 unsigned char carry = 0;
2584 r->ndigits++;
2585 r->digits = xrealloc (r->digits, r->ndigits);
2586 for (p = r->digits + r->ndigits - 1; p > r->digits; p--) {
2587 unsigned char digit = (p[-1]<<1) + carry;
2588 *p = digit % 10;
2589 carry = digit / 10;
2590 }
2591 *p = carry;
2592 arbprec_pack (r);
2593 }
2594
arbprec_div2(arbprec_decimal * r)2595 static void arbprec_div2 (arbprec_decimal *r)
2596 {
2597 unsigned char *p;
2598 unsigned char carry = 0;
2599 r->ndigits++;
2600 r->digits = xrealloc (r->digits, r->ndigits);
2601 for (p = r->digits; p < r->digits + r->ndigits - 1; p++) {
2602 unsigned char digit = (*p>>1) + carry;
2603 carry = (*p&1)*5;
2604 *p = digit;
2605 }
2606 *p = carry;
2607 r->effexp--;
2608 arbprec_pack (r);
2609 }
2610
arbprec_mul16(arbprec_decimal * r)2611 static void arbprec_mul16 (arbprec_decimal *r)
2612 {
2613 arbprec_mul2 (r);
2614 arbprec_mul2 (r);
2615 arbprec_mul2 (r);
2616 arbprec_mul2 (r);
2617 }
2618
arbprec_div16(arbprec_decimal * r)2619 static void arbprec_div16 (arbprec_decimal *r)
2620 {
2621 arbprec_div2 (r);
2622 arbprec_div2 (r);
2623 arbprec_div2 (r);
2624 arbprec_div2 (r);
2625 }
2626
arbprec_add(arbprec_decimal * r1,arbprec_decimal * r2)2627 static void arbprec_add (arbprec_decimal *r1, arbprec_decimal *r2)
2628 {
2629 if (!r2->ndigits)
2630 return;
2631 else if (!r1->ndigits) {
2632 r1->digits = xrealloc (r1->digits, r2->ndigits);
2633 memcpy (r1->digits, r2->digits, r2->ndigits);
2634 r1->ndigits = r2->ndigits;
2635 r1->effexp = r2->effexp;
2636 } else {
2637 if (r1->effexp > r2->effexp) {
2638 long effexpdiff = r1->effexp - r2->effexp;
2639 r1->digits = xrealloc (r1->digits, r1->ndigits + effexpdiff);
2640 memset (r1->digits + r1->ndigits, 0, effexpdiff);
2641 r1->ndigits += effexpdiff;
2642 r1->effexp = r2->effexp;
2643 } else if (r2->effexp > r1->effexp) {
2644 long effexpdiff = r2->effexp - r1->effexp;
2645 r2->digits = xrealloc (r2->digits, r2->ndigits + effexpdiff);
2646 memset (r2->digits + r2->ndigits, 0, effexpdiff);
2647 r2->ndigits += effexpdiff;
2648 r2->effexp = r1->effexp;
2649 }
2650
2651 {
2652 int ndigits = MAX (r1->ndigits, r2->ndigits) + 1;
2653 unsigned char *digits = xmalloc (ndigits);
2654 unsigned char *p, *q1, *q2;
2655 unsigned char carry = 0;
2656 for (p = digits + ndigits - 1, q1 = r1->digits + r1->ndigits - 1,
2657 q2 = r2->digits + r2->ndigits - 1; p > digits; p--, q1--, q2--) {
2658 unsigned char digit = (q1>=r1->digits?*q1:0) + (q2>=r2->digits?*q2:0)
2659 + carry;
2660 *p = digit % 10;
2661 carry = digit / 10;
2662 }
2663 *p = carry;
2664 free (r1->digits);
2665 r1->digits = digits;
2666 r1->ndigits = ndigits;
2667 arbprec_pack (r1);
2668 arbprec_pack (r2);
2669 }
2670 }
2671 }
2672
arbprec_add_n_times(arbprec_decimal * r1,arbprec_decimal * r2,int n)2673 static void arbprec_add_n_times (arbprec_decimal *r1, arbprec_decimal *r2, int n)
2674 {
2675 int i;
2676 for (i = 0; i < n; i++) arbprec_add (r1, r2);
2677 }
2678
arbprec_to_bcd(arbprec_decimal * a,smap_bcd_float * r)2679 static void arbprec_to_bcd (arbprec_decimal *a, smap_bcd_float *r)
2680 {
2681 if (a->ndigits) {
2682 long exponent = a->effexp + a->ndigits - 1;
2683 if (exponent>16383) /* exponent too large, overflow to +infinity */
2684 *r = POSITIVE_INF;
2685 else if (exponent<-16383) /* exponent too small, underflow to +0 */
2686 *r = POSITIVE_ZERO;
2687 else {
2688 int i;
2689 r->exponent = 0x4000 + exponent;
2690 r->mantissa = 0;
2691 for (i = 0; i < a->ndigits && i < 16; i++) {
2692 r->mantissa = (r->mantissa << 4) + a->digits[i];
2693 }
2694 for (; i < 16; i++) {
2695 r->mantissa <<= 4;
2696 }
2697 if (a->ndigits > 16 && a->digits[16] >= 5)
2698 bcdpadd1ulp (r);
2699 }
2700 } else {
2701 *r = UNSIGNED_ZERO;
2702 }
2703 }
2704
real_value_htof(REAL_VALUE_TYPE * res,const char * string)2705 static void real_value_htof (REAL_VALUE_TYPE *res, const char *string)
2706 {
2707 arbprec_decimal r = {0, NULL, 0};
2708 const char *strpart = string;
2709 signed char state = -1;
2710 arbprec_decimal rdiv = {1, NULL, 0};
2711 unsigned int negative = 0;
2712 unsigned short exp = 0;
2713 unsigned short expsign = 1;
2714 rdiv.digits = xmalloc (1);
2715 *(rdiv.digits) = 1;
2716 while (strpart && *strpart)
2717 {
2718 switch (state)
2719 {
2720 case -1:
2721 if (*strpart == 'x' || *strpart == 'X')
2722 state = 0;
2723 else if (*strpart == '+')
2724 negative = 0;
2725 else if (*strpart == '-')
2726 negative = 1;
2727 break;
2728 case 0:
2729 if (*strpart == '.')
2730 {
2731 state = 1;
2732 break;
2733 }
2734 else if ((*strpart >= '0' && *strpart <= '9') || (*strpart >= 'a' && *strpart <= 'f') || (*strpart >= 'A' && *strpart <= 'F'))
2735 arbprec_mul16 (&r);
2736 case 1:
2737 if (state == 1)
2738 arbprec_div16 (&rdiv);
2739 if ((*strpart == 'p') || (*strpart == 'P'))
2740 state = 2;
2741 else if (*strpart >= '0' && *strpart <= '9')
2742 arbprec_add_n_times (&r, &rdiv, (*strpart - '0'));
2743 else if (*strpart >= 'a' && *strpart <= 'f')
2744 arbprec_add_n_times (&r, &rdiv, (*strpart - 'a' + 0xA));
2745 else if (*strpart >= 'A' && *strpart <= 'F')
2746 arbprec_add_n_times (&r, &rdiv, (*strpart - 'A' + 0xA));
2747 break;
2748 case 2:
2749 if (*strpart == '+')
2750 expsign = 1;
2751 else if (*strpart == '-')
2752 expsign = -1;
2753 if ((*strpart) >= '0' && (*strpart) <= '9')
2754 {
2755 exp *= 10;
2756 exp += (*strpart - '0');
2757 }
2758 break;
2759 }
2760 strpart++;
2761 }
2762 free (rdiv.digits);
2763 if (expsign == 1)
2764 {
2765 while (exp)
2766 {
2767 arbprec_mul2 (&r);
2768 exp--;
2769 }
2770 }
2771 else
2772 {
2773 while (exp)
2774 {
2775 arbprec_div2 (&r);
2776 exp--;
2777 }
2778 }
2779 arbprec_to_bcd (&r, res);
2780 free (r.digits);
2781 if (negative)
2782 *res = REAL_VALUE_NEGATE (*res);
2783 }
2784
2785 /* Initialize R from a decimal or hexadecimal string. The string is
2786 assumed to have been syntax checked already. */
2787
2788 void
real_from_string(REAL_VALUE_TYPE * r,const char * str)2789 real_from_string (REAL_VALUE_TYPE *r, const char *str)
2790 {
2791 const char *p = str;
2792
2793 if (*p == '-' || *p == '+')
2794 p++;
2795
2796 if (p[0] == '0' && (p[1] == 'x' || p[1] == 'X'))
2797 real_value_htof (r, str);
2798 else
2799 real_value_dtof (r, str);
2800
2801 #if 0
2802 int exp = 0;
2803 bool sign = false;
2804
2805 get_zero (r, 0);
2806
2807 if (*str == '-')
2808 {
2809 sign = true;
2810 str++;
2811 }
2812 else if (*str == '+')
2813 str++;
2814
2815 if (str[0] == '0' && (str[1] == 'x' || str[1] == 'X'))
2816 {
2817 /* Hexadecimal floating point. */
2818 int pos = SIGNIFICAND_BITS - 4, d;
2819
2820 str += 2;
2821
2822 while (*str == '0')
2823 str++;
2824 while (1)
2825 {
2826 d = hex_value (*str);
2827 if (d == _hex_bad)
2828 break;
2829 if (pos >= 0)
2830 {
2831 r->sig[pos / HOST_BITS_PER_LONG]
2832 |= (unsigned long) d << (pos % HOST_BITS_PER_LONG);
2833 pos -= 4;
2834 }
2835 else if (d)
2836 /* Ensure correct rounding by setting last bit if there is
2837 a subsequent nonzero digit. */
2838 r->sig[0] |= 1;
2839 exp += 4;
2840 str++;
2841 }
2842 if (*str == '.')
2843 {
2844 str++;
2845 if (pos == SIGNIFICAND_BITS - 4)
2846 {
2847 while (*str == '0')
2848 str++, exp -= 4;
2849 }
2850 while (1)
2851 {
2852 d = hex_value (*str);
2853 if (d == _hex_bad)
2854 break;
2855 if (pos >= 0)
2856 {
2857 r->sig[pos / HOST_BITS_PER_LONG]
2858 |= (unsigned long) d << (pos % HOST_BITS_PER_LONG);
2859 pos -= 4;
2860 }
2861 else if (d)
2862 /* Ensure correct rounding by setting last bit if there is
2863 a subsequent nonzero digit. */
2864 r->sig[0] |= 1;
2865 str++;
2866 }
2867 }
2868 if (*str == 'p' || *str == 'P')
2869 {
2870 bool exp_neg = false;
2871
2872 str++;
2873 if (*str == '-')
2874 {
2875 exp_neg = true;
2876 str++;
2877 }
2878 else if (*str == '+')
2879 str++;
2880
2881 d = 0;
2882 while (ISDIGIT (*str))
2883 {
2884 d *= 10;
2885 d += *str - '0';
2886 if (d > MAX_EXP)
2887 {
2888 /* Overflowed the exponent. */
2889 if (exp_neg)
2890 goto underflow;
2891 else
2892 goto overflow;
2893 }
2894 str++;
2895 }
2896 if (exp_neg)
2897 d = -d;
2898
2899 exp += d;
2900 }
2901
2902 r->cl = rvc_normal;
2903 SET_REAL_EXP (r, exp);
2904
2905 normalize (r);
2906 }
2907 else
2908 {
2909 /* Decimal floating point. */
2910 const REAL_VALUE_TYPE *ten = ten_to_ptwo (0);
2911 int d;
2912
2913 while (*str == '0')
2914 str++;
2915 while (ISDIGIT (*str))
2916 {
2917 d = *str++ - '0';
2918 do_multiply (r, r, ten);
2919 if (d)
2920 do_add (r, r, real_digit (d), 0);
2921 }
2922 if (*str == '.')
2923 {
2924 str++;
2925 if (r->cl == rvc_zero)
2926 {
2927 while (*str == '0')
2928 str++, exp--;
2929 }
2930 while (ISDIGIT (*str))
2931 {
2932 d = *str++ - '0';
2933 do_multiply (r, r, ten);
2934 if (d)
2935 do_add (r, r, real_digit (d), 0);
2936 exp--;
2937 }
2938 }
2939
2940 if (*str == 'e' || *str == 'E')
2941 {
2942 bool exp_neg = false;
2943
2944 str++;
2945 if (*str == '-')
2946 {
2947 exp_neg = true;
2948 str++;
2949 }
2950 else if (*str == '+')
2951 str++;
2952
2953 d = 0;
2954 while (ISDIGIT (*str))
2955 {
2956 d *= 10;
2957 d += *str - '0';
2958 if (d > MAX_EXP)
2959 {
2960 /* Overflowed the exponent. */
2961 if (exp_neg)
2962 goto underflow;
2963 else
2964 goto overflow;
2965 }
2966 str++;
2967 }
2968 if (exp_neg)
2969 d = -d;
2970 exp += d;
2971 }
2972
2973 if (exp)
2974 times_pten (r, exp);
2975 }
2976
2977 r->sign = sign;
2978 return;
2979
2980 underflow:
2981 get_zero (r, sign);
2982 return;
2983
2984 overflow:
2985 get_inf (r, sign);
2986 return;
2987 #endif /* 0 */
2988 }
2989
2990 /* Legacy. Similar, but return the result directly. */
2991
2992 REAL_VALUE_TYPE
real_from_string2(const char * s,enum machine_mode mode)2993 real_from_string2 (const char *s, enum machine_mode mode)
2994 {
2995 REAL_VALUE_TYPE r;
2996
2997 real_from_string (&r, s);
2998 if (mode != VOIDmode)
2999 real_convert (&r, mode, &r);
3000
3001 return r;
3002 }
3003
3004 /* Initialize R from the integer pair HIGH+LOW. */
3005
3006 void
real_from_integer(REAL_VALUE_TYPE * r,enum machine_mode mode,unsigned HOST_WIDE_INT low,HOST_WIDE_INT high,int unsigned_p)3007 real_from_integer (REAL_VALUE_TYPE *r, enum machine_mode mode,
3008 unsigned HOST_WIDE_INT low, HOST_WIDE_INT high,
3009 int unsigned_p)
3010 {
3011 int j, e = -1, lastdigit = 0;
3012 bool negative = high < 0 && !unsigned_p;
3013 unsigned HOST_WIDE_INT uhigh = (unsigned HOST_WIDE_INT)high;
3014 unsigned HOST_WIDE_INT ulow = low;
3015 unsigned HOST_WIDE_INT remainderh, remainderl, divisorh, divisorl;
3016 unsigned long long mantissa = 0;
3017
3018 if (negative) {
3019 uhigh = ~uhigh;
3020 if (!ulow) uhigh++; else ulow = -ulow;
3021 }
3022
3023 /* Convert the integer to BCD */
3024 while (uhigh || ulow) {
3025 remainderh = uhigh; remainderl = ulow;
3026 divisorl = uhigh = ulow = 0;
3027 divisorh = (HOST_WIDE_INT)10<<(HOST_BITS_PER_WIDE_INT-4);
3028 /* Divide the pair of HOST_WIDE_INTs by 10 */
3029 for (j=2*HOST_BITS_PER_WIDE_INT-4;j>=0;j--) {
3030 /* Shift the result to the left. */
3031 uhigh = (uhigh<<1) + (ulow>>(HOST_BITS_PER_WIDE_INT-1));
3032 ulow <<= 1;
3033 /* Check if the divisor 10 fits into the dividend. */
3034 if (divisorh<remainderh || (divisorh==remainderh && divisorl<=remainderl)) {
3035 /* Add 1 to the result. */
3036 ulow++;
3037 if (!ulow) uhigh++; /* handle carry */
3038 /* Subtract the divisor from the dividend. */
3039 if (remainderl<divisorl) /* handle carry, use unsigned wraparound overflow */
3040 remainderh--;
3041 remainderl -= divisorl; /* now do the subtraction */
3042 remainderh -= divisorh;
3043 }
3044 if (j) {
3045 /* Shift the divisor to the right. */
3046 divisorl = ((divisorh&1)<<(HOST_BITS_PER_WIDE_INT-1))+(divisorl>>1);
3047 divisorh >>= 1;
3048 }
3049 }
3050 e++;
3051 lastdigit = (mantissa&15);
3052 mantissa = (mantissa>>4) + ((unsigned long long)remainderl<<60);
3053 }
3054
3055 if (e >= 0) {
3056 /* We have a nonnegative exponent. Do the correct rounding for the last
3057 digit. */
3058
3059 r->exponent = e+0x4000;
3060 r->mantissa = mantissa;
3061
3062 if (lastdigit >= 5) bcdpadd1ulp(r);
3063 } else {
3064 /* We don't have any digit, so our number is actually zero. */
3065 *r = UNSIGNED_ZERO;
3066 }
3067
3068 /* Negate *r if the integer was negative */
3069 if (negative) *r = REAL_VALUE_NEGATE (*r);
3070
3071 #if 0
3072 if (low == 0 && high == 0)
3073 get_zero (r, 0);
3074 else
3075 {
3076 memset (r, 0, sizeof (*r));
3077 r->cl = rvc_normal;
3078 r->sign = high < 0 && !unsigned_p;
3079 SET_REAL_EXP (r, 2 * HOST_BITS_PER_WIDE_INT);
3080
3081 if (r->sign)
3082 {
3083 high = ~high;
3084 if (low == 0)
3085 high += 1;
3086 else
3087 low = -low;
3088 }
3089
3090 if (HOST_BITS_PER_LONG == HOST_BITS_PER_WIDE_INT)
3091 {
3092 r->sig[SIGSZ-1] = high;
3093 r->sig[SIGSZ-2] = low;
3094 }
3095 else
3096 {
3097 gcc_assert (HOST_BITS_PER_LONG*2 == HOST_BITS_PER_WIDE_INT);
3098 r->sig[SIGSZ-1] = high >> (HOST_BITS_PER_LONG - 1) >> 1;
3099 r->sig[SIGSZ-2] = high;
3100 r->sig[SIGSZ-3] = low >> (HOST_BITS_PER_LONG - 1) >> 1;
3101 r->sig[SIGSZ-4] = low;
3102 }
3103
3104 normalize (r);
3105 }
3106
3107 if (mode != VOIDmode)
3108 real_convert (r, mode, r);
3109 #endif /* 0 */
3110 }
3111
3112 #if 0
3113 /* Returns 10**2**N. */
3114
3115 static const REAL_VALUE_TYPE *
3116 ten_to_ptwo (int n)
3117 {
3118 static REAL_VALUE_TYPE tens[EXP_BITS];
3119
3120 gcc_assert (n >= 0);
3121 gcc_assert (n < EXP_BITS);
3122
3123 if (tens[n].cl == rvc_zero)
3124 {
3125 if (n < (HOST_BITS_PER_WIDE_INT == 64 ? 5 : 4))
3126 {
3127 HOST_WIDE_INT t = 10;
3128 int i;
3129
3130 for (i = 0; i < n; ++i)
3131 t *= t;
3132
3133 real_from_integer (&tens[n], VOIDmode, t, 0, 1);
3134 }
3135 else
3136 {
3137 const REAL_VALUE_TYPE *t = ten_to_ptwo (n - 1);
3138 do_multiply (&tens[n], t, t);
3139 }
3140 }
3141
3142 return &tens[n];
3143 }
3144
3145 /* Returns 10**(-2**N). */
3146
3147 static const REAL_VALUE_TYPE *
3148 ten_to_mptwo (int n)
3149 {
3150 static REAL_VALUE_TYPE tens[EXP_BITS];
3151
3152 gcc_assert (n >= 0);
3153 gcc_assert (n < EXP_BITS);
3154
3155 if (tens[n].cl == rvc_zero)
3156 do_divide (&tens[n], real_digit (1), ten_to_ptwo (n));
3157
3158 return &tens[n];
3159 }
3160
3161 /* Returns N. */
3162
3163 static const REAL_VALUE_TYPE *
3164 real_digit (int n)
3165 {
3166 static REAL_VALUE_TYPE num[10];
3167
3168 gcc_assert (n >= 0);
3169 gcc_assert (n <= 9);
3170
3171 if (n > 0 && num[n].cl == rvc_zero)
3172 real_from_integer (&num[n], VOIDmode, n, 0, 1);
3173
3174 return &num[n];
3175 }
3176
3177 /* Multiply R by 10**EXP. */
3178
3179 static void
3180 times_pten (REAL_VALUE_TYPE *r, int exp)
3181 {
3182 REAL_VALUE_TYPE pten, *rr;
3183 bool negative = (exp < 0);
3184 int i;
3185
3186 if (negative)
3187 {
3188 exp = -exp;
3189 pten = *real_digit (1);
3190 rr = &pten;
3191 }
3192 else
3193 rr = r;
3194
3195 for (i = 0; exp > 0; ++i, exp >>= 1)
3196 if (exp & 1)
3197 do_multiply (rr, rr, ten_to_ptwo (i));
3198
3199 if (negative)
3200 do_divide (r, r, &pten);
3201 }
3202
3203 /* Fills R with +Inf. */
3204
3205 void
3206 real_inf (REAL_VALUE_TYPE *r)
3207 {
3208 get_inf (r, 0);
3209 }
3210
3211 /* Fills R with a NaN whose significand is described by STR. If QUIET,
3212 we force a QNaN, else we force an SNaN. The string, if not empty,
3213 is parsed as a number and placed in the significand. Return true
3214 if the string was successfully parsed. */
3215
3216 bool
3217 real_nan (REAL_VALUE_TYPE *r, const char *str, int quiet,
3218 enum machine_mode mode)
3219 {
3220 const struct real_format *fmt;
3221
3222 fmt = REAL_MODE_FORMAT (mode);
3223 gcc_assert (fmt);
3224
3225 if (*str == 0)
3226 {
3227 if (quiet)
3228 get_canonical_qnan (r, 0);
3229 else
3230 get_canonical_snan (r, 0);
3231 }
3232 else
3233 {
3234 int base = 10, d;
3235
3236 memset (r, 0, sizeof (*r));
3237 r->cl = rvc_nan;
3238
3239 /* Parse akin to strtol into the significand of R. */
3240
3241 while (ISSPACE (*str))
3242 str++;
3243 if (*str == '-')
3244 str++;
3245 else if (*str == '+')
3246 str++;
3247 if (*str == '0')
3248 {
3249 if (*++str == 'x')
3250 str++, base = 16;
3251 else
3252 base = 8;
3253 }
3254
3255 while ((d = hex_value (*str)) < base)
3256 {
3257 REAL_VALUE_TYPE u;
3258
3259 switch (base)
3260 {
3261 case 8:
3262 lshift_significand (r, r, 3);
3263 break;
3264 case 16:
3265 lshift_significand (r, r, 4);
3266 break;
3267 case 10:
3268 lshift_significand_1 (&u, r);
3269 lshift_significand (r, r, 3);
3270 add_significands (r, r, &u);
3271 break;
3272 default:
3273 gcc_unreachable ();
3274 }
3275
3276 get_zero (&u, 0);
3277 u.sig[0] = d;
3278 add_significands (r, r, &u);
3279
3280 str++;
3281 }
3282
3283 /* Must have consumed the entire string for success. */
3284 if (*str != 0)
3285 return false;
3286
3287 /* Shift the significand into place such that the bits
3288 are in the most significant bits for the format. */
3289 lshift_significand (r, r, SIGNIFICAND_BITS - fmt->pnan);
3290
3291 /* Our MSB is always unset for NaNs. */
3292 r->sig[SIGSZ-1] &= ~SIG_MSB;
3293
3294 /* Force quiet or signalling NaN. */
3295 r->signalling = !quiet;
3296 }
3297
3298 return true;
3299 }
3300
3301 /* Fills R with the largest finite value representable in mode MODE.
3302 If SIGN is nonzero, R is set to the most negative finite value. */
3303
3304 void
3305 real_maxval (REAL_VALUE_TYPE *r, int sign, enum machine_mode mode)
3306 {
3307 const struct real_format *fmt;
3308 int np2;
3309
3310 fmt = REAL_MODE_FORMAT (mode);
3311 gcc_assert (fmt);
3312
3313 r->cl = rvc_normal;
3314 r->sign = sign;
3315 r->signalling = 0;
3316 r->canonical = 0;
3317 SET_REAL_EXP (r, fmt->emax * fmt->log2_b);
3318
3319 np2 = SIGNIFICAND_BITS - fmt->p * fmt->log2_b;
3320 memset (r->sig, -1, SIGSZ * sizeof (unsigned long));
3321 clear_significand_below (r, np2);
3322 }
3323
3324 /* Fills R with 2**N. */
3325
3326 void
3327 real_2expN (REAL_VALUE_TYPE *r, int n)
3328 {
3329 memset (r, 0, sizeof (*r));
3330
3331 n++;
3332 if (n > MAX_EXP)
3333 r->cl = rvc_inf;
3334 else if (n < -MAX_EXP)
3335 ;
3336 else
3337 {
3338 r->cl = rvc_normal;
3339 SET_REAL_EXP (r, n);
3340 r->sig[SIGSZ-1] = SIG_MSB;
3341 }
3342 }
3343
3344
3345 static void
3346 round_for_format (const struct real_format *fmt, REAL_VALUE_TYPE *r)
3347 {
3348 int p2, np2, i, w;
3349 unsigned long sticky;
3350 bool guard, lsb;
3351 int emin2m1, emax2;
3352
3353 p2 = fmt->p * fmt->log2_b;
3354 emin2m1 = (fmt->emin - 1) * fmt->log2_b;
3355 emax2 = fmt->emax * fmt->log2_b;
3356
3357 np2 = SIGNIFICAND_BITS - p2;
3358 switch (r->cl)
3359 {
3360 underflow:
3361 get_zero (r, r->sign);
3362 case rvc_zero:
3363 if (!fmt->has_signed_zero)
3364 r->sign = 0;
3365 return;
3366
3367 overflow:
3368 get_inf (r, r->sign);
3369 case rvc_inf:
3370 return;
3371
3372 case rvc_nan:
3373 clear_significand_below (r, np2);
3374 return;
3375
3376 case rvc_normal:
3377 break;
3378
3379 default:
3380 gcc_unreachable ();
3381 }
3382
3383 /* If we're not base2, normalize the exponent to a multiple of
3384 the true base. */
3385 if (fmt->log2_b != 1)
3386 {
3387 int shift = REAL_EXP (r) & (fmt->log2_b - 1);
3388 if (shift)
3389 {
3390 shift = fmt->log2_b - shift;
3391 r->sig[0] |= sticky_rshift_significand (r, r, shift);
3392 SET_REAL_EXP (r, REAL_EXP (r) + shift);
3393 }
3394 }
3395
3396 /* Check the range of the exponent. If we're out of range,
3397 either underflow or overflow. */
3398 if (REAL_EXP (r) > emax2)
3399 goto overflow;
3400 else if (REAL_EXP (r) <= emin2m1)
3401 {
3402 int diff;
3403
3404 if (!fmt->has_denorm)
3405 {
3406 /* Don't underflow completely until we've had a chance to round. */
3407 if (REAL_EXP (r) < emin2m1)
3408 goto underflow;
3409 }
3410 else
3411 {
3412 diff = emin2m1 - REAL_EXP (r) + 1;
3413 if (diff > p2)
3414 goto underflow;
3415
3416 /* De-normalize the significand. */
3417 r->sig[0] |= sticky_rshift_significand (r, r, diff);
3418 SET_REAL_EXP (r, REAL_EXP (r) + diff);
3419 }
3420 }
3421
3422 /* There are P2 true significand bits, followed by one guard bit,
3423 followed by one sticky bit, followed by stuff. Fold nonzero
3424 stuff into the sticky bit. */
3425
3426 sticky = 0;
3427 for (i = 0, w = (np2 - 1) / HOST_BITS_PER_LONG; i < w; ++i)
3428 sticky |= r->sig[i];
3429 sticky |=
3430 r->sig[w] & (((unsigned long)1 << ((np2 - 1) % HOST_BITS_PER_LONG)) - 1);
3431
3432 guard = test_significand_bit (r, np2 - 1);
3433 lsb = test_significand_bit (r, np2);
3434
3435 /* Round to even. */
3436 if (guard && (sticky || lsb))
3437 {
3438 REAL_VALUE_TYPE u;
3439 get_zero (&u, 0);
3440 set_significand_bit (&u, np2);
3441
3442 if (add_significands (r, r, &u))
3443 {
3444 /* Overflow. Means the significand had been all ones, and
3445 is now all zeros. Need to increase the exponent, and
3446 possibly re-normalize it. */
3447 SET_REAL_EXP (r, REAL_EXP (r) + 1);
3448 if (REAL_EXP (r) > emax2)
3449 goto overflow;
3450 r->sig[SIGSZ-1] = SIG_MSB;
3451
3452 if (fmt->log2_b != 1)
3453 {
3454 int shift = REAL_EXP (r) & (fmt->log2_b - 1);
3455 if (shift)
3456 {
3457 shift = fmt->log2_b - shift;
3458 rshift_significand (r, r, shift);
3459 SET_REAL_EXP (r, REAL_EXP (r) + shift);
3460 if (REAL_EXP (r) > emax2)
3461 goto overflow;
3462 }
3463 }
3464 }
3465 }
3466
3467 /* Catch underflow that we deferred until after rounding. */
3468 if (REAL_EXP (r) <= emin2m1)
3469 goto underflow;
3470
3471 /* Clear out trailing garbage. */
3472 clear_significand_below (r, np2);
3473 }
3474
3475 /* Extend or truncate to a new mode. */
3476
3477 void
3478 real_convert (REAL_VALUE_TYPE *r, enum machine_mode mode,
3479 const REAL_VALUE_TYPE *a)
3480 {
3481 const struct real_format *fmt;
3482
3483 fmt = REAL_MODE_FORMAT (mode);
3484 gcc_assert (fmt);
3485
3486 *r = *a;
3487 round_for_format (fmt, r);
3488
3489 /* round_for_format de-normalizes denormals. Undo just that part. */
3490 if (r->cl == rvc_normal)
3491 normalize (r);
3492 }
3493
3494 /* Legacy. Likewise, except return the struct directly. */
3495
3496 REAL_VALUE_TYPE
3497 real_value_truncate (enum machine_mode mode, REAL_VALUE_TYPE a)
3498 {
3499 REAL_VALUE_TYPE r;
3500 real_convert (&r, mode, &a);
3501 return r;
3502 }
3503
3504 /* Return true if truncating to MODE is exact. */
3505
3506 bool
3507 exact_real_truncate (enum machine_mode mode, const REAL_VALUE_TYPE *a)
3508 {
3509 const struct real_format *fmt;
3510 REAL_VALUE_TYPE t;
3511 int emin2m1;
3512
3513 fmt = REAL_MODE_FORMAT (mode);
3514 gcc_assert (fmt);
3515
3516 /* Don't allow conversion to denormals. */
3517 emin2m1 = (fmt->emin - 1) * fmt->log2_b;
3518 if (REAL_EXP (a) <= emin2m1)
3519 return false;
3520
3521 /* After conversion to the new mode, the value must be identical. */
3522 real_convert (&t, mode, a);
3523 return real_identical (&t, a);
3524 }
3525
3526 /* Write R to the given target format. Place the words of the result
3527 in target word order in BUF. There are always 32 bits in each
3528 long, no matter the size of the host long.
3529
3530 Legacy: return word 0 for implementing REAL_VALUE_TO_TARGET_SINGLE. */
3531
3532 long
3533 real_to_target_fmt (long *buf, const REAL_VALUE_TYPE *r_orig,
3534 const struct real_format *fmt)
3535 {
3536 REAL_VALUE_TYPE r;
3537 long buf1;
3538
3539 r = *r_orig;
3540 round_for_format (fmt, &r);
3541
3542 if (!buf)
3543 buf = &buf1;
3544 (*fmt->encode) (fmt, buf, &r);
3545
3546 return *buf;
3547 }
3548
3549 /* Similar, but look up the format from MODE. */
3550
3551 long
3552 real_to_target (long *buf, const REAL_VALUE_TYPE *r, enum machine_mode mode)
3553 {
3554 const struct real_format *fmt;
3555
3556 fmt = REAL_MODE_FORMAT (mode);
3557 gcc_assert (fmt);
3558
3559 return real_to_target_fmt (buf, r, fmt);
3560 }
3561
3562 /* Read R from the given target format. Read the words of the result
3563 in target word order in BUF. There are always 32 bits in each
3564 long, no matter the size of the host long. */
3565
3566 void
3567 real_from_target_fmt (REAL_VALUE_TYPE *r, const long *buf,
3568 const struct real_format *fmt)
3569 {
3570 (*fmt->decode) (fmt, r, buf);
3571 }
3572
3573 /* Similar, but look up the format from MODE. */
3574
3575 void
3576 real_from_target (REAL_VALUE_TYPE *r, const long *buf, enum machine_mode mode)
3577 {
3578 const struct real_format *fmt;
3579
3580 fmt = REAL_MODE_FORMAT (mode);
3581 gcc_assert (fmt);
3582
3583 (*fmt->decode) (fmt, r, buf);
3584 }
3585
3586 /* Return the number of bits in the significand for MODE. */
3587 /* ??? Legacy. Should get access to real_format directly. */
3588
3589 int
3590 significand_size (enum machine_mode mode)
3591 {
3592 const struct real_format *fmt;
3593
3594 fmt = REAL_MODE_FORMAT (mode);
3595 if (fmt == NULL)
3596 return 0;
3597
3598 return fmt->p * fmt->log2_b;
3599 }
3600 #endif /* 0 */
3601
3602 /* Return a hash value for the given real value. */
3603 /* ??? The "unsigned int" return value is intended to be hashval_t,
3604 but I didn't want to pull hashtab.h into real.h. */
3605
3606 unsigned int
real_hash(const REAL_VALUE_TYPE * r)3607 real_hash (const REAL_VALUE_TYPE *r)
3608 {
3609 /* (TIGCC) Very naive hash for lack of something better. -- Kevin Kofler */
3610 unsigned int h=r->exponent;
3611 return h+(unsigned int)(r->mantissa >>
3612 ((sizeof(unsigned long long)-sizeof(unsigned int))*8));
3613
3614 #if 0
3615 unsigned int h;
3616 size_t i;
3617
3618 h = r->cl | (r->sign << 2);
3619 switch (r->cl)
3620 {
3621 case rvc_zero:
3622 case rvc_inf:
3623 return h;
3624
3625 case rvc_normal:
3626 h |= REAL_EXP (r) << 3;
3627 break;
3628
3629 case rvc_nan:
3630 if (r->signalling)
3631 h ^= (unsigned int)-1;
3632 if (r->canonical)
3633 return h;
3634 break;
3635
3636 default:
3637 gcc_unreachable ();
3638 }
3639
3640 if (sizeof(unsigned long) > sizeof(unsigned int))
3641 for (i = 0; i < SIGSZ; ++i)
3642 {
3643 unsigned long s = r->sig[i];
3644 h ^= s ^ (s >> (HOST_BITS_PER_LONG / 2));
3645 }
3646 else
3647 for (i = 0; i < SIGSZ; ++i)
3648 h ^= r->sig[i];
3649
3650 return h;
3651 #endif /* 0 */
3652 }
3653
3654 #if 0
3655 /* IEEE single-precision format. */
3656
3657 static void encode_ieee_single (const struct real_format *fmt,
3658 long *, const REAL_VALUE_TYPE *);
3659 static void decode_ieee_single (const struct real_format *,
3660 REAL_VALUE_TYPE *, const long *);
3661
3662 static void
3663 encode_ieee_single (const struct real_format *fmt, long *buf,
3664 const REAL_VALUE_TYPE *r)
3665 {
3666 unsigned long image, sig, exp;
3667 unsigned long sign = r->sign;
3668 bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
3669
3670 image = sign << 31;
3671 sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 24)) & 0x7fffff;
3672
3673 switch (r->cl)
3674 {
3675 case rvc_zero:
3676 break;
3677
3678 case rvc_inf:
3679 if (fmt->has_inf)
3680 image |= 255 << 23;
3681 else
3682 image |= 0x7fffffff;
3683 break;
3684
3685 case rvc_nan:
3686 if (fmt->has_nans)
3687 {
3688 if (r->canonical)
3689 sig = 0;
3690 if (r->signalling == fmt->qnan_msb_set)
3691 sig &= ~(1 << 22);
3692 else
3693 sig |= 1 << 22;
3694 /* We overload qnan_msb_set here: it's only clear for
3695 mips_ieee_single, which wants all mantissa bits but the
3696 quiet/signalling one set in canonical NaNs (at least
3697 Quiet ones). */
3698 if (r->canonical && !fmt->qnan_msb_set)
3699 sig |= (1 << 22) - 1;
3700 else if (sig == 0)
3701 sig = 1 << 21;
3702
3703 image |= 255 << 23;
3704 image |= sig;
3705 }
3706 else
3707 image |= 0x7fffffff;
3708 break;
3709
3710 case rvc_normal:
3711 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
3712 whereas the intermediate representation is 0.F x 2**exp.
3713 Which means we're off by one. */
3714 if (denormal)
3715 exp = 0;
3716 else
3717 exp = REAL_EXP (r) + 127 - 1;
3718 image |= exp << 23;
3719 image |= sig;
3720 break;
3721
3722 default:
3723 gcc_unreachable ();
3724 }
3725
3726 buf[0] = image;
3727 }
3728
3729 static void
3730 decode_ieee_single (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3731 const long *buf)
3732 {
3733 unsigned long image = buf[0] & 0xffffffff;
3734 bool sign = (image >> 31) & 1;
3735 int exp = (image >> 23) & 0xff;
3736
3737 memset (r, 0, sizeof (*r));
3738 image <<= HOST_BITS_PER_LONG - 24;
3739 image &= ~SIG_MSB;
3740
3741 if (exp == 0)
3742 {
3743 if (image && fmt->has_denorm)
3744 {
3745 r->cl = rvc_normal;
3746 r->sign = sign;
3747 SET_REAL_EXP (r, -126);
3748 r->sig[SIGSZ-1] = image << 1;
3749 normalize (r);
3750 }
3751 else if (fmt->has_signed_zero)
3752 r->sign = sign;
3753 }
3754 else if (exp == 255 && (fmt->has_nans || fmt->has_inf))
3755 {
3756 if (image)
3757 {
3758 r->cl = rvc_nan;
3759 r->sign = sign;
3760 r->signalling = (((image >> (HOST_BITS_PER_LONG - 2)) & 1)
3761 ^ fmt->qnan_msb_set);
3762 r->sig[SIGSZ-1] = image;
3763 }
3764 else
3765 {
3766 r->cl = rvc_inf;
3767 r->sign = sign;
3768 }
3769 }
3770 else
3771 {
3772 r->cl = rvc_normal;
3773 r->sign = sign;
3774 SET_REAL_EXP (r, exp - 127 + 1);
3775 r->sig[SIGSZ-1] = image | SIG_MSB;
3776 }
3777 }
3778
3779 const struct real_format ieee_single_format =
3780 {
3781 encode_ieee_single,
3782 decode_ieee_single,
3783 2,
3784 1,
3785 24,
3786 24,
3787 -125,
3788 128,
3789 31,
3790 31,
3791 true,
3792 true,
3793 true,
3794 true,
3795 true
3796 };
3797
3798 const struct real_format mips_single_format =
3799 {
3800 encode_ieee_single,
3801 decode_ieee_single,
3802 2,
3803 1,
3804 24,
3805 24,
3806 -125,
3807 128,
3808 31,
3809 31,
3810 true,
3811 true,
3812 true,
3813 true,
3814 false
3815 };
3816
3817
3818 /* IEEE double-precision format. */
3819
3820 static void encode_ieee_double (const struct real_format *fmt,
3821 long *, const REAL_VALUE_TYPE *);
3822 static void decode_ieee_double (const struct real_format *,
3823 REAL_VALUE_TYPE *, const long *);
3824
3825 static void
3826 encode_ieee_double (const struct real_format *fmt, long *buf,
3827 const REAL_VALUE_TYPE *r)
3828 {
3829 unsigned long image_lo, image_hi, sig_lo, sig_hi, exp;
3830 bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
3831
3832 image_hi = r->sign << 31;
3833 image_lo = 0;
3834
3835 if (HOST_BITS_PER_LONG == 64)
3836 {
3837 sig_hi = r->sig[SIGSZ-1];
3838 sig_lo = (sig_hi >> (64 - 53)) & 0xffffffff;
3839 sig_hi = (sig_hi >> (64 - 53 + 1) >> 31) & 0xfffff;
3840 }
3841 else
3842 {
3843 sig_hi = r->sig[SIGSZ-1];
3844 sig_lo = r->sig[SIGSZ-2];
3845 sig_lo = (sig_hi << 21) | (sig_lo >> 11);
3846 sig_hi = (sig_hi >> 11) & 0xfffff;
3847 }
3848
3849 switch (r->cl)
3850 {
3851 case rvc_zero:
3852 break;
3853
3854 case rvc_inf:
3855 if (fmt->has_inf)
3856 image_hi |= 2047 << 20;
3857 else
3858 {
3859 image_hi |= 0x7fffffff;
3860 image_lo = 0xffffffff;
3861 }
3862 break;
3863
3864 case rvc_nan:
3865 if (fmt->has_nans)
3866 {
3867 if (r->canonical)
3868 sig_hi = sig_lo = 0;
3869 if (r->signalling == fmt->qnan_msb_set)
3870 sig_hi &= ~(1 << 19);
3871 else
3872 sig_hi |= 1 << 19;
3873 /* We overload qnan_msb_set here: it's only clear for
3874 mips_ieee_single, which wants all mantissa bits but the
3875 quiet/signalling one set in canonical NaNs (at least
3876 Quiet ones). */
3877 if (r->canonical && !fmt->qnan_msb_set)
3878 {
3879 sig_hi |= (1 << 19) - 1;
3880 sig_lo = 0xffffffff;
3881 }
3882 else if (sig_hi == 0 && sig_lo == 0)
3883 sig_hi = 1 << 18;
3884
3885 image_hi |= 2047 << 20;
3886 image_hi |= sig_hi;
3887 image_lo = sig_lo;
3888 }
3889 else
3890 {
3891 image_hi |= 0x7fffffff;
3892 image_lo = 0xffffffff;
3893 }
3894 break;
3895
3896 case rvc_normal:
3897 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
3898 whereas the intermediate representation is 0.F x 2**exp.
3899 Which means we're off by one. */
3900 if (denormal)
3901 exp = 0;
3902 else
3903 exp = REAL_EXP (r) + 1023 - 1;
3904 image_hi |= exp << 20;
3905 image_hi |= sig_hi;
3906 image_lo = sig_lo;
3907 break;
3908
3909 default:
3910 gcc_unreachable ();
3911 }
3912
3913 if (FLOAT_WORDS_BIG_ENDIAN)
3914 buf[0] = image_hi, buf[1] = image_lo;
3915 else
3916 buf[0] = image_lo, buf[1] = image_hi;
3917 }
3918
3919 static void
3920 decode_ieee_double (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3921 const long *buf)
3922 {
3923 unsigned long image_hi, image_lo;
3924 bool sign;
3925 int exp;
3926
3927 if (FLOAT_WORDS_BIG_ENDIAN)
3928 image_hi = buf[0], image_lo = buf[1];
3929 else
3930 image_lo = buf[0], image_hi = buf[1];
3931 image_lo &= 0xffffffff;
3932 image_hi &= 0xffffffff;
3933
3934 sign = (image_hi >> 31) & 1;
3935 exp = (image_hi >> 20) & 0x7ff;
3936
3937 memset (r, 0, sizeof (*r));
3938
3939 image_hi <<= 32 - 21;
3940 image_hi |= image_lo >> 21;
3941 image_hi &= 0x7fffffff;
3942 image_lo <<= 32 - 21;
3943
3944 if (exp == 0)
3945 {
3946 if ((image_hi || image_lo) && fmt->has_denorm)
3947 {
3948 r->cl = rvc_normal;
3949 r->sign = sign;
3950 SET_REAL_EXP (r, -1022);
3951 if (HOST_BITS_PER_LONG == 32)
3952 {
3953 image_hi = (image_hi << 1) | (image_lo >> 31);
3954 image_lo <<= 1;
3955 r->sig[SIGSZ-1] = image_hi;
3956 r->sig[SIGSZ-2] = image_lo;
3957 }
3958 else
3959 {
3960 image_hi = (image_hi << 31 << 2) | (image_lo << 1);
3961 r->sig[SIGSZ-1] = image_hi;
3962 }
3963 normalize (r);
3964 }
3965 else if (fmt->has_signed_zero)
3966 r->sign = sign;
3967 }
3968 else if (exp == 2047 && (fmt->has_nans || fmt->has_inf))
3969 {
3970 if (image_hi || image_lo)
3971 {
3972 r->cl = rvc_nan;
3973 r->sign = sign;
3974 r->signalling = ((image_hi >> 30) & 1) ^ fmt->qnan_msb_set;
3975 if (HOST_BITS_PER_LONG == 32)
3976 {
3977 r->sig[SIGSZ-1] = image_hi;
3978 r->sig[SIGSZ-2] = image_lo;
3979 }
3980 else
3981 r->sig[SIGSZ-1] = (image_hi << 31 << 1) | image_lo;
3982 }
3983 else
3984 {
3985 r->cl = rvc_inf;
3986 r->sign = sign;
3987 }
3988 }
3989 else
3990 {
3991 r->cl = rvc_normal;
3992 r->sign = sign;
3993 SET_REAL_EXP (r, exp - 1023 + 1);
3994 if (HOST_BITS_PER_LONG == 32)
3995 {
3996 r->sig[SIGSZ-1] = image_hi | SIG_MSB;
3997 r->sig[SIGSZ-2] = image_lo;
3998 }
3999 else
4000 r->sig[SIGSZ-1] = (image_hi << 31 << 1) | image_lo | SIG_MSB;
4001 }
4002 }
4003
4004 const struct real_format ieee_double_format =
4005 {
4006 encode_ieee_double,
4007 decode_ieee_double,
4008 2,
4009 1,
4010 53,
4011 53,
4012 -1021,
4013 1024,
4014 63,
4015 63,
4016 true,
4017 true,
4018 true,
4019 true,
4020 true
4021 };
4022
4023 const struct real_format mips_double_format =
4024 {
4025 encode_ieee_double,
4026 decode_ieee_double,
4027 2,
4028 1,
4029 53,
4030 53,
4031 -1021,
4032 1024,
4033 63,
4034 63,
4035 true,
4036 true,
4037 true,
4038 true,
4039 false
4040 };
4041
4042
4043 /* IEEE extended real format. This comes in three flavors: Intel's as
4044 a 12 byte image, Intel's as a 16 byte image, and Motorola's. Intel
4045 12- and 16-byte images may be big- or little endian; Motorola's is
4046 always big endian. */
4047
4048 /* Helper subroutine which converts from the internal format to the
4049 12-byte little-endian Intel format. Functions below adjust this
4050 for the other possible formats. */
4051 static void
4052 encode_ieee_extended (const struct real_format *fmt, long *buf,
4053 const REAL_VALUE_TYPE *r)
4054 {
4055 unsigned long image_hi, sig_hi, sig_lo;
4056 bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
4057
4058 image_hi = r->sign << 15;
4059 sig_hi = sig_lo = 0;
4060
4061 switch (r->cl)
4062 {
4063 case rvc_zero:
4064 break;
4065
4066 case rvc_inf:
4067 if (fmt->has_inf)
4068 {
4069 image_hi |= 32767;
4070
4071 /* Intel requires the explicit integer bit to be set, otherwise
4072 it considers the value a "pseudo-infinity". Motorola docs
4073 say it doesn't care. */
4074 sig_hi = 0x80000000;
4075 }
4076 else
4077 {
4078 image_hi |= 32767;
4079 sig_lo = sig_hi = 0xffffffff;
4080 }
4081 break;
4082
4083 case rvc_nan:
4084 if (fmt->has_nans)
4085 {
4086 image_hi |= 32767;
4087 if (HOST_BITS_PER_LONG == 32)
4088 {
4089 sig_hi = r->sig[SIGSZ-1];
4090 sig_lo = r->sig[SIGSZ-2];
4091 }
4092 else
4093 {
4094 sig_lo = r->sig[SIGSZ-1];
4095 sig_hi = sig_lo >> 31 >> 1;
4096 sig_lo &= 0xffffffff;
4097 }
4098 if (r->signalling == fmt->qnan_msb_set)
4099 sig_hi &= ~(1 << 30);
4100 else
4101 sig_hi |= 1 << 30;
4102 if ((sig_hi & 0x7fffffff) == 0 && sig_lo == 0)
4103 sig_hi = 1 << 29;
4104
4105 /* Intel requires the explicit integer bit to be set, otherwise
4106 it considers the value a "pseudo-nan". Motorola docs say it
4107 doesn't care. */
4108 sig_hi |= 0x80000000;
4109 }
4110 else
4111 {
4112 image_hi |= 32767;
4113 sig_lo = sig_hi = 0xffffffff;
4114 }
4115 break;
4116
4117 case rvc_normal:
4118 {
4119 int exp = REAL_EXP (r);
4120
4121 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
4122 whereas the intermediate representation is 0.F x 2**exp.
4123 Which means we're off by one.
4124
4125 Except for Motorola, which consider exp=0 and explicit
4126 integer bit set to continue to be normalized. In theory
4127 this discrepancy has been taken care of by the difference
4128 in fmt->emin in round_for_format. */
4129
4130 if (denormal)
4131 exp = 0;
4132 else
4133 {
4134 exp += 16383 - 1;
4135 gcc_assert (exp >= 0);
4136 }
4137 image_hi |= exp;
4138
4139 if (HOST_BITS_PER_LONG == 32)
4140 {
4141 sig_hi = r->sig[SIGSZ-1];
4142 sig_lo = r->sig[SIGSZ-2];
4143 }
4144 else
4145 {
4146 sig_lo = r->sig[SIGSZ-1];
4147 sig_hi = sig_lo >> 31 >> 1;
4148 sig_lo &= 0xffffffff;
4149 }
4150 }
4151 break;
4152
4153 default:
4154 gcc_unreachable ();
4155 }
4156
4157 buf[0] = sig_lo, buf[1] = sig_hi, buf[2] = image_hi;
4158 }
4159
4160 /* Convert from the internal format to the 12-byte Motorola format
4161 for an IEEE extended real. */
4162 static void
4163 encode_ieee_extended_motorola (const struct real_format *fmt, long *buf,
4164 const REAL_VALUE_TYPE *r)
4165 {
4166 long intermed[3];
4167 encode_ieee_extended (fmt, intermed, r);
4168
4169 /* Motorola chips are assumed always to be big-endian. Also, the
4170 padding in a Motorola extended real goes between the exponent and
4171 the mantissa. At this point the mantissa is entirely within
4172 elements 0 and 1 of intermed, and the exponent entirely within
4173 element 2, so all we have to do is swap the order around, and
4174 shift element 2 left 16 bits. */
4175 buf[0] = intermed[2] << 16;
4176 buf[1] = intermed[1];
4177 buf[2] = intermed[0];
4178 }
4179
4180 /* Convert from the internal format to the 12-byte Intel format for
4181 an IEEE extended real. */
4182 static void
4183 encode_ieee_extended_intel_96 (const struct real_format *fmt, long *buf,
4184 const REAL_VALUE_TYPE *r)
4185 {
4186 if (FLOAT_WORDS_BIG_ENDIAN)
4187 {
4188 /* All the padding in an Intel-format extended real goes at the high
4189 end, which in this case is after the mantissa, not the exponent.
4190 Therefore we must shift everything down 16 bits. */
4191 long intermed[3];
4192 encode_ieee_extended (fmt, intermed, r);
4193 buf[0] = ((intermed[2] << 16) | ((unsigned long)(intermed[1] & 0xFFFF0000) >> 16));
4194 buf[1] = ((intermed[1] << 16) | ((unsigned long)(intermed[0] & 0xFFFF0000) >> 16));
4195 buf[2] = (intermed[0] << 16);
4196 }
4197 else
4198 /* encode_ieee_extended produces what we want directly. */
4199 encode_ieee_extended (fmt, buf, r);
4200 }
4201
4202 /* Convert from the internal format to the 16-byte Intel format for
4203 an IEEE extended real. */
4204 static void
4205 encode_ieee_extended_intel_128 (const struct real_format *fmt, long *buf,
4206 const REAL_VALUE_TYPE *r)
4207 {
4208 /* All the padding in an Intel-format extended real goes at the high end. */
4209 encode_ieee_extended_intel_96 (fmt, buf, r);
4210 buf[3] = 0;
4211 }
4212
4213 /* As above, we have a helper function which converts from 12-byte
4214 little-endian Intel format to internal format. Functions below
4215 adjust for the other possible formats. */
4216 static void
4217 decode_ieee_extended (const struct real_format *fmt, REAL_VALUE_TYPE *r,
4218 const long *buf)
4219 {
4220 unsigned long image_hi, sig_hi, sig_lo;
4221 bool sign;
4222 int exp;
4223
4224 sig_lo = buf[0], sig_hi = buf[1], image_hi = buf[2];
4225 sig_lo &= 0xffffffff;
4226 sig_hi &= 0xffffffff;
4227 image_hi &= 0xffffffff;
4228
4229 sign = (image_hi >> 15) & 1;
4230 exp = image_hi & 0x7fff;
4231
4232 memset (r, 0, sizeof (*r));
4233
4234 if (exp == 0)
4235 {
4236 if ((sig_hi || sig_lo) && fmt->has_denorm)
4237 {
4238 r->cl = rvc_normal;
4239 r->sign = sign;
4240
4241 /* When the IEEE format contains a hidden bit, we know that
4242 it's zero at this point, and so shift up the significand
4243 and decrease the exponent to match. In this case, Motorola
4244 defines the explicit integer bit to be valid, so we don't
4245 know whether the msb is set or not. */
4246 SET_REAL_EXP (r, fmt->emin);
4247 if (HOST_BITS_PER_LONG == 32)
4248 {
4249 r->sig[SIGSZ-1] = sig_hi;
4250 r->sig[SIGSZ-2] = sig_lo;
4251 }
4252 else
4253 r->sig[SIGSZ-1] = (sig_hi << 31 << 1) | sig_lo;
4254
4255 normalize (r);
4256 }
4257 else if (fmt->has_signed_zero)
4258 r->sign = sign;
4259 }
4260 else if (exp == 32767 && (fmt->has_nans || fmt->has_inf))
4261 {
4262 /* See above re "pseudo-infinities" and "pseudo-nans".
4263 Short summary is that the MSB will likely always be
4264 set, and that we don't care about it. */
4265 sig_hi &= 0x7fffffff;
4266
4267 if (sig_hi || sig_lo)
4268 {
4269 r->cl = rvc_nan;
4270 r->sign = sign;
4271 r->signalling = ((sig_hi >> 30) & 1) ^ fmt->qnan_msb_set;
4272 if (HOST_BITS_PER_LONG == 32)
4273 {
4274 r->sig[SIGSZ-1] = sig_hi;
4275 r->sig[SIGSZ-2] = sig_lo;
4276 }
4277 else
4278 r->sig[SIGSZ-1] = (sig_hi << 31 << 1) | sig_lo;
4279 }
4280 else
4281 {
4282 r->cl = rvc_inf;
4283 r->sign = sign;
4284 }
4285 }
4286 else
4287 {
4288 r->cl = rvc_normal;
4289 r->sign = sign;
4290 SET_REAL_EXP (r, exp - 16383 + 1);
4291 if (HOST_BITS_PER_LONG == 32)
4292 {
4293 r->sig[SIGSZ-1] = sig_hi;
4294 r->sig[SIGSZ-2] = sig_lo;
4295 }
4296 else
4297 r->sig[SIGSZ-1] = (sig_hi << 31 << 1) | sig_lo;
4298 }
4299 }
4300
4301 /* Convert from the internal format to the 12-byte Motorola format
4302 for an IEEE extended real. */
4303 static void
4304 decode_ieee_extended_motorola (const struct real_format *fmt, REAL_VALUE_TYPE *r,
4305 const long *buf)
4306 {
4307 long intermed[3];
4308
4309 /* Motorola chips are assumed always to be big-endian. Also, the
4310 padding in a Motorola extended real goes between the exponent and
4311 the mantissa; remove it. */
4312 intermed[0] = buf[2];
4313 intermed[1] = buf[1];
4314 intermed[2] = (unsigned long)buf[0] >> 16;
4315
4316 decode_ieee_extended (fmt, r, intermed);
4317 }
4318
4319 /* Convert from the internal format to the 12-byte Intel format for
4320 an IEEE extended real. */
4321 static void
4322 decode_ieee_extended_intel_96 (const struct real_format *fmt, REAL_VALUE_TYPE *r,
4323 const long *buf)
4324 {
4325 if (FLOAT_WORDS_BIG_ENDIAN)
4326 {
4327 /* All the padding in an Intel-format extended real goes at the high
4328 end, which in this case is after the mantissa, not the exponent.
4329 Therefore we must shift everything up 16 bits. */
4330 long intermed[3];
4331
4332 intermed[0] = (((unsigned long)buf[2] >> 16) | (buf[1] << 16));
4333 intermed[1] = (((unsigned long)buf[1] >> 16) | (buf[0] << 16));
4334 intermed[2] = ((unsigned long)buf[0] >> 16);
4335
4336 decode_ieee_extended (fmt, r, intermed);
4337 }
4338 else
4339 /* decode_ieee_extended produces what we want directly. */
4340 decode_ieee_extended (fmt, r, buf);
4341 }
4342
4343 /* Convert from the internal format to the 16-byte Intel format for
4344 an IEEE extended real. */
4345 static void
4346 decode_ieee_extended_intel_128 (const struct real_format *fmt, REAL_VALUE_TYPE *r,
4347 const long *buf)
4348 {
4349 /* All the padding in an Intel-format extended real goes at the high end. */
4350 decode_ieee_extended_intel_96 (fmt, r, buf);
4351 }
4352
4353 const struct real_format ieee_extended_motorola_format =
4354 {
4355 encode_ieee_extended_motorola,
4356 decode_ieee_extended_motorola,
4357 2,
4358 1,
4359 64,
4360 64,
4361 -16382,
4362 16384,
4363 95,
4364 95,
4365 true,
4366 true,
4367 true,
4368 true,
4369 true
4370 };
4371
4372 const struct real_format ieee_extended_intel_96_format =
4373 {
4374 encode_ieee_extended_intel_96,
4375 decode_ieee_extended_intel_96,
4376 2,
4377 1,
4378 64,
4379 64,
4380 -16381,
4381 16384,
4382 79,
4383 79,
4384 true,
4385 true,
4386 true,
4387 true,
4388 true
4389 };
4390
4391 const struct real_format ieee_extended_intel_128_format =
4392 {
4393 encode_ieee_extended_intel_128,
4394 decode_ieee_extended_intel_128,
4395 2,
4396 1,
4397 64,
4398 64,
4399 -16381,
4400 16384,
4401 79,
4402 79,
4403 true,
4404 true,
4405 true,
4406 true,
4407 true
4408 };
4409
4410 /* The following caters to i386 systems that set the rounding precision
4411 to 53 bits instead of 64, e.g. FreeBSD. */
4412 const struct real_format ieee_extended_intel_96_round_53_format =
4413 {
4414 encode_ieee_extended_intel_96,
4415 decode_ieee_extended_intel_96,
4416 2,
4417 1,
4418 53,
4419 53,
4420 -16381,
4421 16384,
4422 79,
4423 79,
4424 true,
4425 true,
4426 true,
4427 true,
4428 true
4429 };
4430
4431 /* IBM 128-bit extended precision format: a pair of IEEE double precision
4432 numbers whose sum is equal to the extended precision value. The number
4433 with greater magnitude is first. This format has the same magnitude
4434 range as an IEEE double precision value, but effectively 106 bits of
4435 significand precision. Infinity and NaN are represented by their IEEE
4436 double precision value stored in the first number, the second number is
4437 +0.0 or -0.0 for Infinity and don't-care for NaN. */
4438
4439 static void encode_ibm_extended (const struct real_format *fmt,
4440 long *, const REAL_VALUE_TYPE *);
4441 static void decode_ibm_extended (const struct real_format *,
4442 REAL_VALUE_TYPE *, const long *);
4443
4444 static void
4445 encode_ibm_extended (const struct real_format *fmt, long *buf,
4446 const REAL_VALUE_TYPE *r)
4447 {
4448 REAL_VALUE_TYPE u, normr, v;
4449 const struct real_format *base_fmt;
4450
4451 base_fmt = fmt->qnan_msb_set ? &ieee_double_format : &mips_double_format;
4452
4453 /* Renormlize R before doing any arithmetic on it. */
4454 normr = *r;
4455 if (normr.cl == rvc_normal)
4456 normalize (&normr);
4457
4458 /* u = IEEE double precision portion of significand. */
4459 u = normr;
4460 round_for_format (base_fmt, &u);
4461 encode_ieee_double (base_fmt, &buf[0], &u);
4462
4463 if (u.cl == rvc_normal)
4464 {
4465 do_add (&v, &normr, &u, 1);
4466 /* Call round_for_format since we might need to denormalize. */
4467 round_for_format (base_fmt, &v);
4468 encode_ieee_double (base_fmt, &buf[2], &v);
4469 }
4470 else
4471 {
4472 /* Inf, NaN, 0 are all representable as doubles, so the
4473 least-significant part can be 0.0. */
4474 buf[2] = 0;
4475 buf[3] = 0;
4476 }
4477 }
4478
4479 static void
4480 decode_ibm_extended (const struct real_format *fmt ATTRIBUTE_UNUSED, REAL_VALUE_TYPE *r,
4481 const long *buf)
4482 {
4483 REAL_VALUE_TYPE u, v;
4484 const struct real_format *base_fmt;
4485
4486 base_fmt = fmt->qnan_msb_set ? &ieee_double_format : &mips_double_format;
4487 decode_ieee_double (base_fmt, &u, &buf[0]);
4488
4489 if (u.cl != rvc_zero && u.cl != rvc_inf && u.cl != rvc_nan)
4490 {
4491 decode_ieee_double (base_fmt, &v, &buf[2]);
4492 do_add (r, &u, &v, 0);
4493 }
4494 else
4495 *r = u;
4496 }
4497
4498 const struct real_format ibm_extended_format =
4499 {
4500 encode_ibm_extended,
4501 decode_ibm_extended,
4502 2,
4503 1,
4504 53 + 53,
4505 53,
4506 -1021 + 53,
4507 1024,
4508 127,
4509 -1,
4510 true,
4511 true,
4512 true,
4513 true,
4514 true
4515 };
4516
4517 const struct real_format mips_extended_format =
4518 {
4519 encode_ibm_extended,
4520 decode_ibm_extended,
4521 2,
4522 1,
4523 53 + 53,
4524 53,
4525 -1021 + 53,
4526 1024,
4527 127,
4528 -1,
4529 true,
4530 true,
4531 true,
4532 true,
4533 false
4534 };
4535
4536
4537 /* IEEE quad precision format. */
4538
4539 static void encode_ieee_quad (const struct real_format *fmt,
4540 long *, const REAL_VALUE_TYPE *);
4541 static void decode_ieee_quad (const struct real_format *,
4542 REAL_VALUE_TYPE *, const long *);
4543
4544 static void
4545 encode_ieee_quad (const struct real_format *fmt, long *buf,
4546 const REAL_VALUE_TYPE *r)
4547 {
4548 unsigned long image3, image2, image1, image0, exp;
4549 bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
4550 REAL_VALUE_TYPE u;
4551
4552 image3 = r->sign << 31;
4553 image2 = 0;
4554 image1 = 0;
4555 image0 = 0;
4556
4557 rshift_significand (&u, r, SIGNIFICAND_BITS - 113);
4558
4559 switch (r->cl)
4560 {
4561 case rvc_zero:
4562 break;
4563
4564 case rvc_inf:
4565 if (fmt->has_inf)
4566 image3 |= 32767 << 16;
4567 else
4568 {
4569 image3 |= 0x7fffffff;
4570 image2 = 0xffffffff;
4571 image1 = 0xffffffff;
4572 image0 = 0xffffffff;
4573 }
4574 break;
4575
4576 case rvc_nan:
4577 if (fmt->has_nans)
4578 {
4579 image3 |= 32767 << 16;
4580
4581 if (r->canonical)
4582 {
4583 /* Don't use bits from the significand. The
4584 initialization above is right. */
4585 }
4586 else if (HOST_BITS_PER_LONG == 32)
4587 {
4588 image0 = u.sig[0];
4589 image1 = u.sig[1];
4590 image2 = u.sig[2];
4591 image3 |= u.sig[3] & 0xffff;
4592 }
4593 else
4594 {
4595 image0 = u.sig[0];
4596 image1 = image0 >> 31 >> 1;
4597 image2 = u.sig[1];
4598 image3 |= (image2 >> 31 >> 1) & 0xffff;
4599 image0 &= 0xffffffff;
4600 image2 &= 0xffffffff;
4601 }
4602 if (r->signalling == fmt->qnan_msb_set)
4603 image3 &= ~0x8000;
4604 else
4605 image3 |= 0x8000;
4606 /* We overload qnan_msb_set here: it's only clear for
4607 mips_ieee_single, which wants all mantissa bits but the
4608 quiet/signalling one set in canonical NaNs (at least
4609 Quiet ones). */
4610 if (r->canonical && !fmt->qnan_msb_set)
4611 {
4612 image3 |= 0x7fff;
4613 image2 = image1 = image0 = 0xffffffff;
4614 }
4615 else if (((image3 & 0xffff) | image2 | image1 | image0) == 0)
4616 image3 |= 0x4000;
4617 }
4618 else
4619 {
4620 image3 |= 0x7fffffff;
4621 image2 = 0xffffffff;
4622 image1 = 0xffffffff;
4623 image0 = 0xffffffff;
4624 }
4625 break;
4626
4627 case rvc_normal:
4628 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
4629 whereas the intermediate representation is 0.F x 2**exp.
4630 Which means we're off by one. */
4631 if (denormal)
4632 exp = 0;
4633 else
4634 exp = REAL_EXP (r) + 16383 - 1;
4635 image3 |= exp << 16;
4636
4637 if (HOST_BITS_PER_LONG == 32)
4638 {
4639 image0 = u.sig[0];
4640 image1 = u.sig[1];
4641 image2 = u.sig[2];
4642 image3 |= u.sig[3] & 0xffff;
4643 }
4644 else
4645 {
4646 image0 = u.sig[0];
4647 image1 = image0 >> 31 >> 1;
4648 image2 = u.sig[1];
4649 image3 |= (image2 >> 31 >> 1) & 0xffff;
4650 image0 &= 0xffffffff;
4651 image2 &= 0xffffffff;
4652 }
4653 break;
4654
4655 default:
4656 gcc_unreachable ();
4657 }
4658
4659 if (FLOAT_WORDS_BIG_ENDIAN)
4660 {
4661 buf[0] = image3;
4662 buf[1] = image2;
4663 buf[2] = image1;
4664 buf[3] = image0;
4665 }
4666 else
4667 {
4668 buf[0] = image0;
4669 buf[1] = image1;
4670 buf[2] = image2;
4671 buf[3] = image3;
4672 }
4673 }
4674
4675 static void
4676 decode_ieee_quad (const struct real_format *fmt, REAL_VALUE_TYPE *r,
4677 const long *buf)
4678 {
4679 unsigned long image3, image2, image1, image0;
4680 bool sign;
4681 int exp;
4682
4683 if (FLOAT_WORDS_BIG_ENDIAN)
4684 {
4685 image3 = buf[0];
4686 image2 = buf[1];
4687 image1 = buf[2];
4688 image0 = buf[3];
4689 }
4690 else
4691 {
4692 image0 = buf[0];
4693 image1 = buf[1];
4694 image2 = buf[2];
4695 image3 = buf[3];
4696 }
4697 image0 &= 0xffffffff;
4698 image1 &= 0xffffffff;
4699 image2 &= 0xffffffff;
4700
4701 sign = (image3 >> 31) & 1;
4702 exp = (image3 >> 16) & 0x7fff;
4703 image3 &= 0xffff;
4704
4705 memset (r, 0, sizeof (*r));
4706
4707 if (exp == 0)
4708 {
4709 if ((image3 | image2 | image1 | image0) && fmt->has_denorm)
4710 {
4711 r->cl = rvc_normal;
4712 r->sign = sign;
4713
4714 SET_REAL_EXP (r, -16382 + (SIGNIFICAND_BITS - 112));
4715 if (HOST_BITS_PER_LONG == 32)
4716 {
4717 r->sig[0] = image0;
4718 r->sig[1] = image1;
4719 r->sig[2] = image2;
4720 r->sig[3] = image3;
4721 }
4722 else
4723 {
4724 r->sig[0] = (image1 << 31 << 1) | image0;
4725 r->sig[1] = (image3 << 31 << 1) | image2;
4726 }
4727
4728 normalize (r);
4729 }
4730 else if (fmt->has_signed_zero)
4731 r->sign = sign;
4732 }
4733 else if (exp == 32767 && (fmt->has_nans || fmt->has_inf))
4734 {
4735 if (image3 | image2 | image1 | image0)
4736 {
4737 r->cl = rvc_nan;
4738 r->sign = sign;
4739 r->signalling = ((image3 >> 15) & 1) ^ fmt->qnan_msb_set;
4740
4741 if (HOST_BITS_PER_LONG == 32)
4742 {
4743 r->sig[0] = image0;
4744 r->sig[1] = image1;
4745 r->sig[2] = image2;
4746 r->sig[3] = image3;
4747 }
4748 else
4749 {
4750 r->sig[0] = (image1 << 31 << 1) | image0;
4751 r->sig[1] = (image3 << 31 << 1) | image2;
4752 }
4753 lshift_significand (r, r, SIGNIFICAND_BITS - 113);
4754 }
4755 else
4756 {
4757 r->cl = rvc_inf;
4758 r->sign = sign;
4759 }
4760 }
4761 else
4762 {
4763 r->cl = rvc_normal;
4764 r->sign = sign;
4765 SET_REAL_EXP (r, exp - 16383 + 1);
4766
4767 if (HOST_BITS_PER_LONG == 32)
4768 {
4769 r->sig[0] = image0;
4770 r->sig[1] = image1;
4771 r->sig[2] = image2;
4772 r->sig[3] = image3;
4773 }
4774 else
4775 {
4776 r->sig[0] = (image1 << 31 << 1) | image0;
4777 r->sig[1] = (image3 << 31 << 1) | image2;
4778 }
4779 lshift_significand (r, r, SIGNIFICAND_BITS - 113);
4780 r->sig[SIGSZ-1] |= SIG_MSB;
4781 }
4782 }
4783
4784 const struct real_format ieee_quad_format =
4785 {
4786 encode_ieee_quad,
4787 decode_ieee_quad,
4788 2,
4789 1,
4790 113,
4791 113,
4792 -16381,
4793 16384,
4794 127,
4795 127,
4796 true,
4797 true,
4798 true,
4799 true,
4800 true
4801 };
4802
4803 const struct real_format mips_quad_format =
4804 {
4805 encode_ieee_quad,
4806 decode_ieee_quad,
4807 2,
4808 1,
4809 113,
4810 113,
4811 -16381,
4812 16384,
4813 127,
4814 127,
4815 true,
4816 true,
4817 true,
4818 true,
4819 false
4820 };
4821
4822 /* Descriptions of VAX floating point formats can be found beginning at
4823
4824 http://h71000.www7.hp.com/doc/73FINAL/4515/4515pro_013.html#f_floating_point_format
4825
4826 The thing to remember is that they're almost IEEE, except for word
4827 order, exponent bias, and the lack of infinities, nans, and denormals.
4828
4829 We don't implement the H_floating format here, simply because neither
4830 the VAX or Alpha ports use it. */
4831
4832 static void encode_vax_f (const struct real_format *fmt,
4833 long *, const REAL_VALUE_TYPE *);
4834 static void decode_vax_f (const struct real_format *,
4835 REAL_VALUE_TYPE *, const long *);
4836 static void encode_vax_d (const struct real_format *fmt,
4837 long *, const REAL_VALUE_TYPE *);
4838 static void decode_vax_d (const struct real_format *,
4839 REAL_VALUE_TYPE *, const long *);
4840 static void encode_vax_g (const struct real_format *fmt,
4841 long *, const REAL_VALUE_TYPE *);
4842 static void decode_vax_g (const struct real_format *,
4843 REAL_VALUE_TYPE *, const long *);
4844
4845 static void
4846 encode_vax_f (const struct real_format *fmt ATTRIBUTE_UNUSED, long *buf,
4847 const REAL_VALUE_TYPE *r)
4848 {
4849 unsigned long sign, exp, sig, image;
4850
4851 sign = r->sign << 15;
4852
4853 switch (r->cl)
4854 {
4855 case rvc_zero:
4856 image = 0;
4857 break;
4858
4859 case rvc_inf:
4860 case rvc_nan:
4861 image = 0xffff7fff | sign;
4862 break;
4863
4864 case rvc_normal:
4865 sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 24)) & 0x7fffff;
4866 exp = REAL_EXP (r) + 128;
4867
4868 image = (sig << 16) & 0xffff0000;
4869 image |= sign;
4870 image |= exp << 7;
4871 image |= sig >> 16;
4872 break;
4873
4874 default:
4875 gcc_unreachable ();
4876 }
4877
4878 buf[0] = image;
4879 }
4880
4881 static void
4882 decode_vax_f (const struct real_format *fmt ATTRIBUTE_UNUSED,
4883 REAL_VALUE_TYPE *r, const long *buf)
4884 {
4885 unsigned long image = buf[0] & 0xffffffff;
4886 int exp = (image >> 7) & 0xff;
4887
4888 memset (r, 0, sizeof (*r));
4889
4890 if (exp != 0)
4891 {
4892 r->cl = rvc_normal;
4893 r->sign = (image >> 15) & 1;
4894 SET_REAL_EXP (r, exp - 128);
4895
4896 image = ((image & 0x7f) << 16) | ((image >> 16) & 0xffff);
4897 r->sig[SIGSZ-1] = (image << (HOST_BITS_PER_LONG - 24)) | SIG_MSB;
4898 }
4899 }
4900
4901 static void
4902 encode_vax_d (const struct real_format *fmt ATTRIBUTE_UNUSED, long *buf,
4903 const REAL_VALUE_TYPE *r)
4904 {
4905 unsigned long image0, image1, sign = r->sign << 15;
4906
4907 switch (r->cl)
4908 {
4909 case rvc_zero:
4910 image0 = image1 = 0;
4911 break;
4912
4913 case rvc_inf:
4914 case rvc_nan:
4915 image0 = 0xffff7fff | sign;
4916 image1 = 0xffffffff;
4917 break;
4918
4919 case rvc_normal:
4920 /* Extract the significand into straight hi:lo. */
4921 if (HOST_BITS_PER_LONG == 64)
4922 {
4923 image0 = r->sig[SIGSZ-1];
4924 image1 = (image0 >> (64 - 56)) & 0xffffffff;
4925 image0 = (image0 >> (64 - 56 + 1) >> 31) & 0x7fffff;
4926 }
4927 else
4928 {
4929 image0 = r->sig[SIGSZ-1];
4930 image1 = r->sig[SIGSZ-2];
4931 image1 = (image0 << 24) | (image1 >> 8);
4932 image0 = (image0 >> 8) & 0xffffff;
4933 }
4934
4935 /* Rearrange the half-words of the significand to match the
4936 external format. */
4937 image0 = ((image0 << 16) | (image0 >> 16)) & 0xffff007f;
4938 image1 = ((image1 << 16) | (image1 >> 16)) & 0xffffffff;
4939
4940 /* Add the sign and exponent. */
4941 image0 |= sign;
4942 image0 |= (REAL_EXP (r) + 128) << 7;
4943 break;
4944
4945 default:
4946 gcc_unreachable ();
4947 }
4948
4949 if (FLOAT_WORDS_BIG_ENDIAN)
4950 buf[0] = image1, buf[1] = image0;
4951 else
4952 buf[0] = image0, buf[1] = image1;
4953 }
4954
4955 static void
4956 decode_vax_d (const struct real_format *fmt ATTRIBUTE_UNUSED,
4957 REAL_VALUE_TYPE *r, const long *buf)
4958 {
4959 unsigned long image0, image1;
4960 int exp;
4961
4962 if (FLOAT_WORDS_BIG_ENDIAN)
4963 image1 = buf[0], image0 = buf[1];
4964 else
4965 image0 = buf[0], image1 = buf[1];
4966 image0 &= 0xffffffff;
4967 image1 &= 0xffffffff;
4968
4969 exp = (image0 >> 7) & 0xff;
4970
4971 memset (r, 0, sizeof (*r));
4972
4973 if (exp != 0)
4974 {
4975 r->cl = rvc_normal;
4976 r->sign = (image0 >> 15) & 1;
4977 SET_REAL_EXP (r, exp - 128);
4978
4979 /* Rearrange the half-words of the external format into
4980 proper ascending order. */
4981 image0 = ((image0 & 0x7f) << 16) | ((image0 >> 16) & 0xffff);
4982 image1 = ((image1 & 0xffff) << 16) | ((image1 >> 16) & 0xffff);
4983
4984 if (HOST_BITS_PER_LONG == 64)
4985 {
4986 image0 = (image0 << 31 << 1) | image1;
4987 image0 <<= 64 - 56;
4988 image0 |= SIG_MSB;
4989 r->sig[SIGSZ-1] = image0;
4990 }
4991 else
4992 {
4993 r->sig[SIGSZ-1] = image0;
4994 r->sig[SIGSZ-2] = image1;
4995 lshift_significand (r, r, 2*HOST_BITS_PER_LONG - 56);
4996 r->sig[SIGSZ-1] |= SIG_MSB;
4997 }
4998 }
4999 }
5000
5001 static void
5002 encode_vax_g (const struct real_format *fmt ATTRIBUTE_UNUSED, long *buf,
5003 const REAL_VALUE_TYPE *r)
5004 {
5005 unsigned long image0, image1, sign = r->sign << 15;
5006
5007 switch (r->cl)
5008 {
5009 case rvc_zero:
5010 image0 = image1 = 0;
5011 break;
5012
5013 case rvc_inf:
5014 case rvc_nan:
5015 image0 = 0xffff7fff | sign;
5016 image1 = 0xffffffff;
5017 break;
5018
5019 case rvc_normal:
5020 /* Extract the significand into straight hi:lo. */
5021 if (HOST_BITS_PER_LONG == 64)
5022 {
5023 image0 = r->sig[SIGSZ-1];
5024 image1 = (image0 >> (64 - 53)) & 0xffffffff;
5025 image0 = (image0 >> (64 - 53 + 1) >> 31) & 0xfffff;
5026 }
5027 else
5028 {
5029 image0 = r->sig[SIGSZ-1];
5030 image1 = r->sig[SIGSZ-2];
5031 image1 = (image0 << 21) | (image1 >> 11);
5032 image0 = (image0 >> 11) & 0xfffff;
5033 }
5034
5035 /* Rearrange the half-words of the significand to match the
5036 external format. */
5037 image0 = ((image0 << 16) | (image0 >> 16)) & 0xffff000f;
5038 image1 = ((image1 << 16) | (image1 >> 16)) & 0xffffffff;
5039
5040 /* Add the sign and exponent. */
5041 image0 |= sign;
5042 image0 |= (REAL_EXP (r) + 1024) << 4;
5043 break;
5044
5045 default:
5046 gcc_unreachable ();
5047 }
5048
5049 if (FLOAT_WORDS_BIG_ENDIAN)
5050 buf[0] = image1, buf[1] = image0;
5051 else
5052 buf[0] = image0, buf[1] = image1;
5053 }
5054
5055 static void
5056 decode_vax_g (const struct real_format *fmt ATTRIBUTE_UNUSED,
5057 REAL_VALUE_TYPE *r, const long *buf)
5058 {
5059 unsigned long image0, image1;
5060 int exp;
5061
5062 if (FLOAT_WORDS_BIG_ENDIAN)
5063 image1 = buf[0], image0 = buf[1];
5064 else
5065 image0 = buf[0], image1 = buf[1];
5066 image0 &= 0xffffffff;
5067 image1 &= 0xffffffff;
5068
5069 exp = (image0 >> 4) & 0x7ff;
5070
5071 memset (r, 0, sizeof (*r));
5072
5073 if (exp != 0)
5074 {
5075 r->cl = rvc_normal;
5076 r->sign = (image0 >> 15) & 1;
5077 SET_REAL_EXP (r, exp - 1024);
5078
5079 /* Rearrange the half-words of the external format into
5080 proper ascending order. */
5081 image0 = ((image0 & 0xf) << 16) | ((image0 >> 16) & 0xffff);
5082 image1 = ((image1 & 0xffff) << 16) | ((image1 >> 16) & 0xffff);
5083
5084 if (HOST_BITS_PER_LONG == 64)
5085 {
5086 image0 = (image0 << 31 << 1) | image1;
5087 image0 <<= 64 - 53;
5088 image0 |= SIG_MSB;
5089 r->sig[SIGSZ-1] = image0;
5090 }
5091 else
5092 {
5093 r->sig[SIGSZ-1] = image0;
5094 r->sig[SIGSZ-2] = image1;
5095 lshift_significand (r, r, 64 - 53);
5096 r->sig[SIGSZ-1] |= SIG_MSB;
5097 }
5098 }
5099 }
5100
5101 const struct real_format vax_f_format =
5102 {
5103 encode_vax_f,
5104 decode_vax_f,
5105 2,
5106 1,
5107 24,
5108 24,
5109 -127,
5110 127,
5111 15,
5112 15,
5113 false,
5114 false,
5115 false,
5116 false,
5117 false
5118 };
5119
5120 const struct real_format vax_d_format =
5121 {
5122 encode_vax_d,
5123 decode_vax_d,
5124 2,
5125 1,
5126 56,
5127 56,
5128 -127,
5129 127,
5130 15,
5131 15,
5132 false,
5133 false,
5134 false,
5135 false,
5136 false
5137 };
5138
5139 const struct real_format vax_g_format =
5140 {
5141 encode_vax_g,
5142 decode_vax_g,
5143 2,
5144 1,
5145 53,
5146 53,
5147 -1023,
5148 1023,
5149 15,
5150 15,
5151 false,
5152 false,
5153 false,
5154 false,
5155 false
5156 };
5157
5158 /* A good reference for these can be found in chapter 9 of
5159 "ESA/390 Principles of Operation", IBM document number SA22-7201-01.
5160 An on-line version can be found here:
5161
5162 http://publibz.boulder.ibm.com/cgi-bin/bookmgr_OS390/BOOKS/DZ9AR001/9.1?DT=19930923083613
5163 */
5164
5165 static void encode_i370_single (const struct real_format *fmt,
5166 long *, const REAL_VALUE_TYPE *);
5167 static void decode_i370_single (const struct real_format *,
5168 REAL_VALUE_TYPE *, const long *);
5169 static void encode_i370_double (const struct real_format *fmt,
5170 long *, const REAL_VALUE_TYPE *);
5171 static void decode_i370_double (const struct real_format *,
5172 REAL_VALUE_TYPE *, const long *);
5173
5174 static void
5175 encode_i370_single (const struct real_format *fmt ATTRIBUTE_UNUSED,
5176 long *buf, const REAL_VALUE_TYPE *r)
5177 {
5178 unsigned long sign, exp, sig, image;
5179
5180 sign = r->sign << 31;
5181
5182 switch (r->cl)
5183 {
5184 case rvc_zero:
5185 image = 0;
5186 break;
5187
5188 case rvc_inf:
5189 case rvc_nan:
5190 image = 0x7fffffff | sign;
5191 break;
5192
5193 case rvc_normal:
5194 sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 24)) & 0xffffff;
5195 exp = ((REAL_EXP (r) / 4) + 64) << 24;
5196 image = sign | exp | sig;
5197 break;
5198
5199 default:
5200 gcc_unreachable ();
5201 }
5202
5203 buf[0] = image;
5204 }
5205
5206 static void
5207 decode_i370_single (const struct real_format *fmt ATTRIBUTE_UNUSED,
5208 REAL_VALUE_TYPE *r, const long *buf)
5209 {
5210 unsigned long sign, sig, image = buf[0];
5211 int exp;
5212
5213 sign = (image >> 31) & 1;
5214 exp = (image >> 24) & 0x7f;
5215 sig = image & 0xffffff;
5216
5217 memset (r, 0, sizeof (*r));
5218
5219 if (exp || sig)
5220 {
5221 r->cl = rvc_normal;
5222 r->sign = sign;
5223 SET_REAL_EXP (r, (exp - 64) * 4);
5224 r->sig[SIGSZ-1] = sig << (HOST_BITS_PER_LONG - 24);
5225 normalize (r);
5226 }
5227 }
5228
5229 static void
5230 encode_i370_double (const struct real_format *fmt ATTRIBUTE_UNUSED,
5231 long *buf, const REAL_VALUE_TYPE *r)
5232 {
5233 unsigned long sign, exp, image_hi, image_lo;
5234
5235 sign = r->sign << 31;
5236
5237 switch (r->cl)
5238 {
5239 case rvc_zero:
5240 image_hi = image_lo = 0;
5241 break;
5242
5243 case rvc_inf:
5244 case rvc_nan:
5245 image_hi = 0x7fffffff | sign;
5246 image_lo = 0xffffffff;
5247 break;
5248
5249 case rvc_normal:
5250 if (HOST_BITS_PER_LONG == 64)
5251 {
5252 image_hi = r->sig[SIGSZ-1];
5253 image_lo = (image_hi >> (64 - 56)) & 0xffffffff;
5254 image_hi = (image_hi >> (64 - 56 + 1) >> 31) & 0xffffff;
5255 }
5256 else
5257 {
5258 image_hi = r->sig[SIGSZ-1];
5259 image_lo = r->sig[SIGSZ-2];
5260 image_lo = (image_lo >> 8) | (image_hi << 24);
5261 image_hi >>= 8;
5262 }
5263
5264 exp = ((REAL_EXP (r) / 4) + 64) << 24;
5265 image_hi |= sign | exp;
5266 break;
5267
5268 default:
5269 gcc_unreachable ();
5270 }
5271
5272 if (FLOAT_WORDS_BIG_ENDIAN)
5273 buf[0] = image_hi, buf[1] = image_lo;
5274 else
5275 buf[0] = image_lo, buf[1] = image_hi;
5276 }
5277
5278 static void
5279 decode_i370_double (const struct real_format *fmt ATTRIBUTE_UNUSED,
5280 REAL_VALUE_TYPE *r, const long *buf)
5281 {
5282 unsigned long sign, image_hi, image_lo;
5283 int exp;
5284
5285 if (FLOAT_WORDS_BIG_ENDIAN)
5286 image_hi = buf[0], image_lo = buf[1];
5287 else
5288 image_lo = buf[0], image_hi = buf[1];
5289
5290 sign = (image_hi >> 31) & 1;
5291 exp = (image_hi >> 24) & 0x7f;
5292 image_hi &= 0xffffff;
5293 image_lo &= 0xffffffff;
5294
5295 memset (r, 0, sizeof (*r));
5296
5297 if (exp || image_hi || image_lo)
5298 {
5299 r->cl = rvc_normal;
5300 r->sign = sign;
5301 SET_REAL_EXP (r, (exp - 64) * 4 + (SIGNIFICAND_BITS - 56));
5302
5303 if (HOST_BITS_PER_LONG == 32)
5304 {
5305 r->sig[0] = image_lo;
5306 r->sig[1] = image_hi;
5307 }
5308 else
5309 r->sig[0] = image_lo | (image_hi << 31 << 1);
5310
5311 normalize (r);
5312 }
5313 }
5314
5315 const struct real_format i370_single_format =
5316 {
5317 encode_i370_single,
5318 decode_i370_single,
5319 16,
5320 4,
5321 6,
5322 6,
5323 -64,
5324 63,
5325 31,
5326 31,
5327 false,
5328 false,
5329 false, /* ??? The encoding does allow for "unnormals". */
5330 false, /* ??? The encoding does allow for "unnormals". */
5331 false
5332 };
5333
5334 const struct real_format i370_double_format =
5335 {
5336 encode_i370_double,
5337 decode_i370_double,
5338 16,
5339 4,
5340 14,
5341 14,
5342 -64,
5343 63,
5344 63,
5345 63,
5346 false,
5347 false,
5348 false, /* ??? The encoding does allow for "unnormals". */
5349 false, /* ??? The encoding does allow for "unnormals". */
5350 false
5351 };
5352
5353 /* The "twos-complement" c4x format is officially defined as
5354
5355 x = s(~s).f * 2**e
5356
5357 This is rather misleading. One must remember that F is signed.
5358 A better description would be
5359
5360 x = -1**s * ((s + 1 + .f) * 2**e
5361
5362 So if we have a (4 bit) fraction of .1000 with a sign bit of 1,
5363 that's -1 * (1+1+(-.5)) == -1.5. I think.
5364
5365 The constructions here are taken from Tables 5-1 and 5-2 of the
5366 TMS320C4x User's Guide wherein step-by-step instructions for
5367 conversion from IEEE are presented. That's close enough to our
5368 internal representation so as to make things easy.
5369
5370 See http://www-s.ti.com/sc/psheets/spru063c/spru063c.pdf */
5371
5372 static void encode_c4x_single (const struct real_format *fmt,
5373 long *, const REAL_VALUE_TYPE *);
5374 static void decode_c4x_single (const struct real_format *,
5375 REAL_VALUE_TYPE *, const long *);
5376 static void encode_c4x_extended (const struct real_format *fmt,
5377 long *, const REAL_VALUE_TYPE *);
5378 static void decode_c4x_extended (const struct real_format *,
5379 REAL_VALUE_TYPE *, const long *);
5380
5381 static void
5382 encode_c4x_single (const struct real_format *fmt ATTRIBUTE_UNUSED,
5383 long *buf, const REAL_VALUE_TYPE *r)
5384 {
5385 unsigned long image, exp, sig;
5386
5387 switch (r->cl)
5388 {
5389 case rvc_zero:
5390 exp = -128;
5391 sig = 0;
5392 break;
5393
5394 case rvc_inf:
5395 case rvc_nan:
5396 exp = 127;
5397 sig = 0x800000 - r->sign;
5398 break;
5399
5400 case rvc_normal:
5401 exp = REAL_EXP (r) - 1;
5402 sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 24)) & 0x7fffff;
5403 if (r->sign)
5404 {
5405 if (sig)
5406 sig = -sig;
5407 else
5408 exp--;
5409 sig |= 0x800000;
5410 }
5411 break;
5412
5413 default:
5414 gcc_unreachable ();
5415 }
5416
5417 image = ((exp & 0xff) << 24) | (sig & 0xffffff);
5418 buf[0] = image;
5419 }
5420
5421 static void
5422 decode_c4x_single (const struct real_format *fmt ATTRIBUTE_UNUSED,
5423 REAL_VALUE_TYPE *r, const long *buf)
5424 {
5425 unsigned long image = buf[0];
5426 unsigned long sig;
5427 int exp, sf;
5428
5429 exp = (((image >> 24) & 0xff) ^ 0x80) - 0x80;
5430 sf = ((image & 0xffffff) ^ 0x800000) - 0x800000;
5431
5432 memset (r, 0, sizeof (*r));
5433
5434 if (exp != -128)
5435 {
5436 r->cl = rvc_normal;
5437
5438 sig = sf & 0x7fffff;
5439 if (sf < 0)
5440 {
5441 r->sign = 1;
5442 if (sig)
5443 sig = -sig;
5444 else
5445 exp++;
5446 }
5447 sig = (sig << (HOST_BITS_PER_LONG - 24)) | SIG_MSB;
5448
5449 SET_REAL_EXP (r, exp + 1);
5450 r->sig[SIGSZ-1] = sig;
5451 }
5452 }
5453
5454 static void
5455 encode_c4x_extended (const struct real_format *fmt ATTRIBUTE_UNUSED,
5456 long *buf, const REAL_VALUE_TYPE *r)
5457 {
5458 unsigned long exp, sig;
5459
5460 switch (r->cl)
5461 {
5462 case rvc_zero:
5463 exp = -128;
5464 sig = 0;
5465 break;
5466
5467 case rvc_inf:
5468 case rvc_nan:
5469 exp = 127;
5470 sig = 0x80000000 - r->sign;
5471 break;
5472
5473 case rvc_normal:
5474 exp = REAL_EXP (r) - 1;
5475
5476 sig = r->sig[SIGSZ-1];
5477 if (HOST_BITS_PER_LONG == 64)
5478 sig = sig >> 1 >> 31;
5479 sig &= 0x7fffffff;
5480
5481 if (r->sign)
5482 {
5483 if (sig)
5484 sig = -sig;
5485 else
5486 exp--;
5487 sig |= 0x80000000;
5488 }
5489 break;
5490
5491 default:
5492 gcc_unreachable ();
5493 }
5494
5495 exp = (exp & 0xff) << 24;
5496 sig &= 0xffffffff;
5497
5498 if (FLOAT_WORDS_BIG_ENDIAN)
5499 buf[0] = exp, buf[1] = sig;
5500 else
5501 buf[0] = sig, buf[0] = exp;
5502 }
5503
5504 static void
5505 decode_c4x_extended (const struct real_format *fmt ATTRIBUTE_UNUSED,
5506 REAL_VALUE_TYPE *r, const long *buf)
5507 {
5508 unsigned long sig;
5509 int exp, sf;
5510
5511 if (FLOAT_WORDS_BIG_ENDIAN)
5512 exp = buf[0], sf = buf[1];
5513 else
5514 sf = buf[0], exp = buf[1];
5515
5516 exp = (((exp >> 24) & 0xff) & 0x80) - 0x80;
5517 sf = ((sf & 0xffffffff) ^ 0x80000000) - 0x80000000;
5518
5519 memset (r, 0, sizeof (*r));
5520
5521 if (exp != -128)
5522 {
5523 r->cl = rvc_normal;
5524
5525 sig = sf & 0x7fffffff;
5526 if (sf < 0)
5527 {
5528 r->sign = 1;
5529 if (sig)
5530 sig = -sig;
5531 else
5532 exp++;
5533 }
5534 if (HOST_BITS_PER_LONG == 64)
5535 sig = sig << 1 << 31;
5536 sig |= SIG_MSB;
5537
5538 SET_REAL_EXP (r, exp + 1);
5539 r->sig[SIGSZ-1] = sig;
5540 }
5541 }
5542
5543 const struct real_format c4x_single_format =
5544 {
5545 encode_c4x_single,
5546 decode_c4x_single,
5547 2,
5548 1,
5549 24,
5550 24,
5551 -126,
5552 128,
5553 23,
5554 -1,
5555 false,
5556 false,
5557 false,
5558 false,
5559 false
5560 };
5561
5562 const struct real_format c4x_extended_format =
5563 {
5564 encode_c4x_extended,
5565 decode_c4x_extended,
5566 2,
5567 1,
5568 32,
5569 32,
5570 -126,
5571 128,
5572 31,
5573 -1,
5574 false,
5575 false,
5576 false,
5577 false,
5578 false
5579 };
5580
5581
5582 /* A synthetic "format" for internal arithmetic. It's the size of the
5583 internal significand minus the two bits needed for proper rounding.
5584 The encode and decode routines exist only to satisfy our paranoia
5585 harness. */
5586
5587 static void encode_internal (const struct real_format *fmt,
5588 long *, const REAL_VALUE_TYPE *);
5589 static void decode_internal (const struct real_format *,
5590 REAL_VALUE_TYPE *, const long *);
5591
5592 static void
5593 encode_internal (const struct real_format *fmt ATTRIBUTE_UNUSED, long *buf,
5594 const REAL_VALUE_TYPE *r)
5595 {
5596 memcpy (buf, r, sizeof (*r));
5597 }
5598
5599 static void
5600 decode_internal (const struct real_format *fmt ATTRIBUTE_UNUSED,
5601 REAL_VALUE_TYPE *r, const long *buf)
5602 {
5603 memcpy (r, buf, sizeof (*r));
5604 }
5605
5606 const struct real_format real_internal_format =
5607 {
5608 encode_internal,
5609 decode_internal,
5610 2,
5611 1,
5612 SIGNIFICAND_BITS - 2,
5613 SIGNIFICAND_BITS - 2,
5614 -MAX_EXP,
5615 MAX_EXP,
5616 -1,
5617 -1,
5618 true,
5619 true,
5620 false,
5621 true,
5622 true
5623 };
5624
5625 /* Calculate the square root of X in mode MODE, and store the result
5626 in R. Return TRUE if the operation does not raise an exception.
5627 For details see "High Precision Division and Square Root",
5628 Alan H. Karp and Peter Markstein, HP Lab Report 93-93-42, June
5629 1993. http://www.hpl.hp.com/techreports/93/HPL-93-42.pdf. */
5630
5631 bool
5632 real_sqrt (REAL_VALUE_TYPE *r, enum machine_mode mode,
5633 const REAL_VALUE_TYPE *x)
5634 {
5635 static REAL_VALUE_TYPE halfthree;
5636 static bool init = false;
5637 REAL_VALUE_TYPE h, t, i;
5638 int iter, exp;
5639
5640 /* sqrt(-0.0) is -0.0. */
5641 if (real_isnegzero (x))
5642 {
5643 *r = *x;
5644 return false;
5645 }
5646
5647 /* Negative arguments return NaN. */
5648 if (real_isneg (x))
5649 {
5650 get_canonical_qnan (r, 0);
5651 return false;
5652 }
5653
5654 /* Infinity and NaN return themselves. */
5655 if (real_isinf (x) || real_isnan (x))
5656 {
5657 *r = *x;
5658 return false;
5659 }
5660
5661 if (!init)
5662 {
5663 do_add (&halfthree, &dconst1, &dconsthalf, 0);
5664 init = true;
5665 }
5666
5667 /* Initial guess for reciprocal sqrt, i. */
5668 exp = real_exponent (x);
5669 real_ldexp (&i, &dconst1, -exp/2);
5670
5671 /* Newton's iteration for reciprocal sqrt, i. */
5672 for (iter = 0; iter < 16; iter++)
5673 {
5674 /* i(n+1) = i(n) * (1.5 - 0.5*i(n)*i(n)*x). */
5675 do_multiply (&t, x, &i);
5676 do_multiply (&h, &t, &i);
5677 do_multiply (&t, &h, &dconsthalf);
5678 do_add (&h, &halfthree, &t, 1);
5679 do_multiply (&t, &i, &h);
5680
5681 /* Check for early convergence. */
5682 if (iter >= 6 && real_identical (&i, &t))
5683 break;
5684
5685 /* ??? Unroll loop to avoid copying. */
5686 i = t;
5687 }
5688
5689 /* Final iteration: r = i*x + 0.5*i*x*(1.0 - i*(i*x)). */
5690 do_multiply (&t, x, &i);
5691 do_multiply (&h, &t, &i);
5692 do_add (&i, &dconst1, &h, 1);
5693 do_multiply (&h, &t, &i);
5694 do_multiply (&i, &dconsthalf, &h);
5695 do_add (&h, &t, &i, 0);
5696
5697 /* ??? We need a Tuckerman test to get the last bit. */
5698
5699 real_convert (r, mode, &h);
5700 return true;
5701 }
5702
5703 /* Calculate X raised to the integer exponent N in mode MODE and store
5704 the result in R. Return true if the result may be inexact due to
5705 loss of precision. The algorithm is the classic "left-to-right binary
5706 method" described in section 4.6.3 of Donald Knuth's "Seminumerical
5707 Algorithms", "The Art of Computer Programming", Volume 2. */
5708
5709 bool
5710 real_powi (REAL_VALUE_TYPE *r, enum machine_mode mode,
5711 const REAL_VALUE_TYPE *x, HOST_WIDE_INT n)
5712 {
5713 unsigned HOST_WIDE_INT bit;
5714 REAL_VALUE_TYPE t;
5715 bool inexact = false;
5716 bool init = false;
5717 bool neg;
5718 int i;
5719
5720 if (n == 0)
5721 {
5722 *r = dconst1;
5723 return false;
5724 }
5725 else if (n < 0)
5726 {
5727 /* Don't worry about overflow, from now on n is unsigned. */
5728 neg = true;
5729 n = -n;
5730 }
5731 else
5732 neg = false;
5733
5734 t = *x;
5735 bit = (unsigned HOST_WIDE_INT) 1 << (HOST_BITS_PER_WIDE_INT - 1);
5736 for (i = 0; i < HOST_BITS_PER_WIDE_INT; i++)
5737 {
5738 if (init)
5739 {
5740 inexact |= do_multiply (&t, &t, &t);
5741 if (n & bit)
5742 inexact |= do_multiply (&t, &t, x);
5743 }
5744 else if (n & bit)
5745 init = true;
5746 bit >>= 1;
5747 }
5748
5749 if (neg)
5750 inexact |= do_divide (&t, &dconst1, &t);
5751
5752 real_convert (r, mode, &t);
5753 return inexact;
5754 }
5755
5756 /* Round X to the nearest integer not larger in absolute value, i.e.
5757 towards zero, placing the result in R in mode MODE. */
5758
5759 void
5760 real_trunc (REAL_VALUE_TYPE *r, enum machine_mode mode,
5761 const REAL_VALUE_TYPE *x)
5762 {
5763 do_fix_trunc (r, x);
5764 if (mode != VOIDmode)
5765 real_convert (r, mode, r);
5766 }
5767
5768 /* Round X to the largest integer not greater in value, i.e. round
5769 down, placing the result in R in mode MODE. */
5770
5771 void
5772 real_floor (REAL_VALUE_TYPE *r, enum machine_mode mode,
5773 const REAL_VALUE_TYPE *x)
5774 {
5775 REAL_VALUE_TYPE t;
5776
5777 do_fix_trunc (&t, x);
5778 if (! real_identical (&t, x) && x->sign)
5779 do_add (&t, &t, &dconstm1, 0);
5780 if (mode != VOIDmode)
5781 real_convert (r, mode, &t);
5782 else
5783 *r = t;
5784 }
5785
5786 /* Round X to the smallest integer not less then argument, i.e. round
5787 up, placing the result in R in mode MODE. */
5788
5789 void
5790 real_ceil (REAL_VALUE_TYPE *r, enum machine_mode mode,
5791 const REAL_VALUE_TYPE *x)
5792 {
5793 REAL_VALUE_TYPE t;
5794
5795 do_fix_trunc (&t, x);
5796 if (! real_identical (&t, x) && ! x->sign)
5797 do_add (&t, &t, &dconst1, 0);
5798 if (mode != VOIDmode)
5799 real_convert (r, mode, &t);
5800 else
5801 *r = t;
5802 }
5803
5804 /* Round X to the nearest integer, but round halfway cases away from
5805 zero. */
5806
5807 void
5808 real_round (REAL_VALUE_TYPE *r, enum machine_mode mode,
5809 const REAL_VALUE_TYPE *x)
5810 {
5811 do_add (r, x, &dconsthalf, x->sign);
5812 do_fix_trunc (r, r);
5813 if (mode != VOIDmode)
5814 real_convert (r, mode, r);
5815 }
5816
5817 /* Set the sign of R to the sign of X. */
5818
5819 void
5820 real_copysign (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *x)
5821 {
5822 r->sign = x->sign;
5823 }
5824 #endif /* 0 */
5825
5826