|
|
1.1 root 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: ** [email protected] 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 1022
63: #define HIDDEND (1L << 20L)
64: #define EXPDBITS 11
65: #define EXPDMASK 0x7FF
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 0xFFFFF /* mask of upper part */
71:
72: /* the following deal with IEEE extended-precision numbers */
73: #define EXCESSX 16382
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 0x7FFFFFFF /* 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: /* convert int to double */
109: double
110: __floatsidf (int a1)
111: {
112: long sign = 0, exp = 31 + EXCESSD;
113: union double_long dl;
114:
115: if (!a1)
116: {
117: dl.l.upper = dl.l.lower = 0;
118: return dl.d;
119: }
120:
121: if (a1 < 0)
122: {
123: sign = SIGNBIT;
124: a1 = -a1;
125: if (a1 < 0)
126: {
127: dl.l.upper = SIGNBIT | ((32 + EXCESSD) << 20L);
128: dl.l.lower = 0;
129: return dl.d;
130: }
131: }
132:
133: while (a1 < 0x1000000)
134: {
135: a1 <<= 4;
136: exp -= 4;
137: }
138:
139: while (a1 < 0x40000000)
140: {
141: a1 <<= 1;
142: exp--;
143: }
144:
145: /* pack up and go home */
146: dl.l.upper = sign;
147: dl.l.upper |= exp << 20L;
148: dl.l.upper |= (a1 >> 10L) & ~HIDDEND;
149: dl.l.lower = a1 << 22L;
150:
151: return dl.d;
152: }
153:
154: /* convert int to float */
155: float
156: __floatsisf (int l)
157: {
158: double foo = __floatsidf (l);
159: return foo;
160: }
161:
162: /* convert float to double */
163: double
164: __extendsfdf2 (float a1)
165: {
166: register union float_long fl1;
167: register union double_long dl;
168: register long exp;
169:
170: fl1.f = a1;
171:
172: if (!fl1.l)
173: {
174: dl.l.upper = dl.l.lower = 0;
175: return dl.d;
176: }
177:
178: dl.l.upper = SIGN (fl1.l);
179: exp = EXP (fl1.l) - EXCESS + EXCESSD;
180: dl.l.upper |= exp << 20;
181: dl.l.upper |= (MANT (fl1.l) & ~HIDDEN) >> 3;
182: dl.l.lower = MANT (fl1.l) << 29;
183:
184: return dl.d;
185: }
186:
187: /* convert double to float */
188: float
189: __truncdfsf2 (double a1)
190: {
191: register long exp;
192: register long mant;
193: register union float_long fl;
194: register union double_long dl1;
195:
196: dl1.d = a1;
197:
198: if (!dl1.l.upper && !dl1.l.lower)
199: return 0;
200:
201: exp = EXPD (dl1) - EXCESSD + EXCESS;
202:
203: /* shift double mantissa 6 bits so we can round */
204: mant = MANTD (dl1) >> 6;
205:
206: /* now round and shift down */
207: mant += 1;
208: mant >>= 1;
209:
210: /* did the round overflow? */
211: if (mant & 0xFF000000)
212: {
213: mant >>= 1;
214: exp++;
215: }
216:
217: mant &= ~HIDDEN;
218:
219: /* pack up and go home */
220: fl.l = PACK (SIGND (dl1), exp, mant);
221: return (fl.f);
222: }
223:
224: /* convert double to int */
225: int
226: __fixdfsi (double a1)
227: {
228: register union double_long dl1;
229: register long exp;
230: register long l;
231:
232: dl1.d = a1;
233:
234: if (!dl1.l.upper && !dl1.l.lower)
235: return 0;
236:
237: exp = EXPD (dl1) - EXCESSD - 31;
238: l = MANTD (dl1);
239:
240: if (exp > 0)
241: {
242: /* Return largest integer. */
243: return SIGND (dl1) ? 0x80000000 : 0x7fffffff;
244: }
245:
246: if (exp <= -32)
247: return 0;
248:
249: /* shift down until exp = 0 */
250: if (exp < 0)
251: l >>= -exp;
252:
253: return (SIGND (dl1) ? -l : l);
254: }
255:
256: /* convert float to int */
257: int
258: __fixsfsi (float a1)
259: {
260: double foo = a1;
261: return __fixdfsi (foo);
262: }
263:
264: #else /* EXTFLOAT */
265:
266: /* Primitive extended precision floating point support.
267:
268: We assume all numbers are normalized, don't do any rounding, etc. */
269:
270: /* Prototypes for the above in case we use them. */
271: double __floatsidf (int);
272: float __floatsisf (int);
273: double __extendsfdf2 (float);
274: float __truncdfsf2 (double);
275: int __fixdfsi (double);
276: int __fixsfsi (float);
277:
278: /* convert double to long double */
279: long double
280: __extenddfxf2 (double d)
281: {
282: register union double_long dl;
283: register union long_double_long ldl;
284: register long exp;
285:
286: dl.d = d;
287: /*printf ("dfxf in: %g\n", d);*/
288:
289: if (!dl.l.upper && !dl.l.lower)
290: return 0;
291:
292: ldl.l.upper = SIGND (dl);
293: exp = EXPD (dl) - EXCESSD + EXCESSX;
294: ldl.l.upper |= exp << 16;
295: ldl.l.middle = HIDDENX;
296: /* 31-20: # mantissa bits in ldl.l.middle - # mantissa bits in dl.l.upper */
297: ldl.l.middle |= (dl.l.upper & MANTDMASK) << (31 - 20);
298: /* 1+20: explicit-integer-bit + # mantissa bits in dl.l.upper */
299: ldl.l.middle |= dl.l.lower >> (1 + 20);
300: /* 32 - 21: # bits of dl.l.lower in ldl.l.middle */
301: ldl.l.lower = dl.l.lower << (32 - 21);
302:
303: /*printf ("dfxf out: %s\n", dumpxf (ldl.ld));*/
304: return ldl.ld;
305: }
306:
307: /* convert long double to double */
308: double
309: __truncxfdf2 (long double ld)
310: {
311: register long exp;
312: register union double_long dl;
313: register union long_double_long ldl;
314:
315: ldl.ld = ld;
316: /*printf ("xfdf in: %s\n", dumpxf (ld));*/
317:
318: if (!ldl.l.upper && !ldl.l.middle && !ldl.l.lower)
319: return 0;
320:
321: exp = EXPX (ldl) - EXCESSX + EXCESSD;
322: /* ??? quick and dirty: keep `exp' sane */
323: if (exp >= EXPDMASK)
324: exp = EXPDMASK - 1;
325: dl.l.upper = SIGNX (ldl);
326: dl.l.upper |= exp << (32 - (EXPDBITS + 1));
327: /* +1-1: add one for sign bit, but take one off for explicit-integer-bit */
328: dl.l.upper |= (ldl.l.middle & MANTXMASK) >> (EXPDBITS + 1 - 1);
329: dl.l.lower = (ldl.l.middle & MANTXMASK) << (32 - (EXPDBITS + 1 - 1));
330: dl.l.lower |= ldl.l.lower >> (EXPDBITS + 1 - 1);
331:
332: /*printf ("xfdf out: %g\n", dl.d);*/
333: return dl.d;
334: }
335:
336: /* convert a float to a long double */
337: long double
338: __extendsfxf2 (float f)
339: {
340: long double foo = __extenddfxf2 (__extendsfdf2 (f));
341: return foo;
342: }
343:
344: /* convert a long double to a float */
345: float
346: __truncxfsf2 (long double ld)
347: {
348: float foo = __truncdfsf2 (__truncxfdf2 (ld));
349: return foo;
350: }
351:
352: /* convert an int to a long double */
353: long double
354: __floatsixf (int l)
355: {
356: double foo = __floatsidf (l);
357: return foo;
358: }
359:
360: /* convert a long double to an int */
361: int
362: __fixxfsi (long double ld)
363: {
364: int foo = __fixdfsi ((double) ld);
365: return foo;
366: }
367:
368: /* The remaining provide crude math support by working in double precision. */
369:
370: long double
371: __addxf3 (long double x1, long double x2)
372: {
373: return (double) x1 + (double) x2;
374: }
375:
376: long double
377: __subxf3 (long double x1, long double x2)
378: {
379: return (double) x1 - (double) x2;
380: }
381:
382: long double
383: __mulxf3 (long double x1, long double x2)
384: {
385: return (double) x1 * (double) x2;
386: }
387:
388: long double
389: __divxf3 (long double x1, long double x2)
390: {
391: return (double) x1 / (double) x2;
392: }
393:
394: long double
395: __negxf2 (long double x1)
396: {
397: return - (double) x1;
398: }
399:
400: long
401: __cmpxf2 (long double x1, long double x2)
402: {
403: return __cmpdf2 ((double) x1, (double) x2);
404: }
405:
406: long
407: __eqxf2 (long double x1, long double x2)
408: {
409: return __cmpdf2 ((double) x1, (double) x2);
410: }
411:
412: long
413: __nexf2 (long double x1, long double x2)
414: {
415: return __cmpdf2 ((double) x1, (double) x2);
416: }
417:
418: long
419: __ltxf2 (long double x1, long double x2)
420: {
421: return __cmpdf2 ((double) x1, (double) x2);
422: }
423:
424: long
425: __lexf2 (long double x1, long double x2)
426: {
427: return __cmpdf2 ((double) x1, (double) x2);
428: }
429:
430: long
431: __gtxf2 (long double x1, long double x2)
432: {
433: return __cmpdf2 ((double) x1, (double) x2);
434: }
435:
436: long
437: __gexf2 (long double x1, long double x2)
438: {
439: return __cmpdf2 ((double) x1, (double) x2);
440: }
441:
442: #endif /* EXTFLOAT */
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.