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