src/float/ffloat/elem/cl_FF_mul.cc

   1 // binary operator *
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/ffloat.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_FF.h"
  13 #include "cl_F.h"
  14 #include "cl_low.h"
  15 #include "cl_ieee.h"
  16
  17 #undef MAYBE_INLINE
  18 #define MAYBE_INLINE inline
  19 #include "cl_FF_zerop.cc"
  20
  21 namespace cln {
  22
  23 NEED_IEEE_FLOATS()
  24
  25 const cl_FF operator* (const cl_FF& x1, const cl_FF& x2)
  26 {
  27 // Methode:
  28 // Falls x1=0.0 oder x2=0.0 -> Ergebnis 0.0
  29 // Sonst: Ergebnis-Vorzeichen = VZ von x1 xor VZ von x2.
  30 //        Ergebnis-Exponent = Summe der Exponenten von x1 und x2.
  31 //        Ergebnis-Mantisse = Produkt der Mantissen von x1 und x2, gerundet:
  32 //          2^-24 * mant1  *  2^-24 * mant2  =  2^-48 * (mant1*mant2),
  33 //          die Klammer ist >=2^46, <=(2^24-1)^2<2^48 .
  34 //          Falls die Klammer >=2^47 ist, um 24 Bit nach rechts schieben und
  35 //            runden: Falls Bit 23 Null, abrunden; falls Bit 23 Eins und
  36 //            Bits 22..0 alle Null, round-to-even; sonst aufrunden.
  37 //          Falls die Klammer <2^47 ist, um 23 Bit nach rechts schieben und
  38 //            runden: Falls Bit 22 Null, abrunden; falls Bit 22 Eins und
  39 //            Bits 21..0 alle Null, round-to-even; sonst aufrunden. Nach
  40 //            Aufrunden: Falls =2^24, um 1 Bit nach rechts schieben. Sonst
  41 //            Exponenten um 1 erniedrigen.
  42   #ifdef FAST_FLOAT
  43       float_to_FF(FF_to_float(x1) * FF_to_float(x2), return ,
  44                   TRUE, TRUE, // Overflow und subnormale Zahl abfangen
  45                   !(zerop(x1) || zerop(x2)), // ein Ergebnis +/- 0.0
  46                               // ist genau dann in Wirklichkeit ein Underflow
  47                   FALSE, FALSE // keine Singularität, kein NaN als Ergebnis möglich
  48                  );
  49   #else
  50       // x1,x2 entpacken:
  51       var cl_signean sign1;
  52       var sintL exp1;
  53       var uintL mant1;
  54       var cl_signean sign2;
  55       var sintL exp2;
  56       var uintL mant2;
  57       FF_decode(x1, { return x1; }, sign1=,exp1=,mant1=);
  58       FF_decode(x2, { return x2; }, sign2=,exp2=,mant2=);
  59       exp1 = exp1 + exp2; // Summe der Exponenten
  60       sign1 = sign1 ^ sign2; // Ergebnis-Vorzeichen
  61       var uintL manthi;
  62       var uintL mantlo;
  63       // Mantissen mant1 und mant2 multiplizieren:
  64       mulu24(mant1,mant2, manthi=,mantlo=);
  65       manthi = (manthi << (32-FF_mant_len)) | (mantlo >> FF_mant_len);
  66       mantlo = mantlo & (bit(FF_mant_len)-1);
  67       // Nun ist 2^FF_mant_len * manthi + mantlo = mant1 * mant2.
  68       if (manthi >= bit(FF_mant_len+1))
  69         // mant1*mant2 >= 2^(2*FF_mant_len+1)
  70         { if ( ((manthi & bit(0)) ==0) // Bit FF_mant_len =0 -> abrunden
  71                || ( (mantlo ==0) // Bit FF_mant_len =1 und Bits FF_mant_len-1..0 >0 -> aufrunden
  72                     // round-to-even, je nach Bit FF_mant_len+1 :
  73                     && ((manthi & bit(1)) ==0)
  74              )    )
  75             // abrunden
  76             { manthi = manthi >> 1; goto ab; }
  77             else
  78             // aufrunden
  79             { manthi = manthi >> 1; goto auf; }
  80         }
  81         else
  82         // mant1*mant2 < 2^(2*FF_mant_len+1)
  83         { exp1 = exp1-1; // Exponenten decrementieren
  84           if ( ((mantlo & bit(FF_mant_len-1)) ==0) // Bit FF_mant_len-1 =0 -> abrunden
  85                || ( ((mantlo & (bit(FF_mant_len-1)-1)) ==0) // Bit FF_mant_len-1 =1 und Bits FF_mant_len-2..0 >0 -> aufrunden
  86                     // round-to-even, je nach Bit FF_mant_len :
  87                     && ((manthi & bit(0)) ==0)
  88              )    )
  89             // abrunden
  90             goto ab;
  91             else
  92             // aufrunden
  93             goto auf;
  94         }
  95       auf:
  96       manthi = manthi+1;
  97       // Hier ist 2^FF_mant_len <= manthi <= 2^(FF_mant_len+1)
  98       if (manthi >= bit(FF_mant_len+1)) // rounding overflow?
  99         { manthi = manthi>>1; exp1 = exp1+1; } // Shift nach rechts
 100       ab:
 101       // Runden fertig, 2^FF_mant_len <= manthi < 2^(FF_mant_len+1)
 102       return encode_FF(sign1,exp1,manthi);
 103   #endif
 104 }
 105
 106 }  // namespace cln