workbench/libs/mathieeedoubtrans/ieeedpsqrt.c

   1 /*
   2     Copyright © 1995-2003, The AROS Development Team. All rights reserved.
   3     $Id$
   4 */
   5
   6 #include "mathieeedoubtrans_intern.h"
   7
   8 /*
   9     FUNCTION
  10       Calculate square root of IEEE double precision floating point 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       ALGORITHM:
  30         <p>First check for a zero and a negative argument and take
  31         appropriate action.</p>
  32
  33         <code>fnum = M * 2^E</code>
  34
  35         <p>Exponent is an odd number:</p>
  36
  37         <code>
  38         fnum = ( M*2 ) * 2^ (E-1)
  39         Now E' = E-1 is an even number and
  40            -> sqrt(fnum) = sqrt(M)   * sqrt(2)   * sqrt (2^E')
  41                          = sqrt(M)   * sqrt(2)   * 2^(E'/2)
  42         (with sqrt(M*2)>1)
  43                          = sqrt(M)   * sqrt(2)   * 2^(E'/2)
  44                                                                                      = sqrt(M)   * 1/sqrt(2) * 2^(1+(E'/2))
  45                          = sqrt(M/2)             * 2^(1+(E'/2))
  46
  47         </code>
  48
  49         <p>Exponent is an even number:</p>
  50
  51         <code>-> sqrt(fnum) = sqrt(M) * sqrt (2^E) =
  52                       = sqrt(M) * 2^(E/2)
  53         </code>
  54
  55         <p>Now calculate the square root of the mantisse.
  56         The following algorithm calculates the square of a number + delta
  57         and compares it to the mantisse. If the square of that  number +
  58         delta is less than the mantisse then keep that number + delta.
  59         Otherwise calculate a lower offset and try again.
  60         Start out with number = 0;</p>
  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         </code>
  72
  73         until you`re happy with the accuracy
  74
  75     HISTORY
  76 */
  77
  78 AROS_LHQUAD1(double, IEEEDPSqrt,
  79     AROS_LHAQUAD(double, y, D0, D1),
  80     struct MathIeeeDoubTransBase *, MathIeeeDoubTransBase, 16, MathIeeeDoubTrans
  81 )
  82 {
  83     AROS_LIBFUNC_INIT
  84
  85     QUAD Res, ResSquared, Delta, X, TargetMantisse, y2, tmp;
  86     int z;
  87     ULONG Exponent;
  88
  89     if (is_eqC(y, 0, 0)) /* 0 == y */
  90     {
  91         SetSR(Zero_Bit, Zero_Bit | Negative_Bit | Overflow_Bit);
  92         Set_Value64C(Res, 0, 0);
  93         return Res;
  94     }
  95     /* is fnum negative? */
  96     if (is_lessSC(y, 0, 0)) /* y < 0 */
  97     {
  98         SetSR(Overflow_Bit, Zero_Bit | Negative_Bit | Overflow_Bit);
  99         Set_Value64C(Res, IEEEDPNAN_Hi, IEEEDPNAN_Lo);
 100         return Res;
 101     }
 102
 103     /* we can calculate the square-root now! */
 104
 105     Set_Value64(y2, y);
 106     AND64QC(y2, IEEEDPMantisse_Mask_Hi, IEEEDPMantisse_Mask_Lo);
 107     OR64QC(y2, 0x00100000, 0x00000000 );
 108     SHL64(TargetMantisse, y2, 11);
 109
 110     Exponent = (Get_High32of64(y) >> 1) & IEEEDPExponent_Mask_Hi;
 111
 112     /* do we have an odd exponent? */
 113     if (Get_High32of64(y) & 0x00100000)
 114     {
 115         /* TargetMantisse = TargetMantisse / 2; */
 116         SHRU64(TargetMantisse, TargetMantisse, 1); /* TargetMantisse >>= 1; */
 117         Exponent += 0x20000000;
 118     }
 119     else
 120     {
 121         Exponent += 0x1ff00000;
 122     }
 123
 124     Set_Value64C(Res,        0x0, 0x0);
 125     Set_Value64C(ResSquared, 0x0, 0x0);
 126     z = 0;
 127
 128     /*
 129         this calculates the sqrt of the mantisse. It`s short, isn`t it?
 130         Delta starts out with 0.5, then 0.25, 0.125 etc.
 131     */
 132     while ( is_neq(ResSquared, TargetMantisse) && z < 53)
 133     {
 134         QUAD Restmp, Deltatmp;
 135         Set_Value64C(Delta, 0x80000000, 0x00000000);
 136         SHRU64(Delta, Delta, z); /* Delta >> = z */
 137
 138         /*
 139             X = (Res+Delta)^2 = Res^2 + 2*Res*Delta + Delta^2
 140         */
 141         SHRU64(Restmp, Res, z);     /* Restmp   = Res   >> z */
 142         z++;
 143         SHRU64(Deltatmp, Delta, z); /* Deltatmp = Delta >> (z+1) */
 144
 145         /* X = ResSquared + (Res >> z)  + (Delta >> ++z); */
 146         Set_Value64(X, ResSquared);
 147         ADD64Q(X, Restmp);
 148         ADD64Q(X, Deltatmp);
 149
 150         if (is_leq(X, TargetMantisse)) /* X <= TargetMantisse */
 151         {
 152             //printf("setting a bit!\n");
 153             //OUTPUT(X);OUTPUT(TargetMantisse);
 154             ADD64(Res, Res, Delta );     /* Res += Delta */
 155             Set_Value64(ResSquared, X);   /* ResSquared = X; */
 156         }
 157     }
 158
 159     SHRU64(Res, Res, 11);
 160     AND64QC(Res, IEEEDPMantisse_Mask_Hi, IEEEDPMantisse_Mask_Lo );
 161     SHL32(tmp, Exponent, 32);
 162     OR64Q(Res, tmp);
 163
 164     return Res;
 165
 166     AROS_LIBFUNC_EXIT
 167 }