// 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<m<900 (definitly negligible)
+// 800 never 850<m<1500
+// 1000 never 1030<m<2000
+// 1200 1400<m<1700 1230<m<2700
+// 1500 1590<m<2500 1530<m<5100, 6600<m (ridge negligible)
+// 2000 2060<m<3600 2030<m
+// 3000 3040<m 3030<m
+// 4000 4030<m 4030<m
+// 5000 5030<m 5030<m
+// 8000 8030<m 8030<m
+// Break-even-point, should be acceptable for both architectures.
+// When in doubt, prefer to choose the standard algorithm.
+#if CL_USE_GMP
+ static inline bool cl_recip_suitable (uintC m, uintC n) // m > 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
// 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);
// 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.
//
// 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<intDsize.
+// s:=intDsize-1-(Nummer des höchsten Bits in b[n-1]), 0<=s<intDsize.
// Schiebe a und b um s Bits nach links und kopiere sie dabei. r:=a.
// r=[r[m],...,r[0]], b=[b[n-1],...,b[0]] mit b[n-1]>=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:
// 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,
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;
}
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; }
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);
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);
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
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
#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
#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)
{ // 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...
}
// Bit complexity (N = a_len): O(M(N)).
+} // namespace cln
git://www.ginac.de / cln.git / blobdiff