// pi().
// General includes.
-#include "cl_sysdep.h"
+#include "base/cl_sysdep.h"
// Specification.
-#include "cl_F_tran.h"
+#include "float/transcendental/cl_F_tran.h"
// Implementation.
#include "cln/lfloat.h"
-#include "cl_LF_tran.h"
-#include "cl_LF.h"
+#include "float/transcendental/cl_LF_tran.h"
+#include "float/lfloat/cl_LF.h"
#include "cln/integer.h"
-#include "cl_alloca.h"
+#include "base/cl_alloca.h"
namespace cln {
ALL_cl_LF_OPERATIONS_SAME_PRECISION()
// For the next algorithms, I warmly recommend
-// [Jörg Arndt: hfloat documentation, august 1996,
+// [Jörg Arndt: hfloat documentation, august 1996,
// http://www.tu-chemnitz.de/~arndt/ hfdoc.dvi
// But beware of the typos in his formulas! ]
// functions. J. ACM 23(1976), 242-251.]
// [Jonathan M. Borwein, Peter B. Borwein: Pi and the AGM.
// Wiley 1987. Algorithm 2.2, p. 48.]
- // [Jörg Arndt, formula (4.51)-(4.52).]
+ // [Jörg Arndt, formula (4.51)-(4.52).]
// pi = AGM(1,1/sqrt(2))^2 * 2/(1 - sum(k=0..infty, 2^k c_k^2)).
// where the AGM iteration reads
// a_0 := 1, b_0 := 1/sqrt(2).
// )
// (/ (expt a 2) t)
// )
- var uintC actuallen = len + 1; // 1 Schutz-Digit
- var uintC uexp_limit = LF_exp_mid - intDsize*len;
- // Ein Long-Float ist genau dann betragsmäßig <2^-n, wenn
+ var uintC actuallen = len + 1; // 1 guard digit
+ var uintE uexp_limit = LF_exp_mid - intDsize*len;
+ // Ein Long-Float ist genau dann betragsmäßig <2^-n, wenn
// sein Exponent < LF_exp_mid-n = uexp_limit ist.
var cl_LF a = cl_I_to_LF(1,actuallen);
var cl_LF b = sqrt(scale_float(a,-1));
k++;
}
var cl_LF pires = square(a)/t; // a^2/t
- return shorten(pires,len); // verkürzen und fertig
+ return shorten(pires,len); // verkürzen und fertig
}
// Bit complexity (N := len): O(log(N)*M(N)).
const cl_LF compute_pi_brent_salamin_quartic (uintC len)
{
// See [Borwein, Borwein, section 1.4, exercise 3, p. 17].
- // See also [Jörg Arndt], formulas (4.52) and (3.30)[wrong!].
+ // See also [Jörg Arndt], formulas (4.52) and (3.30)[wrong!].
// As above, we are using the formula
// pi = AGM(1,1/sqrt(2))^2 * 2/(1 - sum(k=0..infty, 2^k c_k^2)).
// where the AGM iteration reads
// = 2^(k+1) * [wa_k^4 - ((wa_k^2+wb_k^2)/2)^2].
// Hence,
// pi = AGM(1,1/sqrt(2))^2 * 1/(1/2 - sum(k even, 2^k*[....])).
- var uintC actuallen = len + 1; // 1 Schutz-Digit
- var uintC uexp_limit = LF_exp_mid - intDsize*len;
+ var uintC actuallen = len + 1; // 1 guard digit
+ var uintE uexp_limit = LF_exp_mid - intDsize*len;
var cl_LF one = cl_I_to_LF(1,actuallen);
var cl_LF a = one;
var cl_LF wa = one;
k += 2;
}
var cl_LF pires = square(a)/t;
- return shorten(pires,len); // verkürzen und fertig
+ return shorten(pires,len); // verkürzen und fertig
}
// Bit complexity (N := len): O(log(N)*M(N)).
// mit J = -53360^3 = - (2^4 5 23 29)^3
// A = 163096908 = 2^2 3 13 1045493
// B = 6541681608 = 2^3 3^3 7 11 19 127 163
- // See [Jörg Arndt], formulas (4.27)-(4.30).
+ // See [Jörg Arndt], formulas (4.27)-(4.30).
// This is also the formula used in Pari.
// The absolute value of the n-th summand is approximately
// |J|^-n * n^(-1/2) * B/(2*pi^(3/2)),
// in precision.
// The sum is best evaluated using fixed-point arithmetic,
// so that the precision is reduced for the later summands.
- var uintC actuallen = len + 4; // 4 Schutz-Digits
+ var uintC actuallen = len + 4; // 4 guard digits
var sintC scale = intDsize*actuallen;
static const cl_I A = "163096908";
static const cl_I B = "6541681608";
var cl_LF fsum = scale_float(cl_I_to_LF(sum,actuallen),-scale);
static const cl_I J3 = "262537412640768000"; // -1728*J
var cl_LF pires = sqrt(cl_I_to_LF(J3,actuallen)) / fsum;
- return shorten(pires,len); // verkürzen und fertig
+ return shorten(pires,len); // verkürzen und fertig
}
// Bit complexity (N := len): O(N^2).
-#if defined(__mips__) && !defined(__GNUC__) // workaround SGI CC bug
-#define A A_fast
-#define B B_fast
-#define J3 J3_fast
-#endif
-
const cl_LF compute_pi_ramanujan_163_fast (uintC len)
{
// Same formula as above, using a binary splitting evaluation.
// See [Borwein, Borwein, section 10.2.3].
- var uintC actuallen = len + 4; // 4 Schutz-Digits
+ struct rational_series_stream : cl_pqa_series_stream {
+ uintC n;
+ static cl_pqa_series_term computenext (cl_pqa_series_stream& thisss)
+ {
+ static const cl_I A = "163096908";
+ static const cl_I B = "6541681608";
+ static const cl_I J1 = "10939058860032000"; // 72*abs(J)
+ var rational_series_stream& thiss = (rational_series_stream&)thisss;
+ var uintC n = thiss.n;
+ var cl_pqa_series_term result;
+ if (n==0) {
+ result.p = 1;
+ result.q = 1;
+ } else {
+ result.p = -((cl_I)(6*n-5)*(cl_I)(2*n-1)*(cl_I)(6*n-1));
+ result.q = (cl_I)n*(cl_I)n*(cl_I)n*J1;
+ }
+ result.a = A+n*B;
+ thiss.n = n+1;
+ return result;
+ }
+ rational_series_stream ()
+ : cl_pqa_series_stream (rational_series_stream::computenext),
+ n (0) {}
+ } series;
+ var uintC actuallen = len + 4; // 4 guard digits
static const cl_I A = "163096908";
static const cl_I B = "6541681608";
static const cl_I J1 = "10939058860032000"; // 72*abs(J)
var uintC N = (n_slope*actuallen)/32 + 1;
// N > intDsize*log(2)/log(|J|) * actuallen, hence
// |J|^-N < 2^(-intDsize*actuallen).
- CL_ALLOCA_STACK;
- var cl_I* av = (cl_I*) cl_alloca(N*sizeof(cl_I));
- var cl_I* pv = (cl_I*) cl_alloca(N*sizeof(cl_I));
- var cl_I* qv = (cl_I*) cl_alloca(N*sizeof(cl_I));
- var uintC* qsv = (uintC*) cl_alloca(N*sizeof(uintC));
- var uintC n;
- for (n = 0; n < N; n++) {
- init1(cl_I, av[n]) (A+n*B);
- if (n==0) {
- init1(cl_I, pv[n]) (1);
- init1(cl_I, qv[n]) (1);
- } else {
- init1(cl_I, pv[n]) (-((cl_I)(6*n-5)*(cl_I)(2*n-1)*(cl_I)(6*n-1)));
- init1(cl_I, qv[n]) ((cl_I)n*(cl_I)n*(cl_I)n*J1);
- }
- }
- var cl_pqa_series series;
- series.av = av;
- series.pv = pv; series.qv = qv;
- series.qsv = (len >= 35 ? qsv : 0); // 5% speedup for large len's
- var cl_LF fsum = eval_rational_series(N,series,actuallen);
- for (n = 0; n < N; n++) {
- av[n].~cl_I();
- pv[n].~cl_I();
- qv[n].~cl_I();
- }
+ var cl_LF fsum = eval_rational_series<false>(N,series,actuallen,actuallen);
static const cl_I J3 = "262537412640768000"; // -1728*J
var cl_LF pires = sqrt(cl_I_to_LF(J3,actuallen)) / fsum;
- return shorten(pires,len); // verkürzen und fertig
+ return shorten(pires,len); // verkürzen und fertig
}
// Bit complexity (N := len): O(log(N)^2*M(N)).
const cl_LF pi (uintC len)
{
- var uintC oldlen = TheLfloat(cl_LF_pi)->len; // vorhandene Länge
+ var uintC oldlen = TheLfloat(cl_LF_pi())->len; // vorhandene Länge
if (len < oldlen)
- return shorten(cl_LF_pi,len);
+ return shorten(cl_LF_pi(),len);
if (len == oldlen)
- return cl_LF_pi;
+ return cl_LF_pi();
- // TheLfloat(cl_LF_pi)->len um mindestens einen konstanten Faktor
- // > 1 wachsen lassen, damit es nicht zu häufig nachberechnet wird:
+ // TheLfloat(cl_LF_pi())->len um mindestens einen konstanten Faktor
+ // > 1 wachsen lassen, damit es nicht zu häufig nachberechnet wird:
var uintC newlen = len;
oldlen += floor(oldlen,2); // oldlen * 3/2
if (newlen < oldlen)
newlen = oldlen;
- // gewünschte > vorhandene Länge -> muß nachberechnen:
- cl_LF_pi = compute_pi_ramanujan_163_fast(newlen);
- return (len < newlen ? shorten(cl_LF_pi,len) : cl_LF_pi);
+ // gewünschte > vorhandene Länge -> muß nachberechnen:
+ cl_LF_pi() = compute_pi_ramanujan_163_fast(newlen);
+ return (len < newlen ? shorten(cl_LF_pi(),len) : cl_LF_pi());
}
} // namespace cln
git://www.ginac.de / cln.git / blobdiff