src/float/sfloat/elem/cl_SF_mul.cc

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