yacas 1.3.6-2 / usr / share / yacas / scripts / specfunc.rep / bernou.ys

/// Simple implementation of the recurrence relation: create an array of Bernoulli numbers
// special cases: n=0 or n=1
10 # Internal'BernoulliArray(n_IsInteger)_(n=0 Or n=1) <-- [
	Local(B);
	B:=Array'Create(n+1,0);
	B[1] := 1;
	If(n=1, B[2] := -1/2);
	B;
];
/// Assume n>=2
20 # Internal'BernoulliArray(n_IsInteger) <-- [
	Local(B, i, k, k2, bin);
	If (InVerboseMode(), Echo({"Internal'BernoulliArray: using direct recursion, n = ", n}));
	B:=Array'Create(n+1, 0);	// array of B[k], k=1,2,... where B[1] is the 0th Bernoulli number
	// it would be better not to store the odd elements but let's optimize this later
	// we could also maintain a global cache of Bernoulli numbers computed so far, but it won't really speed up things at large n
	// all odd elements after B[2] are zero
	B[1] := 1;
	B[2] := -1/2;
	B[3] := 1/6;
	For(i:=4, i<=n, i := i+2)	// compute and store B[i]
	[	// maintain binomial coefficient
		bin := 1;	// Bin(i+1,0)
		// do not sum over odd elements that are zero anyway - cuts time in half
		B[i+1] := 1/2-1/(i+1)*(1 + Sum(k, 1, i/2-1,
			[
				bin := bin * (i+3-2*k) * (i+2-2*k)/ (2*k-1) / (2*k);
				B[2*k+1]*bin;	// *Bin(i+1, 2*k)
			]
		) );
	];
	B;
];

/// Find the fractional part of Bernoulli number with even index >=2
/// return negative if the sign of the Bernoulli number is negative
BernoulliFracPart(n_IsEven)_(n>=2) <-- [
	Local(p, sum);
	// always 2 and 3
	sum := 1/2+1/3;
	// check whether n+1 and n/2+1 are prime
	If(IsPrime(n+1), sum := sum+1/(n+1));
	If(IsPrime(n/2+1), sum := sum+1/(n/2+1));
	// sum over all primes p such that n / p-1 is integer
	// enough to check up to n/3 now
	For(p:=5, p<=n/3+1, p:=NextPrime(p))
		If(Mod(n, p-1)=0, sum := sum + 1/p);
	// for negative Bernoulli numbers, let's change sign
	// Mod(n/2, 2) is 0 for negative Bernoulli numbers and 1 for positive ones
	Div(Numer(sum), Denom(sum)) - sum
		 + Mod(n/2,2);	// we'll return a negative number if the Bernoulli itself is negative -- slightly against our definitions in the manual
		//+ 1;	// this would be exactly like the manual says
];

/// Find one Bernoulli number for large index
/// compute Riemann's zeta function and combine with the fractional part
Bernoulli1(n_IsEven)_(n>=2) <-- [
	Local(B, prec);
	prec := Builtin'Precision'Get();
	// estimate the size of B[n] using Stirling formula
	// and compute Ln(B[n])/Ln(10) to find the number of digits
	Builtin'Precision'Set(10);
	Builtin'Precision'Set(
		Ceil(N((1/2*Ln(8*Pi*n)-n+n*Ln(n/2/Pi))/Ln(10)))+3	// 3 guard digits
	);
	If (InVerboseMode(), Echo({"Bernoulli: using zeta funcion, precision ", Builtin'Precision'Set(), ", n = ", n}));
	B := Floor(N(	// compute integer part of B
		If(	// use different methods to compute Zeta function
			n>250,	// threshold is roughly right for internal math
			Internal'ZetaNum2(n, n/17+1),	// with this method, a single Bernoulli number n is computed in O(n*M(P)) operations where P = O(n*Ln(n)) is the required precision
			// Brent's method requires n^2*P+n*M(P)
			// simple array method requires 
			Internal'ZetaNum1(n, n/17+1)	// this gives O(n*Ln(n)*M(P))
		)
		*N(2*n! /(2*Pi)^n)))
		// 2*Pi*e is approx. 17, add 1 to guard precision
		* (2*Mod(n/2,2)-1)	// sign of B
		+ BernoulliFracPart(n);	// this already has the right sign
	Builtin'Precision'Set(prec);	// restore old precision
	B;
];

/// Bernoulli numbers; algorithm from: R. P. Brent, "A FORTRAN multiple-precision arithmetic package", ACM TOMS vol. 4, no. 1, p. 57 (1978).
/// this may be good for floating-point (not exact) evaluation of B[n] at large n
/// but is not good at all for exact evaluation! (too slow)
/// Brent claims that the usual recurrence is numerically unstable
/// but we can't check this because Yacas internal math is fixed-point and Brent's algorithm needs real floating point (C[k] are very small and then multiplied by (2*k)! )
Internal'BernoulliArray1(n_IsEven) _ (n>=2) <--
[
	Local(C, f, k, j, denom, sum);
	C := Array'Create(n+1, 0);
	f := Array'Create(n/2, 0);
	C[1] := 1;
	C[2] := -1/2;
	C[3] := 1/12;	// C[2*k+1] = B[2*k]/(2*k)!
	f[1] := 2;	// f[k] = (2k)!
	For(k:=2, k<=n/2, k++)	// we could start with k=1 but it would be awkward to compute f[] recursively
	[
		// compute f[k]
		f[k] := f[k-1] * (2*k)*(2*k-1);
		// compute C[k]
		C[2*k+1] := 1/(1-4^(-k))/2*(
			[
				denom := 4;	// = 4^1
				sum := 0;
				For(j:=1, j<=k-1, j++)
				[
					sum := sum + C[2*(k-j)+1]/denom/f[j];	// + C[k-j]/(2*j)! /4^j
					denom := denom * 4;
				];
				(2*k-1)/denom/f[k] - sum;
			]
		);
//	Echo({n, k, denom, C[k]});
	];
	// multiply C's with factorials to get B's
	For(k:=1, k<=n/2, k++)
		C[2*k+1] := C[2*k+1] * f[k];
	// return array object
	C;
];