1 | /* This is a software floating point library which can be used instead of |
---|
2 | the floating point routines in libgcc1.c for targets without hardware |
---|
3 | floating point. |
---|
4 | Copyright (C) 1994, 1995, 1996, 1997 Free Software Foundation, Inc. |
---|
5 | |
---|
6 | This file is free software; you can redistribute it and/or modify it |
---|
7 | under the terms of the GNU General Public License as published by the |
---|
8 | Free Software Foundation; either version 2, or (at your option) any |
---|
9 | later version. |
---|
10 | |
---|
11 | In addition to the permissions in the GNU General Public License, the |
---|
12 | Free Software Foundation gives you unlimited permission to link the |
---|
13 | compiled version of this file with other programs, and to distribute |
---|
14 | those programs without any restriction coming from the use of this |
---|
15 | file. (The General Public License restrictions do apply in other |
---|
16 | respects; for example, they cover modification of the file, and |
---|
17 | distribution when not linked into another program.) |
---|
18 | |
---|
19 | This file is distributed in the hope that it will be useful, but |
---|
20 | WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
21 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
---|
22 | General Public License for more details. |
---|
23 | |
---|
24 | You should have received a copy of the GNU General Public License |
---|
25 | along with this program; see the file COPYING. If not, write to |
---|
26 | the Free Software Foundation, 59 Temple Place - Suite 330, |
---|
27 | Boston, MA 02111-1307, USA. */ |
---|
28 | |
---|
29 | /* As a special exception, if you link this library with other files, |
---|
30 | some of which are compiled with GCC, to produce an executable, |
---|
31 | this library does not by itself cause the resulting executable |
---|
32 | to be covered by the GNU General Public License. |
---|
33 | This exception does not however invalidate any other reasons why |
---|
34 | the executable file might be covered by the GNU General Public License. */ |
---|
35 | |
---|
36 | /* This implements IEEE 754 format arithmetic, but does not provide a |
---|
37 | mechanism for setting the rounding mode, or for generating or handling |
---|
38 | exceptions. |
---|
39 | |
---|
40 | The original code by Steve Chamberlain, hacked by Mark Eichin and Jim |
---|
41 | Wilson, all of Cygnus Support. */ |
---|
42 | |
---|
43 | /* The intended way to use this file is to make two copies, add `#define FLOAT' |
---|
44 | to one copy, then compile both copies and add them to libgcc.a. */ |
---|
45 | |
---|
46 | /* The following macros can be defined to change the behaviour of this file: |
---|
47 | FLOAT: Implement a `float', aka SFmode, fp library. If this is not |
---|
48 | defined, then this file implements a `double', aka DFmode, fp library. |
---|
49 | FLOAT_ONLY: Used with FLOAT, to implement a `float' only library, i.e. |
---|
50 | don't include float->double conversion which requires the double library. |
---|
51 | This is useful only for machines which can't support doubles, e.g. some |
---|
52 | 8-bit processors. |
---|
53 | CMPtype: Specify the type that floating point compares should return. |
---|
54 | This defaults to SItype, aka int. |
---|
55 | US_SOFTWARE_GOFAST: This makes all entry points use the same names as the |
---|
56 | US Software goFast library. If this is not defined, the entry points use |
---|
57 | the same names as libgcc1.c. |
---|
58 | _DEBUG_BITFLOAT: This makes debugging the code a little easier, by adding |
---|
59 | two integers to the FLO_union_type. |
---|
60 | NO_NANS: Disable nan and infinity handling |
---|
61 | SMALL_MACHINE: Useful when operations on QIs and HIs are faster |
---|
62 | than on an SI */ |
---|
63 | |
---|
64 | /* We don't currently support extended floats (long doubles) on machines |
---|
65 | without hardware to deal with them. |
---|
66 | |
---|
67 | These stubs are just to keep the linker from complaining about unresolved |
---|
68 | references which can be pulled in from libio & libstdc++, even if the |
---|
69 | user isn't using long doubles. However, they may generate an unresolved |
---|
70 | external to abort if abort is not used by the function, and the stubs |
---|
71 | are referenced from within libc, since libgcc goes before and after the |
---|
72 | system library. */ |
---|
73 | |
---|
74 | #ifdef EXTENDED_FLOAT_STUBS |
---|
75 | __truncxfsf2 (){ abort(); } |
---|
76 | __extendsfxf2 (){ abort(); } |
---|
77 | __addxf3 (){ abort(); } |
---|
78 | __divxf3 (){ abort(); } |
---|
79 | __eqxf2 (){ abort(); } |
---|
80 | __extenddfxf2 (){ abort(); } |
---|
81 | __gtxf2 (){ abort(); } |
---|
82 | __lexf2 (){ abort(); } |
---|
83 | __ltxf2 (){ abort(); } |
---|
84 | __mulxf3 (){ abort(); } |
---|
85 | __negxf2 (){ abort(); } |
---|
86 | __nexf2 (){ abort(); } |
---|
87 | __subxf3 (){ abort(); } |
---|
88 | __truncxfdf2 (){ abort(); } |
---|
89 | |
---|
90 | __trunctfsf2 (){ abort(); } |
---|
91 | __extendsftf2 (){ abort(); } |
---|
92 | __addtf3 (){ abort(); } |
---|
93 | __divtf3 (){ abort(); } |
---|
94 | __eqtf2 (){ abort(); } |
---|
95 | __extenddftf2 (){ abort(); } |
---|
96 | __gttf2 (){ abort(); } |
---|
97 | __letf2 (){ abort(); } |
---|
98 | __lttf2 (){ abort(); } |
---|
99 | __multf3 (){ abort(); } |
---|
100 | __negtf2 (){ abort(); } |
---|
101 | __netf2 (){ abort(); } |
---|
102 | __subtf3 (){ abort(); } |
---|
103 | __trunctfdf2 (){ abort(); } |
---|
104 | #else /* !EXTENDED_FLOAT_STUBS, rest of file */ |
---|
105 | |
---|
106 | |
---|
107 | typedef SFtype __attribute__ ((mode (SF))); |
---|
108 | typedef DFtype __attribute__ ((mode (DF))); |
---|
109 | |
---|
110 | typedef int HItype __attribute__ ((mode (HI))); |
---|
111 | typedef int SItype __attribute__ ((mode (SI))); |
---|
112 | typedef int DItype __attribute__ ((mode (DI))); |
---|
113 | |
---|
114 | /* The type of the result of a fp compare */ |
---|
115 | #ifndef CMPtype |
---|
116 | #define CMPtype SItype |
---|
117 | #endif |
---|
118 | |
---|
119 | typedef unsigned int UHItype __attribute__ ((mode (HI))); |
---|
120 | typedef unsigned int USItype __attribute__ ((mode (SI))); |
---|
121 | typedef unsigned int UDItype __attribute__ ((mode (DI))); |
---|
122 | |
---|
123 | #define MAX_SI_INT ((SItype) ((unsigned) (~0)>>1)) |
---|
124 | #define MAX_USI_INT ((USItype) ~0) |
---|
125 | |
---|
126 | |
---|
127 | #ifdef FLOAT_ONLY |
---|
128 | #define NO_DI_MODE |
---|
129 | #endif |
---|
130 | |
---|
131 | #ifdef FLOAT |
---|
132 | # define NGARDS 7L |
---|
133 | # define GARDROUND 0x3f |
---|
134 | # define GARDMASK 0x7f |
---|
135 | # define GARDMSB 0x40 |
---|
136 | # define EXPBITS 8 |
---|
137 | # define EXPBIAS 127 |
---|
138 | # define FRACBITS 23 |
---|
139 | # define EXPMAX (0xff) |
---|
140 | # define QUIET_NAN 0x100000L |
---|
141 | # define FRAC_NBITS 32 |
---|
142 | # define FRACHIGH 0x80000000L |
---|
143 | # define FRACHIGH2 0xc0000000L |
---|
144 | # define pack_d pack_f |
---|
145 | # define unpack_d unpack_f |
---|
146 | typedef USItype fractype; |
---|
147 | typedef UHItype halffractype; |
---|
148 | typedef SFtype FLO_type; |
---|
149 | typedef SItype intfrac; |
---|
150 | |
---|
151 | #else |
---|
152 | # define PREFIXFPDP dp |
---|
153 | # define PREFIXSFDF df |
---|
154 | # define NGARDS 8L |
---|
155 | # define GARDROUND 0x7f |
---|
156 | # define GARDMASK 0xff |
---|
157 | # define GARDMSB 0x80 |
---|
158 | # define EXPBITS 11 |
---|
159 | # define EXPBIAS 1023 |
---|
160 | # define FRACBITS 52 |
---|
161 | # define EXPMAX (0x7ff) |
---|
162 | # define QUIET_NAN 0x8000000000000LL |
---|
163 | # define FRAC_NBITS 64 |
---|
164 | # define FRACHIGH 0x8000000000000000LL |
---|
165 | # define FRACHIGH2 0xc000000000000000LL |
---|
166 | typedef UDItype fractype; |
---|
167 | typedef USItype halffractype; |
---|
168 | typedef DFtype FLO_type; |
---|
169 | typedef DItype intfrac; |
---|
170 | #endif |
---|
171 | |
---|
172 | #ifdef US_SOFTWARE_GOFAST |
---|
173 | # ifdef FLOAT |
---|
174 | # define add fpadd |
---|
175 | # define sub fpsub |
---|
176 | # define multiply fpmul |
---|
177 | # define divide fpdiv |
---|
178 | # define compare fpcmp |
---|
179 | # define si_to_float sitofp |
---|
180 | # define float_to_si fptosi |
---|
181 | # define float_to_usi fptoui |
---|
182 | # define negate __negsf2 |
---|
183 | # define sf_to_df fptodp |
---|
184 | # define dptofp dptofp |
---|
185 | #else |
---|
186 | # define add dpadd |
---|
187 | # define sub dpsub |
---|
188 | # define multiply dpmul |
---|
189 | # define divide dpdiv |
---|
190 | # define compare dpcmp |
---|
191 | # define si_to_float litodp |
---|
192 | # define float_to_si dptoli |
---|
193 | # define float_to_usi dptoul |
---|
194 | # define negate __negdf2 |
---|
195 | # define df_to_sf dptofp |
---|
196 | #endif |
---|
197 | #else |
---|
198 | # ifdef FLOAT |
---|
199 | # define add __addsf3 |
---|
200 | # define sub __subsf3 |
---|
201 | # define multiply __mulsf3 |
---|
202 | # define divide __divsf3 |
---|
203 | # define compare __cmpsf2 |
---|
204 | # define _eq_f2 __eqsf2 |
---|
205 | # define _ne_f2 __nesf2 |
---|
206 | # define _gt_f2 __gtsf2 |
---|
207 | # define _ge_f2 __gesf2 |
---|
208 | # define _lt_f2 __ltsf2 |
---|
209 | # define _le_f2 __lesf2 |
---|
210 | # define si_to_float __floatsisf |
---|
211 | # define float_to_si __fixsfsi |
---|
212 | # define float_to_usi __fixunssfsi |
---|
213 | # define negate __negsf2 |
---|
214 | # define sf_to_df __extendsfdf2 |
---|
215 | #else |
---|
216 | # define add __adddf3 |
---|
217 | # define sub __subdf3 |
---|
218 | # define multiply __muldf3 |
---|
219 | # define divide __divdf3 |
---|
220 | # define compare __cmpdf2 |
---|
221 | # define _eq_f2 __eqdf2 |
---|
222 | # define _ne_f2 __nedf2 |
---|
223 | # define _gt_f2 __gtdf2 |
---|
224 | # define _ge_f2 __gedf2 |
---|
225 | # define _lt_f2 __ltdf2 |
---|
226 | # define _le_f2 __ledf2 |
---|
227 | # define si_to_float __floatsidf |
---|
228 | # define float_to_si __fixdfsi |
---|
229 | # define float_to_usi __fixunsdfsi |
---|
230 | # define negate __negdf2 |
---|
231 | # define df_to_sf __truncdfsf2 |
---|
232 | # endif |
---|
233 | #endif |
---|
234 | |
---|
235 | |
---|
236 | #define INLINE __inline__ |
---|
237 | |
---|
238 | /* Preserve the sticky-bit when shifting fractions to the right. */ |
---|
239 | #define LSHIFT(a) { a = (a & 1) | (a >> 1); } |
---|
240 | |
---|
241 | /* numeric parameters */ |
---|
242 | /* F_D_BITOFF is the number of bits offset between the MSB of the mantissa |
---|
243 | of a float and of a double. Assumes there are only two float types. |
---|
244 | (double::FRAC_BITS+double::NGARDS-(float::FRAC_BITS-float::NGARDS)) |
---|
245 | */ |
---|
246 | #define F_D_BITOFF (52+8-(23+7)) |
---|
247 | |
---|
248 | |
---|
249 | #define NORMAL_EXPMIN (-(EXPBIAS)+1) |
---|
250 | #define IMPLICIT_1 (1LL<<(FRACBITS+NGARDS)) |
---|
251 | #define IMPLICIT_2 (1LL<<(FRACBITS+1+NGARDS)) |
---|
252 | |
---|
253 | /* common types */ |
---|
254 | |
---|
255 | typedef enum |
---|
256 | { |
---|
257 | CLASS_SNAN, |
---|
258 | CLASS_QNAN, |
---|
259 | CLASS_ZERO, |
---|
260 | CLASS_NUMBER, |
---|
261 | CLASS_INFINITY |
---|
262 | } fp_class_type; |
---|
263 | |
---|
264 | typedef struct |
---|
265 | { |
---|
266 | #ifdef SMALL_MACHINE |
---|
267 | char class; |
---|
268 | unsigned char sign; |
---|
269 | short normal_exp; |
---|
270 | #else |
---|
271 | fp_class_type class; |
---|
272 | unsigned int sign; |
---|
273 | int normal_exp; |
---|
274 | #endif |
---|
275 | |
---|
276 | union |
---|
277 | { |
---|
278 | fractype ll; |
---|
279 | halffractype l[2]; |
---|
280 | } fraction; |
---|
281 | } fp_number_type; |
---|
282 | |
---|
283 | typedef union |
---|
284 | { |
---|
285 | FLO_type value; |
---|
286 | fractype value_raw; |
---|
287 | |
---|
288 | #ifndef FLOAT |
---|
289 | halffractype words[2]; |
---|
290 | #endif |
---|
291 | |
---|
292 | #ifdef FLOAT_BIT_ORDER_MISMATCH |
---|
293 | struct |
---|
294 | { |
---|
295 | fractype fraction:FRACBITS __attribute__ ((packed)); |
---|
296 | unsigned int exp:EXPBITS __attribute__ ((packed)); |
---|
297 | unsigned int sign:1 __attribute__ ((packed)); |
---|
298 | } |
---|
299 | bits; |
---|
300 | #endif |
---|
301 | |
---|
302 | #ifdef _DEBUG_BITFLOAT |
---|
303 | struct |
---|
304 | { |
---|
305 | unsigned int sign:1 __attribute__ ((packed)); |
---|
306 | unsigned int exp:EXPBITS __attribute__ ((packed)); |
---|
307 | fractype fraction:FRACBITS __attribute__ ((packed)); |
---|
308 | } |
---|
309 | bits_big_endian; |
---|
310 | |
---|
311 | struct |
---|
312 | { |
---|
313 | fractype fraction:FRACBITS __attribute__ ((packed)); |
---|
314 | unsigned int exp:EXPBITS __attribute__ ((packed)); |
---|
315 | unsigned int sign:1 __attribute__ ((packed)); |
---|
316 | } |
---|
317 | bits_little_endian; |
---|
318 | #endif |
---|
319 | } |
---|
320 | FLO_union_type; |
---|
321 | |
---|
322 | |
---|
323 | /* end of header */ |
---|
324 | |
---|
325 | /* IEEE "special" number predicates */ |
---|
326 | |
---|
327 | #ifdef NO_NANS |
---|
328 | |
---|
329 | #define nan() 0 |
---|
330 | #define isnan(x) 0 |
---|
331 | #define isinf(x) 0 |
---|
332 | #else |
---|
333 | |
---|
334 | INLINE |
---|
335 | static fp_number_type * |
---|
336 | nan () |
---|
337 | { |
---|
338 | static fp_number_type thenan; |
---|
339 | |
---|
340 | return &thenan; |
---|
341 | } |
---|
342 | |
---|
343 | INLINE |
---|
344 | static int |
---|
345 | isnan ( fp_number_type * x) |
---|
346 | { |
---|
347 | return x->class == CLASS_SNAN || x->class == CLASS_QNAN; |
---|
348 | } |
---|
349 | |
---|
350 | INLINE |
---|
351 | static int |
---|
352 | isinf ( fp_number_type * x) |
---|
353 | { |
---|
354 | return x->class == CLASS_INFINITY; |
---|
355 | } |
---|
356 | |
---|
357 | #endif |
---|
358 | |
---|
359 | INLINE |
---|
360 | static int |
---|
361 | iszero ( fp_number_type * x) |
---|
362 | { |
---|
363 | return x->class == CLASS_ZERO; |
---|
364 | } |
---|
365 | |
---|
366 | INLINE |
---|
367 | static void |
---|
368 | flip_sign ( fp_number_type * x) |
---|
369 | { |
---|
370 | x->sign = !x->sign; |
---|
371 | } |
---|
372 | |
---|
373 | static FLO_type |
---|
374 | pack_d ( fp_number_type * src) |
---|
375 | { |
---|
376 | FLO_union_type dst; |
---|
377 | fractype fraction = src->fraction.ll; /* wasn't unsigned before? */ |
---|
378 | int sign = src->sign; |
---|
379 | int exp = 0; |
---|
380 | |
---|
381 | if (isnan (src)) |
---|
382 | { |
---|
383 | exp = EXPMAX; |
---|
384 | if (src->class == CLASS_QNAN || 1) |
---|
385 | { |
---|
386 | fraction |= QUIET_NAN; |
---|
387 | } |
---|
388 | } |
---|
389 | else if (isinf (src)) |
---|
390 | { |
---|
391 | exp = EXPMAX; |
---|
392 | fraction = 0; |
---|
393 | } |
---|
394 | else if (iszero (src)) |
---|
395 | { |
---|
396 | exp = 0; |
---|
397 | fraction = 0; |
---|
398 | } |
---|
399 | else if (fraction == 0) |
---|
400 | { |
---|
401 | exp = 0; |
---|
402 | sign = 0; |
---|
403 | } |
---|
404 | else |
---|
405 | { |
---|
406 | if (src->normal_exp < NORMAL_EXPMIN) |
---|
407 | { |
---|
408 | /* This number's exponent is too low to fit into the bits |
---|
409 | available in the number, so we'll store 0 in the exponent and |
---|
410 | shift the fraction to the right to make up for it. */ |
---|
411 | |
---|
412 | int shift = NORMAL_EXPMIN - src->normal_exp; |
---|
413 | |
---|
414 | exp = 0; |
---|
415 | |
---|
416 | if (shift > FRAC_NBITS - NGARDS) |
---|
417 | { |
---|
418 | /* No point shifting, since it's more that 64 out. */ |
---|
419 | fraction = 0; |
---|
420 | } |
---|
421 | else |
---|
422 | { |
---|
423 | /* Shift by the value */ |
---|
424 | fraction >>= shift; |
---|
425 | } |
---|
426 | fraction >>= NGARDS; |
---|
427 | } |
---|
428 | else if (src->normal_exp > EXPBIAS) |
---|
429 | { |
---|
430 | exp = EXPMAX; |
---|
431 | fraction = 0; |
---|
432 | } |
---|
433 | else |
---|
434 | { |
---|
435 | exp = src->normal_exp + EXPBIAS; |
---|
436 | /* IF the gard bits are the all zero, but the first, then we're |
---|
437 | half way between two numbers, choose the one which makes the |
---|
438 | lsb of the answer 0. */ |
---|
439 | if ((fraction & GARDMASK) == GARDMSB) |
---|
440 | { |
---|
441 | if (fraction & (1 << NGARDS)) |
---|
442 | fraction += GARDROUND + 1; |
---|
443 | } |
---|
444 | else |
---|
445 | { |
---|
446 | /* Add a one to the guards to round up */ |
---|
447 | fraction += GARDROUND; |
---|
448 | } |
---|
449 | if (fraction >= IMPLICIT_2) |
---|
450 | { |
---|
451 | fraction >>= 1; |
---|
452 | exp += 1; |
---|
453 | } |
---|
454 | fraction >>= NGARDS; |
---|
455 | } |
---|
456 | } |
---|
457 | |
---|
458 | /* We previously used bitfields to store the number, but this doesn't |
---|
459 | handle little/big endian systems conveniently, so use shifts and |
---|
460 | masks */ |
---|
461 | #ifdef FLOAT_BIT_ORDER_MISMATCH |
---|
462 | dst.bits.fraction = fraction; |
---|
463 | dst.bits.exp = exp; |
---|
464 | dst.bits.sign = sign; |
---|
465 | #else |
---|
466 | dst.value_raw = fraction & ((((fractype)1) << FRACBITS) - (fractype)1); |
---|
467 | dst.value_raw |= ((fractype) (exp & ((1 << EXPBITS) - 1))) << FRACBITS; |
---|
468 | dst.value_raw |= ((fractype) (sign & 1)) << (FRACBITS | EXPBITS); |
---|
469 | #endif |
---|
470 | |
---|
471 | #if defined(FLOAT_WORD_ORDER_MISMATCH) && !defined(FLOAT) |
---|
472 | { |
---|
473 | halffractype tmp = dst.words[0]; |
---|
474 | dst.words[0] = dst.words[1]; |
---|
475 | dst.words[1] = tmp; |
---|
476 | } |
---|
477 | #endif |
---|
478 | |
---|
479 | return dst.value; |
---|
480 | } |
---|
481 | |
---|
482 | static void |
---|
483 | unpack_d (FLO_union_type * src, fp_number_type * dst) |
---|
484 | { |
---|
485 | /* We previously used bitfields to store the number, but this doesn't |
---|
486 | handle little/big endian systems conveniently, so use shifts and |
---|
487 | masks */ |
---|
488 | fractype fraction; |
---|
489 | int exp; |
---|
490 | int sign; |
---|
491 | |
---|
492 | #if defined(FLOAT_WORD_ORDER_MISMATCH) && !defined(FLOAT) |
---|
493 | FLO_union_type swapped; |
---|
494 | |
---|
495 | swapped.words[0] = src->words[1]; |
---|
496 | swapped.words[1] = src->words[0]; |
---|
497 | src = &swapped; |
---|
498 | #endif |
---|
499 | |
---|
500 | #ifdef FLOAT_BIT_ORDER_MISMATCH |
---|
501 | fraction = src->bits.fraction; |
---|
502 | exp = src->bits.exp; |
---|
503 | sign = src->bits.sign; |
---|
504 | #else |
---|
505 | fraction = src->value_raw & ((((fractype)1) << FRACBITS) - (fractype)1); |
---|
506 | exp = ((int)(src->value_raw >> FRACBITS)) & ((1 << EXPBITS) - 1); |
---|
507 | sign = ((int)(src->value_raw >> (FRACBITS + EXPBITS))) & 1; |
---|
508 | #endif |
---|
509 | |
---|
510 | dst->sign = sign; |
---|
511 | if (exp == 0) |
---|
512 | { |
---|
513 | /* Hmm. Looks like 0 */ |
---|
514 | if (fraction == 0) |
---|
515 | { |
---|
516 | /* tastes like zero */ |
---|
517 | dst->class = CLASS_ZERO; |
---|
518 | } |
---|
519 | else |
---|
520 | { |
---|
521 | /* Zero exponent with non zero fraction - it's denormalized, |
---|
522 | so there isn't a leading implicit one - we'll shift it so |
---|
523 | it gets one. */ |
---|
524 | dst->normal_exp = exp - EXPBIAS + 1; |
---|
525 | fraction <<= NGARDS; |
---|
526 | |
---|
527 | dst->class = CLASS_NUMBER; |
---|
528 | #if 1 |
---|
529 | while (fraction < IMPLICIT_1) |
---|
530 | { |
---|
531 | fraction <<= 1; |
---|
532 | dst->normal_exp--; |
---|
533 | } |
---|
534 | #endif |
---|
535 | dst->fraction.ll = fraction; |
---|
536 | } |
---|
537 | } |
---|
538 | else if (exp == EXPMAX) |
---|
539 | { |
---|
540 | /* Huge exponent*/ |
---|
541 | if (fraction == 0) |
---|
542 | { |
---|
543 | /* Attached to a zero fraction - means infinity */ |
---|
544 | dst->class = CLASS_INFINITY; |
---|
545 | } |
---|
546 | else |
---|
547 | { |
---|
548 | /* Non zero fraction, means nan */ |
---|
549 | if (fraction & QUIET_NAN) |
---|
550 | { |
---|
551 | dst->class = CLASS_QNAN; |
---|
552 | } |
---|
553 | else |
---|
554 | { |
---|
555 | dst->class = CLASS_SNAN; |
---|
556 | } |
---|
557 | /* Keep the fraction part as the nan number */ |
---|
558 | dst->fraction.ll = fraction; |
---|
559 | } |
---|
560 | } |
---|
561 | else |
---|
562 | { |
---|
563 | /* Nothing strange about this number */ |
---|
564 | dst->normal_exp = exp - EXPBIAS; |
---|
565 | dst->class = CLASS_NUMBER; |
---|
566 | dst->fraction.ll = (fraction << NGARDS) | IMPLICIT_1; |
---|
567 | } |
---|
568 | } |
---|
569 | |
---|
570 | static fp_number_type * |
---|
571 | _fpadd_parts (fp_number_type * a, |
---|
572 | fp_number_type * b, |
---|
573 | fp_number_type * tmp) |
---|
574 | { |
---|
575 | intfrac tfraction; |
---|
576 | |
---|
577 | /* Put commonly used fields in local variables. */ |
---|
578 | int a_normal_exp; |
---|
579 | int b_normal_exp; |
---|
580 | fractype a_fraction; |
---|
581 | fractype b_fraction; |
---|
582 | |
---|
583 | if (isnan (a)) |
---|
584 | { |
---|
585 | return a; |
---|
586 | } |
---|
587 | if (isnan (b)) |
---|
588 | { |
---|
589 | return b; |
---|
590 | } |
---|
591 | if (isinf (a)) |
---|
592 | { |
---|
593 | /* Adding infinities with opposite signs yields a NaN. */ |
---|
594 | if (isinf (b) && a->sign != b->sign) |
---|
595 | return nan (); |
---|
596 | return a; |
---|
597 | } |
---|
598 | if (isinf (b)) |
---|
599 | { |
---|
600 | return b; |
---|
601 | } |
---|
602 | if (iszero (b)) |
---|
603 | { |
---|
604 | return a; |
---|
605 | } |
---|
606 | if (iszero (a)) |
---|
607 | { |
---|
608 | return b; |
---|
609 | } |
---|
610 | |
---|
611 | /* Got two numbers. shift the smaller and increment the exponent till |
---|
612 | they're the same */ |
---|
613 | { |
---|
614 | int diff; |
---|
615 | |
---|
616 | a_normal_exp = a->normal_exp; |
---|
617 | b_normal_exp = b->normal_exp; |
---|
618 | a_fraction = a->fraction.ll; |
---|
619 | b_fraction = b->fraction.ll; |
---|
620 | |
---|
621 | diff = a_normal_exp - b_normal_exp; |
---|
622 | |
---|
623 | if (diff < 0) |
---|
624 | diff = -diff; |
---|
625 | if (diff < FRAC_NBITS) |
---|
626 | { |
---|
627 | /* ??? This does shifts one bit at a time. Optimize. */ |
---|
628 | while (a_normal_exp > b_normal_exp) |
---|
629 | { |
---|
630 | b_normal_exp++; |
---|
631 | LSHIFT (b_fraction); |
---|
632 | } |
---|
633 | while (b_normal_exp > a_normal_exp) |
---|
634 | { |
---|
635 | a_normal_exp++; |
---|
636 | LSHIFT (a_fraction); |
---|
637 | } |
---|
638 | } |
---|
639 | else |
---|
640 | { |
---|
641 | /* Somethings's up.. choose the biggest */ |
---|
642 | if (a_normal_exp > b_normal_exp) |
---|
643 | { |
---|
644 | b_normal_exp = a_normal_exp; |
---|
645 | b_fraction = 0; |
---|
646 | } |
---|
647 | else |
---|
648 | { |
---|
649 | a_normal_exp = b_normal_exp; |
---|
650 | a_fraction = 0; |
---|
651 | } |
---|
652 | } |
---|
653 | } |
---|
654 | |
---|
655 | if (a->sign != b->sign) |
---|
656 | { |
---|
657 | if (a->sign) |
---|
658 | { |
---|
659 | tfraction = -a_fraction + b_fraction; |
---|
660 | } |
---|
661 | else |
---|
662 | { |
---|
663 | tfraction = a_fraction - b_fraction; |
---|
664 | } |
---|
665 | if (tfraction > 0) |
---|
666 | { |
---|
667 | tmp->sign = 0; |
---|
668 | tmp->normal_exp = a_normal_exp; |
---|
669 | tmp->fraction.ll = tfraction; |
---|
670 | } |
---|
671 | else |
---|
672 | { |
---|
673 | tmp->sign = 1; |
---|
674 | tmp->normal_exp = a_normal_exp; |
---|
675 | tmp->fraction.ll = -tfraction; |
---|
676 | } |
---|
677 | /* and renormalize it */ |
---|
678 | |
---|
679 | while (tmp->fraction.ll < IMPLICIT_1 && tmp->fraction.ll) |
---|
680 | { |
---|
681 | tmp->fraction.ll <<= 1; |
---|
682 | tmp->normal_exp--; |
---|
683 | } |
---|
684 | } |
---|
685 | else |
---|
686 | { |
---|
687 | tmp->sign = a->sign; |
---|
688 | tmp->normal_exp = a_normal_exp; |
---|
689 | tmp->fraction.ll = a_fraction + b_fraction; |
---|
690 | } |
---|
691 | tmp->class = CLASS_NUMBER; |
---|
692 | /* Now the fraction is added, we have to shift down to renormalize the |
---|
693 | number */ |
---|
694 | |
---|
695 | if (tmp->fraction.ll >= IMPLICIT_2) |
---|
696 | { |
---|
697 | LSHIFT (tmp->fraction.ll); |
---|
698 | tmp->normal_exp++; |
---|
699 | } |
---|
700 | return tmp; |
---|
701 | |
---|
702 | } |
---|
703 | |
---|
704 | FLO_type |
---|
705 | add (FLO_type arg_a, FLO_type arg_b) |
---|
706 | { |
---|
707 | fp_number_type a; |
---|
708 | fp_number_type b; |
---|
709 | fp_number_type tmp; |
---|
710 | fp_number_type *res; |
---|
711 | |
---|
712 | unpack_d ((FLO_union_type *) & arg_a, &a); |
---|
713 | unpack_d ((FLO_union_type *) & arg_b, &b); |
---|
714 | |
---|
715 | res = _fpadd_parts (&a, &b, &tmp); |
---|
716 | |
---|
717 | return pack_d (res); |
---|
718 | } |
---|
719 | |
---|
720 | FLO_type |
---|
721 | sub (FLO_type arg_a, FLO_type arg_b) |
---|
722 | { |
---|
723 | fp_number_type a; |
---|
724 | fp_number_type b; |
---|
725 | fp_number_type tmp; |
---|
726 | fp_number_type *res; |
---|
727 | |
---|
728 | unpack_d ((FLO_union_type *) & arg_a, &a); |
---|
729 | unpack_d ((FLO_union_type *) & arg_b, &b); |
---|
730 | |
---|
731 | b.sign ^= 1; |
---|
732 | |
---|
733 | res = _fpadd_parts (&a, &b, &tmp); |
---|
734 | |
---|
735 | return pack_d (res); |
---|
736 | } |
---|
737 | |
---|
738 | static fp_number_type * |
---|
739 | _fpmul_parts ( fp_number_type * a, |
---|
740 | fp_number_type * b, |
---|
741 | fp_number_type * tmp) |
---|
742 | { |
---|
743 | fractype low = 0; |
---|
744 | fractype high = 0; |
---|
745 | |
---|
746 | if (isnan (a)) |
---|
747 | { |
---|
748 | a->sign = a->sign != b->sign; |
---|
749 | return a; |
---|
750 | } |
---|
751 | if (isnan (b)) |
---|
752 | { |
---|
753 | b->sign = a->sign != b->sign; |
---|
754 | return b; |
---|
755 | } |
---|
756 | if (isinf (a)) |
---|
757 | { |
---|
758 | if (iszero (b)) |
---|
759 | return nan (); |
---|
760 | a->sign = a->sign != b->sign; |
---|
761 | return a; |
---|
762 | } |
---|
763 | if (isinf (b)) |
---|
764 | { |
---|
765 | if (iszero (a)) |
---|
766 | { |
---|
767 | return nan (); |
---|
768 | } |
---|
769 | b->sign = a->sign != b->sign; |
---|
770 | return b; |
---|
771 | } |
---|
772 | if (iszero (a)) |
---|
773 | { |
---|
774 | a->sign = a->sign != b->sign; |
---|
775 | return a; |
---|
776 | } |
---|
777 | if (iszero (b)) |
---|
778 | { |
---|
779 | b->sign = a->sign != b->sign; |
---|
780 | return b; |
---|
781 | } |
---|
782 | |
---|
783 | /* Calculate the mantissa by multiplying both 64bit numbers to get a |
---|
784 | 128 bit number */ |
---|
785 | { |
---|
786 | fractype x = a->fraction.ll; |
---|
787 | fractype ylow = b->fraction.ll; |
---|
788 | fractype yhigh = 0; |
---|
789 | int bit; |
---|
790 | |
---|
791 | #if defined(NO_DI_MODE) |
---|
792 | { |
---|
793 | /* ??? This does multiplies one bit at a time. Optimize. */ |
---|
794 | for (bit = 0; bit < FRAC_NBITS; bit++) |
---|
795 | { |
---|
796 | int carry; |
---|
797 | |
---|
798 | if (x & 1) |
---|
799 | { |
---|
800 | carry = (low += ylow) < ylow; |
---|
801 | high += yhigh + carry; |
---|
802 | } |
---|
803 | yhigh <<= 1; |
---|
804 | if (ylow & FRACHIGH) |
---|
805 | { |
---|
806 | yhigh |= 1; |
---|
807 | } |
---|
808 | ylow <<= 1; |
---|
809 | x >>= 1; |
---|
810 | } |
---|
811 | } |
---|
812 | #elif defined(FLOAT) |
---|
813 | { |
---|
814 | /* Multiplying two 32 bit numbers to get a 64 bit number on |
---|
815 | a machine with DI, so we're safe */ |
---|
816 | |
---|
817 | DItype answer = (DItype)(a->fraction.ll) * (DItype)(b->fraction.ll); |
---|
818 | |
---|
819 | high = answer >> 32; |
---|
820 | low = answer; |
---|
821 | } |
---|
822 | #else |
---|
823 | /* Doing a 64*64 to 128 */ |
---|
824 | { |
---|
825 | UDItype nl = a->fraction.ll & 0xffffffff; |
---|
826 | UDItype nh = a->fraction.ll >> 32; |
---|
827 | UDItype ml = b->fraction.ll & 0xffffffff; |
---|
828 | UDItype mh = b->fraction.ll >>32; |
---|
829 | UDItype pp_ll = ml * nl; |
---|
830 | UDItype pp_hl = mh * nl; |
---|
831 | UDItype pp_lh = ml * nh; |
---|
832 | UDItype pp_hh = mh * nh; |
---|
833 | UDItype res2 = 0; |
---|
834 | UDItype res0 = 0; |
---|
835 | UDItype ps_hh__ = pp_hl + pp_lh; |
---|
836 | if (ps_hh__ < pp_hl) |
---|
837 | res2 += 0x100000000LL; |
---|
838 | pp_hl = (ps_hh__ << 32) & 0xffffffff00000000LL; |
---|
839 | res0 = pp_ll + pp_hl; |
---|
840 | if (res0 < pp_ll) |
---|
841 | res2++; |
---|
842 | res2 += ((ps_hh__ >> 32) & 0xffffffffL) + pp_hh; |
---|
843 | high = res2; |
---|
844 | low = res0; |
---|
845 | } |
---|
846 | #endif |
---|
847 | } |
---|
848 | |
---|
849 | tmp->normal_exp = a->normal_exp + b->normal_exp; |
---|
850 | tmp->sign = a->sign != b->sign; |
---|
851 | #ifdef FLOAT |
---|
852 | tmp->normal_exp += 2; /* ??????????????? */ |
---|
853 | #else |
---|
854 | tmp->normal_exp += 4; /* ??????????????? */ |
---|
855 | #endif |
---|
856 | while (high >= IMPLICIT_2) |
---|
857 | { |
---|
858 | tmp->normal_exp++; |
---|
859 | if (high & 1) |
---|
860 | { |
---|
861 | low >>= 1; |
---|
862 | low |= FRACHIGH; |
---|
863 | } |
---|
864 | high >>= 1; |
---|
865 | } |
---|
866 | while (high < IMPLICIT_1) |
---|
867 | { |
---|
868 | tmp->normal_exp--; |
---|
869 | |
---|
870 | high <<= 1; |
---|
871 | if (low & FRACHIGH) |
---|
872 | high |= 1; |
---|
873 | low <<= 1; |
---|
874 | } |
---|
875 | /* rounding is tricky. if we only round if it won't make us round later. */ |
---|
876 | #if 0 |
---|
877 | if (low & FRACHIGH2) |
---|
878 | { |
---|
879 | if (((high & GARDMASK) != GARDMSB) |
---|
880 | && (((high + 1) & GARDMASK) == GARDMSB)) |
---|
881 | { |
---|
882 | /* don't round, it gets done again later. */ |
---|
883 | } |
---|
884 | else |
---|
885 | { |
---|
886 | high++; |
---|
887 | } |
---|
888 | } |
---|
889 | #endif |
---|
890 | if ((high & GARDMASK) == GARDMSB) |
---|
891 | { |
---|
892 | if (high & (1 << NGARDS)) |
---|
893 | { |
---|
894 | /* half way, so round to even */ |
---|
895 | high += GARDROUND + 1; |
---|
896 | } |
---|
897 | else if (low) |
---|
898 | { |
---|
899 | /* but we really weren't half way */ |
---|
900 | high += GARDROUND + 1; |
---|
901 | } |
---|
902 | } |
---|
903 | tmp->fraction.ll = high; |
---|
904 | tmp->class = CLASS_NUMBER; |
---|
905 | return tmp; |
---|
906 | } |
---|
907 | |
---|
908 | FLO_type |
---|
909 | multiply (FLO_type arg_a, FLO_type arg_b) |
---|
910 | { |
---|
911 | fp_number_type a; |
---|
912 | fp_number_type b; |
---|
913 | fp_number_type tmp; |
---|
914 | fp_number_type *res; |
---|
915 | |
---|
916 | unpack_d ((FLO_union_type *) & arg_a, &a); |
---|
917 | unpack_d ((FLO_union_type *) & arg_b, &b); |
---|
918 | |
---|
919 | res = _fpmul_parts (&a, &b, &tmp); |
---|
920 | |
---|
921 | return pack_d (res); |
---|
922 | } |
---|
923 | |
---|
924 | static fp_number_type * |
---|
925 | _fpdiv_parts (fp_number_type * a, |
---|
926 | fp_number_type * b, |
---|
927 | fp_number_type * tmp) |
---|
928 | { |
---|
929 | fractype low = 0; |
---|
930 | fractype high = 0; |
---|
931 | fractype r0, r1, y0, y1, bit; |
---|
932 | fractype q; |
---|
933 | fractype numerator; |
---|
934 | fractype denominator; |
---|
935 | fractype quotient; |
---|
936 | fractype remainder; |
---|
937 | |
---|
938 | if (isnan (a)) |
---|
939 | { |
---|
940 | return a; |
---|
941 | } |
---|
942 | if (isnan (b)) |
---|
943 | { |
---|
944 | return b; |
---|
945 | } |
---|
946 | |
---|
947 | a->sign = a->sign ^ b->sign; |
---|
948 | |
---|
949 | if (isinf (a) || iszero (a)) |
---|
950 | { |
---|
951 | if (a->class == b->class) |
---|
952 | return nan (); |
---|
953 | return a; |
---|
954 | } |
---|
955 | |
---|
956 | if (isinf (b)) |
---|
957 | { |
---|
958 | a->fraction.ll = 0; |
---|
959 | a->normal_exp = 0; |
---|
960 | return a; |
---|
961 | } |
---|
962 | if (iszero (b)) |
---|
963 | { |
---|
964 | a->class = CLASS_INFINITY; |
---|
965 | return b; |
---|
966 | } |
---|
967 | |
---|
968 | /* Calculate the mantissa by multiplying both 64bit numbers to get a |
---|
969 | 128 bit number */ |
---|
970 | { |
---|
971 | int carry; |
---|
972 | intfrac d0, d1; /* weren't unsigned before ??? */ |
---|
973 | |
---|
974 | /* quotient = |
---|
975 | ( numerator / denominator) * 2^(numerator exponent - denominator exponent) |
---|
976 | */ |
---|
977 | |
---|
978 | a->normal_exp = a->normal_exp - b->normal_exp; |
---|
979 | numerator = a->fraction.ll; |
---|
980 | denominator = b->fraction.ll; |
---|
981 | |
---|
982 | if (numerator < denominator) |
---|
983 | { |
---|
984 | /* Fraction will be less than 1.0 */ |
---|
985 | numerator *= 2; |
---|
986 | a->normal_exp--; |
---|
987 | } |
---|
988 | bit = IMPLICIT_1; |
---|
989 | quotient = 0; |
---|
990 | /* ??? Does divide one bit at a time. Optimize. */ |
---|
991 | while (bit) |
---|
992 | { |
---|
993 | if (numerator >= denominator) |
---|
994 | { |
---|
995 | quotient |= bit; |
---|
996 | numerator -= denominator; |
---|
997 | } |
---|
998 | bit >>= 1; |
---|
999 | numerator *= 2; |
---|
1000 | } |
---|
1001 | |
---|
1002 | if ((quotient & GARDMASK) == GARDMSB) |
---|
1003 | { |
---|
1004 | if (quotient & (1 << NGARDS)) |
---|
1005 | { |
---|
1006 | /* half way, so round to even */ |
---|
1007 | quotient += GARDROUND + 1; |
---|
1008 | } |
---|
1009 | else if (numerator) |
---|
1010 | { |
---|
1011 | /* but we really weren't half way, more bits exist */ |
---|
1012 | quotient += GARDROUND + 1; |
---|
1013 | } |
---|
1014 | } |
---|
1015 | |
---|
1016 | a->fraction.ll = quotient; |
---|
1017 | return (a); |
---|
1018 | } |
---|
1019 | } |
---|
1020 | |
---|
1021 | FLO_type |
---|
1022 | divide (FLO_type arg_a, FLO_type arg_b) |
---|
1023 | { |
---|
1024 | fp_number_type a; |
---|
1025 | fp_number_type b; |
---|
1026 | fp_number_type tmp; |
---|
1027 | fp_number_type *res; |
---|
1028 | |
---|
1029 | unpack_d ((FLO_union_type *) & arg_a, &a); |
---|
1030 | unpack_d ((FLO_union_type *) & arg_b, &b); |
---|
1031 | |
---|
1032 | res = _fpdiv_parts (&a, &b, &tmp); |
---|
1033 | |
---|
1034 | return pack_d (res); |
---|
1035 | } |
---|
1036 | |
---|
1037 | /* according to the demo, fpcmp returns a comparison with 0... thus |
---|
1038 | a<b -> -1 |
---|
1039 | a==b -> 0 |
---|
1040 | a>b -> +1 |
---|
1041 | */ |
---|
1042 | |
---|
1043 | static int |
---|
1044 | _fpcmp_parts (fp_number_type * a, fp_number_type * b) |
---|
1045 | { |
---|
1046 | #if 0 |
---|
1047 | /* either nan -> unordered. Must be checked outside of this routine. */ |
---|
1048 | if (isnan (a) && isnan (b)) |
---|
1049 | { |
---|
1050 | return 1; /* still unordered! */ |
---|
1051 | } |
---|
1052 | #endif |
---|
1053 | |
---|
1054 | if (isnan (a) || isnan (b)) |
---|
1055 | { |
---|
1056 | return 1; /* how to indicate unordered compare? */ |
---|
1057 | } |
---|
1058 | if (isinf (a) && isinf (b)) |
---|
1059 | { |
---|
1060 | /* +inf > -inf, but +inf != +inf */ |
---|
1061 | /* b \a| +inf(0)| -inf(1) |
---|
1062 | ______\+--------+-------- |
---|
1063 | +inf(0)| a==b(0)| a<b(-1) |
---|
1064 | -------+--------+-------- |
---|
1065 | -inf(1)| a>b(1) | a==b(0) |
---|
1066 | -------+--------+-------- |
---|
1067 | So since unordered must be non zero, just line up the columns... |
---|
1068 | */ |
---|
1069 | return b->sign - a->sign; |
---|
1070 | } |
---|
1071 | /* but not both... */ |
---|
1072 | if (isinf (a)) |
---|
1073 | { |
---|
1074 | return a->sign ? -1 : 1; |
---|
1075 | } |
---|
1076 | if (isinf (b)) |
---|
1077 | { |
---|
1078 | return b->sign ? 1 : -1; |
---|
1079 | } |
---|
1080 | if (iszero (a) && iszero (b)) |
---|
1081 | { |
---|
1082 | return 0; |
---|
1083 | } |
---|
1084 | if (iszero (a)) |
---|
1085 | { |
---|
1086 | return b->sign ? 1 : -1; |
---|
1087 | } |
---|
1088 | if (iszero (b)) |
---|
1089 | { |
---|
1090 | return a->sign ? -1 : 1; |
---|
1091 | } |
---|
1092 | /* now both are "normal". */ |
---|
1093 | if (a->sign != b->sign) |
---|
1094 | { |
---|
1095 | /* opposite signs */ |
---|
1096 | return a->sign ? -1 : 1; |
---|
1097 | } |
---|
1098 | /* same sign; exponents? */ |
---|
1099 | if (a->normal_exp > b->normal_exp) |
---|
1100 | { |
---|
1101 | return a->sign ? -1 : 1; |
---|
1102 | } |
---|
1103 | if (a->normal_exp < b->normal_exp) |
---|
1104 | { |
---|
1105 | return a->sign ? 1 : -1; |
---|
1106 | } |
---|
1107 | /* same exponents; check size. */ |
---|
1108 | if (a->fraction.ll > b->fraction.ll) |
---|
1109 | { |
---|
1110 | return a->sign ? -1 : 1; |
---|
1111 | } |
---|
1112 | if (a->fraction.ll < b->fraction.ll) |
---|
1113 | { |
---|
1114 | return a->sign ? 1 : -1; |
---|
1115 | } |
---|
1116 | /* after all that, they're equal. */ |
---|
1117 | return 0; |
---|
1118 | } |
---|
1119 | |
---|
1120 | CMPtype |
---|
1121 | compare (FLO_type arg_a, FLO_type arg_b) |
---|
1122 | { |
---|
1123 | fp_number_type a; |
---|
1124 | fp_number_type b; |
---|
1125 | |
---|
1126 | unpack_d ((FLO_union_type *) & arg_a, &a); |
---|
1127 | unpack_d ((FLO_union_type *) & arg_b, &b); |
---|
1128 | |
---|
1129 | return _fpcmp_parts (&a, &b); |
---|
1130 | } |
---|
1131 | |
---|
1132 | #ifndef US_SOFTWARE_GOFAST |
---|
1133 | |
---|
1134 | /* These should be optimized for their specific tasks someday. */ |
---|
1135 | |
---|
1136 | CMPtype |
---|
1137 | _eq_f2 (FLO_type arg_a, FLO_type arg_b) |
---|
1138 | { |
---|
1139 | fp_number_type a; |
---|
1140 | fp_number_type b; |
---|
1141 | |
---|
1142 | unpack_d ((FLO_union_type *) & arg_a, &a); |
---|
1143 | unpack_d ((FLO_union_type *) & arg_b, &b); |
---|
1144 | |
---|
1145 | if (isnan (&a) || isnan (&b)) |
---|
1146 | return 1; /* false, truth == 0 */ |
---|
1147 | |
---|
1148 | return _fpcmp_parts (&a, &b) ; |
---|
1149 | } |
---|
1150 | |
---|
1151 | CMPtype |
---|
1152 | _ne_f2 (FLO_type arg_a, FLO_type arg_b) |
---|
1153 | { |
---|
1154 | fp_number_type a; |
---|
1155 | fp_number_type b; |
---|
1156 | |
---|
1157 | unpack_d ((FLO_union_type *) & arg_a, &a); |
---|
1158 | unpack_d ((FLO_union_type *) & arg_b, &b); |
---|
1159 | |
---|
1160 | if (isnan (&a) || isnan (&b)) |
---|
1161 | return 1; /* true, truth != 0 */ |
---|
1162 | |
---|
1163 | return _fpcmp_parts (&a, &b) ; |
---|
1164 | } |
---|
1165 | |
---|
1166 | CMPtype |
---|
1167 | _gt_f2 (FLO_type arg_a, FLO_type arg_b) |
---|
1168 | { |
---|
1169 | fp_number_type a; |
---|
1170 | fp_number_type b; |
---|
1171 | |
---|
1172 | unpack_d ((FLO_union_type *) & arg_a, &a); |
---|
1173 | unpack_d ((FLO_union_type *) & arg_b, &b); |
---|
1174 | |
---|
1175 | if (isnan (&a) || isnan (&b)) |
---|
1176 | return -1; /* false, truth > 0 */ |
---|
1177 | |
---|
1178 | return _fpcmp_parts (&a, &b); |
---|
1179 | } |
---|
1180 | |
---|
1181 | CMPtype |
---|
1182 | _ge_f2 (FLO_type arg_a, FLO_type arg_b) |
---|
1183 | { |
---|
1184 | fp_number_type a; |
---|
1185 | fp_number_type b; |
---|
1186 | |
---|
1187 | unpack_d ((FLO_union_type *) & arg_a, &a); |
---|
1188 | unpack_d ((FLO_union_type *) & arg_b, &b); |
---|
1189 | |
---|
1190 | if (isnan (&a) || isnan (&b)) |
---|
1191 | return -1; /* false, truth >= 0 */ |
---|
1192 | return _fpcmp_parts (&a, &b) ; |
---|
1193 | } |
---|
1194 | |
---|
1195 | CMPtype |
---|
1196 | _lt_f2 (FLO_type arg_a, FLO_type arg_b) |
---|
1197 | { |
---|
1198 | fp_number_type a; |
---|
1199 | fp_number_type b; |
---|
1200 | |
---|
1201 | unpack_d ((FLO_union_type *) & arg_a, &a); |
---|
1202 | unpack_d ((FLO_union_type *) & arg_b, &b); |
---|
1203 | |
---|
1204 | if (isnan (&a) || isnan (&b)) |
---|
1205 | return 1; /* false, truth < 0 */ |
---|
1206 | |
---|
1207 | return _fpcmp_parts (&a, &b); |
---|
1208 | } |
---|
1209 | |
---|
1210 | CMPtype |
---|
1211 | _le_f2 (FLO_type arg_a, FLO_type arg_b) |
---|
1212 | { |
---|
1213 | fp_number_type a; |
---|
1214 | fp_number_type b; |
---|
1215 | |
---|
1216 | unpack_d ((FLO_union_type *) & arg_a, &a); |
---|
1217 | unpack_d ((FLO_union_type *) & arg_b, &b); |
---|
1218 | |
---|
1219 | if (isnan (&a) || isnan (&b)) |
---|
1220 | return 1; /* false, truth <= 0 */ |
---|
1221 | |
---|
1222 | return _fpcmp_parts (&a, &b) ; |
---|
1223 | } |
---|
1224 | |
---|
1225 | #endif /* ! US_SOFTWARE_GOFAST */ |
---|
1226 | |
---|
1227 | FLO_type |
---|
1228 | si_to_float (SItype arg_a) |
---|
1229 | { |
---|
1230 | fp_number_type in; |
---|
1231 | |
---|
1232 | in.class = CLASS_NUMBER; |
---|
1233 | in.sign = arg_a < 0; |
---|
1234 | if (!arg_a) |
---|
1235 | { |
---|
1236 | in.class = CLASS_ZERO; |
---|
1237 | } |
---|
1238 | else |
---|
1239 | { |
---|
1240 | in.normal_exp = FRACBITS + NGARDS; |
---|
1241 | if (in.sign) |
---|
1242 | { |
---|
1243 | /* Special case for minint, since there is no +ve integer |
---|
1244 | representation for it */ |
---|
1245 | if (arg_a == 0x80000000) |
---|
1246 | { |
---|
1247 | return -2147483648.0; |
---|
1248 | } |
---|
1249 | in.fraction.ll = (-arg_a); |
---|
1250 | } |
---|
1251 | else |
---|
1252 | in.fraction.ll = arg_a; |
---|
1253 | |
---|
1254 | while (in.fraction.ll < (1LL << (FRACBITS + NGARDS))) |
---|
1255 | { |
---|
1256 | in.fraction.ll <<= 1; |
---|
1257 | in.normal_exp -= 1; |
---|
1258 | } |
---|
1259 | } |
---|
1260 | return pack_d (&in); |
---|
1261 | } |
---|
1262 | |
---|
1263 | SItype |
---|
1264 | float_to_si (FLO_type arg_a) |
---|
1265 | { |
---|
1266 | fp_number_type a; |
---|
1267 | SItype tmp; |
---|
1268 | |
---|
1269 | unpack_d ((FLO_union_type *) & arg_a, &a); |
---|
1270 | if (iszero (&a)) |
---|
1271 | return 0; |
---|
1272 | if (isnan (&a)) |
---|
1273 | return 0; |
---|
1274 | /* get reasonable MAX_SI_INT... */ |
---|
1275 | if (isinf (&a)) |
---|
1276 | return a.sign ? (-MAX_SI_INT)-1 : MAX_SI_INT; |
---|
1277 | /* it is a number, but a small one */ |
---|
1278 | if (a.normal_exp < 0) |
---|
1279 | return 0; |
---|
1280 | if (a.normal_exp > 30) |
---|
1281 | return a.sign ? (-MAX_SI_INT)-1 : MAX_SI_INT; |
---|
1282 | tmp = a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp); |
---|
1283 | return a.sign ? (-tmp) : (tmp); |
---|
1284 | } |
---|
1285 | |
---|
1286 | #ifdef US_SOFTWARE_GOFAST |
---|
1287 | /* While libgcc2.c defines its own __fixunssfsi and __fixunsdfsi routines, |
---|
1288 | we also define them for GOFAST because the ones in libgcc2.c have the |
---|
1289 | wrong names and I'd rather define these here and keep GOFAST CYG-LOC's |
---|
1290 | out of libgcc2.c. We can't define these here if not GOFAST because then |
---|
1291 | there'd be duplicate copies. */ |
---|
1292 | |
---|
1293 | USItype |
---|
1294 | float_to_usi (FLO_type arg_a) |
---|
1295 | { |
---|
1296 | fp_number_type a; |
---|
1297 | |
---|
1298 | unpack_d ((FLO_union_type *) & arg_a, &a); |
---|
1299 | if (iszero (&a)) |
---|
1300 | return 0; |
---|
1301 | if (isnan (&a)) |
---|
1302 | return 0; |
---|
1303 | /* it is a negative number */ |
---|
1304 | if (a.sign) |
---|
1305 | return 0; |
---|
1306 | /* get reasonable MAX_USI_INT... */ |
---|
1307 | if (isinf (&a)) |
---|
1308 | return MAX_USI_INT; |
---|
1309 | /* it is a number, but a small one */ |
---|
1310 | if (a.normal_exp < 0) |
---|
1311 | return 0; |
---|
1312 | if (a.normal_exp > 31) |
---|
1313 | return MAX_USI_INT; |
---|
1314 | else if (a.normal_exp > (FRACBITS + NGARDS)) |
---|
1315 | return a.fraction.ll << (a.normal_exp - (FRACBITS + NGARDS)); |
---|
1316 | else |
---|
1317 | return a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp); |
---|
1318 | } |
---|
1319 | #endif |
---|
1320 | |
---|
1321 | FLO_type |
---|
1322 | negate (FLO_type arg_a) |
---|
1323 | { |
---|
1324 | fp_number_type a; |
---|
1325 | |
---|
1326 | unpack_d ((FLO_union_type *) & arg_a, &a); |
---|
1327 | flip_sign (&a); |
---|
1328 | return pack_d (&a); |
---|
1329 | } |
---|
1330 | |
---|
1331 | #ifdef FLOAT |
---|
1332 | |
---|
1333 | SFtype |
---|
1334 | __make_fp(fp_class_type class, |
---|
1335 | unsigned int sign, |
---|
1336 | int exp, |
---|
1337 | USItype frac) |
---|
1338 | { |
---|
1339 | fp_number_type in; |
---|
1340 | |
---|
1341 | in.class = class; |
---|
1342 | in.sign = sign; |
---|
1343 | in.normal_exp = exp; |
---|
1344 | in.fraction.ll = frac; |
---|
1345 | return pack_d (&in); |
---|
1346 | } |
---|
1347 | |
---|
1348 | #ifndef FLOAT_ONLY |
---|
1349 | |
---|
1350 | /* This enables one to build an fp library that supports float but not double. |
---|
1351 | Otherwise, we would get an undefined reference to __make_dp. |
---|
1352 | This is needed for some 8-bit ports that can't handle well values that |
---|
1353 | are 8-bytes in size, so we just don't support double for them at all. */ |
---|
1354 | |
---|
1355 | extern DFtype __make_dp (fp_class_type, unsigned int, int, UDItype frac); |
---|
1356 | |
---|
1357 | DFtype |
---|
1358 | sf_to_df (SFtype arg_a) |
---|
1359 | { |
---|
1360 | fp_number_type in; |
---|
1361 | |
---|
1362 | unpack_d ((FLO_union_type *) & arg_a, &in); |
---|
1363 | return __make_dp (in.class, in.sign, in.normal_exp, |
---|
1364 | ((UDItype) in.fraction.ll) << F_D_BITOFF); |
---|
1365 | } |
---|
1366 | |
---|
1367 | #endif |
---|
1368 | #endif |
---|
1369 | |
---|
1370 | #ifndef FLOAT |
---|
1371 | |
---|
1372 | extern SFtype __make_fp (fp_class_type, unsigned int, int, USItype); |
---|
1373 | |
---|
1374 | DFtype |
---|
1375 | __make_dp (fp_class_type class, unsigned int sign, int exp, UDItype frac) |
---|
1376 | { |
---|
1377 | fp_number_type in; |
---|
1378 | |
---|
1379 | in.class = class; |
---|
1380 | in.sign = sign; |
---|
1381 | in.normal_exp = exp; |
---|
1382 | in.fraction.ll = frac; |
---|
1383 | return pack_d (&in); |
---|
1384 | } |
---|
1385 | |
---|
1386 | SFtype |
---|
1387 | df_to_sf (DFtype arg_a) |
---|
1388 | { |
---|
1389 | fp_number_type in; |
---|
1390 | USItype sffrac; |
---|
1391 | |
---|
1392 | unpack_d ((FLO_union_type *) & arg_a, &in); |
---|
1393 | |
---|
1394 | sffrac = in.fraction.ll >> F_D_BITOFF; |
---|
1395 | |
---|
1396 | /* We set the lowest guard bit in SFFRAC if we discarded any non |
---|
1397 | zero bits. */ |
---|
1398 | if ((in.fraction.ll & (((USItype) 1 << F_D_BITOFF) - 1)) != 0) |
---|
1399 | sffrac |= 1; |
---|
1400 | |
---|
1401 | return __make_fp (in.class, in.sign, in.normal_exp, sffrac); |
---|
1402 | } |
---|
1403 | |
---|
1404 | #endif |
---|
1405 | #endif /* !EXTENDED_FLOAT_STUBS */ |
---|