src/float/sfloat/elem/cl_SF_plus.cc

   1 // binary operator +
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/sfloat.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_SF.h"
  13 #include "cl_xmacros.h"
  14
  15 namespace cln {
  16
  17 const cl_SF operator+ (const cl_SF& x1, const cl_SF& x2)
  18 {
  19 // Methode (nach [Knuth, II, Seminumerical Algorithms, Abschnitt 4.2.1., S.200]):
  20 // x1=0.0 -> Ergebnis x2.
  21 // x2=0.0 -> Ergebnis x1.
  22 // Falls e1<e2, vertausche x1 und x2.
  23 // Also e1 >= e2.
  24 // Falls e1 - e2 >= 16 + 3, Ergebnis x1.
  25 // Schiebe beide Mantissen um 3 Bits nach links (Vorbereitung der Rundung:
  26 //   Bei e1-e2=0,1 ist keine Rundung nötig, bei e1-e2>1 ist der Exponent des
  27 //   Ergebnisses =e1-1, =e1 oder =e1+1. Brauche daher 1 Schutzbit und zwei
  28 //   Rundungsbits: 00 exakt, 01 1.Hälfte, 10 exakte Mitte, 11 2.Hälfte.)
  29 // Schiebe die Mantisse von x2 um e0-e1 Bits nach rechts. (Dabei die Rundung
  30 // ausführen: Bit 0 ist das logische Oder der Bits 0,-1,-2,...)
  31 // Falls x1,x2 selbes Vorzeichen haben: Addiere dieses zur Mantisse von x1.
  32 // Falls x1,x2 verschiedenes Vorzeichen haben: Subtrahiere dieses von der
  33 //   Mantisse von x1. <0 -> (Es war e1=e2) Vertausche die Vorzeichen, negiere.
  34 //                    =0 -> Ergebnis 0.0
  35 // Exponent ist e1.
  36 // Normalisiere, fertig.
  37       // x1,x2 entpacken:
  38       var cl_signean sign1;
  39       var sintL exp1;
  40       var uintL mant1;
  41       var cl_signean sign2;
  42       var sintL exp2;
  43       var uintL mant2;
  44       SF_decode(x1, { return x2; }, sign1=,exp1=,mant1=);
  45       SF_decode(x2, { return x1; }, sign2=,exp2=,mant2=);
  46       var cl_uint max_x1_x2 = x1.word;
  47       if (exp1 < exp2)
  48         { max_x1_x2 = x2.word;
  49           swap(cl_signean, sign1,sign2);
  50           swap(sintL,      exp1 ,exp2 );
  51           swap(uintL,      mant1,mant2);
  52         }
  53       // Nun ist exp1>=exp2.
  54      {var uintL expdiff = exp1 - exp2; // Exponentendifferenz
  55       if (expdiff >= SF_mant_len+3) // >= 16+3 ?
  56         { return cl_SF_from_word(max_x1_x2); }
  57       mant1 = mant1 << 3; mant2 = mant2 << 3;
  58       // Nun 2^(SF_mant_len+3) <= mant1,mant2 < 2^(SF_mant_len+4).
  59       {var uintL mant2_last = mant2 & (bit(expdiff)-1); // letzte expdiff Bits von mant2
  60        mant2 = mant2 >> expdiff; if (!(mant2_last==0)) { mant2 |= bit(0); }
  61       }
  62       // mant2 = um expdiff Bits nach rechts geschobene und gerundete Mantisse
  63       // von x2.
  64       if ((x1.word ^ x2.word) & bit(SF_sign_shift))
  65         // verschiedene Vorzeichen -> Mantissen subtrahieren
  66         { if (mant1 > mant2) { mant1 = mant1 - mant2; goto norm_2; }
  67           if (mant1 == mant2) // Ergebnis 0 ?
  68             { return SF_0; }
  69           // negatives Subtraktionsergebnis
  70           mant1 = mant2 - mant1; sign1 = sign2; goto norm_2;
  71         }
  72         else
  73         // gleiche Vorzeichen -> Mantissen addieren
  74         { mant1 = mant1 + mant2; }
  75       // mant1 = Ergebnis-Mantisse >0, sign1 = Ergebnis-Vorzeichen,
  76       // exp1 = Ergebnis-Exponent.
  77       // Außerdem: Bei expdiff=0,1 sind die zwei letzten Bits von mant1 Null,
  78       // bei expdiff>=2 ist mant1 >= 2^(SF_mant_len+2).
  79       // Stets ist mant1 < 2^(SF_mant_len+5). (Daher werden die 2 Rundungsbits
  80       // nachher um höchstens eine Position nach links geschoben werden.)
  81       // [Knuth, S.201, leicht modifiziert:
  82       //   N1. m>=1 -> goto N4.
  83       //   N2. [Hier m<1] m>=1/2 -> goto N5.
  84       //       N3. m:=2*m, e:=e-1, goto N2.
  85       //   N4. [Hier 1<=m<2] m:=m/2, e:=e+1.
  86       //   N5. [Hier 1/2<=m<1] Runde m auf 17 Bits hinterm Komma.
  87       //       Falls hierdurch m=1 geworden, setze m:=m/2, e:=e+1.
  88       // ]
  89       // Bei uns ist m=mant1/2^(SF_mant_len+4),
  90       // ab Schritt N5 ist m=mant1/2^(SF_mant_len+1).
  91       norm_1: // [Knuth, S.201, Schritt N1]
  92       if (mant1 >= bit(SF_mant_len+4)) goto norm_4;
  93       norm_2: // [Knuth, S.201, Schritt N2]
  94               // Hier ist mant1 < 2^(SF_mant_len+4)
  95       if (mant1 >= bit(SF_mant_len+3)) goto norm_5;
  96       // [Knuth, S.201, Schritt N3]
  97       mant1 = mant1 << 1; exp1 = exp1-1; // Mantisse links schieben
  98       goto norm_2;
  99       norm_4: // [Knuth, S.201, Schritt N4]
 100               // Hier ist 2^(SF_mant_len+4) <= mant1 < 2^(SF_mant_len+5)
 101       exp1 = exp1+1;
 102       mant1 = (mant1>>1) | (mant1 & bit(0)); // Mantisse rechts schieben
 103       norm_5: // [Knuth, S.201, Schritt N5]
 104               // Hier ist 2^(SF_mant_len+3) <= mant1 < 2^(SF_mant_len+4)
 105       // Auf SF_mant_len echte Mantissenbits runden, d.h. rechte 3 Bits
 106       // wegrunden, und dabei mant1 um 3 Bits nach rechts schieben:
 107       {var uintL rounding_bits = mant1 & (bit(3)-1);
 108        mant1 = mant1 >> 3;
 109        if ( (rounding_bits < bit(2)) // 000,001,010,011 werden abgerundet
 110             || ( (rounding_bits == bit(2)) // 100 (genau halbzahlig)
 111                  && ((mant1 & bit(0)) ==0) // -> round-to-even
 112           )    )
 113          // abrunden
 114          {}
 115          else
 116          // aufrunden
 117          { mant1 = mant1+1;
 118            if (mant1 >= bit(SF_mant_len+1))
 119              // Bei Überlauf während der Rundung nochmals rechts schieben
 120              // (Runden ist hier überflüssig):
 121              { mant1 = mant1>>1; exp1 = exp1+1; } // Mantisse rechts schieben
 122          }
 123       }// Runden fertig
 124       return encode_SF(sign1,exp1,mant1);
 125      }
 126 }
 127
 128 }  // namespace cln