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     FUNCTION
  10       Calculate square root of ffp number
  11
  12     RESULT
  13       Motorola fast floating point 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
  30       ALGORITHM:
  31         First check for a zero and a negative argument and take
  32         appropriate action.
  33         fnum1 = M * 2^E
  34
  35         If exponent is an odd number:
  36      <code>
  37           fnum = ( M*2 ) * 2^ (E-1)
  38           Now E' = E-1 is an even number and
  39                -> sqrt(fnum) = sqrt(M)   * sqrt(2)   * sqrt (2^E')
  40                              = sqrt(M)   * sqrt(2)   * 2^(E'/2)
  41            (with sqrt(M*2)>1)
  42                              = sqrt(M)   * sqrt(2)   * 2^(E'/2)
  43                              = sqrt(M)   * 1/sqrt(2) * 2^(1+(E'/2))
  44                              = sqrt(M/2)             * 2^(1+(E'/2))
  45      </code>
  46
  47
  48        If Exponent is an even number:
  49      <code>
  50            -> sqrt(fnum) = sqrt(M) * sqrt (2^E) =
  51                          = sqrt(M) * 2^(E/2)
  52      </code>
  53
  54         Now calculate the square root of the mantisse.
  55         The following algorithm calculates the square of a number + delta
  56         and compares it to the mantisse. If the square of that  number +
  57         delta is less than the mantisse then keep that number + delta.
  58         Otherwise calculate a lower delta and try again.
  59         Start out with number = 0;
  60
  61
  62      <code>
  63         Exponent = -1;
  64         Root = 0;
  65         repeat
  66         {
  67           if ( ( Root + 2^Exponent ) ^2 < Mantisse)
  68             Root += 2^Exponent
  69           Exponent --;
  70         }
  71         until you`re happy with the accuracy
  72      </code>
  73
  74
  75     HISTORY
  76 */
  77
  78 AROS_LH1(float, SPSqrt,
  79     AROS_LHA(float, fnum1, D0),
  80     struct Library *, MathTransBase, 16, MathTrans
  81 )
  82 {
  83     AROS_LIBFUNC_INIT
  84
  85     ULONG Res, ResSquared, Delta, X, TargetMantisse;
  86     int z;
  87     BYTE Exponent;
  88
  89     if (0 == fnum1)
  90     {
  91         SetSR(Zero_Bit, Zero_Bit | Negative_Bit | Overflow_Bit);
  92         return 0;
  93     }
  94     /* is fnum negative */
  95     if ((char)fnum1 < 0)
  96     {
  97         SetSR(Overflow_Bit, Zero_Bit | Negative_Bit | Overflow_Bit);
  98         return fnum1;
  99     }
 100
 101     /* let us calculate the square-root now! */
 102
 103     TargetMantisse = ((ULONG)fnum1 & FFPMantisse_Mask);
 104
 105     if (fnum1 & 1)
 106     {
 107         /* TargetMantisse = TargetMantisse / 2; */
 108         TargetMantisse >>= 1;
 109         Exponent = ((BYTE)fnum1 >> 1) + 0x21;
 110     }
 111     else
 112     {
 113         Exponent = ((BYTE)fnum1 >> 1) + 0x20;
 114     }
 115
 116     Res = 0;
 117     ResSquared = 0;
 118     z = 0;
 119
 120     /*
 121         this calculates the sqrt of the mantisse. It`s short, isn`t it?
 122         Delta starts out with 0.5, then 0.25, 0.125 etc.
 123     */
 124     while (ResSquared != TargetMantisse && z < 25)
 125     {
 126         Delta = (0x80000000 >> z);
 127         /*
 128             X = (Res+Delta)^2 =
 129               = Res^2      + 2*Res*Delta + Delta^2
 130         */
 131         X = ResSquared + (Res >> z)  + (Delta >> (z+1));
 132         z++;
 133         if (X <= TargetMantisse)
 134         {
 135             Res += Delta;
 136             ResSquared = X;
 137         }
 138     }
 139
 140     /* an adjustment for precision */
 141     if ((char) Res < 0) Res += 0x100;
 142
 143     Res &= FFPMantisse_Mask;
 144     Res |= Exponent;
 145
 146     return Res;
 147
 148     AROS_LIBFUNC_EXIT
 149 }