nd4j-common/src/main/java/org/nd4j/linalg/util/BigDecimalMath.java
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);


    }
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.

@param a  The vector of 8 integer arguments
@param mc Specification of the accuracy of the result
@return S_(n, p)(a)
@see \protect\vrule width0pt\protect\href{http://arxiv.org/abs/math/9803067}{arXiv:math/9803067}