src/float/lfloat/elem/cl_LF_I_mul.cc

   1 // cl_LF_I_mul().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_LF.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_LF_impl.h"
  13 #include "cln/integer.h"
  14 #include "cl_I.h"
  15 #include "cl_DS.h"
  16 #include "cl_F.h"
  17
  18 namespace cln {
  19
  20 const cl_R cl_LF_I_mul (const cl_LF& x, const cl_I& y)
  21 {
  22 // Method:
  23 // If y=0, return 0.
  24 // If x=0.0, return x.
  25 // If y is longer than x, convert y to a float and multiply.
  26 // Else multiply the mantissa of x with the absolute value of y, then round.
  27         if (eq(y,0)) { return 0; }
  28         if (TheLfloat(x)->expo == 0) { return x; }
  29         var cl_signean sign = -(cl_signean)minusp(y); // Vorzeichen von y
  30         var cl_I abs_y = (sign==0 ? y : -y);
  31         var uintC y_exp = integer_length(abs_y);
  32         var uintC len = TheLfloat(x)->len;
  33 #ifndef CL_LF_PEDANTIC
  34         if (ceiling(y_exp,intDsize) > len)
  35                 return x * cl_I_to_LF(y,len);
  36 #endif
  37         // x länger als y, direkt multiplizieren.
  38         CL_ALLOCA_STACK;
  39         var const uintD* y_MSDptr;
  40         var uintC y_len;
  41         var const uintD* y_LSDptr;
  42         I_to_NDS_nocopy(abs_y, y_MSDptr=,y_len=,y_LSDptr=,false,); // NDS zu y bilden, y_len>0
  43         if (mspref(y_MSDptr,0)==0) y_len--; // NUDS zu y bilden, y_len>0
  44         // Multiplizieren.
  45         var uintD* prodMSDptr;
  46         var uintC prodlen;
  47         UDS_UDS_mul_UDS(len,arrayLSDptr(TheLfloat(x)->data,len),
  48                         y_len,y_LSDptr,
  49                         prodMSDptr=,prodlen=,);
  50         // x fing mit 0 Nullbits an, y mit maximal intDsize-1 Nullbits,
  51         // daher fängt das Produkt mit maximal intDsize Nullbits an.
  52         var uintL shiftcount;
  53         if (mspref(prodMSDptr,0)==0) {
  54                 shiftcount = intDsize;
  55                 msshrink(prodMSDptr); prodlen--;
  56         } else {
  57                 integerlengthD(mspref(prodMSDptr,0), shiftcount = intDsize -);
  58                 if (shiftcount > 0)
  59                         shiftleft_loop_lsp(prodMSDptr mspop (len+1),len+1,shiftcount,0);
  60         }
  61         // Produkt ist nun normalisiert: höchstes Bit =1.
  62         // exponent := exponent(x) + intDsize*y_len - shiftcount
  63         var uintE uexp = TheLfloat(x)->expo;
  64         var uintE iexp = intDsize*y_len - shiftcount; // >= 0 !
  65         uexp = uexp + iexp;
  66         if ((uexp < iexp) || (uexp > LF_exp_high))
  67                 throw floating_point_overflow_exception();
  68         // Runden:
  69         var uintD* midptr = prodMSDptr mspop len;
  70         var uintC restlen = prodlen - len;
  71         if ( (restlen==0)
  72              || ((sintD)mspref(midptr,0) >= 0) // nächstes Bit =0 -> abrunden
  73              || ( ((mspref(midptr,0) & ((uintD)bit(intDsize-1)-1)) ==0) // Bit =1, weitere Bits >0 -> aufrunden
  74                   && !test_loop_msp(midptr mspop 1,restlen-1)
  75                   // round-to-even
  76                   && ((lspref(midptr,0) & bit(0)) ==0)
  77            )    )
  78           // abrunden
  79           {}
  80           else
  81           // aufrunden
  82           { if ( inc_loop_lsp(midptr,len) )
  83               // Übertrag durchs Aufrunden
  84               { mspref(prodMSDptr,0) = bit(intDsize-1); // Mantisse := 10...0
  85                 if (++uexp == LF_exp_high+1) { throw floating_point_overflow_exception(); }
  86           }   }
  87         return encode_LFu(TheLfloat(x)->sign ^ sign, uexp, prodMSDptr, len);
  88 }
  89 // Bit complexity (N = max(length(x),length(y))): O(M(N)).
  90
  91 }  // namespace cln