src/float/lfloat/elem/cl_LF_div.cc

   1 // binary operator /
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_lfloat.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_LF.h"
  13 #include "cl_LF_impl.h"
  14 #include "cl_DS.h"
  15 #include "cl_F.h"
  16 #include "cl_N.h"
  17
  18 // Workaround gcc-2.7.0 bug on i386.
  19 #if defined(__GNUC__)
  20   #if (__GNUC__ == 2)
  21     #if (__GNUC_MINOR__ == 7)
  22       #define workaround_gcc_bug()  *&uexp1 = *&uexp1;
  23     #endif
  24   #endif
  25 #endif
  26 #ifndef workaround_gcc_bug
  27   #define workaround_gcc_bug()
  28 #endif
  29
  30 const cl_LF operator/ (const cl_LF& x1, const cl_LF& x2)
  31 {
  32 // Methode:
  33 // x2 = 0.0 -> Error
  34 // x1 = 0.0 -> Ergebnis 0.0
  35 // Sonst:
  36 // Ergebnis-Vorzeichen = xor der beiden Vorzeichen von x1 und x2
  37 // Ergebnis-Exponent = Differenz der beiden Exponenten von x1 und x2
  38 // Ergebnis-Mantisse = Mantisse mant1 / Mantisse mant2, gerundet.
  39 //   mant1/mant2 > 1/2, mant1/mant2 < 2;
  40 //   nach Rundung mant1/mant2 >=1/2, <=2*mant1<2.
  41 //   Bei mant1/mant2 >=1 brauche 16n-1 Nachkommabits,
  42 //   bei mant1/mant2 <1 brauche 16n Nachkommabits.
  43 //   Fürs Runden: brauche ein Rundungsbit (Rest gibt an, ob exakt).
  44 //   Brauche daher insgesamt 16n+1 Nachkommabits von mant1/mant2.
  45 //   Dividiere daher (als Unsigned Integers)
  46 //     2^16(n+1)*(2^16n*m0) durch (2^16n*m1).
  47 //   Falls der Quotient >=2^16(n+1) ist, schiebe ihn um 1 Bit nach rechts,
  48 //     erhöhe den Exponenten um 1 und runde das letzte Digit weg.
  49 //   Falls der Quotient <2^16(n+1) ist, runde das letzte Digit weg. Bei rounding
  50 //     overflow schiebe um 1 Bit nach rechts und erhöhe den Exponenten um 1.
  51       var uintC len1 = TheLfloat(x1)->len;
  52       var uintC len2 = TheLfloat(x2)->len;
  53       var uintC len = (len1 < len2 ? len1 : len2); // min. Länge n von x1 und x2
  54       var uintL uexp2 = TheLfloat(x2)->expo;
  55       if (uexp2==0) { cl_error_division_by_0(); } // x2=0.0 -> Error
  56       var uintL uexp1 = TheLfloat(x1)->expo;
  57       if (uexp1==0) // x1=0.0 -> Ergebnis 0.0
  58         { if (len < len1) return shorten(x1,len); else return x1; }
  59       // Exponenten subtrahieren:
  60       // (uexp1-LF_exp_mid) - (uexp2-LF_exp_mid) = (uexp1-uexp2+LF_exp_mid)-LF_exp_mid
  61       if (uexp1 >= uexp2)
  62         { uexp1 = uexp1 - uexp2; // kein Carry
  63           workaround_gcc_bug();
  64           if (uexp1 > LF_exp_high-LF_exp_mid) { cl_error_floating_point_overflow(); }
  65           uexp1 = uexp1 + LF_exp_mid;
  66         }
  67         else
  68         { uexp1 = uexp1 - uexp2; // Carry
  69           workaround_gcc_bug();
  70           if (uexp1 < (uintL)(LF_exp_low-1-LF_exp_mid))
  71             { if (underflow_allowed())
  72                 { cl_error_floating_point_underflow(); }
  73                 else
  74                 { return encode_LF0(len); } // Ergebnis 0.0
  75             }
  76           uexp1 = uexp1 + LF_exp_mid;
  77         }
  78       // Nun ist LF_exp_low-1 <= uexp1 <= LF_exp_high.
  79       // neues Long-Float allozieren:
  80       var Lfloat y = allocate_lfloat(len,uexp1,
  81                                      TheLfloat(x1)->sign ^ TheLfloat(x2)->sign // Vorzeichen kombinieren
  82                                     );
  83       // Nenner bilden:
  84       var uintL n_len;
  85       n_len = len2;
  86 #ifndef CL_LF_PEDANTIC
  87       if (n_len > len) { n_len = len+1; }
  88 #endif
  89       // Zähler bilden:
  90       CL_ALLOCA_STACK;
  91       var uintD* z_MSDptr;
  92       var uintL z_len;
  93       var uintD* z_LSDptr;
  94       z_len = n_len + len + 1;
  95       num_stack_alloc(z_len, z_MSDptr=,z_LSDptr=);
  96       if (z_len > len1)
  97         { var uintD* ptr =
  98             copy_loop_msp(arrayMSDptr(TheLfloat(x1)->data,len1),z_MSDptr,len1); // n Digits kopieren
  99           clear_loop_msp(ptr,z_len-len1); // und n+1 Null-Digits
 100         }
 101         else
 102         { copy_loop_msp(arrayMSDptr(TheLfloat(x1)->data,len1),z_MSDptr,z_len); }
 103       // Quotienten bilden: 2n+1-Digit-Zahl durch n-Digit-Zahl dividieren
 104       {var DS q;
 105        var DS r;
 106        {var uintD* x2_mantMSDptr = arrayMSDptr(TheLfloat(x2)->data,len2);
 107         UDS_divide(z_MSDptr,z_len,z_LSDptr,
 108                    x2_mantMSDptr,n_len,x2_mantMSDptr mspop n_len,
 109                    &q, &r
 110                   );
 111        }
 112        // q ist der Quotient mit n+1 oder n+2 Digits, r der Rest.
 113        if (q.len > len+1)
 114          // Quotient hat n+2 Digits -> um 1 Bit nach rechts schieben:
 115          { var uintD* y_mantMSDptr = arrayMSDptr(TheLfloat(y)->data,len);
 116            var uintD carry_rechts =
 117              shiftrightcopy_loop_msp(q.MSDptr mspop 1,y_mantMSDptr,len,1,
 118                                      /* carry links = mspref(q.MSDptr,0) = 1 */ 1 );
 119            // Exponenten incrementieren:
 120            if (++(TheLfloat(y)->expo) == LF_exp_high+1) { cl_error_floating_point_overflow(); }
 121            // Runden:
 122            if ( (carry_rechts == 0) // herausgeschobenes Bit =0 -> abrunden
 123                 || ( (lspref(q.LSDptr,0)==0) // =1 und weitere Bits >0 oder Rest >0 -> aufrunden
 124                      && (r.len==0)
 125                      // round-to-even
 126                      && ((lspref(q.LSDptr,1) & bit(1)) ==0)
 127               )    )
 128              // abrunden
 129              {}
 130              else
 131              // aufrunden
 132              { inc_loop_lsp(y_mantMSDptr mspop len,len); }
 133          }
 134          else
 135          // Quotient hat n+1 Digits -> nur kopieren:
 136          { var uintD* y_mantMSDptr = arrayMSDptr(TheLfloat(y)->data,len);
 137            copy_loop_msp(q.MSDptr,y_mantMSDptr,len);
 138            // Runden:
 139            if ( ((sintD)lspref(q.LSDptr,0) >= 0) // nächstes Bit =0 -> abrunden
 140                 || ( ((lspref(q.LSDptr,0) & ((uintD)bit(intDsize-1)-1)) ==0) // =1 und weitere Bits >0 oder Rest >0 -> aufrunden
 141                      && (r.len==0)
 142                      // round-to-even
 143                      && ((lspref(q.LSDptr,1) & bit(0)) ==0)
 144               )    )
 145              // abrunden
 146              {}
 147              else
 148              // aufrunden
 149              { if ( inc_loop_lsp(y_mantMSDptr mspop len,len) )
 150                  // Übertrag durchs Aufrunden
 151                  { mspref(y_mantMSDptr,0) = bit(intDsize-1); // Mantisse := 10...0
 152                    // Exponenten incrementieren:
 153                    if (++(TheLfloat(y)->expo) == LF_exp_high+1) { cl_error_floating_point_overflow(); }
 154              }   }
 155          }
 156       }
 157       // LF_exp_low <= exp <= LF_exp_high sicherstellen:
 158       if (TheLfloat(y)->expo == LF_exp_low-1)
 159         { if (underflow_allowed())
 160             { cl_error_floating_point_underflow(); }
 161             else
 162             { return encode_LF0(len); } // Ergebnis 0.0
 163         }
 164       return y;
 165 }
 166 // Bit complexity (N := max(length(x1),length(x2))): O(M(N)).
 167