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