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