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 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
__floatsidf(int a1)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
__floatsisf(int l)156 __floatsisf (int l)
157 {
158 double foo = __floatsidf (l);
159 return foo;
160 }
161
162 /* convert float to double */
163 double
__extendsfdf2(float a1)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
__truncdfsf2(double a1)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
__fixdfsi(double a1)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
__fixsfsi(float a1)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
__extenddfxf2(double d)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
__truncxfdf2(long double ld)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
__extendsfxf2(float f)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
__truncxfsf2(long double ld)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
__floatsixf(int l)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
__fixxfsi(long double ld)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
__addxf3(long double x1,long double x2)371 __addxf3 (long double x1, long double x2)
372 {
373 return (double) x1 + (double) x2;
374 }
375
376 long double
__subxf3(long double x1,long double x2)377 __subxf3 (long double x1, long double x2)
378 {
379 return (double) x1 - (double) x2;
380 }
381
382 long double
__mulxf3(long double x1,long double x2)383 __mulxf3 (long double x1, long double x2)
384 {
385 return (double) x1 * (double) x2;
386 }
387
388 long double
__divxf3(long double x1,long double x2)389 __divxf3 (long double x1, long double x2)
390 {
391 return (double) x1 / (double) x2;
392 }
393
394 long double
__negxf2(long double x1)395 __negxf2 (long double x1)
396 {
397 return - (double) x1;
398 }
399
400 long
__cmpxf2(long double x1,long double x2)401 __cmpxf2 (long double x1, long double x2)
402 {
403 return __cmpdf2 ((double) x1, (double) x2);
404 }
405
406 long
__eqxf2(long double x1,long double x2)407 __eqxf2 (long double x1, long double x2)
408 {
409 return __cmpdf2 ((double) x1, (double) x2);
410 }
411
412 long
__nexf2(long double x1,long double x2)413 __nexf2 (long double x1, long double x2)
414 {
415 return __cmpdf2 ((double) x1, (double) x2);
416 }
417
418 long
__ltxf2(long double x1,long double x2)419 __ltxf2 (long double x1, long double x2)
420 {
421 return __cmpdf2 ((double) x1, (double) x2);
422 }
423
424 long
__lexf2(long double x1,long double x2)425 __lexf2 (long double x1, long double x2)
426 {
427 return __cmpdf2 ((double) x1, (double) x2);
428 }
429
430 long
__gtxf2(long double x1,long double x2)431 __gtxf2 (long double x1, long double x2)
432 {
433 return __cmpdf2 ((double) x1, (double) x2);
434 }
435
436 long
__gexf2(long double x1,long double x2)437 __gexf2 (long double x1, long double x2)
438 {
439 return __cmpdf2 ((double) x1, (double) x2);
440 }
441
442 #endif /* EXTFLOAT */
443