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.S.
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
81union double_long
82{
83  double d;
84  struct {
85      long upper;
86      unsigned long lower;
87    } l;
88};
89
90union float_long {
91  float f;
92  long l;
93};
94
95union 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
108int
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
122int
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 */
139double
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 */
172double
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 */
218float
219__floatunsisf (unsigned long l)
220{
221  double foo = __floatunsidf (l);
222  return foo;
223}
224
225/* convert int to float */
226float
227__floatsisf (long l)
228{
229  double foo = __floatsidf (l);
230  return foo;
231}
232
233/* convert float to double */
234double
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 */
273float
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 */
341long
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 */
373long
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.  */
391double __floatunsidf (unsigned long);
392double __floatsidf (long);
393float __floatsisf (long);
394double __extendsfdf2 (float);
395float __truncdfsf2 (double);
396long __fixdfsi (double);
397long __fixsfsi (float);
398
399int
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 */
416long 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 */
449double
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 */
481long 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 */
489float
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 */
497long 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 */
505long 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 */
513long
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
522long double
523__addxf3 (long double x1, long double x2)
524{
525  return (double) x1 + (double) x2;
526}
527
528long double
529__subxf3 (long double x1, long double x2)
530{
531  return (double) x1 - (double) x2;
532}
533
534long double
535__mulxf3 (long double x1, long double x2)
536{
537  return (double) x1 * (double) x2;
538}
539
540long double
541__divxf3 (long double x1, long double x2)
542{
543  return (double) x1 / (double) x2;
544}
545
546long double
547__negxf2 (long double x1)
548{
549  return - (double) x1;
550}
551
552long
553__cmpxf2 (long double x1, long double x2)
554{
555  return __cmpdf2 ((double) x1, (double) x2);
556}
557
558long
559__eqxf2 (long double x1, long double x2)
560{
561  return __cmpdf2 ((double) x1, (double) x2);
562}
563
564long
565__nexf2 (long double x1, long double x2)
566{
567  return __cmpdf2 ((double) x1, (double) x2);
568}
569
570long
571__ltxf2 (long double x1, long double x2)
572{
573  return __cmpdf2 ((double) x1, (double) x2);
574}
575
576long
577__lexf2 (long double x1, long double x2)
578{
579  return __cmpdf2 ((double) x1, (double) x2);
580}
581
582long
583__gtxf2 (long double x1, long double x2)
584{
585  return __cmpdf2 ((double) x1, (double) x2);
586}
587
588long
589__gexf2 (long double x1, long double x2)
590{
591  return __cmpdf2 ((double) x1, (double) x2);
592}
593
594#endif /* !__mcoldfire__ */
595#endif /* EXTFLOAT */
596