home *** CD-ROM | disk | FTP | other *** search
/ OpenStep (Enterprise) / OpenStepENTCD.toast / OEDEV / GNUSRC.Z / floatlib.c < prev    next >
C/C++ Source or Header  |  1992-01-04  |  11KB  |  582 lines
  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. }
  582.