src/float/sfloat/elem/cl_SF_mul.cc

   1 // binary operator *
   2
   3 // General includes.
   4 #include "base/cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/sfloat.h"
   8
   9
  10 // Implementation.
  11
  12 #include "float/sfloat/cl_SF.h"
  13 #include "base/cl_low.h"
  14
  15 namespace cln {
  16
  17 const cl_SF operator* (const cl_SF& x1, const cl_SF& x2)
  18 {
  19 // Methode:
  20 // Falls x1=0.0 oder x2=0.0 -> Ergebnis 0.0
  21 // Sonst: Ergebnis-Vorzeichen = VZ von x1 xor VZ von x2.
  22 //        Ergebnis-Exponent = Summe der Exponenten von x1 und x2.
  23 //        Ergebnis-Mantisse = Produkt der Mantissen von x1 und x2, gerundet:
  24 //          2^-17 * (2^16 + m1)  *  2^-17 * (2^16 + m2)
  25 //          = 2^-34 * (2^32 + 2^16*m1 + 2^16*m2 + m1*m2),
  26 //          die Klammer ist >=2^32, <=(2^17-1)^2<2^34 .
  27 //          Falls die Klammer >=2^33 ist, um 17 Bit nach rechts schieben und
  28 //            runden: Falls Bit 16 Null, abrunden; falls Bit 16 Eins und
  29 //            Bits 15..0 alle Null, round-to-even; sonst aufrunden.
  30 //          Falls die Klammer <2^33 ist, um 16 Bit nach rechts schieben und
  31 //            runden: Falls Bit 15 Null, abrunden; falls Bit 15 Eins und
  32 //            Bits 14..0 alle Null, round-to-even; sonst aufrunden. Nach
  33 //            Aufrunden: Falls =2^17, um 1 Bit nach rechts schieben. Sonst
  34 //            Exponenten um 1 erniedrigen.
  35       // x1,x2 entpacken:
  36       var cl_signean sign1;
  37       var sintL exp1;
  38       var uintL mant1;
  39       var cl_signean sign2;
  40       var sintL exp2;
  41       var uintL mant2;
  42       SF_decode(x1, { return x1; }, sign1=,exp1=,mant1=);
  43       SF_decode(x2, { return x2; }, sign2=,exp2=,mant2=);
  44       exp1 = exp1 + exp2; // Summe der Exponenten
  45       sign1 = sign1 ^ sign2; // Ergebnis-Vorzeichen
  46       var uintL manthi;
  47       var uintL mantlo;
  48       // Mantissen mant1 und mant2 multiplizieren:
  49       #if (SF_mant_len<16)
  50       mantlo = mulu16(mant1,mant2);
  51       manthi = mantlo >> SF_mant_len;
  52       mantlo = mantlo & (bit(SF_mant_len)-1);
  53       #elif (SF_mant_len==16)
  54       manthi = mulu16(low16(mant1),low16(mant2));
  55       mantlo = low16(manthi);
  56       manthi = (uint32)(high16(manthi)) + (uint32)(low16(mant1)) + mant2;
  57       #else // (SF_mant_len>16)
  58       mulu24(mant1,mant2, manthi=,mantlo=);
  59       manthi = (manthi << (32-SF_mant_len)) | (mantlo >> SF_mant_len);
  60       mantlo = mantlo & (bit(SF_mant_len)-1);
  61       #endif
  62       // Nun ist 2^SF_mant_len * manthi + mantlo = mant1 * mant2.
  63       if (manthi >= bit(SF_mant_len+1))
  64         // mant1*mant2 >= 2^(2*SF_mant_len+1)
  65         { if ( ((manthi & bit(0)) ==0) // Bit SF_mant_len =0 -> abrunden
  66                || ( (mantlo ==0) // Bit SF_mant_len =1 und Bits SF_mant_len-1..0 >0 -> aufrunden
  67                     // round-to-even, je nach Bit SF_mant_len+1 :
  68                     && ((manthi & bit(1)) ==0)
  69              )    )
  70             // abrunden
  71             { manthi = manthi >> 1; goto ab; }
  72             else
  73             // aufrunden
  74             { manthi = manthi >> 1; goto auf; }
  75         }
  76         else
  77         // mant1*mant2 < 2^(2*SF_mant_len+1)
  78         { exp1 = exp1-1; // Exponenten decrementieren
  79           if ( ((mantlo & bit(SF_mant_len-1)) ==0) // Bit SF_mant_len-1 =0 -> abrunden
  80                || ( ((mantlo & (bit(SF_mant_len-1)-1)) ==0) // Bit SF_mant_len-1 =1 und Bits SF_mant_len-2..0 >0 -> aufrunden
  81                     // round-to-even, je nach Bit SF_mant_len :
  82                     && ((manthi & bit(0)) ==0)
  83              )    )
  84             // abrunden
  85             goto ab;
  86             else
  87             // aufrunden
  88             goto auf;
  89         }
  90       auf:
  91       manthi = manthi+1;
  92       // Hier ist 2^SF_mant_len <= manthi <= 2^(SF_mant_len+1)
  93       if (manthi >= bit(SF_mant_len+1)) // rounding overflow?
  94         { manthi = manthi>>1; exp1 = exp1+1; } // Shift nach rechts
  95       ab:
  96       // Runden fertig, 2^SF_mant_len <= manthi < 2^(SF_mant_len+1)
  97       return encode_SF(sign1,exp1,manthi);
  98 }
  99
 100 }  // namespace cln