src/float/lfloat/elem/cl_LF_square.cc

   1 // square().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_lfloat.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_LF.h"
  13 #include "cl_LF_impl.h"
  14 #include "cl_DS.h"
  15 #include "cl_F.h"
  16
  17 const cl_LF square (const cl_LF& x)
  18 {
  19 // Methode: wie operator*(x,x).
  20       var uintC len = TheLfloat(x)->len;
  21       var uintL uexp = TheLfloat(x)->expo;
  22       if (uexp==0) // x=0.0 -> Ergebnis 0.0
  23         { return x; }
  24       // Exponenten addieren:
  25       // (uexp-LF_exp_mid) + (uexp-LF_exp_mid) = (2*uexp-LF_exp_mid)-LF_exp_mid
  26       if ((sintL)uexp >= 0)
  27         // kein Carry
  28         { uexp = 2*uexp;
  29           if (uexp < LF_exp_mid+LF_exp_low)
  30             { if (underflow_allowed())
  31                 { cl_error_floating_point_underflow(); }
  32                 else
  33                 { return encode_LF0(len); } // Ergebnis 0.0
  34         }   }
  35         else
  36         // Carry
  37         { uexp = 2*uexp;
  38           if (uexp > (uintL)(LF_exp_mid+LF_exp_high+1)) { cl_error_floating_point_overflow(); }
  39         }
  40       uexp = uexp - LF_exp_mid;
  41       // Nun ist LF_exp_low <= uexp <= LF_exp_high+1.
  42       // neues Long-Float allozieren:
  43       var Lfloat y = allocate_lfloat(len,uexp,0);
  44       // Produkt bilden:
  45       var const uintD* x_LSDptr = arrayLSDptr(TheLfloat(x)->data,len);
  46       var uintD* MSDptr;
  47       var uintD* LSDptr;
  48       CL_ALLOCA_STACK;
  49       num_stack_alloc(2*len,MSDptr=,LSDptr=);
  50       cl_UDS_mul_square(x_LSDptr,len,LSDptr);
  51       {var uintD* midptr = MSDptr mspop len; // Pointer in die Mitte der 2*len Digits
  52        if ((sintD)mspref(MSDptr,0) >= 0) // führendes Bit abtesten
  53          { // erste n+1 Digits um 1 Bit nach links schieben:
  54            shift1left_loop_lsp(midptr mspop 1,len+1);
  55            // Exponenten decrementieren:
  56            if ((TheLfloat(y)->expo)-- == LF_exp_low-1)
  57              { if (underflow_allowed())
  58                  { cl_error_floating_point_underflow(); }
  59                  else
  60                  { return encode_LF0(len); } // Ergebnis 0.0
  61              }
  62          }
  63        // erste Hälfte des Mantissenprodukts übertragen:
  64        {var uintD* y_mantMSDptr = arrayMSDptr(TheLfloat(y)->data,len);
  65         var uintD* y_mantLSDptr = copy_loop_msp(MSDptr,y_mantMSDptr,len);
  66         // Runden:
  67         if ( ((sintD)mspref(midptr,0) >= 0) // nächstes Bit =0 -> abrunden
  68              || ( ((mspref(midptr,0) & ((uintD)bit(intDsize-1)-1)) ==0) // Bit =1, weitere Bits >0 -> aufrunden
  69                   && !test_loop_msp(midptr mspop 1,len-1)
  70                   // round-to-even
  71                   && ((lspref(midptr,0) & bit(0)) ==0)
  72            )    )
  73           // abrunden
  74           {}
  75           else
  76           // aufrunden
  77           { if ( inc_loop_lsp(y_mantLSDptr,len) )
  78               { // Übertrag durchs Aufrunden (kann nur auftreten,
  79                 // wenn vorhin um 1 Bit nach links geschoben wurde)
  80                 mspref(y_mantMSDptr,0) = bit(intDsize-1); // Mantisse := 10...0
  81                 (TheLfloat(y)->expo)++; // Exponent wieder zurück-erhöhen
  82           }   }
  83         // LF_exp_low <= exp <= LF_exp_high sicherstellen:
  84         if (TheLfloat(y)->expo == LF_exp_high+1) { cl_error_floating_point_overflow(); }
  85       }}
  86       return y;
  87 }
  88 // Bit complexity (N = length(x)): O(M(N)).
  89