Commit | Line | Data |
---|---|---|
9221ab2b C |
1 | /* |
2 | ** libgcc support for software floating point. | |
3 | ** Copyright (C) 1991 by Pipeline Associates, Inc. All rights reserved. | |
4 | ** Permission is granted to do *anything* you want with this file, | |
5 | ** commercial or otherwise, provided this message remains intact. So there! | |
6 | ** I would appreciate receiving any updates/patches/changes that anyone | |
7 | ** makes, and am willing to be the repository for said changes (am I | |
8 | ** making a big mistake?). | |
9 | ||
10 | Warning! Only single-precision is actually implemented. This file | |
11 | won't really be much use until double-precision is supported. | |
12 | ||
13 | However, once that is done, this file might eventually become a | |
14 | replacement for libgcc1.c. It might also make possible | |
15 | cross-compilation for an IEEE target machine from a non-IEEE | |
16 | host such as a VAX. | |
17 | ||
18 | If you'd like to work on completing this, please talk to rms@gnu.ai.mit.edu. | |
19 | ||
20 | ||
21 | ** | |
22 | ** Pat Wood | |
23 | ** Pipeline Associates, Inc. | |
24 | ** pipeline!phw@motown.com or | |
25 | ** sun!pipeline!phw or | |
26 | ** uunet!motown!pipeline!phw | |
27 | ** | |
28 | ** 05/01/91 -- V1.0 -- first release to gcc mailing lists | |
29 | ** 05/04/91 -- V1.1 -- added float and double prototypes and return values | |
30 | ** -- fixed problems with adding and subtracting zero | |
31 | ** -- fixed rounding in truncdfsf2 | |
32 | ** -- fixed SWAP define and tested on 386 | |
33 | */ | |
34 | ||
35 | /* | |
36 | ** The following are routines that replace the libgcc soft floating point | |
37 | ** routines that are called automatically when -msoft-float is selected. | |
38 | ** The support single and double precision IEEE format, with provisions | |
39 | ** for byte-swapped machines (tested on 386). Some of the double-precision | |
40 | ** routines work at full precision, but most of the hard ones simply punt | |
41 | ** and call the single precision routines, producing a loss of accuracy. | |
42 | ** long long support is not assumed or included. | |
43 | ** Overall accuracy is close to IEEE (actually 68882) for single-precision | |
44 | ** arithmetic. I think there may still be a 1 in 1000 chance of a bit | |
45 | ** being rounded the wrong way during a multiply. I'm not fussy enough to | |
46 | ** bother with it, but if anyone is, knock yourself out. | |
47 | ** | |
48 | ** Efficiency has only been addressed where it was obvious that something | |
49 | ** would make a big difference. Anyone who wants to do this right for | |
50 | ** best speed should go in and rewrite in assembler. | |
51 | ** | |
52 | ** I have tested this only on a 68030 workstation and 386/ix integrated | |
53 | ** in with -msoft-float. | |
54 | */ | |
55 | ||
56 | /* the following deal with IEEE single-precision numbers */ | |
57 | #define EXCESS 126 | |
58 | #define SIGNBIT 0x80000000 | |
59 | #define HIDDEN (1 << 23) | |
60 | #define SIGN(fp) ((fp) & SIGNBIT) | |
61 | #define EXP(fp) (((fp) >> 23) & 0xFF) | |
62 | #define MANT(fp) (((fp) & 0x7FFFFF) | HIDDEN) | |
63 | #define PACK(s,e,m) ((s) | ((e) << 23) | (m)) | |
64 | ||
65 | /* the following deal with IEEE double-precision numbers */ | |
66 | #define EXCESSD 1022 | |
67 | #define HIDDEND (1 << 20) | |
68 | #define EXPD(fp) (((fp.l.upper) >> 20) & 0x7FF) | |
69 | #define SIGND(fp) ((fp.l.upper) & SIGNBIT) | |
70 | #define MANTD(fp) (((((fp.l.upper) & 0xFFFFF) | HIDDEND) << 10) | \ | |
71 | (fp.l.lower >> 22)) | |
72 | ||
73 | /* define SWAP for 386/960 reverse-byte-order brain-damaged CPUs */ | |
74 | union double_long | |
75 | { | |
76 | double d; | |
77 | #ifdef SWAP | |
78 | struct { | |
79 | unsigned long lower; | |
80 | long upper; | |
81 | } l; | |
82 | #else | |
83 | struct { | |
84 | long upper; | |
85 | unsigned long lower; | |
86 | } l; | |
87 | #endif | |
88 | }; | |
89 | ||
90 | union float_long | |
91 | { | |
92 | float f; | |
93 | long l; | |
94 | }; | |
95 | ||
96 | /* add two floats */ | |
97 | float | |
98 | __addsf3 (float a1, float a2) | |
99 | { | |
100 | register long mant1, mant2; | |
101 | register union float_long fl1, fl2; | |
102 | register int exp1, exp2; | |
103 | int sign = 0; | |
104 | ||
105 | fl1.f = a1; | |
106 | fl2.f = a2; | |
107 | ||
108 | /* check for zero args */ | |
109 | if (!fl1.l) | |
110 | return (fl2.f); | |
111 | if (!fl2.l) | |
112 | return (fl1.f); | |
113 | ||
114 | exp1 = EXP (fl1.l); | |
115 | exp2 = EXP (fl2.l); | |
116 | ||
117 | if (exp1 > exp2 + 25) | |
118 | return (fl1.l); | |
119 | if (exp2 > exp1 + 25) | |
120 | return (fl2.l); | |
121 | ||
122 | /* do everything in excess precision so's we can round later */ | |
123 | mant1 = MANT (fl1.l) << 6; | |
124 | mant2 = MANT (fl2.l) << 6; | |
125 | ||
126 | if (SIGN (fl1.l)) | |
127 | mant1 = -mant1; | |
128 | if (SIGN (fl2.l)) | |
129 | mant2 = -mant2; | |
130 | ||
131 | if (exp1 > exp2) | |
132 | { | |
133 | mant2 >>= exp1 - exp2; | |
134 | } | |
135 | else | |
136 | { | |
137 | mant1 >>= exp2 - exp1; | |
138 | exp1 = exp2; | |
139 | } | |
140 | mant1 += mant2; | |
141 | ||
142 | if (mant1 < 0) | |
143 | { | |
144 | mant1 = -mant1; | |
145 | sign = SIGNBIT; | |
146 | } | |
147 | else if (!mant1) | |
148 | return (0); | |
149 | ||
150 | /* normalize up */ | |
151 | while (!(mant1 & 0xE0000000)) | |
152 | { | |
153 | mant1 <<= 1; | |
154 | exp1--; | |
155 | } | |
156 | ||
157 | /* normalize down? */ | |
158 | if (mant1 & (1 << 30)) | |
159 | { | |
160 | mant1 >>= 1; | |
161 | exp1++; | |
162 | } | |
163 | ||
164 | /* round to even */ | |
165 | mant1 += (mant1 & 0x40) ? 0x20 : 0x1F; | |
166 | ||
167 | /* normalize down? */ | |
168 | if (mant1 & (1 << 30)) | |
169 | { | |
170 | mant1 >>= 1; | |
171 | exp1++; | |
172 | } | |
173 | ||
174 | /* lose extra precision */ | |
175 | mant1 >>= 6; | |
176 | ||
177 | /* turn off hidden bit */ | |
178 | mant1 &= ~HIDDEN; | |
179 | ||
180 | /* pack up and go home */ | |
181 | fl1.l = PACK (sign, exp1, mant1); | |
182 | return (fl1.f); | |
183 | } | |
184 | ||
185 | /* subtract two floats */ | |
186 | float | |
187 | __subsf3 (float a1, float a2) | |
188 | { | |
189 | register union float_long fl1, fl2; | |
190 | ||
191 | fl1.f = a1; | |
192 | fl2.f = a2; | |
193 | ||
194 | /* check for zero args */ | |
195 | if (!fl2.l) | |
196 | return (fl1.f); | |
197 | if (!fl1.l) | |
198 | return (-fl2.f); | |
199 | ||
200 | /* twiddle sign bit and add */ | |
201 | fl2.l ^= SIGNBIT; | |
202 | return __addsf3 (a1, fl2.f); | |
203 | } | |
204 | ||
205 | /* compare two floats */ | |
206 | long | |
207 | __cmpsf2 (float a1, float a2) | |
208 | { | |
209 | register union float_long fl1, fl2; | |
210 | ||
211 | fl1.f = a1; | |
212 | fl2.f = a2; | |
213 | ||
214 | if (SIGN (fl1.l) && SIGN (fl2.l)) | |
215 | { | |
216 | fl1.l ^= SIGNBIT; | |
217 | fl2.l ^= SIGNBIT; | |
218 | } | |
219 | if (fl1.l < fl2.l) | |
220 | return (-1); | |
221 | if (fl1.l > fl2.l) | |
222 | return (1); | |
223 | return (0); | |
224 | } | |
225 | ||
226 | /* multiply two floats */ | |
227 | float | |
228 | __mulsf3 (float a1, float a2) | |
229 | { | |
230 | register union float_long fl1, fl2; | |
231 | register unsigned long result; | |
232 | register int exp; | |
233 | int sign; | |
234 | ||
235 | fl1.f = a1; | |
236 | fl2.f = a2; | |
237 | ||
238 | if (!fl1.l || !fl2.l) | |
239 | return (0); | |
240 | ||
241 | /* compute sign and exponent */ | |
242 | sign = SIGN (fl1.l) ^ SIGN (fl2.l); | |
243 | exp = EXP (fl1.l) - EXCESS; | |
244 | exp += EXP (fl2.l); | |
245 | ||
246 | fl1.l = MANT (fl1.l); | |
247 | fl2.l = MANT (fl2.l); | |
248 | ||
249 | /* the multiply is done as one 16x16 multiply and two 16x8 multiples */ | |
250 | result = (fl1.l >> 8) * (fl2.l >> 8); | |
251 | result += ((fl1.l & 0xFF) * (fl2.l >> 8)) >> 8; | |
252 | result += ((fl2.l & 0xFF) * (fl1.l >> 8)) >> 8; | |
253 | ||
254 | if (result & 0x80000000) | |
255 | { | |
256 | /* round */ | |
257 | result += 0x80; | |
258 | result >>= 8; | |
259 | } | |
260 | else | |
261 | { | |
262 | /* round */ | |
263 | result += 0x40; | |
264 | result >>= 7; | |
265 | exp--; | |
266 | } | |
267 | ||
268 | result &= ~HIDDEN; | |
269 | ||
270 | /* pack up and go home */ | |
271 | fl1.l = PACK (sign, exp, result); | |
272 | return (fl1.f); | |
273 | } | |
274 | ||
275 | /* divide two floats */ | |
276 | float | |
277 | __divsf3 (float a1, float a2) | |
278 | { | |
279 | register union float_long fl1, fl2; | |
280 | register int result; | |
281 | register int mask; | |
282 | register int exp, sign; | |
283 | ||
284 | fl1.f = a1; | |
285 | fl2.f = a2; | |
286 | ||
287 | /* subtract exponents */ | |
288 | exp = EXP (fl1.l) - EXP (fl2.l) + EXCESS; | |
289 | ||
290 | /* compute sign */ | |
291 | sign = SIGN (fl1.l) ^ SIGN (fl2.l); | |
292 | ||
293 | /* divide by zero??? */ | |
294 | if (!fl2.l) | |
295 | /* return NaN or -NaN */ | |
296 | return (sign ? 0xFFFFFFFF : 0x7FFFFFFF); | |
297 | ||
298 | /* numerator zero??? */ | |
299 | if (!fl1.l) | |
300 | return (0); | |
301 | ||
302 | /* now get mantissas */ | |
303 | fl1.l = MANT (fl1.l); | |
304 | fl2.l = MANT (fl2.l); | |
305 | ||
306 | /* this assures we have 25 bits of precision in the end */ | |
307 | if (fl1.l < fl2.l) | |
308 | { | |
309 | fl1.l <<= 1; | |
310 | exp--; | |
311 | } | |
312 | ||
313 | /* now we perform repeated subtraction of fl2.l from fl1.l */ | |
314 | mask = 0x1000000; | |
315 | result = 0; | |
316 | while (mask) | |
317 | { | |
318 | if (fl1.l >= fl2.l) | |
319 | { | |
320 | result |= mask; | |
321 | fl1.l -= fl2.l; | |
322 | } | |
323 | fl1.l <<= 1; | |
324 | mask >>= 1; | |
325 | } | |
326 | ||
327 | /* round */ | |
328 | result += 1; | |
329 | ||
330 | /* normalize down */ | |
331 | exp++; | |
332 | result >>= 1; | |
333 | ||
334 | result &= ~HIDDEN; | |
335 | ||
336 | /* pack up and go home */ | |
337 | fl1.l = PACK (sign, exp, result); | |
338 | return (fl1.f); | |
339 | } | |
340 | ||
341 | /* convert int to double */ | |
342 | double | |
343 | __floatsidf (register long a1) | |
344 | { | |
345 | register int sign = 0, exp = 31 + EXCESSD; | |
346 | union double_long dl; | |
347 | ||
348 | if (!a1) | |
349 | { | |
350 | dl.l.upper = dl.l.lower = 0; | |
351 | return (dl.d); | |
352 | } | |
353 | ||
354 | if (a1 < 0) | |
355 | { | |
356 | sign = SIGNBIT; | |
357 | a1 = -a1; | |
358 | } | |
359 | ||
360 | while (a1 < 0x1000000) | |
361 | { | |
362 | a1 <<= 4; | |
363 | exp -= 4; | |
364 | } | |
365 | ||
366 | while (a1 < 0x40000000) | |
367 | { | |
368 | a1 <<= 1; | |
369 | exp--; | |
370 | } | |
371 | ||
372 | /* pack up and go home */ | |
373 | dl.l.upper = sign; | |
374 | dl.l.upper |= exp << 20; | |
375 | dl.l.upper |= (a1 >> 10) & ~HIDDEND; | |
376 | dl.l.lower = a1 << 22; | |
377 | ||
378 | return (dl.d); | |
379 | } | |
380 | ||
381 | /* negate a float */ | |
382 | float | |
383 | __negsf2 (float a1) | |
384 | { | |
385 | register union float_long fl1; | |
386 | ||
387 | fl1.f = a1; | |
388 | if (!fl1.l) | |
389 | return (0); | |
390 | ||
391 | fl1.l ^= SIGNBIT; | |
392 | return (fl1.f); | |
393 | } | |
394 | ||
395 | /* negate a double */ | |
396 | double | |
397 | __negdf2 (double a1) | |
398 | { | |
399 | register union double_long dl1; | |
400 | ||
401 | dl1.d = a1; | |
402 | ||
403 | if (!dl1.l.upper && !dl1.l.lower) | |
404 | return (dl1.d); | |
405 | ||
406 | dl1.l.upper ^= SIGNBIT; | |
407 | return (dl1.d); | |
408 | } | |
409 | ||
410 | /* convert float to double */ | |
411 | double | |
412 | __extendsfdf2 (float a1) | |
413 | { | |
414 | register union float_long fl1; | |
415 | register union double_long dl; | |
416 | register int exp; | |
417 | ||
418 | fl1.f = a1; | |
419 | ||
420 | if (!fl1.l) | |
421 | { | |
422 | dl.l.upper = dl.l.lower = 0; | |
423 | return (dl.d); | |
424 | } | |
425 | ||
426 | dl.l.upper = SIGN (fl1.l); | |
427 | exp = EXP (fl1.l) - EXCESS + EXCESSD; | |
428 | dl.l.upper |= exp << 20; | |
429 | dl.l.upper |= (MANT (fl1.l) & ~HIDDEN) >> 3; | |
430 | dl.l.lower = MANT (fl1.l) << 29; | |
431 | ||
432 | return (dl.d); | |
433 | } | |
434 | ||
435 | /* convert double to float */ | |
436 | float | |
437 | __truncdfsf2 (double a1) | |
438 | { | |
439 | register int exp; | |
440 | register long mant; | |
441 | register union float_long fl; | |
442 | register union double_long dl1; | |
443 | ||
444 | dl1.d = a1; | |
445 | ||
446 | if (!dl1.l.upper && !dl1.l.lower) | |
447 | return (0); | |
448 | ||
449 | exp = EXPD (dl1) - EXCESSD + EXCESS; | |
450 | ||
451 | /* shift double mantissa 6 bits so we can round */ | |
452 | mant = MANTD (dl1) >> 6; | |
453 | ||
454 | /* now round and shift down */ | |
455 | mant += 1; | |
456 | mant >>= 1; | |
457 | ||
458 | /* did the round overflow? */ | |
459 | if (mant & 0xFF000000) | |
460 | { | |
461 | mant >>= 1; | |
462 | exp++; | |
463 | } | |
464 | ||
465 | mant &= ~HIDDEN; | |
466 | ||
467 | /* pack up and go home */ | |
468 | fl.l = PACK (SIGND (dl1), exp, mant); | |
469 | return (fl.f); | |
470 | } | |
471 | ||
472 | /* compare two doubles */ | |
473 | long | |
474 | __cmpdf2 (double a1, double a2) | |
475 | { | |
476 | register union double_long dl1, dl2; | |
477 | ||
478 | dl1.d = a1; | |
479 | dl2.d = a2; | |
480 | ||
481 | if (SIGND (dl1) && SIGND (dl2)) | |
482 | { | |
483 | dl1.l.upper ^= SIGNBIT; | |
484 | dl2.l.upper ^= SIGNBIT; | |
485 | } | |
486 | if (dl1.l.upper < dl2.l.upper) | |
487 | return (-1); | |
488 | if (dl1.l.upper > dl2.l.upper) | |
489 | return (1); | |
490 | if (dl1.l.lower < dl2.l.lower) | |
491 | return (-1); | |
492 | if (dl1.l.lower > dl2.l.lower) | |
493 | return (1); | |
494 | return (0); | |
495 | } | |
496 | ||
497 | /* convert double to int */ | |
498 | long | |
499 | __fixdfsi (double a1) | |
500 | { | |
501 | register union double_long dl1; | |
502 | register int exp; | |
503 | register long l; | |
504 | ||
505 | dl1.d = a1; | |
506 | ||
507 | if (!dl1.l.upper && !dl1.l.lower) | |
508 | return (0); | |
509 | ||
510 | exp = EXPD (dl1) - EXCESSD - 31; | |
511 | l = MANTD (dl1); | |
512 | ||
513 | if (exp > 0) | |
514 | return (0x7FFFFFFF | SIGND (dl1)); /* largest integer */ | |
515 | ||
516 | /* shift down until exp = 0 or l = 0 */ | |
517 | if (exp < 0 && exp > -32 && l) | |
518 | l >>= -exp; | |
519 | else | |
520 | return (0); | |
521 | ||
522 | return (SIGND (dl1) ? -l : l); | |
523 | } | |
524 | ||
525 | /* convert double to unsigned int */ | |
526 | unsigned | |
527 | long __fixunsdfsi (double a1) | |
528 | { | |
529 | register union double_long dl1; | |
530 | register int exp; | |
531 | register unsigned long l; | |
532 | ||
533 | dl1.d = a1; | |
534 | ||
535 | if (!dl1.l.upper && !dl1.l.lower) | |
536 | return (0); | |
537 | ||
538 | exp = EXPD (dl1) - EXCESSD - 32; | |
539 | l = (((((dl1.l.upper) & 0xFFFFF) | HIDDEND) << 11) | (dl1.l.lower >> 21)); | |
540 | ||
541 | if (exp > 0) | |
542 | return (0xFFFFFFFF); /* largest integer */ | |
543 | ||
544 | /* shift down until exp = 0 or l = 0 */ | |
545 | if (exp < 0 && exp > -32 && l) | |
546 | l >>= -exp; | |
547 | else | |
548 | return (0); | |
549 | ||
550 | return (l); | |
551 | } | |
552 | ||
553 | /* For now, the hard double-precision routines simply | |
554 | punt and do it in single */ | |
555 | /* addtwo doubles */ | |
556 | double | |
557 | __adddf3 (double a1, double a2) | |
558 | { | |
559 | return ((float) a1 + (float) a2); | |
560 | } | |
561 | ||
562 | /* subtract two doubles */ | |
563 | double | |
564 | __subdf3 (double a1, double a2) | |
565 | { | |
566 | return ((float) a1 - (float) a2); | |
567 | } | |
568 | ||
569 | /* multiply two doubles */ | |
570 | double | |
571 | __muldf3 (double a1, double a2) | |
572 | { | |
573 | return ((float) a1 * (float) a2); | |
574 | } | |
575 | ||
576 | /* divide two doubles */ | |
577 | double | |
578 | __divdf3 (double a1, double a2) | |
579 | { | |
580 | return ((float) a1 / (float) a2); | |
581 | } |