12 #include "cln/lfloat.h"
13 #include "cl_LF_impl.h"
14 #include "cln/integer.h"
22 const cl_LF cl_LF_I_div (const cl_LF& x, const cl_I& y)
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) {
30 throw division_by_0_exception();
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);
42 // x länger als y, direkt dividieren.
44 var const uintD* y_MSDptr;
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)
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
61 UDS_divide(z_MSDptr,z_len,z_LSDptr,
62 y_MSDptr,y_len,y_LSDptr,
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;
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);
78 carry_rechts = shiftright_loop_msp(MSDptr,len+1,shiftcount%intDsize);
80 shiftcount += intDsize;
81 carry_rechts |= (mspref(MSDptr,len+1)==0 ? 0 : 1);
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(); }
93 { return encode_LF0(len); }
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
100 && ((mspref(MSDptr,len-1) & bit(0)) ==0)
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(); }
112 return encode_LFu(TheLfloat(x)->sign ^ sign, uexp, MSDptr, len);
114 // Bit complexity (N := max(length(x),length(y))): O(M(N)).