* Optimizations needed... too long!!!... (memoization???)
*
*/
static public BigDecimal zeta(final int n, final MathContext mc, Bernoulli bern_cache, Factorial fact_cache) {
MathContext nmc=new MathContext(3+mc.getPrecision(),mc.getRoundingMode());
if (n <= 0)
throw new ProviderException("Not implemented: zeta at negative argument " + n);
if (n == 1)
throw new ArithmeticException("Pole at zeta(1) ");
if (n % 2 == 0) {
/*
* Even indices. Abramowitz-Stegun 23.2.16. Start with 2^(n-1)*B(n)/n!
*/
//System.out.print("'");
Rational b = (bern_cache!=null ? bern_cache : new Bernoulli()).at(n).abs();
b = b.divide((fact_cache!=null ? fact_cache : new Factorial()).at(n));
b = b.multiply(BigInteger.ONE.shiftLeft(n - 1));
//System.out.print("'");
/* to be multiplied by pi^n. Absolute error in the result of pi^n is n times
* error in pi times pi^(n-1). Relative error is n*error(pi)/pi, requested by mc.
* Need one more digit in pi if n=10, two digits if n=100 etc, and add one extra digit.
*/
MathContext mcpi = new MathContext(mc.getPrecision() + (int) (Math.log10(10.0 * n)));
BigDecimal piton = pi(mcpi).pow(n, mc);
return piton.multiply(b.BigDecimalValue(mc),mc);
}
else if (n == 3) {
/* Broadhurst BBP <a href="http://arxiv.org/abs/math/9803067">arXiv:math/9803067</a>
* Error propagation: S31 is roughly 0.087, S33 roughly 0.131
*/
int[] a31 = { 1,
-7,
-1,
10,
-1,
-7,
1,
0 };
int[] a33 = { 1,
1,
-1,
-2,
-1,
1,
1,
0 };
BigDecimal S31 = broadhurstBBP(3, 1, a31, mc);
BigDecimal S33 = broadhurstBBP(3, 3, a33, mc);
S31 = S31.multiply(new BigDecimal(48),nmc);
S33 = S33.multiply(new BigDecimal(32),nmc);
return S31.add(S33).divide(new BigDecimal(7), mc);
}
else if (n == 5) {
/* Broadhurst BBP <a href=http://arxiv.org/abs/math/9803067">arXiv:math/9803067</a>
* Error propagation: S51 is roughly -11.15, S53 roughly 22.165, S55 is roughly 0.031
* 9*2048*S51/6265 = -3.28. 7*2038*S53/61651= 5.07. 738*2048*S55/61651= 0.747.
* The result is of the order 1.03, so we add 2 digits to S51 and S52 and one digit to S55.
*/
int[] a51 = { 31,
-1614,
-31,
-6212,
-31,
-1614,
31,
74552 };
int[] a53 = { 173,
284,
-173,
-457,
-173,
284,
173,
-111 };
int[] a55 = { 1,
0,
-1,
-1,
-1,
0,
1,
1 };
BigDecimal S51 = broadhurstBBP(5, 1, a51, new MathContext(2 + mc.getPrecision()));
BigDecimal S53 = broadhurstBBP(5, 3, a53, new MathContext(2 + mc.getPrecision()));
BigDecimal S55 = broadhurstBBP(5, 5, a55, new MathContext(1 + mc.getPrecision()));
S51 = S51.multiply(new BigDecimal(18432),nmc);
S53 = S53.multiply(new BigDecimal(14336),nmc);
S55 = S55.multiply(new BigDecimal(1511424),nmc);
return S51.add(S53).subtract(S55).divide(new BigDecimal(62651), mc);
}
else {
/* Cohen et al Exp Math 1 (1) (1992) 25
*/
//System.out.print("'");
Rational betsum = new Rational();
Bernoulli bern = bern_cache!=null ? bern_cache : new Bernoulli();
Factorial fact = fact_cache!=null ? fact_cache : new Factorial();
for (int npr = 0; npr <= (n + 1) / 2; npr++) {
Rational b = bern.at(2 * npr).multiply(bern.at(n + 1 - 2 * npr));
b = b.divide(fact.at(2 * npr)).divide(fact.at(n + 1 - 2 * npr));
b = b.multiply(1 - 2 * npr);
if (npr % 2 == 0)
betsum = betsum.add(b);
else
betsum = betsum.subtract(b);
}
//System.out.print("'");
betsum = betsum.divide(n - 1);
/* The first term, including the facor (2pi)^n, is essentially most
* of the result, near one. The second term below is roughly in the range 0.003 to 0.009.
* So the precision here is matching the precisionn requested by mc, and the precision
* requested for 2*pi is in absolute terms adjusted.
*/
MathContext mcloc = new MathContext(2 + mc.getPrecision() + (int) (Math.log10((double) (n))));
BigDecimal ftrm = pi(mcloc).multiply(TWO,mcloc);
ftrm = ftrm.pow(n,mcloc);
ftrm = ftrm.multiply(betsum.BigDecimalValue(mcloc),mcloc);
BigDecimal exps = new BigDecimal(0);