src/float/dfloat/algebraic/cl_DF_sqrt.cc

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