src/float/lfloat/elem/cl_LF_1plus.cc

   1 // LF_LF_plus_LF().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_LF.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_LF_impl.h"
  13 #include "cl_DS.h"
  14 #include "cl_F.h"
  15 #include "cl_xmacros.h"
  16
  17 namespace cln {
  18
  19 const cl_LF LF_LF_plus_LF (const cl_LF& arg1, const cl_LF& arg2)
  20 {
  21 // Methode (nach [Knuth, II, Seminumerical Algorithms, Abschnitt 4.2.1., S.200]):
  22 // Falls e1<e2, vertausche x1 und x2.
  23 // Also e1 >= e2.
  24 // Falls e2=0, also x2=0.0, Ergebnis x1.
  25 // Falls e1 - e2 >= 16n+2, Ergebnis x1.
  26 // Erweitere die Mantissen rechts um 3 Bits (Bit -1 als Schutzbit, Bits -2,-3
  27 //   als Rundungsbits: 00 exakt, 01 1.Hälfte, 10 exakte Mitte, 11 2.Hälfte.)
  28 // Schiebe die Mantisse von x2 um e0-e1 Bits nach rechts. (Dabei die Rundung
  29 // ausführen: Bit -3 ist das logische Oder der Bits -3,-4,-5,...)
  30 // Falls x1,x2 selbes Vorzeichen haben: Addiere dieses zur Mantisse von x1.
  31 // Falls x1,x2 verschiedenes Vorzeichen haben: Subtrahiere dieses von der
  32 //   Mantisse von x1. <0 -> (Es war e1=e2) Vertausche die Vorzeichen, negiere.
  33 //                    =0 -> Ergebnis 0.0
  34 // Exponent ist e1.
  35 // Normalisiere, fertig.
  36       var cl_LF x1 = arg1;
  37       var cl_LF x2 = arg2;
  38       var uintL uexp1 = TheLfloat(arg1)->expo;
  39       var uintL uexp2 = TheLfloat(arg2)->expo;
  40       if (uexp1 < uexp2)
  41         // x1 und x2 vertauschen
  42         { x1 = arg2; x2 = arg1; swap(uintL, uexp1,uexp2); }
  43       // uexp1 >= uexp2
  44       if (uexp2==0) { return x1; } // x2=0.0 -> x1 als Ergebnis
  45       var uintC len = TheLfloat(x1)->len; // Länge n von x1 und x2
  46       var uintL expdiff = uexp1-uexp2; // e1-e2
  47       if ((expdiff == 0) && (TheLfloat(x1)->sign != TheLfloat(x2)->sign))
  48         // verschiedene Vorzeichen, aber gleicher Exponent
  49         { // Vorzeichen des Ergebnisses festlegen:
  50           var cl_signean erg = // Mantissen (je len Digits) vergleichen
  51             compare_loop_msp(arrayMSDptr(TheLfloat(x1)->data,len),arrayMSDptr(TheLfloat(x2)->data,len),len);
  52           if (erg==0) // Mantissen gleich
  53             { return encode_LF0(len); } // Ergebnis 0.0
  54           if (erg<0) // |x1| < |x2|
  55             // x1 und x2 vertauschen, expdiff bleibt =0
  56             { x1.pointer = arg2.pointer; x2.pointer = arg1.pointer;
  57               swap(uintL, uexp1,uexp2);
  58             }
  59         }
  60       if (expdiff >= intDsize * (uintL)len + 2) // e1-e2 >= 16n+2 ?
  61         { return x1; } // ja -> x1 als Ergebnis
  62       // neues Long-Float allozieren:
  63       var Lfloat y = allocate_lfloat(len,uexp1,TheLfloat(x1)->sign);
  64       var uintL i = floor(expdiff,intDsize); // e1-e2 div 16 (>=0, <=n)
  65       var uintL j = expdiff % intDsize; // e1-e2 mod 16 (>=0, <16)
  66       // Mantisse von x2 muß um intDsize*i+j Bits nach rechts geschoben werden.
  67       var uintC x2_len = len - i; // n-i Digits von x2 gebraucht
  68       // x2_len Digits um j Bits nach rechts schieben und dabei kopieren:
  69       CL_ALLOCA_STACK;
  70       var uintD* x2_MSDptr;
  71       var uintD* x2_LSDptr;
  72       var uintD rounding_bits;
  73       num_stack_alloc(x2_len, x2_MSDptr=,x2_LSDptr=); // x2_len Digits Platz
  74       if (j==0)
  75         { copy_loop_msp(arrayMSDptr(TheLfloat(x2)->data,len),x2_MSDptr,x2_len); rounding_bits = 0; }
  76         else
  77         { rounding_bits = shiftrightcopy_loop_msp(arrayMSDptr(TheLfloat(x2)->data,len),x2_MSDptr,x2_len,j,0); }
  78       // x2_MSDptr/x2_len/x2_LSDptr sind die essentiellen Digits von x2.
  79       // rounding_bits enthält die letzten j herausgeschobenen Bits.
  80       // Aus rounding_bits und den nächsten i Digits die 3 Rundungsbits
  81       // (als Bits intDsize-1..intDsize-3 von rounding_bits) aufbauen:
  82       if (j>=2)
  83         // j>=2 -> Bits -1,-2 sind OK, Bit -3 bestimmen:
  84         { if ((rounding_bits & (bit(intDsize-3)-1)) ==0)
  85             { if (test_loop_msp(arrayMSDptr(TheLfloat(x2)->data,len) mspop x2_len,i))
  86                 { rounding_bits |= bit(intDsize-3); } // Rundungsbit -3 setzen
  87             }
  88             else
  89             { rounding_bits |= bit(intDsize-3); // Rundungsbit -3 setzen
  90               rounding_bits &= bitm(intDsize)-bit(intDsize-3); // andere Bits löschen
  91         }   }
  92         else
  93         // j<=3 -> Bits intDsize-4..0 von rounding_bits sind bereits Null.
  94         // nächstes und weitere i-1 Digits heranziehen:
  95         { if (i > 0) // i=0 -> Bits -1,-2,-3 sind OK.
  96             { var uintD* ptr = arrayMSDptr(TheLfloat(x2)->data,len) mspop x2_len;
  97               rounding_bits |= (mspref(ptr,0) >> j); // weitere relevante Bits des nächsten Digit dazu
  98               if ((rounding_bits & (bit(intDsize-3)-1)) ==0) // Alle Bits -3,-4,... =0 ?
  99                 { if (   (!((mspref(ptr,0) & (bit(3)-1)) ==0)) // j (<=3) untere Bits von ptr[0] alle =0 ?
 100                       || test_loop_msp(ptr mspop 1,i-1)
 101                      )
 102                     { rounding_bits |= bit(intDsize-3); } // Rundungsbit -3 setzen
 103                 }
 104                 else
 105                 { rounding_bits |= bit(intDsize-3); // Rundungsbit -3 setzen
 106                   rounding_bits &= bitm(intDsize)-bit(intDsize-3); // andere Bits löschen
 107         }   }   }
 108       // x2 liegt in verschobener Form in der UDS x2_MSDptr/x2_len/x2_LSDptr
 109       // vor, mit Rundungsbits in Bit intDsize-1..intDsize-3 von rounding_bits.
 110       {var uintD* y_mantMSDptr = arrayMSDptr(TheLfloat(y)->data,len);
 111        var uintD* y_mantLSDptr = arrayLSDptr(TheLfloat(y)->data,len);
 112        if (TheLfloat(x1)->sign == TheLfloat(x2)->sign)
 113          // gleiche Vorzeichen -> Mantissen addieren
 114          { // erst rechten Mantissenteil (x2_len Digits) durch Addition:
 115            var uintD carry =
 116              add_loop_lsp(arrayLSDptr(TheLfloat(x1)->data,len),x2_LSDptr,
 117                           y_mantLSDptr, x2_len
 118                          );
 119            // dann linken Mantissenteil (i Digits) direkt kopieren:
 120            var uintD* ptr =
 121              copy_loop_msp(arrayMSDptr(TheLfloat(x1)->data,len),y_mantMSDptr,i);
 122            // dann Übertrag vom rechten zum linken Mantissenteil addieren:
 123            if (!(carry==0))
 124              { if ( inc_loop_lsp(ptr,i) )
 125                  // Übertrag über das erste Digit hinaus
 126                  { // Exponent von y incrementieren:
 127                    if ( ++(TheLfloat(y)->expo) == LF_exp_high+1 ) { cl_error_floating_point_overflow(); }
 128                    // normalisiere durch Schieben um 1 Bit nach rechts:
 129                   {var uintD carry_rechts =
 130                      shift1right_loop_msp(y_mantMSDptr,len,~(uintD)0);
 131                    rounding_bits = rounding_bits>>1; // Rundungsbits mitschieben
 132                    if (!(carry_rechts==0)) { rounding_bits |= bit(intDsize-1); }
 133              }   }}
 134          }
 135          else
 136          // verschiedene Vorzeichen -> Mantissen subtrahieren
 137          { // erst rechten Mantissenteil (x2_len Digits) durch Subtraktion:
 138            rounding_bits = -rounding_bits;
 139            {var uintD carry =
 140               subx_loop_lsp(arrayLSDptr(TheLfloat(x1)->data,len),x2_LSDptr,
 141                             y_mantLSDptr, x2_len,
 142                             (rounding_bits==0 ? 0 : ~(uintD)0)
 143                            );
 144             // dann linken Mantissenteil (i Digits) direkt kopieren:
 145             var uintD* ptr =
 146               copy_loop_msp(arrayMSDptr(TheLfloat(x1)->data,len),y_mantMSDptr,i);
 147             // dann Übertrag des rechten vom linken Mantissenteil subtrahieren:
 148             if (!(carry==0))
 149               { if ( dec_loop_lsp(ptr,i) )
 150                   // Übertrag über das erste Digit hinaus, also e1=e2
 151                   { NOTREACHED } // diesen Fall haben wir schon behandelt
 152               }
 153            }
 154            // UDS y_mantMSDptr/len/y_mantLSDptr/rounding_bits normalisieren:
 155            {var uintD* ptr = y_mantMSDptr;
 156             var uintL k = 0;
 157             var uintC count;
 158             dotimesC(count,len,
 159               { if (!(mspref(ptr,0)==0)) goto nonzero_found;
 160                 ptr = ptr mspop 1; k++;
 161               });
 162             if (!(rounding_bits==0)) goto nonzero_found;
 163             // Die UDS ist ganz Null. Also war e1=e2, keine Rundungsbits.
 164             { NOTREACHED } // diesen Fall haben wir schon behandelt
 165             nonzero_found: // Digit /=0 gefunden
 166             // UDS von ptr nach y_mantMSDptr um k Digits nach unten kopieren:
 167             if (k>0)
 168               // mindestens ein führendes Nulldigit. Also war e1-e2 = 0 oder 1.
 169               { ptr = copy_loop_msp(ptr,y_mantMSDptr,len-k); // len-k Digits verschieben
 170                 msprefnext(ptr) = rounding_bits; // Rundungsbits als weiteres Digit
 171                 clear_loop_msp(ptr,k-1); // dann k-1 Nulldigits
 172                 rounding_bits = 0; // und keine weiteren Rundungsbits
 173                 // Exponenten um intDsize*k erniedrigen:
 174                 k = intDsize*k;
 175                {var uintL uexp = TheLfloat(y)->expo;
 176                 #if !(LF_exp_low==1)
 177                 if (uexp < k+LF_exp_low)
 178                 #else
 179                 if (uexp <= k)
 180                 #endif
 181                   { if (underflow_allowed())
 182                       { cl_error_floating_point_underflow(); }
 183                       else
 184                       { return encode_LF0(len); } // Ergebnis 0.0
 185                   }
 186                 TheLfloat(y)->expo = uexp - k;
 187               }}
 188            }
 189            // NUDS y_mantMSDptr/len/y_mantLSDptr/rounding_bits normalisieren:
 190            {var uintL s;
 191             integerlengthD(mspref(y_mantMSDptr,0), s = intDsize - );
 192             // s = Anzahl der führenden Nullbits im ersten Word (>=0, <intDsize)
 193             if (s > 0)
 194               { // Muß die NUDS y_mantMSDptr/len/y_mantLSDptr/rounding_bits
 195                 // um s Bits nach links schieben.
 196                 // (Bei e1-e2>1 ist dabei zwangsläufig s=1.)
 197                 if (s==1)
 198                   { shift1left_loop_lsp(y_mantLSDptr,len);
 199                     if (rounding_bits & bit(intDsize-1))
 200                       { lspref(y_mantLSDptr,0) |= bit(0); }
 201                     rounding_bits = rounding_bits << 1;
 202                   }
 203                   else // s>1, also e1-e2 <= 1 <= s.
 204                   { shiftleft_loop_lsp(y_mantLSDptr,len,s,rounding_bits>>(intDsize-s));
 205                     rounding_bits = 0; // = rounding_bits << s;
 206                   }
 207                 // Exponenten um s erniedrigen:
 208                {var uintL uexp = TheLfloat(y)->expo;
 209                 #if !(LF_exp_low==1)
 210                 if (uexp < s+LF_exp_low)
 211                 #else
 212                 if (uexp <= s)
 213                 #endif
 214                   { if (underflow_allowed())
 215                       { cl_error_floating_point_underflow(); }
 216                       else
 217                       { return encode_LF0(len); } // Ergebnis 0.0
 218                   }
 219                 TheLfloat(y)->expo = uexp - s;
 220               }}
 221          } }
 222        // Hier enthält rounding_bits Bit -1 als Bit intDsize-1, Bit -2 als
 223        // Bit intDsize-2, Bit -3 als Oder(Bits intDsize-3..0) !
 224        // Runden. Dazu rounding_bits inspizieren:
 225        if ((rounding_bits & bit(intDsize-1)) ==0) goto ab; // Bit -1 gelöscht -> abrunden
 226        rounding_bits = rounding_bits<<1; // Bits -2,-3
 227        if (!(rounding_bits==0)) goto auf; // Bit -2 oder Bit -3 gesetzt -> aufrunden
 228        // round-to-even:
 229        if ((lspref(y_mantLSDptr,0) & bit(0)) ==0) goto ab;
 230        auf: // aufrunden
 231          if ( inc_loop_lsp(y_mantLSDptr,len) )
 232            { // Übertrag durchs Aufrunden
 233              mspref(y_mantMSDptr,0) = bit(intDsize-1); // Mantisse := 10...0
 234              // Exponent erhöhen:
 235              if (++(TheLfloat(y)->expo) == LF_exp_high+1) { cl_error_floating_point_overflow(); }
 236            }
 237        ab: // abrunden
 238          ;
 239       }
 240       // y fertig.
 241       return y;
 242 }
 243
 244 }  // namespace cln