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