src/float/lfloat/elem/cl_LF_I_div.cc

   1 // cl_LF_I_div().
   2
   3 // General includes.
   4 #include "base/cl_sysdep.h"
   5
   6 // Specification.
   7 #include "float/lfloat/cl_LF.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cln/lfloat.h"
  13 #include "float/lfloat/cl_LF_impl.h"
  14 #include "cln/integer.h"
  15 #include "integer/cl_I.h"
  16 #include "base/digitseq/cl_DS.h"
  17 #include "float/cl_F.h"
  18 #include "base/cl_N.h"
  19
  20 namespace cln {
  21
  22 const cl_LF cl_LF_I_div (const cl_LF& x, const cl_I& y)
  23 {
  24 // Method:
  25 // If x = 0.0, return x.
  26 // If y is longer than x, convert y to a float and divide.
  27 // Else divide the mantissa of x by the absolute value of y, then round.
  28         if (TheLfloat(x)->expo == 0) {
  29                 if (zerop(y))
  30                         throw division_by_0_exception();
  31                 else
  32                         return x;
  33         }
  34         var cl_signean sign = -(cl_signean)minusp(y); // Vorzeichen von y
  35         var cl_I abs_y = (sign==0 ? y : -y);
  36         var uintC y_exp = integer_length(abs_y);
  37         var uintC len = TheLfloat(x)->len;
  38 #ifndef CL_LF_PEDANTIC
  39         if (ceiling(y_exp,intDsize) > len)
  40                 return x / cl_I_to_LF(y,len);
  41 #endif
  42         // x länger als y, direkt dividieren.
  43         CL_ALLOCA_STACK;
  44         var const uintD* y_MSDptr;
  45         var uintC y_len;
  46         var const uintD* y_LSDptr;
  47         I_to_NDS_nocopy(abs_y, y_MSDptr=,y_len=,y_LSDptr=,false,); // NDS zu y bilden, y_len>0
  48         // y nicht zu einer NUDS normalisieren! (Damit ein Bit Spielraum ist.)
  49         // Zähler bilden: x * 2^(intDsize*y_len)
  50         var uintD* z_MSDptr;
  51         var uintC z_len;
  52         var uintD* z_LSDptr;
  53         z_len = len + y_len;
  54         num_stack_alloc(z_len, z_MSDptr=,z_LSDptr=);
  55         { var uintD* ptr = copy_loop_msp(arrayMSDptr(TheLfloat(x)->data,len),z_MSDptr,len); // len Digits
  56           clear_loop_msp(ptr,y_len); // und y_len Null-Digits
  57         }
  58         // Quotienten bilden:
  59         var DS q;
  60         var DS r;
  61         UDS_divide(z_MSDptr,z_len,z_LSDptr,
  62                    y_MSDptr,y_len,y_LSDptr,
  63                    &q, &r
  64                   );
  65         // q ist der Quotient,
  66         // q = floor(x/y*2^(intDsize*y_len)) >= 2^(intDsize*len),
  67         // q <= x/y*2^(intDsize*y_len) < 2^(1+intDsize+intDsize*len),
  68         // also mit len+1 oder len+2 Digits. r der Rest.
  69         var uintD* MSDptr = q.MSDptr;
  70         var uintL shiftcount;
  71         integerlengthD(mspref(MSDptr,0),shiftcount=);
  72         // Bei q.len = len+2 : shiftcount=1,
  73         // bei q.len = len+1 : 1<=shiftcount<=intDsize.
  74         var uintD carry_rechts;
  75         if (shiftcount==intDsize) {
  76                 carry_rechts = mspref(MSDptr,len);
  77         } else {
  78                 carry_rechts = shiftright_loop_msp(MSDptr,len+1,shiftcount%intDsize);
  79                 if (q.len > len+1) {
  80                         shiftcount += intDsize;
  81                         carry_rechts |= (mspref(MSDptr,len+1)==0 ? 0 : 1);
  82                 }
  83                 msshrink(MSDptr);
  84         }
  85         // Quotient MSDptr/len/.. ist nun normalisiert: höchstes Bit =1.
  86         // exponent := exponent(x) - intDsize*y_len + shiftcount
  87         var uintE uexp = TheLfloat(x)->expo;
  88         var uintE dexp = intDsize*y_len - shiftcount; // >= 0 !
  89         if ((uexp < dexp) || ((uexp = uexp - dexp) < LF_exp_low)) {
  90                 if (underflow_allowed())
  91                         { throw floating_point_underflow_exception(); }
  92                 else
  93                         { return encode_LF0(len); }
  94         }
  95         // Runden:
  96         if ( ((sintD)carry_rechts >= 0) // herausgeschobenes Bit =0 -> abrunden
  97              || ( (carry_rechts == (uintD)bit(intDsize-1)) // =1 und weitere Bits >0 oder Rest >0 -> aufrunden
  98                   && (r.len==0)
  99                   // round-to-even
 100                   && ((mspref(MSDptr,len-1) & bit(0)) ==0)
 101            )    )
 102           // abrunden
 103           {}
 104           else
 105           // aufrunden
 106           { if ( inc_loop_lsp(MSDptr mspop len,len) )
 107               // Übertrag durchs Aufrunden
 108               { mspref(MSDptr,0) = bit(intDsize-1); // Mantisse := 10...0
 109                 // Exponenten incrementieren:
 110                 if (++uexp ==  LF_exp_high+1) { throw floating_point_overflow_exception(); }
 111           }   }
 112         return encode_LFu(TheLfloat(x)->sign ^ sign, uexp, MSDptr, len);
 113 }
 114 // Bit complexity (N := max(length(x),length(y))): O(M(N)).
 115
 116 }  // namespace cln