0
|
1 /* This is a stripped down version of floatlib.c. It supplies only those
|
|
2 functions which exist in libgcc, but for which there is not assembly
|
|
3 language versions in m68k/lb1sf68.asm.
|
|
4
|
|
5 It also includes simplistic support for extended floats (by working in
|
|
6 double precision). You must compile this file again with -DEXTFLOAT
|
|
7 to get this support. */
|
|
8
|
|
9 /*
|
|
10 ** gnulib support for software floating point.
|
|
11 ** Copyright (C) 1991 by Pipeline Associates, Inc. All rights reserved.
|
|
12 ** Permission is granted to do *anything* you want with this file,
|
|
13 ** commercial or otherwise, provided this message remains intact. So there!
|
|
14 ** I would appreciate receiving any updates/patches/changes that anyone
|
|
15 ** makes, and am willing to be the repository for said changes (am I
|
|
16 ** making a big mistake?).
|
|
17 **
|
|
18 ** Pat Wood
|
|
19 ** Pipeline Associates, Inc.
|
|
20 ** pipeline!phw@motown.com or
|
|
21 ** sun!pipeline!phw or
|
|
22 ** uunet!motown!pipeline!phw
|
|
23 **
|
|
24 ** 05/01/91 -- V1.0 -- first release to gcc mailing lists
|
|
25 ** 05/04/91 -- V1.1 -- added float and double prototypes and return values
|
|
26 ** -- fixed problems with adding and subtracting zero
|
|
27 ** -- fixed rounding in truncdfsf2
|
|
28 ** -- fixed SWAP define and tested on 386
|
|
29 */
|
|
30
|
|
31 /*
|
|
32 ** The following are routines that replace the gnulib soft floating point
|
|
33 ** routines that are called automatically when -msoft-float is selected.
|
|
34 ** The support single and double precision IEEE format, with provisions
|
|
35 ** for byte-swapped machines (tested on 386). Some of the double-precision
|
|
36 ** routines work at full precision, but most of the hard ones simply punt
|
|
37 ** and call the single precision routines, producing a loss of accuracy.
|
|
38 ** long long support is not assumed or included.
|
|
39 ** Overall accuracy is close to IEEE (actually 68882) for single-precision
|
|
40 ** arithmetic. I think there may still be a 1 in 1000 chance of a bit
|
|
41 ** being rounded the wrong way during a multiply. I'm not fussy enough to
|
|
42 ** bother with it, but if anyone is, knock yourself out.
|
|
43 **
|
|
44 ** Efficiency has only been addressed where it was obvious that something
|
|
45 ** would make a big difference. Anyone who wants to do this right for
|
|
46 ** best speed should go in and rewrite in assembler.
|
|
47 **
|
|
48 ** I have tested this only on a 68030 workstation and 386/ix integrated
|
|
49 ** in with -msoft-float.
|
|
50 */
|
|
51
|
|
52 /* the following deal with IEEE single-precision numbers */
|
|
53 #define EXCESS 126L
|
|
54 #define SIGNBIT 0x80000000L
|
|
55 #define HIDDEN (1L << 23L)
|
|
56 #define SIGN(fp) ((fp) & SIGNBIT)
|
|
57 #define EXP(fp) (((fp) >> 23L) & 0xFF)
|
|
58 #define MANT(fp) (((fp) & 0x7FFFFFL) | HIDDEN)
|
|
59 #define PACK(s,e,m) ((s) | ((e) << 23L) | (m))
|
|
60
|
|
61 /* the following deal with IEEE double-precision numbers */
|
|
62 #define EXCESSD 1022L
|
|
63 #define HIDDEND (1L << 20L)
|
|
64 #define EXPDBITS 11
|
|
65 #define EXPDMASK 0x7FFL
|
|
66 #define EXPD(fp) (((fp.l.upper) >> 20L) & 0x7FFL)
|
|
67 #define SIGND(fp) ((fp.l.upper) & SIGNBIT)
|
|
68 #define MANTD(fp) (((((fp.l.upper) & 0xFFFFF) | HIDDEND) << 10) | \
|
|
69 (fp.l.lower >> 22))
|
|
70 #define MANTDMASK 0xFFFFFL /* mask of upper part */
|
|
71
|
|
72 /* the following deal with IEEE extended-precision numbers */
|
|
73 #define EXCESSX 16382L
|
|
74 #define HIDDENX (1L << 31L)
|
|
75 #define EXPXBITS 15
|
|
76 #define EXPXMASK 0x7FFF
|
|
77 #define EXPX(fp) (((fp.l.upper) >> 16) & EXPXMASK)
|
|
78 #define SIGNX(fp) ((fp.l.upper) & SIGNBIT)
|
|
79 #define MANTXMASK 0x7FFFFFFFL /* mask of upper part */
|
|
80
|
|
81 union double_long
|
|
82 {
|
|
83 double d;
|
|
84 struct {
|
|
85 long upper;
|
|
86 unsigned long lower;
|
|
87 } l;
|
|
88 };
|
|
89
|
|
90 union float_long {
|
|
91 float f;
|
|
92 long l;
|
|
93 };
|
|
94
|
|
95 union long_double_long
|
|
96 {
|
|
97 long double ld;
|
|
98 struct
|
|
99 {
|
|
100 long upper;
|
|
101 unsigned long middle;
|
|
102 unsigned long lower;
|
|
103 } l;
|
|
104 };
|
|
105
|
|
106 #ifndef EXTFLOAT
|
|
107
|
|
108 int
|
|
109 __unordsf2(float a, float b)
|
|
110 {
|
|
111 union float_long fl;
|
|
112
|
|
113 fl.f = a;
|
|
114 if (EXP(fl.l) == EXP(~0u) && (MANT(fl.l) & ~HIDDEN) != 0)
|
|
115 return 1;
|
|
116 fl.f = b;
|
|
117 if (EXP(fl.l) == EXP(~0u) && (MANT(fl.l) & ~HIDDEN) != 0)
|
|
118 return 1;
|
|
119 return 0;
|
|
120 }
|
|
121
|
|
122 int
|
|
123 __unorddf2(double a, double b)
|
|
124 {
|
|
125 union double_long dl;
|
|
126
|
|
127 dl.d = a;
|
|
128 if (EXPD(dl) == EXPDMASK
|
|
129 && ((dl.l.upper & MANTDMASK) != 0 || dl.l.lower != 0))
|
|
130 return 1;
|
|
131 dl.d = b;
|
|
132 if (EXPD(dl) == EXPDMASK
|
|
133 && ((dl.l.upper & MANTDMASK) != 0 || dl.l.lower != 0))
|
|
134 return 1;
|
|
135 return 0;
|
|
136 }
|
|
137
|
|
138 /* convert unsigned int to double */
|
|
139 double
|
|
140 __floatunsidf (unsigned long a1)
|
|
141 {
|
|
142 long exp = 32 + EXCESSD;
|
|
143 union double_long dl;
|
|
144
|
|
145 if (!a1)
|
|
146 {
|
|
147 dl.l.upper = dl.l.lower = 0;
|
|
148 return dl.d;
|
|
149 }
|
|
150
|
|
151 while (a1 < 0x2000000L)
|
|
152 {
|
|
153 a1 <<= 4;
|
|
154 exp -= 4;
|
|
155 }
|
|
156
|
|
157 while (a1 < 0x80000000L)
|
|
158 {
|
|
159 a1 <<= 1;
|
|
160 exp--;
|
|
161 }
|
|
162
|
|
163 /* pack up and go home */
|
|
164 dl.l.upper = exp << 20L;
|
|
165 dl.l.upper |= (a1 >> 11L) & ~HIDDEND;
|
|
166 dl.l.lower = a1 << 21L;
|
|
167
|
|
168 return dl.d;
|
|
169 }
|
|
170
|
|
171 /* convert int to double */
|
|
172 double
|
|
173 __floatsidf (long a1)
|
|
174 {
|
|
175 long sign = 0, exp = 31 + EXCESSD;
|
|
176 union double_long dl;
|
|
177
|
|
178 if (!a1)
|
|
179 {
|
|
180 dl.l.upper = dl.l.lower = 0;
|
|
181 return dl.d;
|
|
182 }
|
|
183
|
|
184 if (a1 < 0)
|
|
185 {
|
|
186 sign = SIGNBIT;
|
|
187 a1 = (long)-(unsigned long)a1;
|
|
188 if (a1 < 0)
|
|
189 {
|
|
190 dl.l.upper = SIGNBIT | ((32 + EXCESSD) << 20L);
|
|
191 dl.l.lower = 0;
|
|
192 return dl.d;
|
|
193 }
|
|
194 }
|
|
195
|
|
196 while (a1 < 0x1000000L)
|
|
197 {
|
|
198 a1 <<= 4;
|
|
199 exp -= 4;
|
|
200 }
|
|
201
|
|
202 while (a1 < 0x40000000L)
|
|
203 {
|
|
204 a1 <<= 1;
|
|
205 exp--;
|
|
206 }
|
|
207
|
|
208 /* pack up and go home */
|
|
209 dl.l.upper = sign;
|
|
210 dl.l.upper |= exp << 20L;
|
|
211 dl.l.upper |= (a1 >> 10L) & ~HIDDEND;
|
|
212 dl.l.lower = a1 << 22L;
|
|
213
|
|
214 return dl.d;
|
|
215 }
|
|
216
|
|
217 /* convert unsigned int to float */
|
|
218 float
|
|
219 __floatunsisf (unsigned long l)
|
|
220 {
|
|
221 double foo = __floatunsidf (l);
|
|
222 return foo;
|
|
223 }
|
|
224
|
|
225 /* convert int to float */
|
|
226 float
|
|
227 __floatsisf (long l)
|
|
228 {
|
|
229 double foo = __floatsidf (l);
|
|
230 return foo;
|
|
231 }
|
|
232
|
|
233 /* convert float to double */
|
|
234 double
|
|
235 __extendsfdf2 (float a1)
|
|
236 {
|
|
237 register union float_long fl1;
|
|
238 register union double_long dl;
|
|
239 register long exp;
|
|
240 register long mant;
|
|
241
|
|
242 fl1.f = a1;
|
|
243
|
|
244 dl.l.upper = SIGN (fl1.l);
|
|
245 if ((fl1.l & ~SIGNBIT) == 0)
|
|
246 {
|
|
247 dl.l.lower = 0;
|
|
248 return dl.d;
|
|
249 }
|
|
250
|
|
251 exp = EXP(fl1.l);
|
|
252 mant = MANT (fl1.l) & ~HIDDEN;
|
|
253 if (exp == 0)
|
|
254 {
|
|
255 /* Denormal. */
|
|
256 exp = 1;
|
|
257 while (!(mant & HIDDEN))
|
|
258 {
|
|
259 mant <<= 1;
|
|
260 exp--;
|
|
261 }
|
|
262 mant &= ~HIDDEN;
|
|
263 }
|
|
264 exp = exp - EXCESS + EXCESSD;
|
|
265 dl.l.upper |= exp << 20;
|
|
266 dl.l.upper |= mant >> 3;
|
|
267 dl.l.lower = mant << 29;
|
|
268
|
|
269 return dl.d;
|
|
270 }
|
|
271
|
|
272 /* convert double to float */
|
|
273 float
|
|
274 __truncdfsf2 (double a1)
|
|
275 {
|
|
276 register long exp;
|
|
277 register long mant;
|
|
278 register union float_long fl;
|
|
279 register union double_long dl1;
|
|
280 int sticky;
|
|
281 int shift;
|
|
282
|
|
283 dl1.d = a1;
|
|
284
|
|
285 if ((dl1.l.upper & ~SIGNBIT) == 0 && !dl1.l.lower)
|
|
286 {
|
|
287 fl.l = SIGND(dl1);
|
|
288 return fl.f;
|
|
289 }
|
|
290
|
|
291 exp = EXPD (dl1) - EXCESSD + EXCESS;
|
|
292
|
|
293 sticky = dl1.l.lower & ((1 << 22) - 1);
|
|
294 mant = MANTD (dl1);
|
|
295 /* shift double mantissa 6 bits so we can round */
|
|
296 sticky |= mant & ((1 << 6) - 1);
|
|
297 mant >>= 6;
|
|
298
|
|
299 /* Check for underflow and denormals. */
|
|
300 if (exp <= 0)
|
|
301 {
|
|
302 if (exp < -24)
|
|
303 {
|
|
304 sticky |= mant;
|
|
305 mant = 0;
|
|
306 }
|
|
307 else
|
|
308 {
|
|
309 sticky |= mant & ((1 << (1 - exp)) - 1);
|
|
310 mant >>= 1 - exp;
|
|
311 }
|
|
312 exp = 0;
|
|
313 }
|
|
314
|
|
315 /* now round */
|
|
316 shift = 1;
|
|
317 if ((mant & 1) && (sticky || (mant & 2)))
|
|
318 {
|
|
319 int rounding = exp ? 2 : 1;
|
|
320
|
|
321 mant += 1;
|
|
322
|
|
323 /* did the round overflow? */
|
|
324 if (mant >= (HIDDEN << rounding))
|
|
325 {
|
|
326 exp++;
|
|
327 shift = rounding;
|
|
328 }
|
|
329 }
|
|
330 /* shift down */
|
|
331 mant >>= shift;
|
|
332
|
|
333 mant &= ~HIDDEN;
|
|
334
|
|
335 /* pack up and go home */
|
|
336 fl.l = PACK (SIGND (dl1), exp, mant);
|
|
337 return (fl.f);
|
|
338 }
|
|
339
|
|
340 /* convert double to int */
|
|
341 long
|
|
342 __fixdfsi (double a1)
|
|
343 {
|
|
344 register union double_long dl1;
|
|
345 register long exp;
|
|
346 register long l;
|
|
347
|
|
348 dl1.d = a1;
|
|
349
|
|
350 if (!dl1.l.upper && !dl1.l.lower)
|
|
351 return 0;
|
|
352
|
|
353 exp = EXPD (dl1) - EXCESSD - 31;
|
|
354 l = MANTD (dl1);
|
|
355
|
|
356 if (exp > 0)
|
|
357 {
|
|
358 /* Return largest integer. */
|
|
359 return SIGND (dl1) ? 0x80000000L : 0x7fffffffL;
|
|
360 }
|
|
361
|
|
362 if (exp <= -32)
|
|
363 return 0;
|
|
364
|
|
365 /* shift down until exp = 0 */
|
|
366 if (exp < 0)
|
|
367 l >>= -exp;
|
|
368
|
|
369 return (SIGND (dl1) ? -l : l);
|
|
370 }
|
|
371
|
|
372 /* convert float to int */
|
|
373 long
|
|
374 __fixsfsi (float a1)
|
|
375 {
|
|
376 double foo = a1;
|
|
377 return __fixdfsi (foo);
|
|
378 }
|
|
379
|
|
380 #else /* EXTFLOAT */
|
|
381
|
|
382 /* We do not need these routines for coldfire, as it has no extended
|
|
383 float format. */
|
|
384 #if !defined (__mcoldfire__)
|
|
385
|
|
386 /* Primitive extended precision floating point support.
|
|
387
|
|
388 We assume all numbers are normalized, don't do any rounding, etc. */
|
|
389
|
|
390 /* Prototypes for the above in case we use them. */
|
|
391 double __floatunsidf (unsigned long);
|
|
392 double __floatsidf (long);
|
|
393 float __floatsisf (long);
|
|
394 double __extendsfdf2 (float);
|
|
395 float __truncdfsf2 (double);
|
|
396 long __fixdfsi (double);
|
|
397 long __fixsfsi (float);
|
|
398
|
|
399 int
|
|
400 __unordxf2(long double a, long double b)
|
|
401 {
|
|
402 union long_double_long ldl;
|
|
403
|
|
404 ldl.ld = a;
|
|
405 if (EXPX(ldl) == EXPXMASK
|
|
406 && ((ldl.l.middle & MANTXMASK) != 0 || ldl.l.lower != 0))
|
|
407 return 1;
|
|
408 ldl.ld = b;
|
|
409 if (EXPX(ldl) == EXPXMASK
|
|
410 && ((ldl.l.middle & MANTXMASK) != 0 || ldl.l.lower != 0))
|
|
411 return 1;
|
|
412 return 0;
|
|
413 }
|
|
414
|
|
415 /* convert double to long double */
|
|
416 long double
|
|
417 __extenddfxf2 (double d)
|
|
418 {
|
|
419 register union double_long dl;
|
|
420 register union long_double_long ldl;
|
|
421 register long exp;
|
|
422
|
|
423 dl.d = d;
|
|
424 /*printf ("dfxf in: %g\n", d);*/
|
|
425
|
|
426 ldl.l.upper = SIGND (dl);
|
|
427 if ((dl.l.upper & ~SIGNBIT) == 0 && !dl.l.lower)
|
|
428 {
|
|
429 ldl.l.middle = 0;
|
|
430 ldl.l.lower = 0;
|
|
431 return ldl.ld;
|
|
432 }
|
|
433
|
|
434 exp = EXPD (dl) - EXCESSD + EXCESSX;
|
|
435 ldl.l.upper |= exp << 16;
|
|
436 ldl.l.middle = HIDDENX;
|
|
437 /* 31-20: # mantissa bits in ldl.l.middle - # mantissa bits in dl.l.upper */
|
|
438 ldl.l.middle |= (dl.l.upper & MANTDMASK) << (31 - 20);
|
|
439 /* 1+20: explicit-integer-bit + # mantissa bits in dl.l.upper */
|
|
440 ldl.l.middle |= dl.l.lower >> (1 + 20);
|
|
441 /* 32 - 21: # bits of dl.l.lower in ldl.l.middle */
|
|
442 ldl.l.lower = dl.l.lower << (32 - 21);
|
|
443
|
|
444 /*printf ("dfxf out: %s\n", dumpxf (ldl.ld));*/
|
|
445 return ldl.ld;
|
|
446 }
|
|
447
|
|
448 /* convert long double to double */
|
|
449 double
|
|
450 __truncxfdf2 (long double ld)
|
|
451 {
|
|
452 register long exp;
|
|
453 register union double_long dl;
|
|
454 register union long_double_long ldl;
|
|
455
|
|
456 ldl.ld = ld;
|
|
457 /*printf ("xfdf in: %s\n", dumpxf (ld));*/
|
|
458
|
|
459 dl.l.upper = SIGNX (ldl);
|
|
460 if ((ldl.l.upper & ~SIGNBIT) == 0 && !ldl.l.middle && !ldl.l.lower)
|
|
461 {
|
|
462 dl.l.lower = 0;
|
|
463 return dl.d;
|
|
464 }
|
|
465
|
|
466 exp = EXPX (ldl) - EXCESSX + EXCESSD;
|
|
467 /* ??? quick and dirty: keep `exp' sane */
|
|
468 if (exp >= EXPDMASK)
|
|
469 exp = EXPDMASK - 1;
|
|
470 dl.l.upper |= exp << (32 - (EXPDBITS + 1));
|
|
471 /* +1-1: add one for sign bit, but take one off for explicit-integer-bit */
|
|
472 dl.l.upper |= (ldl.l.middle & MANTXMASK) >> (EXPDBITS + 1 - 1);
|
|
473 dl.l.lower = (ldl.l.middle & MANTXMASK) << (32 - (EXPDBITS + 1 - 1));
|
|
474 dl.l.lower |= ldl.l.lower >> (EXPDBITS + 1 - 1);
|
|
475
|
|
476 /*printf ("xfdf out: %g\n", dl.d);*/
|
|
477 return dl.d;
|
|
478 }
|
|
479
|
|
480 /* convert a float to a long double */
|
|
481 long double
|
|
482 __extendsfxf2 (float f)
|
|
483 {
|
|
484 long double foo = __extenddfxf2 (__extendsfdf2 (f));
|
|
485 return foo;
|
|
486 }
|
|
487
|
|
488 /* convert a long double to a float */
|
|
489 float
|
|
490 __truncxfsf2 (long double ld)
|
|
491 {
|
|
492 float foo = __truncdfsf2 (__truncxfdf2 (ld));
|
|
493 return foo;
|
|
494 }
|
|
495
|
|
496 /* convert an int to a long double */
|
|
497 long double
|
|
498 __floatsixf (long l)
|
|
499 {
|
|
500 double foo = __floatsidf (l);
|
|
501 return foo;
|
|
502 }
|
|
503
|
|
504 /* convert an unsigned int to a long double */
|
|
505 long double
|
|
506 __floatunsixf (unsigned long l)
|
|
507 {
|
|
508 double foo = __floatunsidf (l);
|
|
509 return foo;
|
|
510 }
|
|
511
|
|
512 /* convert a long double to an int */
|
|
513 long
|
|
514 __fixxfsi (long double ld)
|
|
515 {
|
|
516 long foo = __fixdfsi ((double) ld);
|
|
517 return foo;
|
|
518 }
|
|
519
|
|
520 /* The remaining provide crude math support by working in double precision. */
|
|
521
|
|
522 long double
|
|
523 __addxf3 (long double x1, long double x2)
|
|
524 {
|
|
525 return (double) x1 + (double) x2;
|
|
526 }
|
|
527
|
|
528 long double
|
|
529 __subxf3 (long double x1, long double x2)
|
|
530 {
|
|
531 return (double) x1 - (double) x2;
|
|
532 }
|
|
533
|
|
534 long double
|
|
535 __mulxf3 (long double x1, long double x2)
|
|
536 {
|
|
537 return (double) x1 * (double) x2;
|
|
538 }
|
|
539
|
|
540 long double
|
|
541 __divxf3 (long double x1, long double x2)
|
|
542 {
|
|
543 return (double) x1 / (double) x2;
|
|
544 }
|
|
545
|
|
546 long double
|
|
547 __negxf2 (long double x1)
|
|
548 {
|
|
549 return - (double) x1;
|
|
550 }
|
|
551
|
|
552 long
|
|
553 __cmpxf2 (long double x1, long double x2)
|
|
554 {
|
|
555 return __cmpdf2 ((double) x1, (double) x2);
|
|
556 }
|
|
557
|
|
558 long
|
|
559 __eqxf2 (long double x1, long double x2)
|
|
560 {
|
|
561 return __cmpdf2 ((double) x1, (double) x2);
|
|
562 }
|
|
563
|
|
564 long
|
|
565 __nexf2 (long double x1, long double x2)
|
|
566 {
|
|
567 return __cmpdf2 ((double) x1, (double) x2);
|
|
568 }
|
|
569
|
|
570 long
|
|
571 __ltxf2 (long double x1, long double x2)
|
|
572 {
|
|
573 return __cmpdf2 ((double) x1, (double) x2);
|
|
574 }
|
|
575
|
|
576 long
|
|
577 __lexf2 (long double x1, long double x2)
|
|
578 {
|
|
579 return __cmpdf2 ((double) x1, (double) x2);
|
|
580 }
|
|
581
|
|
582 long
|
|
583 __gtxf2 (long double x1, long double x2)
|
|
584 {
|
|
585 return __cmpdf2 ((double) x1, (double) x2);
|
|
586 }
|
|
587
|
|
588 long
|
|
589 __gexf2 (long double x1, long double x2)
|
|
590 {
|
|
591 return __cmpdf2 ((double) x1, (double) x2);
|
|
592 }
|
|
593
|
|
594 #endif /* !__mcoldfire__ */
|
|
595 #endif /* EXTFLOAT */
|