src/float/dfloat/elem/cl_DF_mul.cc

   1 // binary operator *
   2
   3 // General includes.
   4 #include "base/cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/dfloat.h"
   8
   9
  10 // Implementation.
  11
  12 #include "float/dfloat/cl_DF.h"
  13 #include "float/cl_F.h"
  14 #include "base/cl_low.h"
  15 #include "base/digitseq/cl_DS.h"
  16
  17 #include "base/cl_inline.h"
  18 #include "float/dfloat/elem/cl_DF_zerop.cc"
  19
  20 namespace cln {
  21
  22
  23 const cl_DF operator* (const cl_DF& x1, const cl_DF& x2)
  24 {
  25 // Methode:
  26 // Falls x1=0.0 oder x2=0.0 -> Ergebnis 0.0
  27 // Sonst: Ergebnis-Vorzeichen = VZ von x1 xor VZ von x2.
  28 //        Ergebnis-Exponent = Summe der Exponenten von x1 und x2.
  29 //        Ergebnis-Mantisse = Produkt der Mantissen von x1 und x2, gerundet:
  30 //          2^-53 * mant1  *  2^-53 * mant2  =  2^-106 * (mant1*mant2),
  31 //          die Klammer ist >=2^104, <=(2^53-1)^2<2^106 .
  32 //          Falls die Klammer >=2^105 ist, um 53 Bit nach rechts schieben und
  33 //            runden: Falls Bit 52 Null, abrunden; falls Bit 52 Eins und
  34 //            Bits 51..0 alle Null, round-to-even; sonst aufrunden.
  35 //          Falls die Klammer <2^105 ist, um 52 Bit nach rechts schieben und
  36 //            runden: Falls Bit 51 Null, abrunden; falls Bit 51 Eins und
  37 //            Bits 50..0 alle Null, round-to-even; sonst aufrunden. Nach
  38 //            Aufrunden: Falls =2^53, um 1 Bit nach rechts schieben. Sonst
  39 //            Exponenten um 1 erniedrigen.
  40 #ifdef FAST_DOUBLE
  41       double_to_DF(DF_to_double(x1) * DF_to_double(x2), return ,
  42                    TRUE, TRUE, // Overflow und subnormale Zahl abfangen
  43                    !(zerop_inline(x1) || zerop_inline(x2)), // ein Ergebnis +/- 0.0
  44                                // ist genau dann in Wirklichkeit ein Underflow
  45                    FALSE, FALSE // keine Singularität, kein NaN als Ergebnis möglich
  46                   );
  47 #else
  48       // x1,x2 entpacken:
  49       var cl_signean sign1;
  50       var sintL exp1;
  51       #if (intDsize<=32)
  52       var uintL manthi1;
  53       var uintL mantlo1;
  54       #endif
  55       var cl_signean sign2;
  56       var sintL exp2;
  57       #if (intDsize<=32)
  58       var uintL manthi2;
  59       var uintL mantlo2;
  60       #endif
  61       #if (cl_word_size==64)
  62       var uint64 mantx1;
  63       DF_decode(x1, { return x1; }, sign1=,exp1=,mantx1=);
  64       #if (intDsize<=32)
  65       manthi1 = high32(mantx1); mantlo1 = low32(mantx1);
  66       #endif
  67       var uint64 mantx2;
  68       DF_decode(x2, { return x2; }, sign2=,exp2=,mantx2=);
  69       #if (intDsize<=32)
  70       manthi2 = high32(mantx2); mantlo2 = low32(mantx2);
  71       #endif
  72       #else
  73       DF_decode2(x1, { return x1; }, sign1=,exp1=,manthi1=,mantlo1=);
  74       DF_decode2(x2, { return x2; }, sign2=,exp2=,manthi2=,mantlo2=);
  75       #endif
  76       exp1 = exp1 + exp2; // Summe der Exponenten
  77       sign1 = sign1 ^ sign2; // Ergebnis-Vorzeichen
  78       // Mantissen mant1 und mant2 multiplizieren (64x64-Bit-Multiplikation):
  79       var uintD mant1 [64/intDsize];
  80       var uintD mant2 [64/intDsize];
  81       var uintD mant [128/intDsize];
  82       #if (intDsize==64)
  83       arrayLSref(mant1,64/intDsize,0) = mantx1;
  84       arrayLSref(mant2,64/intDsize,0) = mantx2;
  85       #elif (intDsize==32) || (intDsize==16) || (intDsize==8)
  86       set_32_Dptr(arrayMSDptr(mant1,64/intDsize),manthi1);
  87       set_32_Dptr(arrayMSDptr(mant1,64/intDsize) mspop 32/intDsize,mantlo1);
  88       set_32_Dptr(arrayMSDptr(mant2,64/intDsize),manthi2);
  89       set_32_Dptr(arrayMSDptr(mant2,64/intDsize) mspop 32/intDsize,mantlo2);
  90       #else
  91       {var uintD* ptr;
  92        ptr = arrayLSDptr(mant1,64/intDsize);
  93        doconsttimes(32/intDsize, { lsprefnext(ptr) = (uintD)mantlo1; mantlo1 = mantlo1>>intDsize; } );
  94        doconsttimes(32/intDsize, { lsprefnext(ptr) = (uintD)manthi1; manthi1 = manthi1>>intDsize; } );
  95       }
  96       {var uintD* ptr;
  97        ptr = arrayLSDptr(mant2,64/intDsize);
  98        doconsttimes(32/intDsize, { lsprefnext(ptr) = (uintD)mantlo2; mantlo2 = mantlo2>>intDsize; } );
  99        doconsttimes(32/intDsize, { lsprefnext(ptr) = (uintD)manthi2; manthi2 = manthi2>>intDsize; } );
 100       }
 101       #endif
 102       cl_UDS_mul(arrayLSDptr(mant1,64/intDsize),64/intDsize,
 103                  arrayLSDptr(mant2,64/intDsize),64/intDsize,
 104                  arrayLSDptr(mant,128/intDsize)
 105                 );
 106       {
 107         #if (cl_word_size==64)
 108         var uint64 manterg;
 109         #else
 110         var uintL manthi;
 111         var uintL mantlo;
 112         #endif
 113         // Produkt mant = mant1 * mant2 ist >= 2^104, < 2^106. Bit 105 abtesten:
 114         #define mant_bit(k)  (arrayLSref(mant,128/intDsize,floor(k,intDsize)) & bit((k)%intDsize))
 115         if (mant_bit(2*DF_mant_len+1))
 116           // mant>=2^(2*DF_mant_len+1), um DF_mant_len+1 Bits nach rechts schieben:
 117           { // Bits 105..53 holen:
 118             #if (cl_word_size==64) && (intDsize==64)
 119               manterg = ((uint64)arrayLSref(mant,2,1) << 11) | ((uint64)arrayLSref(mant,2,0) >> 53); // Bits 116..53
 120               #define mantrest() (arrayLSref(mant,2,0) & (bit(53)-1))
 121             #elif (cl_word_size==64) && (intDsize==32)
 122               manterg = ((uint64)arrayLSref(mant,4,3) << 43) | ((uint64)arrayLSref(mant,4,2) << 11) | ((uint64)arrayLSref(mant,4,1) >> 21); // Bits 116..53
 123               #define mantrest() ((arrayLSref(mant,4,1) & (bit(21)-1)) || arrayLSref(mant,4,0))
 124             #elif (intDsize==32)
 125               manthi = ((uint32)arrayLSref(mant,4,3) << 11) | ((uint32)arrayLSref(mant,4,2) >> 21); // Bits 116..85
 126               mantlo = ((uint32)arrayLSref(mant,4,2) << 11) | ((uint32)arrayLSref(mant,4,1) >> 21); // Bits 84..53
 127               #define mantrest() ((arrayLSref(mant,4,1) & (bit(21)-1)) || arrayLSref(mant,4,0))
 128             #elif (intDsize==16)
 129               manthi = ((uint32)arrayLSref(mant,8,7) << 27) | ((uint32)arrayLSref(mant,8,6) << 11) | ((uint32)arrayLSref(mant,8,5) >> 5); // Bits 116..85
 130               mantlo = ((uint32)arrayLSref(mant,8,5) << 27) | ((uint32)arrayLSref(mant,8,4) << 11) | ((uint32)arrayLSref(mant,8,3) >> 5); // Bits 84..53
 131               #define mantrest() ((arrayLSref(mant,8,3) & (bit(5)-1)) || arrayLSref(mant,8,2) || arrayLSref(mant,8,1) || arrayLSref(mant,8,0))
 132             #elif (intDsize==8)
 133               manthi = ((uint32)arrayLSref(mant,16,14) << 27) | ((uint32)arrayLSref(mant,16,13) << 19) | ((uint32)arrayLSref(mant,16,12) << 11) | ((uint32)arrayLSref(mant,16,11) << 3) | ((uint32)arrayLSref(mant,16,10) >> 5); // Bits 116..85
 134               mantlo = ((uint32)arrayLSref(mant,16,10) << 27) | ((uint32)arrayLSref(mant,16,9) << 19) | ((uint32)arrayLSref(mant,16,8) << 11) | ((uint32)arrayLSref(mant,16,7) << 3) | ((uint32)arrayLSref(mant,16,6) >> 5); // Bits 84..53
 135               #define mantrest() ((arrayLSref(mant,16,6) & (bit(5)-1)) || arrayLSref(mant,16,5) || arrayLSref(mant,16,4) || arrayLSref(mant,16,3) || arrayLSref(mant,16,2) || arrayLSref(mant,16,1) || arrayLSref(mant,16,0))
 136             #endif
 137             if ( (mant_bit(DF_mant_len) ==0) // Bit DF_mant_len =0 -> abrunden
 138                  || ( !mantrest() // Bit DF_mant_len =1 und Bits DF_mant_len-1..0 >0 -> aufrunden
 139                       // round-to-even, je nach Bit DF_mant_len+1 :
 140                       && (mant_bit(DF_mant_len+1) ==0)
 141                )    )
 142               // abrunden
 143               goto ab;
 144               else
 145               // aufrunden
 146               goto auf;
 147             #undef mantrest
 148           }
 149           else
 150           // mant<2^(2*DF_mant_len+1), um DF_mant_len Bits nach rechts schieben:
 151           { exp1 = exp1-1; // Exponenten decrementieren
 152             // Bits 104..52 holen:
 153             #if (cl_word_size==64) && (intDsize==64)
 154               manterg = ((uint64)arrayLSref(mant,2,1) << 12) | ((uint64)arrayLSref(mant,2,0) >> 52); // Bits 115..52
 155               #define mantrest() (arrayLSref(mant,2,0) & (bit(52)-1))
 156             #elif (cl_word_size==64) && (intDsize==32)
 157               manterg = ((uint64)arrayLSref(mant,4,3) << 44) | ((uint64)arrayLSref(mant,4,2) << 12) | ((uint64)arrayLSref(mant,4,1) >> 20); // Bits 115..52
 158               #define mantrest() ((arrayLSref(mant,4,1) & (bit(20)-1)) || arrayLSref(mant,4,0))
 159             #elif (intDsize==32)
 160               manthi = ((uint32)arrayLSref(mant,4,3) << 12) | ((uint32)arrayLSref(mant,4,2) >> 20); // Bits 115..84
 161               mantlo = ((uint32)arrayLSref(mant,4,2) << 12) | ((uint32)arrayLSref(mant,4,1) >> 20); // Bits 83..52
 162               #define mantrest() ((arrayLSref(mant,4,1) & (bit(20)-1)) || arrayLSref(mant,4,0))
 163             #elif (intDsize==16)
 164               manthi = // ((uint32)arrayLSref(mant,8,7) << 28) | ((uint32)arrayLSref(mant,8,6) << 12) | ((uint32)arrayLSref(mant,8,5) >> 4); // Bits 115..84
 165               mantlo = // ((uint32)arrayLSref(mant,8,5) << 28) | ((uint32)arrayLSref(mant,8,4) << 12) | ((uint32)arrayLSref(mant,8,3) >> 4); // Bits 83..52
 166               #define mantrest() ((arrayLSref(mant,8,3) & (bit(4)-1)) || arrayLSref(mant,8,2) || arrayLSref(mant,8,1) || arrayLSref(mant,8,0))
 167             #elif (intDsize==8)
 168               manthi = ((uint32)arrayLSref(mant,16,14) << 28) | ((uint32)arrayLSref(mant,16,13) << 20) | ((uint32)arrayLSref(mant,16,12) << 12) | ((uint32)arrayLSref(mant,16,11) << 4) | ((uint32)arrayLSref(mant,16,10) >> 4); // Bits 115..84
 169               mantlo = ((uint32)arrayLSref(mant,16,10) << 28) | ((uint32)arrayLSref(mant,16,9) << 20) | ((uint32)arrayLSref(mant,16,8) << 12) | ((uint32)arrayLSref(mant,16,7) << 4) | ((uint32)arrayLSref(mant,16,6) >> 4); // Bits 83..52
 170               #define mantrest() ((arrayLSref(mant,16,6) & (bit(4)-1)) || arrayLSref(mant,16,5) || arrayLSref(mant,16,4) || arrayLSref(mant,16,3) || arrayLSref(mant,16,2) || arrayLSref(mant,16,1) || arrayLSref(mant,16,0))
 171             #endif
 172             if ( (mant_bit(DF_mant_len-1) ==0) // Bit DF_mant_len-1 =0 -> abrunden
 173                  || ( !mantrest() // Bit DF_mant_len-1 =1 und Bits DF_mant_len-2..0 >0 -> aufrunden
 174                       // round-to-even, je nach Bit DF_mant_len :
 175                       && (mant_bit(DF_mant_len) ==0)
 176                )    )
 177               // abrunden
 178               goto ab;
 179               else
 180               // aufrunden
 181               goto auf;
 182             #undef mantrest
 183           }
 184         #undef mant_bit
 185         auf:
 186         #if (cl_word_size==64)
 187         manterg = manterg+1;
 188         // Hier ist 2^DF_mant_len <= manterg <= 2^(DF_mant_len+1)
 189         if (manterg >= bit(DF_mant_len+1)) // rounding overflow?
 190           { manterg = manterg>>1; exp1 = exp1+1; } // Shift nach rechts
 191         #else
 192         mantlo = mantlo+1;
 193         if (mantlo==0)
 194           { manthi = manthi+1;
 195             // Hier ist 2^(DF_mant_len-32) <= manthi <= 2^(DF_mant_len-32+1)
 196             if (manthi >= bit(DF_mant_len-32+1)) // rounding overflow?
 197               { manthi = manthi>>1; exp1 = exp1+1; } // Shift nach rechts
 198           }
 199         #endif
 200         ab:
 201         // Runden fertig, 2^DF_mant_len <= manterg < 2^(DF_mant_len+1)
 202         #if (cl_word_size==64)
 203         return encode_DF(sign1,exp1,manterg);
 204         #else
 205         return encode_DF(sign1,exp1,manthi,mantlo);
 206         #endif
 207       }
 208 #endif
 209 }
 210
 211 }  // namespace cln