workbench/libs/mathtrans/spsqrt.c

   1 /*
   2     Copyright © 1995-2003, The AROS Development Team. All rights reserved.
   3     $Id$
   4 */
   5
   6 #include "mathtrans_intern.h"
   7
   8 /*****************************************************************************
   9
  10     NAME */
  11
  12         AROS_LH1(float, SPSqrt,
  13
  14 /*  SYNOPSIS */
  15         AROS_LHA(float, fnum1, D0),
  16
  17 /*  LOCATION */
  18         struct Library *, MathTransBase, 16, MathTrans)
  19
  20 /*  FUNCTION
  21         Calculate square root of ffp number
  22
  23     INPUTS
  24
  25     RESULT
  26         Motorola fast floating point 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         fnum1 = M * 2^E
  40
  41         If exponent is an odd number:
  42
  43         fnum = ( M*2 ) * 2^ (E-1)
  44          Now E' = E-1 is an even number and
  45               -> sqrt(fnum) = sqrt(M)   * sqrt(2)   * sqrt (2^E')
  46                             = sqrt(M)   * sqrt(2)   * 2^(E'/2)
  47           (with sqrt(M*2)>1)
  48                             = sqrt(M)   * sqrt(2)   * 2^(E'/2)
  49                             = sqrt(M)   * 1/sqrt(2) * 2^(1+(E'/2))
  50                             = sqrt(M/2)             * 2^(1+(E'/2))
  51
  52
  53         If Exponent is an even number:
  54            -> sqrt(fnum) = sqrt(M) * sqrt (2^E) =
  55                          = sqrt(M) * 2^(E/2)
  56
  57         Now calculate the square root of the mantisse.
  58         The following algorithm calculates the square of a number + delta
  59         and compares it to the mantisse. If the square of that  number +
  60         delta is less than the mantisse then keep that number + delta.
  61         Otherwise calculate a lower delta and try again.
  62         Start out with number = 0;
  63
  64
  65         Exponent = -1;
  66         Root = 0;
  67         repeat
  68         {
  69           if ( ( Root + 2^Exponent ) ^2 < Mantisse)
  70             Root += 2^Exponent
  71           Exponent --;
  72         }
  73         until you`re happy with the accuracy
  74
  75 *****************************************************************************/
  76 {
  77     AROS_LIBFUNC_INIT
  78
  79     ULONG Res, ResSquared, Delta, X, TargetMantisse;
  80     int z;
  81     BYTE Exponent;
  82
  83     if (0 == fnum1)
  84     {
  85         SetSR(Zero_Bit, Zero_Bit | Negative_Bit | Overflow_Bit);
  86         return 0;
  87     }
  88     /* is fnum negative */
  89     if ((char)fnum1 < 0)
  90     {
  91         SetSR(Overflow_Bit, Zero_Bit | Negative_Bit | Overflow_Bit);
  92         return fnum1;
  93     }
  94
  95     /* let us calculate the square-root now! */
  96
  97     TargetMantisse = ((ULONG)fnum1 & FFPMantisse_Mask);
  98
  99     if (fnum1 & 1)
 100     {
 101         /* TargetMantisse = TargetMantisse / 2; */
 102         TargetMantisse >>= 1;
 103         Exponent = ((BYTE)fnum1 >> 1) + 0x21;
 104     }
 105     else
 106     {
 107         Exponent = ((BYTE)fnum1 >> 1) + 0x20;
 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 =
 123               = Res^2      + 2*Res*Delta + Delta^2
 124         */
 125         X = ResSquared + (Res >> z)  + (Delta >> (z+1));
 126         z++;
 127         if (X <= TargetMantisse)
 128         {
 129             Res += Delta;
 130             ResSquared = X;
 131         }
 132     }
 133
 134     /* an adjustment for precision */
 135     if ((char) Res < 0) Res += 0x100;
 136
 137     Res &= FFPMantisse_Mask;
 138     Res |= Exponent;
 139
 140     return Res;
 141
 142     AROS_LIBFUNC_EXIT
 143 }