src/float/dfloat/algebraic/cl_DF_sqrt.cc

   1 // sqrt().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_dfloat.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_DF.h"
  13 #include "cl_F.h"
  14 #include "cl_low.h"
  15 #include "cl_DS.h"
  16
  17 const cl_DF sqrt (const cl_DF& x)
  18 {
  19 // Methode:
  20 // x = 0.0 -> Ergebnis 0.0
  21 // Ergebnis-Vorzeichen := positiv,
  22 // Ergebnis-Exponent := ceiling(e/2),
  23 // Ergebnis-Mantisse:
  24 //   Bilde aus [1,m51,...,m0,(55 Nullbits)] bei geradem e,
  25 //         aus [0,1,m51,...,m0,(54 Nullbits)] bei ungeradem e
  26 //   die Ganzzahl-Wurzel, eine 54-Bit-Zahl mit einer führenden 1.
  27 //   Runde das letzte Bit weg:
  28 //     Bit 0 = 0 -> abrunden,
  29 //     Bit 0 = 1 und Wurzel exakt -> round-to-even,
  30 //     Bit 0 = 1 und Rest >0 -> aufrunden.
  31 //   Dabei um ein Bit nach rechts schieben.
  32 //   Bei Aufrundung auf 2^53 (rounding overflow) Mantisse um 1 Bit nach rechts
  33 //     schieben und Exponent incrementieren.
  34 #if (cl_word_size==64)
  35       // x entpacken:
  36       var sintL exp;
  37       var uint64 mantx;
  38       DF_decode(x, { return x; }, ,exp=,mantx=);
  39       // Um die 128-Bit-Ganzzahl-Wurzel ausnutzen zu können, fügen wir beim
  40       // Radikanden 74 bzw. 75 statt 54 bzw. 55 Nullbits an.
  41       if (exp & bit(0))
  42         // e ungerade
  43         { mantx = mantx << (63-(DF_mant_len+1)); exp = exp+1; }
  44         else
  45         // e gerade
  46         { mantx = mantx << (64-(DF_mant_len+1)); }
  47       exp = exp >> 1; // exp := exp/2
  48       var uintD mant [128/intDsize];
  49       #if (intDsize==64)
  50       arrayLSref(mant,128/intDsize,1) = mantx;
  51       arrayLSref(mant,128/intDsize,0) = 0;
  52       #else // (intDsize<=32)
  53       set_32_Dptr(arrayMSDptr(mant,128/intDsize),(uint32)(mantx>>32));
  54       set_32_Dptr(arrayMSDptr(mant,128/intDsize) mspop 32/intDsize,(uint32)mantx);
  55       set_32_Dptr(arrayMSDptr(mant,128/intDsize) mspop 2*32/intDsize,0);
  56       set_32_Dptr(arrayMSDptr(mant,128/intDsize) mspop 3*32/intDsize,0);
  57       #endif
  58       {CL_ALLOCA_STACK;
  59        var DS wurzel;
  60        var cl_boolean exactp;
  61        UDS_sqrt(arrayMSDptr(mant,128/intDsize),128/intDsize,arrayLSDptr(mant,128/intDsize), &wurzel, exactp=);
  62        // wurzel = isqrt(2^74_75 * mant), eine 64-Bit-Zahl.
  63        mantx = get_64_Dptr(wurzel.MSDptr);
  64        // Die hinteren 63-DF_mant_len Bits wegrunden:
  65        if ( ((mantx & bit(62-DF_mant_len)) ==0) // Bit 10 =0 -> abrunden
  66             || ( ((mantx & (bit(62-DF_mant_len)-1)) ==0) // Bit 10 =1 und Bits 9..0 >0 -> aufrunden
  67                  && exactp                   // Bit 10 =1 und Bits 9..0 =0, aber Rest -> aufrunden
  68                  // round-to-even, je nach Bit 11 :
  69                  && ((mantx & bit(63-DF_mant_len)) ==0)
  70           )    )
  71          // abrunden
  72          { mantx = mantx >> (63-DF_mant_len); }
  73          else
  74          // aufrunden
  75          { mantx = mantx >> (63-DF_mant_len);
  76            mantx += 1;
  77            if (mantx >= bit(DF_mant_len+1)) // rounding overflow?
  78              { mantx = mantx>>1; exp = exp+1; }
  79          }
  80       }
  81       return encode_DF(0,exp,mantx);
  82 #else
  83       // x entpacken:
  84       var sintL exp;
  85       var uint32 manthi;
  86       var uint32 mantlo;
  87       DF_decode2(x, { return x; }, ,exp=,manthi=,mantlo=);
  88       // Um die 128-Bit-Ganzzahl-Wurzel ausnutzen zu können, fügen wir beim
  89       // Radikanden 74 bzw. 75 statt 54 bzw. 55 Nullbits an.
  90       if (exp & bit(0))
  91         // e ungerade
  92         { manthi = (manthi << (63-(DF_mant_len+1))) | (mantlo >> ((DF_mant_len+1)-31));
  93           mantlo = mantlo << (63-(DF_mant_len+1));
  94           exp = exp+1;
  95         }
  96         else
  97         // e gerade
  98         { manthi = (manthi << (64-(DF_mant_len+1))) | (mantlo >> ((DF_mant_len+1)-32));
  99           mantlo = mantlo << (64-(DF_mant_len+1));
 100         }
 101       exp = exp >> 1; // exp := exp/2
 102       var uintD mant [128/intDsize];
 103       #if (intDsize==32) || (intDsize==16) || (intDsize==8)
 104       set_32_Dptr(arrayMSDptr(mant,128/intDsize),manthi);
 105       set_32_Dptr(arrayMSDptr(mant,128/intDsize) mspop 32/intDsize,mantlo);
 106       set_32_Dptr(arrayMSDptr(mant,128/intDsize) mspop 2*32/intDsize,0);
 107       set_32_Dptr(arrayMSDptr(mant,128/intDsize) mspop 3*32/intDsize,0);
 108       #else
 109       {var uintD* ptr;
 110        ptr = arrayLSDptr(mant,128/intDsize);
 111        doconsttimes(64/intDsize, { lsprefnext(ptr) = 0; } );
 112        doconsttimes(32/intDsize, { lsprefnext(ptr) = (uintD)mantlo; mantlo = mantlo>>intDsize; } );
 113        doconsttimes(32/intDsize, { lsprefnext(ptr) = (uintD)manthi; manthi = manthi>>intDsize; } );
 114       }
 115       #endif
 116       {CL_ALLOCA_STACK;
 117        var DS wurzel;
 118        var cl_boolean exactp;
 119        UDS_sqrt(arrayMSDptr(mant,128/intDsize),128/intDsize,arrayLSDptr(mant,128/intDsize), &wurzel, exactp=);
 120        // wurzel = isqrt(2^74_75 * mant), eine 64-Bit-Zahl.
 121        {var uintD* ptr = wurzel.MSDptr;
 122         manthi = get_32_Dptr(ptr); mantlo = get_32_Dptr(ptr mspop 32/intDsize);
 123        }
 124        // Die hinteren 63-DF_mant_len Bits wegrunden:
 125        if ( ((mantlo & bit(62-DF_mant_len)) ==0) // Bit 10 =0 -> abrunden
 126             || ( ((mantlo & (bit(62-DF_mant_len)-1)) ==0) // Bit 10 =1 und Bits 9..0 >0 -> aufrunden
 127                  && exactp                   // Bit 10 =1 und Bits 9..0 =0, aber Rest -> aufrunden
 128                  // round-to-even, je nach Bit 11 :
 129                  && ((mantlo & bit(63-DF_mant_len)) ==0)
 130           )    )
 131          // abrunden
 132          { mantlo = (mantlo >> (63-DF_mant_len)) | (manthi << (DF_mant_len-32+1));
 133            manthi = manthi >> (63-DF_mant_len);
 134          }
 135          else
 136          // aufrunden
 137          { mantlo = (mantlo >> (63-DF_mant_len)) | (manthi << (DF_mant_len-32+1));
 138            manthi = manthi >> (63-DF_mant_len);
 139            mantlo += 1;
 140            if (mantlo==0)
 141              { manthi += 1;
 142                if (manthi >= bit(DF_mant_len-32+1)) // rounding overflow?
 143                  { manthi = manthi>>1; exp = exp+1; }
 144          }   }
 145       }
 146       return encode_DF(0,exp,manthi,mantlo);
 147 #endif
 148 }