X-Git-Url: https://ginac.de/CLN/cln.git//cln.git?a=blobdiff_plain;f=src%2Fbase%2Fdigitseq%2Fcl_DS_div.cc;h=7393dc85cc6a6a44f84bc6a4672c8c7f51f90d8a;hb=22549ef70fee95faab1e9f2adaf710ba9e0bdabf;hp=de88316f9c907ba91d933ea7bf8e0c77dc470c3f;hpb=dd9e0f894eec7e2a8cf85078330ddc0a6639090b;p=cln.git diff --git a/src/base/digitseq/cl_DS_div.cc b/src/base/digitseq/cl_DS_div.cc index de88316..7393dc8 100644 --- a/src/base/digitseq/cl_DS_div.cc +++ b/src/base/digitseq/cl_DS_div.cc @@ -9,36 +9,79 @@ // Implementation. -#include "cl_N.h" -#include "cl_abort.h" +#include "cln/exception.h" -// We observe the following timings: -// Time for dividing 2*N digits by N digits, on a i486 33 MHz running Linux: -// N standard Newton -// 10 0.0003 0.0012 -// 25 0.0013 0.0044 -// 50 0.0047 0.0125 -// 100 0.017 0.037 -// 250 0.108 0.146 -// 500 0.43 0.44 -// 1000 1.72 1.32 -// 2500 11.2 4.1 -// 5000 44.3 9.5 -// 10000 187 20.6 -// -----> Newton faster for N >= 550. -// Time for dividing 3*N digits by N digits, on a i486 33 MHz running Linux: -// N standard Newton -// 10 0.0006 0.0025 -// 25 0.0026 0.0103 -// 50 0.0092 0.030 -// 100 0.035 0.089 -// 250 0.215 0.362 -// 500 0.85 1.10 -// 1000 3.44 3.21 -// 2500 23.3 7.9 -// 5000 89.0 15.6 -// 10000 362 33.1 -// -----> Newton faster for N >= 740. +namespace cln { + +// We observe the following timings in seconds: +// Time for dividing a 2*n word number by a n word number: +// OS: Linux 2.2, intDsize==32, OS: TRU64/4.0, intDsize==64, +// Machine: P-III/450MHz Machine: EV5/300MHz: +// n standard Newton standard Newton +// 10 0.000010 0.000024 0.000036 0.000058 +// 30 0.000026 0.000080 0.00012 0.00027 +// 100 0.00018 0.00048 0.00084 0.0016 +// 300 0.0013 0.0028 0.0062 0.0090 +// 1000 0.014 0.019 0.064 0.066 <-(~2200) +// 2000 0.058 0.058 <-(~2000) 0.26 0.20 +// 3000 0.20 0.11 0.57 0.24 +// 10000 2.3 0.50 6.7 1.2 +// 30000 24.4 1.2 62.0 2.8 +// Time for dividing a 3*n word number by a n word number: +// OS: Linux 2.2, intDsize==32, OS: TRU64/4.0, intDsize==64, +// Machine: P-III/450MHz Machine: EV5/300MHz: +// n standard Newton standard Newton +// 10 0.000013 0.000040 0.000063 0.00011 +// 30 0.000046 0.00018 0.00024 0.00062 +// 100 0.00035 0.0012 0.0016 0.0040 +// 300 0.0027 0.0071 0.012 0.021 +// 1000 0.029 0.047 0.13 0.16 +// 2000 0.12 0.14 <-(~2200) 0.51 0.45 <-(~1600) +// 3000 0.40 0.22 1.1 0.52 +// 10000 4.5 0.76 13.2 2.0 +// 30000 42.0 2.8 123.0 6.0 +// Time for dividing m digits by n digits: +// OS: Linux 2.2, intDsize==32, OS: TRU64/4.0, intDsize==64, +// Machine: P-III/450MHz Machine: EV5/300MHz: +// n Newton faster for: Newton faster for: +// 2-400 never never +// 600 never 670 n + { if (n < 900) + return false; + else + if (n < 2200) + return (m >= n+50) && (m < 2*n-600); + else + return m >= n+30; + } +#else +// Use the old default values from CLN version <= 1.0.3 as a crude estimate. +// They came from the timings for dividing m digits by n digits on an i486/33: +// Dividing 2*N digits by N digits: Dividing 3*N digits by N digits: +// N standard Newton standard Newton +// 10 0.0003 0.0012 0.0006 0.0025 +// 25 0.0013 0.0044 0.0026 0.0103 +// 50 0.0047 0.0125 0.0092 0.030 +// 100 0.017 0.037 0.035 0.089 +// 250 0.108 0.146 0.215 0.362 +// 500 0.43 0.44 <-(~550) 0.85 1.10 +// 1000 1.72 1.32 3.44 3.21 <-(~740) +// 2500 11.2 4.1 23.3 7.9 +// 5000 44.3 9.5 89.0 15.6 +// 10000 187 20.6 362 33.1 // Time for dividing m digits by n digits: // n = 2,3,5,10,25,50,100,250: Newton never faster. // n = 400: Newton faster for m >= 440, m < 600 @@ -52,16 +95,16 @@ // n = 2000: Newton faster for m >= 2020 // n = 2500: Newton faster for m >= 2520 // n = 5000: Newton faster for m >= 5020 -// Break-even-point. When in doubt, prefer to choose the standard algorithm. - static inline cl_boolean cl_recip_suitable (uintL m, uintL n) // m > n + static inline bool cl_recip_suitable (uintC m, uintC n) // m > n { if (n < 500) - return cl_false; + return false; else if (n < 1000) - return (cl_boolean)((m >= n+30) && (m < 3*n-600)); + return (m >= n+30) && (m < 3*n-600); else - return (cl_boolean)(m >= n+20); + return m >= n+20; } +#endif // Dividiert zwei Unsigned Digit sequences durcheinander. // UDS_divide(a_MSDptr,a_len,a_LSDptr, b_MSDptr,b_len,b_LSDptr, &q,&r); @@ -72,7 +115,7 @@ // q = q_MSDptr/q_len/q_LSDptr, r = r_MSDptr/r_len/r_LSDptr beides // Normalized Unsigned Digit sequences. // Vorsicht: q_LSDptr <= r_MSDptr, -// Vorzeichenerweiterung von r kann q zerstören! +// Vorzeichenerweiterung von r kann q zerstören! // Vorzeichenerweiterung von q ist erlaubt. // a und b werden nicht modifiziert. // @@ -89,19 +132,19 @@ // Normalisiere [q[m-1],...,q[0]], liefert q. // Falls m>=n>1, Multiple-Precision-Division: // Es gilt a/b < beta^(m-n+1). -// s:=intDsize-1-(Nummer des höchsten Bits in b[n-1]), 0<=s=beta/2. -// Für j=m-n,...,0: {Hier 0 <= r < b*beta^(j+1).} +// Für j=m-n,...,0: {Hier 0 <= r < b*beta^(j+1).} // Berechne q* : // q* := floor((r[j+n]*beta+r[j+n-1])/b[n-1]). -// Bei Überlauf (q* >= beta) setze q* := beta-1. +// Bei Überlauf (q* >= beta) setze q* := beta-1. // Berechne c2 := ((r[j+n]*beta+r[j+n-1]) - q* * b[n-1])*beta + r[j+n-2] // und c3 := b[n-2] * q*. // {Es ist 0 <= c2 < 2*beta^2, sogar 0 <= c2 < beta^2 falls kein -// Überlauf aufgetreten war. Ferner 0 <= c3 < beta^2. -// Bei Überlauf und r[j+n]*beta+r[j+n-1] - q* * b[n-1] >= beta, -// das heißt c2 >= beta^2, kann man die nächste Abfrage überspringen.} +// Überlauf aufgetreten war. Ferner 0 <= c3 < beta^2. +// Bei Überlauf und r[j+n]*beta+r[j+n-1] - q* * b[n-1] >= beta, +// das heißt c2 >= beta^2, kann man die nächste Abfrage überspringen.} // Solange c3 > c2, {hier 0 <= c2 < c3 < beta^2} setze // q* := q* - 1, c2 := c2 + b[n-1]*beta, c3 := c3 - b[n-2]. // Falls q* > 0: @@ -113,8 +156,8 @@ // u:=u div beta (+ 1, falls bei der Subtraktion Carry) // r[n+j]:=r[n+j]-u. // {Da stets u = (q* * [b[i-1],...,b[0]] div beta^i) + 1 -// < q* + 1 <= beta, läuft der Übertrag u nicht über.} -// Tritt dabei ein negativer Übertrag auf, so setze q* := q* - 1 +// < q* + 1 <= beta, läuft der Übertrag u nicht über.} +// Tritt dabei ein negativer Übertrag auf, so setze q* := q* - 1 // und [r[n+j],...,r[j]] := [r[n+j],...,r[j]] + [0,b[n-1],...,b[0]]. // Setze q[j] := q*. // Normalisiere [q[m-n],..,q[0]] und erhalte den Quotienten q, @@ -125,11 +168,11 @@ const uintD* b_MSDptr, uintC b_len, const uintD* b_LSDptr, uintD* roomptr, // ab roomptr kommen a_len+1 freie Digits DS* q_, DS* r_) - { // a normalisieren (a_MSDptr erhöhen, a_len erniedrigen): + { // a normalisieren (a_MSDptr erhöhen, a_len erniedrigen): while ((a_len>0) && (mspref(a_MSDptr,0)==0)) { msshrink(a_MSDptr); a_len--; } - // b normalisieren (b_MSDptr erhöhen, b_len erniedrigen): + // b normalisieren (b_MSDptr erhöhen, b_len erniedrigen): loop - { if (b_len==0) { cl_error_division_by_0(); } + { if (b_len==0) { throw division_by_0_exception(); } if (mspref(b_MSDptr,0)==0) { msshrink(b_MSDptr); b_len--; } else break; } @@ -201,7 +244,7 @@ var uintD* r_LSDptr = roomptr mspop (a_len+1); // Speicheraufbau: r_MSDptr/ a_len+1 /r_LSDptr // | r | - // später: q_MSDptr/a_len-b_len+1/r_MSDptr/b_len/r_LSDptr + // später: q_MSDptr/a_len-b_len+1/r_MSDptr/b_len/r_LSDptr // | q | r | if (s==0) { copy_loop_lsp(a_LSDptr,r_LSDptr,a_len); mspref(r_MSDptr,0) = 0; } @@ -238,7 +281,7 @@ var uintD* p_MSDptr; var uintD* p_LSDptr; UDS_UDS_mul_UDS(j+2,d_LSDptr, b_len,b_LSDptr, p_MSDptr=,,p_LSDptr=); - // d ist um höchstens 2 zu groß, muß also evtl. zweimal um 1 + // d ist um höchstens 2 zu groß, muß also evtl. zweimal um 1 // decrementieren, bis das Produkt <= a wird. if ((mspref(p_MSDptr,0) > 0) || (compare_loop_msp(p_MSDptr mspop 1,r_MSDptr,a_len+1) > 0)) { dec_loop_lsp(d_LSDptr,j+2); @@ -249,13 +292,13 @@ if (subfrom_loop_lsp(b_LSDptr,p_LSDptr,b_len)) dec_loop_lsp(p_LSDptr lspop b_len,j+2); if ((mspref(p_MSDptr,0) > 0) || (compare_loop_msp(p_MSDptr mspop 1,r_MSDptr,a_len+1) > 0)) - cl_abort(); + throw runtime_exception(); } } // Rest bestimmen: subfrom_loop_lsp(p_LSDptr,r_LSDptr,a_len+1); - if (test_loop_msp(r_MSDptr,j)) cl_abort(); + if (test_loop_msp(r_MSDptr,j)) throw runtime_exception(); r_MSDptr = r_LSDptr lspop b_len; // = r_MSDptr mspop (j+1); - // d ist um höchstens 2 zu klein, muß also evtl. zweimal um 1 + // d ist um höchstens 2 zu klein, muß also evtl. zweimal um 1 // incrementieren, bis der Rest < b wird. if ((lspref(r_MSDptr,0) > 0) || (compare_loop_msp(r_MSDptr,b_MSDptr,b_len) >= 0)) { inc_loop_lsp(d_LSDptr,j+2); @@ -266,10 +309,10 @@ if (subfrom_loop_lsp(b_LSDptr,r_LSDptr,b_len)) lspref(r_LSDptr,b_len) -= 1; if ((lspref(r_MSDptr,0) > 0) || (compare_loop_msp(r_MSDptr,b_MSDptr,b_len) >= 0)) - cl_abort(); + throw runtime_exception(); } } // r ist fertig, q := d. - if (mspref(d_MSDptr,0) > 0) cl_abort(); + if (mspref(d_MSDptr,0) > 0) throw runtime_exception(); q_len = j+1; copy_loop_msp(d_MSDptr mspop 1,q_MSDptr,q_len); } else @@ -281,13 +324,13 @@ var uintDD b_msdd = highlowDD(b_msd,b_2msd); // b[n-1]*beta+b[n-2] #endif // Divisions-Schleife: (wird m-n+1 mal durchlaufen) - // j = Herabzähler, b_MSDptr/b_len/b_LSDptr = [b[n-1],...,b[0]], b_len=n, + // j = Herabzähler, b_MSDptr/b_len/b_LSDptr = [b[n-1],...,b[0]], b_len=n, // r_MSDptr = Pointer auf r[n+j] = Pointer auf q[j], // r_ptr = Pointer oberhalb von r[j]. do { var uintD q_stern; var uintD c1; if (mspref(r_MSDptr,0) < b_msd) // r[j+n] < b[n-1] ? - { // Dividiere r[j+n]*beta+r[j+n-1] durch b[n-1], ohne Überlauf: + { // Dividiere r[j+n]*beta+r[j+n-1] durch b[n-1], ohne Überlauf: #if HAVE_DD divuD(highlowDD(mspref(r_MSDptr,0),mspref(r_MSDptr,1)),b_msd, q_stern=,c1=); #else @@ -295,7 +338,7 @@ #endif } else - { // Überlauf, also r[j+n]*beta+r[j+n-1] >= beta*b[n-1] + { // Überlauf, also r[j+n]*beta+r[j+n-1] >= beta*b[n-1] q_stern = bitm(intDsize)-1; // q* = beta-1 // Teste ob r[j+n]*beta+r[j+n-1] - (beta-1)*b[n-1] >= beta // <==> r[j+n]*beta+r[j+n-1] + b[n-1] >= beta*b[n-1]+beta @@ -313,11 +356,11 @@ #if HAVE_DD { var uintDD c2 = highlowDD(c1,mspref(r_MSDptr,2)); // c1*beta+r[j+n-2] var uintDD c3 = muluD(b_2msd,q_stern); // b[n-2] * q* - // Solange c2 < c3, c2 erhöhen, c3 erniedrigen: + // Solange c2 < c3, c2 erhöhen, c3 erniedrigen: // Rechne dabei mit c3-c2: // solange >0, um b[n-1]*beta+b[n-2] erniedrigen. // Dies kann wegen b[n-1]*beta+b[n-2] >= beta^2/2 - // höchstens zwei mal auftreten. + // höchstens zwei mal auftreten. if (c3 > c2) { q_stern = q_stern-1; // q* := q* - 1 if (c3-c2 > b_msdd) @@ -343,9 +386,9 @@ { // Subtraktionsschleife: r := r - b * q* * beta^j var uintD carry = mulusub_loop_lsp(q_stern,b_LSDptr,r_ptr,b_len); // Noch r_ptr[-b_len-1] -= carry, d.h. r_MSDptr[0] -= carry - // durchführen und danach r_MSDptr[0] vergessen: + // durchführen und danach r_MSDptr[0] vergessen: if (carry > mspref(r_MSDptr,0)) - // Subtraktion ergab Übertrag + // Subtraktion ergab Übertrag { q_stern = q_stern-1; // q* := q* - 1 addto_loop_lsp(b_LSDptr,r_ptr,b_len); // Additionsschleife // r[n+j] samt Carry kann vergessen werden... @@ -374,3 +417,4 @@ } // Bit complexity (N = a_len): O(M(N)). +} // namespace cln