workbench/libs/mathieeesingtrans/ieeespsqrt.c

   1 /*
   2     Copyright © 1995-2003, The AROS Development Team. All rights reserved.
   3     $Id$
   4 */
   5
   6 #include "mathieeesingtrans_intern.h"
   7
   8 /*
   9     FUNCTION
  10       Calculate square root of IEEE single precision number
  11
  12     RESULT
  13       IEEE single precision number
  14
  15       flags:
  16         zero     : result is zero
  17         negative : 0
  18         overflow : square root could not be calculated
  19
  20     NOTES
  21
  22     EXAMPLE
  23
  24     BUGS
  25
  26     SEE ALSO
  27
  28     INTERNALS
  29       ALGORITHM:
  30       <code>
  31         First check for a zero and a negative argument and take
  32         appropriate action.
  33         y = M * 2^E
  34
  35         Exponent is an odd number:
  36         y = ( M*2 ) * 2^ (E-1)
  37         Now E' = E-1 is an even number and
  38               -> sqrt(y) = sqrt(M)   * sqrt(2)   * sqrt (2^E')
  39                          = sqrt(M)   * sqrt(2)   * 2^(E'/2)
  40         (with sqrt(M*2)>1)
  41                          = sqrt(M)   * sqrt(2)   * 2^(E'/2)
  42                                                                                      = sqrt(M)   * 1/sqrt(2) * 2^(1+(E'/2))
  43                          = sqrt(M/2)             * 2^(1+(E'/2))
  44
  45
  46         Exponent is an even number:
  47         -> sqrt(y) = sqrt(M) * sqrt (2^E) =
  48                    = sqrt(M) * 2^(E/2)
  49
  50         Now calculate the square root of the mantisse.
  51         The following algorithm calculates the square of a number + delta
  52         and compares it to the mantisse. If the square of that  number +
  53         delta is less than the mantisse then keep that number + delta.
  54         Otherwise calculate a lower offset and try again.
  55         Start out with number = 0;
  56
  57
  58         Exponent = -1;
  59         Root = 0;
  60         repeat
  61         {
  62           if ( (Root + 2^Exponent)^2 < Mantisse)
  63           Root += 2^Exponent
  64           Exponent --;
  65         }
  66         until you`re happy with the accuracy
  67       </code>
  68
  69     HISTORY
  70 */
  71
  72 AROS_LH1(float, IEEESPSqrt,
  73     AROS_LHA(float, y, D0),
  74     struct Library *, MathIeeeSingTransBase, 16, MathIeeeSingTrans
  75 )
  76 {
  77     AROS_LIBFUNC_INIT
  78
  79     ULONG Res, ResSquared, Delta, X, TargetMantisse;
  80     int z;
  81     ULONG Exponent;
  82
  83     if (0 == y)
  84     {
  85         SetSR(Zero_Bit, Zero_Bit | Negative_Bit | Overflow_Bit);
  86         return 0;
  87     }
  88     /* is fnum negative? */
  89     if (y < 0)
  90     {
  91         SetSR(Overflow_Bit, Zero_Bit | Negative_Bit | Overflow_Bit);
  92         return y;
  93     }
  94
  95     /* we can calculate the square-root now! */
  96
  97     TargetMantisse = ((y & IEEESPMantisse_Mask) | 0x00800000) << 8;
  98
  99     if (y & 0x00800000)
 100     {
 101         /* TargetMantisse = TargetMantisse / 2; */
 102         TargetMantisse >>= 1;
 103         Exponent = ((y >> 1) & IEEESPExponent_Mask) + 0x20000000;
 104     }
 105     else
 106     {
 107         Exponent = ((y >> 1) & IEEESPExponent_Mask) + 0x1f800000;
 108     }
 109
 110     Res = 0;
 111     ResSquared = 0;
 112     z = 0;
 113
 114     /*
 115         this calculates the sqrt of the mantisse. It`s short, isn`t it?
 116         Delta starts out with 0.5, then 0.25, 0.125 etc.
 117     */
 118     while (ResSquared != TargetMantisse && z < 25)
 119     {
 120         Delta = (0x80000000 >> z);
 121         /*
 122             X = (Res+Delta)^2 = Res^2 + 2*Res*Delta + Delta^2
 123         */
 124         X = ResSquared + (Res >> z)  + (Delta >> (z+1));
 125         z++;
 126         if (X <= TargetMantisse)
 127         {
 128             Res += Delta;
 129             ResSquared = X;
 130         }
 131     }
 132
 133     /* an adjustment for precision
 134     */
 135     if ((char) Res < 0) Res += 0x100;
 136
 137     Res >>=8;
 138     Res &= IEEESPMantisse_Mask;
 139     Res |= Exponent;
 140
 141     return Res;
 142
 143     AROS_LIBFUNC_EXIT
 144 }