src/float/dfloat/elem/cl_DF_plus.cc

   1 // binary operator +
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/dfloat.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_DF.h"
  13 #include "cl_F.h"
  14 #include "cl_ieee.h"
  15 #include "cl_xmacros.h"
  16
  17 namespace cln {
  18
  19 NEED_IEEE_FLOATS()
  20
  21 const cl_DF operator+ (const cl_DF& x1, const cl_DF& x2)
  22 {
  23 // Methode (nach [Knuth, II, Seminumerical Algorithms, Abschnitt 4.2.1., S.200]):
  24 // x1=0.0 -> Ergebnis x2.
  25 // x2=0.0 -> Ergebnis x1.
  26 // Falls e1<e2, vertausche x1 und x2.
  27 // Also e1 >= e2.
  28 // Falls e1 - e2 >= 52 + 3, Ergebnis x1.
  29 // Schiebe beide Mantissen um 3 Bits nach links (Vorbereitung der Rundung:
  30 //   Bei e1-e2=0,1 ist keine Rundung nötig, bei e1-e2>1 ist der Exponent des
  31 //   Ergebnisses =e1-1, =e1 oder =e1+1. Brauche daher 1 Schutzbit und zwei
  32 //   Rundungsbits: 00 exakt, 01 1.Hälfte, 10 exakte Mitte, 11 2.Hälfte.)
  33 // Schiebe die Mantisse von x2 um e0-e1 Bits nach rechts. (Dabei die Rundung
  34 // ausführen: Bit 0 ist das logische Oder der Bits 0,-1,-2,...)
  35 // Falls x1,x2 selbes Vorzeichen haben: Addiere dieses zur Mantisse von x1.
  36 // Falls x1,x2 verschiedenes Vorzeichen haben: Subtrahiere dieses von der
  37 //   Mantisse von x1. <0 -> (Es war e1=e2) Vertausche die Vorzeichen, negiere.
  38 //                    =0 -> Ergebnis 0.0
  39 // Exponent ist e1.
  40 // Normalisiere, fertig.
  41 #ifdef FAST_DOUBLE
  42       double_to_DF(DF_to_double(x1) + DF_to_double(x2), return ,
  43                    TRUE, TRUE, // Overflow und subnormale Zahl abfangen
  44                    FALSE, // kein Underflow mit Ergebnis +/- 0.0 möglich
  45                           // (nach Definition der subnormalen Zahlen)
  46                    FALSE, FALSE // keine Singularität, kein NaN als Ergebnis möglich
  47                   );
  48 #else
  49 #if (cl_word_size==64)
  50       // x1,x2 entpacken:
  51       var cl_signean sign1;
  52       var sintL exp1;
  53       var uint64 mant1;
  54       var cl_signean sign2;
  55       var sintL exp2;
  56       var uint64 mant2;
  57       DF_decode(x1, { return x2; }, sign1=,exp1=,mant1=);
  58       DF_decode(x2, { return x1; }, sign2=,exp2=,mant2=);
  59       var cl_DF max_x1_x2 = x1;
  60       if (exp1 < exp2)
  61         { max_x1_x2 = x2;
  62           swap(cl_signean, sign1,sign2);
  63           swap(sintL,      exp1 ,exp2 );
  64           swap(uint64,     mant1,mant2);
  65         }
  66       // Nun ist exp1>=exp2.
  67       var uintL expdiff = exp1 - exp2; // Exponentendifferenz
  68       if (expdiff >= DF_mant_len+3) // >= 52+3 ?
  69         { return max_x1_x2; }
  70       mant1 = mant1 << 3; mant2 = mant2 << 3;
  71       // Nun 2^(DF_mant_len+3) <= mant1,mant2 < 2^(DF_mant_len+4).
  72       {var uint64 mant2_last = mant2 & (bit(expdiff)-1); // letzte expdiff Bits von mant2
  73        mant2 = mant2 >> expdiff; if (!(mant2_last==0)) { mant2 |= bit(0); }
  74       }
  75       // mant2 = um expdiff Bits nach rechts geschobene und gerundete Mantisse
  76       // von x2.
  77       if (!(sign1==sign2))
  78         // verschiedene Vorzeichen -> Mantissen subtrahieren
  79         { if (mant1 > mant2) { mant1 = mant1 - mant2; goto norm_2; }
  80           if (mant1 == mant2) // Ergebnis 0 ?
  81             { return cl_DF_0; }
  82           // negatives Subtraktionsergebnis
  83           mant1 = mant2 - mant1; sign1 = sign2; goto norm_2;
  84         }
  85         else
  86         // gleiche Vorzeichen -> Mantissen addieren
  87         { mant1 = mant1 + mant2; }
  88       // mant1 = Ergebnis-Mantisse >0, sign1 = Ergebnis-Vorzeichen,
  89       // exp1 = Ergebnis-Exponent.
  90       // Außerdem: Bei expdiff=0,1 sind die zwei letzten Bits von mant1 Null,
  91       // bei expdiff>=2 ist mant1 >= 2^(DF_mant_len+2).
  92       // Stets ist mant1 < 2^(DF_mant_len+5). (Daher werden die 2 Rundungsbits
  93       // nachher um höchstens eine Position nach links geschoben werden.)
  94       // [Knuth, S.201, leicht modifiziert:
  95       //   N1. m>=1 -> goto N4.
  96       //   N2. [Hier m<1] m>=1/2 -> goto N5.
  97       //       N3. m:=2*m, e:=e-1, goto N2.
  98       //   N4. [Hier 1<=m<2] m:=m/2, e:=e+1.
  99       //   N5. [Hier 1/2<=m<1] Runde m auf 53 Bits hinterm Komma.
 100       //       Falls hierdurch m=1 geworden, setze m:=m/2, e:=e+1.
 101       // ]
 102       // Bei uns ist m=mant1/2^(DF_mant_len+4),
 103       // ab Schritt N5 ist m=mant1/2^(DF_mant_len+1).
 104       norm_1: // [Knuth, S.201, Schritt N1]
 105       if (mant1 >= bit(DF_mant_len+4)) goto norm_4;
 106       norm_2: // [Knuth, S.201, Schritt N2]
 107               // Hier ist mant1 < 2^(DF_mant_len+4)
 108       if (mant1 >= bit(DF_mant_len+3)) goto norm_5;
 109       // [Knuth, S.201, Schritt N3]
 110       mant1 = mant1 << 1; exp1 = exp1-1; // Mantisse links schieben
 111       goto norm_2;
 112       norm_4: // [Knuth, S.201, Schritt N4]
 113               // Hier ist 2^(DF_mant_len+4) <= mant1 < 2^(DF_mant_len+5)
 114       exp1 = exp1+1;
 115       mant1 = (mant1>>1) | (mant1 & bit(0)); // Mantisse rechts schieben
 116       norm_5: // [Knuth, S.201, Schritt N5]
 117               // Hier ist 2^(DF_mant_len+3) <= mant1 < 2^(DF_mant_len+4)
 118       // Auf DF_mant_len echte Mantissenbits runden, d.h. rechte 3 Bits
 119       // wegrunden, und dabei mant1 um 3 Bits nach rechts schieben:
 120       {var uint64 rounding_bits = mant1 & (bit(3)-1);
 121        mant1 = mant1 >> 3;
 122        if ( (rounding_bits < bit(2)) // 000,001,010,011 werden abgerundet
 123             || ( (rounding_bits == bit(2)) // 100 (genau halbzahlig)
 124                  && ((mant1 & bit(0)) ==0) // -> round-to-even
 125           )    )
 126          // abrunden
 127          {}
 128          else
 129          // aufrunden
 130          { mant1 = mant1+1;
 131            if (mant1 >= bit(DF_mant_len+1))
 132              // Bei Überlauf während der Rundung nochmals rechts schieben
 133              // (Runden ist hier überflüssig):
 134              { mant1 = mant1>>1; exp1 = exp1+1; } // Mantisse rechts schieben
 135          }
 136       }// Runden fertig
 137       return encode_DF(sign1,exp1,mant1);
 138 #else
 139       // x1,x2 entpacken:
 140       var cl_signean sign1;
 141       var sintL exp1;
 142       var uintL manthi1;
 143       var uintL mantlo1;
 144       var cl_signean sign2;
 145       var sintL exp2;
 146       var uintL manthi2;
 147       var uintL mantlo2;
 148       DF_decode2(x1, { return x2; }, sign1=,exp1=,manthi1=,mantlo1=);
 149       DF_decode2(x2, { return x1; }, sign2=,exp2=,manthi2=,mantlo2=);
 150       var cl_DF max_x1_x2 = x1;
 151       if (exp1 < exp2)
 152         { max_x1_x2 = x2;
 153           swap(cl_signean, sign1,sign2);
 154           swap(sintL,      exp1 ,exp2 );
 155           swap(uintL,      manthi1,manthi2);
 156           swap(uintL,      mantlo1,mantlo2);
 157         }
 158       // Nun ist exp1>=exp2.
 159       var uintL expdiff = exp1 - exp2; // Exponentendifferenz
 160       if (expdiff >= DF_mant_len+3) // >= 52+3 ?
 161         { return max_x1_x2; }
 162       manthi1 = (manthi1 << 3) | (mantlo1 >> (32-3)); mantlo1 = mantlo1 << 3;
 163       manthi2 = (manthi2 << 3) | (mantlo2 >> (32-3)); mantlo2 = mantlo2 << 3;
 164       // Nun 2^(DF_mant_len+3) <= mant1,mant2 < 2^(DF_mant_len+4).
 165       if (expdiff<32)
 166         {if (!(expdiff==0))
 167            {var uintL mant2_last = mantlo2 & (bit(expdiff)-1); // letzte expdiff Bits von mant2
 168             mantlo2 = (mantlo2 >> expdiff) | (manthi2 << (32-expdiff));
 169             manthi2 = manthi2 >> expdiff;
 170             if (!(mant2_last==0)) { mantlo2 |= bit(0); }
 171         }  }
 172         else
 173         {var uintL mant2_last = (manthi2 & (bit(expdiff-32)-1)) | mantlo2; // letzte expdiff Bits von mant2
 174          mantlo2 = manthi2 >> (expdiff-32); manthi2 = 0;
 175          if (!(mant2_last==0)) { mantlo2 |= bit(0); }
 176         }
 177       // mant2 = um expdiff Bits nach rechts geschobene und gerundete Mantisse
 178       // von x2.
 179       if (!(sign1==sign2))
 180         // verschiedene Vorzeichen -> Mantissen subtrahieren
 181         { if (manthi1 > manthi2)
 182             { manthi1 = manthi1 - manthi2;
 183               if (mantlo1 < mantlo2) { manthi1 -= 1; }
 184               mantlo1 = mantlo1 - mantlo2;
 185               goto norm_2;
 186             }
 187           if (manthi1 == manthi2)
 188             { if (mantlo1 > mantlo2)
 189                 { manthi1 = 0; mantlo1 = mantlo1 - mantlo2; goto norm_2; }
 190               if (mantlo1 == mantlo2) // Ergebnis 0 ?
 191                 { return cl_DF_0; }
 192             }
 193           // Hier ((manthi1 < manthi2) || ((manthi1 == manthi2) && (mantlo1 < mantlo2))).
 194           // negatives Subtraktionsergebnis
 195           manthi1 = manthi2 - manthi1;
 196           if (mantlo2 < mantlo1) { manthi1 -= 1; }
 197           mantlo1 = mantlo2 - mantlo1;
 198           sign1 = sign2;
 199           goto norm_2;
 200         }
 201         else
 202         // gleiche Vorzeichen -> Mantissen addieren
 203         { manthi1 = manthi1 + manthi2;
 204           if ((mantlo1 = mantlo1 + mantlo2) < mantlo2) { manthi1 += 1; }
 205         }
 206       // mant1 = Ergebnis-Mantisse >0, sign1 = Ergebnis-Vorzeichen,
 207       // exp1 = Ergebnis-Exponent.
 208       // Außerdem: Bei expdiff=0,1 sind die zwei letzten Bits von mant1 Null,
 209       // bei expdiff>=2 ist mant1 >= 2^(DF_mant_len+2).
 210       // Stets ist mant1 < 2^(DF_mant_len+5). (Daher werden die 2 Rundungsbits
 211       // nachher um höchstens eine Position nach links geschoben werden.)
 212       // [Knuth, S.201, leicht modifiziert:
 213       //   N1. m>=1 -> goto N4.
 214       //   N2. [Hier m<1] m>=1/2 -> goto N5.
 215       //       N3. m:=2*m, e:=e-1, goto N2.
 216       //   N4. [Hier 1<=m<2] m:=m/2, e:=e+1.
 217       //   N5. [Hier 1/2<=m<1] Runde m auf 53 Bits hinterm Komma.
 218       //       Falls hierdurch m=1 geworden, setze m:=m/2, e:=e+1.
 219       // ]
 220       // Bei uns ist m=mant1/2^(DF_mant_len+4),
 221       // ab Schritt N5 ist m=mant1/2^(DF_mant_len+1).
 222       norm_1: // [Knuth, S.201, Schritt N1]
 223       if (manthi1 >= bit(DF_mant_len-32+4)) goto norm_4;
 224       norm_2: // [Knuth, S.201, Schritt N2]
 225               // Hier ist mant1 < 2^(DF_mant_len+4)
 226       if (manthi1 >= bit(DF_mant_len-32+3)) goto norm_5;
 227       // [Knuth, S.201, Schritt N3]
 228       manthi1 = (manthi1 << 1) | (mantlo1 >> 31); // Mantisse links schieben
 229       mantlo1 = mantlo1 << 1;
 230       exp1 = exp1-1;
 231       goto norm_2;
 232       norm_4: // [Knuth, S.201, Schritt N4]
 233               // Hier ist 2^(DF_mant_len+4) <= mant1 < 2^(DF_mant_len+5)
 234       exp1 = exp1+1;
 235       mantlo1 = (mantlo1 >> 1) | (manthi1 << 31) | (mantlo1 & bit(0)); // Mantisse rechts schieben
 236       manthi1 = (manthi1 >> 1);
 237       norm_5: // [Knuth, S.201, Schritt N5]
 238               // Hier ist 2^(DF_mant_len+3) <= mant1 < 2^(DF_mant_len+4)
 239       // Auf DF_mant_len echte Mantissenbits runden, d.h. rechte 3 Bits
 240       // wegrunden, und dabei mant1 um 3 Bits nach rechts schieben:
 241       {var uintL rounding_bits = mantlo1 & (bit(3)-1);
 242        mantlo1 = (mantlo1 >> 3) | (manthi1 << (32-3)); manthi1 = manthi1 >> 3;
 243        if ( (rounding_bits < bit(2)) // 000,001,010,011 werden abgerundet
 244             || ( (rounding_bits == bit(2)) // 100 (genau halbzahlig)
 245                  && ((mantlo1 & bit(0)) ==0) // -> round-to-even
 246           )    )
 247          // abrunden
 248          {}
 249          else
 250          // aufrunden
 251          { mantlo1 = mantlo1+1;
 252            if (mantlo1==0)
 253              { manthi1 = manthi1+1;
 254                if (manthi1 >= bit(DF_mant_len-32+1))
 255                  // Bei Überlauf während der Rundung nochmals rechts schieben
 256                  // (Runden ist hier überflüssig):
 257                  { manthi1 = manthi1>>1; exp1 = exp1+1; } // Mantisse rechts schieben
 258          }   }
 259       }// Runden fertig
 260       return encode_DF(sign1,exp1,manthi1,mantlo1);
 261 #endif
 262 #endif
 263 }
 264
 265 }  // namespace cln