src/float/lfloat/elem/cl_LF_mul.cc

   1 // binary operator *
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/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 namespace cln {
  18
  19 const cl_LF operator* (const cl_LF& x1, const cl_LF& x2)
  20 {
  21 // Methode:
  22 // Falls x1=0.0 oder x2=0.0 -> Ergebnis 0.0
  23 // Sonst: Ergebnis-Vorzeichen = VZ von x1 xor VZ von x2.
  24 //        Ergebnis-Exponent = Summe der Exponenten von x1 und x2.
  25 //        Produkt der Mantissen bilden (2n Digits).
  26 //        Falls das fhrende Bit =0 ist: Mantissenprodukt um 1 Bit nach links
  27 //          schieben (die vorderen n+1 Digits gengen)
  28 //          und Exponent decrementieren.
  29 //        Runden auf n Digits liefert die Ergebnis-Mantisse.
  30       var uintC len1 = TheLfloat(x1)->len;
  31       var uintC len2 = TheLfloat(x2)->len;
  32       var uintC len = (len1 < len2 ? len1 : len2); // min. L�ge n von x1 und x2
  33       var uintE uexp1 = TheLfloat(x1)->expo;
  34       if (uexp1==0) // x1=0.0 -> Ergebnis 0.0
  35         { if (len < len1) return shorten(x1,len); else return x1; }
  36       var uintE uexp2 = TheLfloat(x2)->expo;
  37       if (uexp2==0) // x2=0.0 -> Ergebnis 0.0
  38         { if (len < len2) return shorten(x2,len); else return x2; }
  39       // Exponenten addieren:
  40       // (uexp1-LF_exp_mid) + (uexp2-LF_exp_mid) = (uexp1+uexp2-LF_exp_mid)-LF_exp_mid
  41       uexp1 = uexp1 + uexp2;
  42       if (uexp1 >= uexp2)
  43         // kein Carry
  44         { if (uexp1 < LF_exp_mid+LF_exp_low)
  45             { if (underflow_allowed())
  46                 { throw floating_point_underflow_exception(); }
  47                 else
  48                 { return encode_LF0(len); } // Ergebnis 0.0
  49         }   }
  50         else
  51         // Carry
  52         { if (uexp1 > (uintE)(LF_exp_mid+LF_exp_high+1)) { throw floating_point_overflow_exception(); } }
  53       uexp1 = uexp1 - LF_exp_mid;
  54       // Nun ist LF_exp_low <= uexp1 <= LF_exp_high+1.
  55       // neues Long-Float allozieren:
  56       var Lfloat y = allocate_lfloat(len,uexp1,
  57                                      TheLfloat(x1)->sign ^ TheLfloat(x2)->sign // Vorzeichen kombinieren
  58                                     );
  59       // Produkt bilden:
  60       var const uintD* x1_LSDptr = arrayLSDptr(TheLfloat(x1)->data,len1);
  61       var const uintD* x2_LSDptr = arrayLSDptr(TheLfloat(x2)->data,len2);
  62 #ifndef CL_LF_PEDANTIC
  63       if (len1 > len2)
  64         { x1_LSDptr = x1_LSDptr lspop (len1-(len2+1)); len1 = len2+1; }
  65       else if (len1 < len2)
  66         { x2_LSDptr = x2_LSDptr lspop (len2-(len1+1)); len2 = len1+1; }
  67 #endif
  68       var uintD* MSDptr;
  69       CL_ALLOCA_STACK;
  70       UDS_UDS_mul_UDS(len1,x1_LSDptr,
  71                       len2,x2_LSDptr,
  72                       MSDptr=,,);
  73       {var uintD* midptr = MSDptr mspop len; // Pointer in die Mitte der len1+len2 Digits
  74        if ((sintD)mspref(MSDptr,0) >= 0) // fhrendes Bit abtesten
  75          { // erste n+1 Digits um 1 Bit nach links schieben:
  76            shift1left_loop_lsp(midptr mspop 1,len+1);
  77            // Exponenten decrementieren:
  78            if (--(TheLfloat(y)->expo) == LF_exp_low-1)
  79              { if (underflow_allowed())
  80                  { throw floating_point_underflow_exception(); }
  81                  else
  82                  { return encode_LF0(len); } // Ergebnis 0.0
  83              }
  84          }
  85        // erste H�fte des Mantissenprodukts bertragen:
  86        {var uintD* y_mantMSDptr = arrayMSDptr(TheLfloat(y)->data,len);
  87         var uintD* y_mantLSDptr = copy_loop_msp(MSDptr,y_mantMSDptr,len);
  88         // Runden:
  89         if ( ((sintD)mspref(midptr,0) >= 0) // n�hstes Bit =0 -> abrunden
  90              || ( ((mspref(midptr,0) & ((uintD)bit(intDsize-1)-1)) ==0) // Bit =1, weitere Bits >0 -> aufrunden
  91                   && !test_loop_msp(midptr mspop 1,len1+len2-len-1)
  92                   // round-to-even
  93                   && ((lspref(midptr,0) & bit(0)) ==0)
  94            )    )
  95           // abrunden
  96           {}
  97           else
  98           // aufrunden
  99           { if ( inc_loop_lsp(y_mantLSDptr,len) )
 100               { // �ertrag durchs Aufrunden (kann nur auftreten,
 101                 // wenn vorhin um 1 Bit nach links geschoben wurde)
 102                 mspref(y_mantMSDptr,0) = bit(intDsize-1); // Mantisse := 10...0
 103                 (TheLfloat(y)->expo)++; // Exponent wieder zurck-erhhen
 104           }   }
 105         // LF_exp_low <= exp <= LF_exp_high sicherstellen:
 106         if (TheLfloat(y)->expo == LF_exp_high+1) { throw floating_point_overflow_exception(); }
 107       }}
 108       return y;
 109 }
 110 // Bit complexity (N = max(length(x1),length(x2))): O(M(N)).
 111
 112 }  // namespace cln