![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to floatlib.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Apple XNU] / GNUtools / cc|
Return to floatlib.c CVS log | Up to [Apple XNU] / GNUtools / cc |
1.1 root 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 [email protected].
19:
20:
21: **
22: ** Pat Wood
23: ** Pipeline Associates, Inc.
24: ** [email protected] 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: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.