static protected BigDecimal broadhurstBBP ( final int n , final int p , final int a [ ] , MathContext mc ) { /* Explore the actual magnitude of the result first with a quick estimate. */ double x = 0.0 ; for ( int k = 1 ; k < 10 ; k ++ ) { x += a [ ( k - 1 ) % 8 ] / Math . pow ( 2. , p * ( k + 1 ) / 2 ) / Math . pow ( ( double ) k , n ) ; } /* Convert the relative precision and estimate of the result into an absolute precision. */ double eps = prec2err ( x , mc . getPrecision ( ) ) ; /* Divide this through the number of terms in the sum to account for error accumulation * The divisor 2^(p(k+1)/2) means that on the average each 8th term in k has shrunk by * relative to the 8th predecessor by 1/2^(4p). 1/2^(4pc) = 10^(-precision) with c the 8term * cycles yields c=log_2( 10^precision)/4p = 3.3*precision/4p with k=8c */ int kmax = ( int ) ( 6.6 * mc . getPrecision ( ) / p ) ; /* Now eps is the absolute error in each term */ eps /= kmax ; BigDecimal res = BigDecimal . ZERO ; for ( int c = 0 ; ; c ++ ) { Rational r = new Rational ( ) ; for ( int k = 0 ; k < 8 ; k ++ ) { Rational tmp = new Rational ( BigInteger . valueOf ( a [ k ] ) , BigInteger . valueOf ( ( 1 + 8 * c + k ) ) . pow ( n ) ) ; /* floor( (pk+p)/2) */ int pk1h = p * ( 2 + 8 * c + k ) / 2 ; tmp = tmp . divide ( BigInteger . ONE . shiftLeft ( pk1h ) ) ; r = r . add ( tmp ) ; } if ( Math . abs ( r . doubleValue ( ) ) < eps ) { break ; } MathContext mcloc = new MathContext ( 1 + err2prec ( r . doubleValue ( ) , eps ) ) ; res = res . add ( r . BigDecimalValue ( mcloc ) ) ; } return res . round ( mc ) ; }
Broadhurst ladder sequence .