nd4j-common/src/main/java/org/nd4j/linalg/util/BigDecimalMath.java
static public BigDecimal zeta(final int n, final MathContext mc) {
        if (n <= 0) {
            throw new ProviderException("Unimplemented 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!
             */
            Rational b = (new Bernoulli()).at(n).abs();
            b = b.divide((new Factorial()).at(n));
            b = b.multiply(BigInteger.ONE.shiftLeft(n - 1));
            /* 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)));
            final BigDecimal piton = pi(mcpi).pow(n, mc);


            return multiplyRound(piton, b);


        } else if (n == 3) {
            /* Broadhurst BBP \protect\vrule width0pt\protect\href{http://arxiv.org/abs/math/9803067}{arXiv:math/9803067}
             * 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));
            S33 = S33.multiply(new BigDecimal(32));


            return S31.add(S33).divide(new BigDecimal(7), mc);


        } else if (n == 5) {
            /* Broadhurst BBP \protect\vrule width0pt\protect\href{http://arxiv.org/abs/math/9803067}{arXiv:math/9803067}
             * 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));
            S53 = S53.multiply(new BigDecimal(14336));
            S55 = S55.multiply(new BigDecimal(1511424));


            return S51.add(S53).subtract(S55).divide(new BigDecimal(62651), mc);


        } else {
            /* Cohen et al Exp Math 1 (1) (1992) 25
             */
            Rational betsum = new Rational();
            Bernoulli bern = new Bernoulli();
            Factorial fact = 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);
                }


            }
            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(new BigDecimal(2));
            ftrm = ftrm.pow(n);


            ftrm = multiplyRound(ftrm, betsum.BigDecimalValue(mcloc));
            BigDecimal exps = new BigDecimal(0);
            /* the basic accuracy of the accumulated terms before multiplication with 2
             */


            double eps = Math.pow(10., -mc.getPrecision());


            if (n % 4 == 3) {
                /* since the argument n is at least 7 here, the drop
                 * of the terms is at rather constant pace at least 10^-3, for example
                 * 0.0018, 0.2e-7, 0.29e-11, 0.74e-15 etc for npr=1,2,3.... We want 2 times these terms
                 * fall below eps/10.
                 */
                int kmax = mc.getPrecision() / 3;
                eps /= kmax;
                /* need an error of eps for 2/(exp(2pi)-1) = 0.0037
                 * The absolute error is 4*exp(2pi)*err(pi)/(exp(2pi)-1)^2=0.0075*err(pi)
                 */
                BigDecimal exp2p = pi(new MathContext(3 + err2prec(3.14, eps / 0.0075)));
                exp2p = exp(exp2p.multiply(new BigDecimal(2)));
                BigDecimal c = exp2p.subtract(BigDecimal.ONE);
                exps = divideRound(1, c);


                for (int npr = 2; npr <= kmax; npr++) {
                    /* the error estimate above for npr=1 is the worst case of
                     * the absolute error created by an error in 2pi. So we can
                     * safely re-use the exp2p value computed above without
                     * reassessment of its error.
                     */
                    c = powRound(exp2p, npr).subtract(BigDecimal.ONE);
                    c = multiplyRound(c, (BigInteger.valueOf(npr)).pow(n));
                    c = divideRound(1, c);
                    exps = exps.add(c);


                }
            } else {
                /* since the argument n is at least 9 here, the drop
                 * of the terms is at rather constant pace at least 10^-3, for example
                 * 0.0096, 0.5e-7, 0.3e-11, 0.6e-15 etc. We want these terms
                 * fall below eps/10.
                 */
                int kmax = (1 + mc.getPrecision()) / 3;
                eps /= kmax;
                /* need an error of eps for 2/(exp(2pi)-1)*(1+4*Pi/8/(1-exp(-2pi)) = 0.0096
                 * at k=7 or = 0.00766 at k=13 for example.
                 * The absolute error is 0.017*err(pi) at k=9, 0.013*err(pi) at k=13, 0.012 at k=17
                 */
                BigDecimal twop = pi(new MathContext(3 + err2prec(3.14, eps / 0.017)));
                twop = twop.multiply(new BigDecimal(2));
                BigDecimal exp2p = exp(twop);
                BigDecimal c = exp2p.subtract(BigDecimal.ONE);
                exps = divideRound(1, c);
                c = BigDecimal.ONE.subtract(divideRound(1, exp2p));
                c = divideRound(twop, c).multiply(new BigDecimal(2));
                c = divideRound(c, n - 1).add(BigDecimal.ONE);
                exps = multiplyRound(exps, c);


                for (int npr = 2; npr <= kmax; npr++) {
                    c = powRound(exp2p, npr).subtract(BigDecimal.ONE);
                    c = multiplyRound(c, (BigInteger.valueOf(npr)).pow(n));
                    BigDecimal d = divideRound(1, exp2p.pow(npr));
                    d = BigDecimal.ONE.subtract(d);
                    d = divideRound(twop, d).multiply(new BigDecimal(2 * npr));
                    d = divideRound(d, n - 1).add(BigDecimal.ONE);
                    d = divideRound(d, c);
                    exps = exps.add(d);


                }
            }
            exps = exps.multiply(new BigDecimal(2));


            return ftrm.subtract(exps, mc);


        }
    }
static public BigDecimal zeta(final int n, final MathContext mc) {
        if (n <= 0) {
            throw new ProviderException("Unimplemented 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!
             */
            Rational b = (new Bernoulli()).at(n).abs();
            b = b.divide((new Factorial()).at(n));
            b = b.multiply(BigInteger.ONE.shiftLeft(n - 1));
            /* 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)));
            final BigDecimal piton = pi(mcpi).pow(n, mc);


            return multiplyRound(piton, b);


        } else if (n == 3) {
            /* Broadhurst BBP \protect\vrule width0pt\protect\href{http://arxiv.org/abs/math/9803067}{arXiv:math/9803067}
             * 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));
            S33 = S33.multiply(new BigDecimal(32));


            return S31.add(S33).divide(new BigDecimal(7), mc);


        } else if (n == 5) {
            /* Broadhurst BBP \protect\vrule width0pt\protect\href{http://arxiv.org/abs/math/9803067}{arXiv:math/9803067}
             * 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));
            S53 = S53.multiply(new BigDecimal(14336));
            S55 = S55.multiply(new BigDecimal(1511424));


            return S51.add(S53).subtract(S55).divide(new BigDecimal(62651), mc);


        } else {
            /* Cohen et al Exp Math 1 (1) (1992) 25
             */
            Rational betsum = new Rational();
            Bernoulli bern = new Bernoulli();
            Factorial fact = 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);
                }


            }
            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(new BigDecimal(2));
            ftrm = ftrm.pow(n);


            ftrm = multiplyRound(ftrm, betsum.BigDecimalValue(mcloc));
            BigDecimal exps = new BigDecimal(0);
            /* the basic accuracy of the accumulated terms before multiplication with 2
             */


            double eps = Math.pow(10., -mc.getPrecision());


            if (n % 4 == 3) {
                /* since the argument n is at least 7 here, the drop
                 * of the terms is at rather constant pace at least 10^-3, for example
                 * 0.0018, 0.2e-7, 0.29e-11, 0.74e-15 etc for npr=1,2,3.... We want 2 times these terms
                 * fall below eps/10.
                 */
                int kmax = mc.getPrecision() / 3;
                eps /= kmax;
                /* need an error of eps for 2/(exp(2pi)-1) = 0.0037
                 * The absolute error is 4*exp(2pi)*err(pi)/(exp(2pi)-1)^2=0.0075*err(pi)
                 */
                BigDecimal exp2p = pi(new MathContext(3 + err2prec(3.14, eps / 0.0075)));
                exp2p = exp(exp2p.multiply(new BigDecimal(2)));
                BigDecimal c = exp2p.subtract(BigDecimal.ONE);
                exps = divideRound(1, c);


                for (int npr = 2; npr <= kmax; npr++) {
                    /* the error estimate above for npr=1 is the worst case of
                     * the absolute error created by an error in 2pi. So we can
                     * safely re-use the exp2p value computed above without
                     * reassessment of its error.
                     */
                    c = powRound(exp2p, npr).subtract(BigDecimal.ONE);
                    c = multiplyRound(c, (BigInteger.valueOf(npr)).pow(n));
                    c = divideRound(1, c);
                    exps = exps.add(c);


                }
            } else {
                /* since the argument n is at least 9 here, the drop
                 * of the terms is at rather constant pace at least 10^-3, for example
                 * 0.0096, 0.5e-7, 0.3e-11, 0.6e-15 etc. We want these terms
                 * fall below eps/10.
                 */
                int kmax = (1 + mc.getPrecision()) / 3;
                eps /= kmax;
                /* need an error of eps for 2/(exp(2pi)-1)*(1+4*Pi/8/(1-exp(-2pi)) = 0.0096
                 * at k=7 or = 0.00766 at k=13 for example.
                 * The absolute error is 0.017*err(pi) at k=9, 0.013*err(pi) at k=13, 0.012 at k=17
                 */
                BigDecimal twop = pi(new MathContext(3 + err2prec(3.14, eps / 0.017)));
                twop = twop.multiply(new BigDecimal(2));
                BigDecimal exp2p = exp(twop);
                BigDecimal c = exp2p.subtract(BigDecimal.ONE);
                exps = divideRound(1, c);
                c = BigDecimal.ONE.subtract(divideRound(1, exp2p));
                c = divideRound(twop, c).multiply(new BigDecimal(2));
                c = divideRound(c, n - 1).add(BigDecimal.ONE);
                exps = multiplyRound(exps, c);


                for (int npr = 2; npr <= kmax; npr++) {
                    c = powRound(exp2p, npr).subtract(BigDecimal.ONE);
                    c = multiplyRound(c, (BigInteger.valueOf(npr)).pow(n));
                    BigDecimal d = divideRound(1, exp2p.pow(npr));
                    d = BigDecimal.ONE.subtract(d);
                    d = divideRound(twop, d).multiply(new BigDecimal(2 * npr));
                    d = divideRound(d, n - 1).add(BigDecimal.ONE);
                    d = divideRound(d, c);
                    exps = exps.add(d);


                }
            }
            exps = exps.multiply(new BigDecimal(2));


            return ftrm.subtract(exps, mc);


        }
    }
static public BigDecimal zeta ( final int n , final MathContext mc ) { if ( n <= 0 ) { throw new ProviderException ( "Unimplemented 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!
             */ Rational b = ( new Bernoulli ( ) ) . at ( n ) . abs ( ) ; b = b . divide ( ( new Factorial ( ) ) . at ( n ) ) ; b = b . multiply ( BigInteger . ONE . shiftLeft ( n - 1 ) ) ; /* 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 ) ) ) ; final BigDecimal piton = pi ( mcpi ) . pow ( n , mc ) ; return multiplyRound ( piton , b ) ; } else if ( n == 3 ) { /* Broadhurst BBP \protect\vrule width0pt\protect\href{http://arxiv.org/abs/math/9803067}{arXiv:math/9803067}
             * 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 ) ) ; S33 = S33 . multiply ( new BigDecimal ( 32 ) ) ; return S31 . add ( S33 ) . divide ( new BigDecimal ( 7 ) , mc ) ; } else if ( n == 5 ) { /* Broadhurst BBP \protect\vrule width0pt\protect\href{http://arxiv.org/abs/math/9803067}{arXiv:math/9803067}
             * 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 ) ) ; S53 = S53 . multiply ( new BigDecimal ( 14336 ) ) ; S55 = S55 . multiply ( new BigDecimal ( 1511424 ) ) ; return S51 . add ( S53 ) . subtract ( S55 ) . divide ( new BigDecimal ( 62651 ) , mc ) ; } else { /* Cohen et al Exp Math 1 (1) (1992) 25
             */ Rational betsum = new Rational ( ) ; Bernoulli bern = new Bernoulli ( ) ; Factorial fact = 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 ) ; } } 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 ( new BigDecimal ( 2 ) ) ; ftrm = ftrm . pow ( n ) ; ftrm = multiplyRound ( ftrm , betsum . BigDecimalValue ( mcloc ) ) ; BigDecimal exps = new BigDecimal ( 0 ) ; /* the basic accuracy of the accumulated terms before multiplication with 2
             */ double eps = Math . pow ( 10. , - mc . getPrecision ( ) ) ; if ( n % 4 == 3 ) { /* since the argument n is at least 7 here, the drop
                 * of the terms is at rather constant pace at least 10^-3, for example
                 * 0.0018, 0.2e-7, 0.29e-11, 0.74e-15 etc for npr=1,2,3.... We want 2 times these terms
                 * fall below eps/10.
                 */ int kmax = mc . getPrecision ( ) / 3 ; eps /= kmax ; /* need an error of eps for 2/(exp(2pi)-1) = 0.0037
                 * The absolute error is 4*exp(2pi)*err(pi)/(exp(2pi)-1)^2=0.0075*err(pi)
                 */ BigDecimal exp2p = pi ( new MathContext ( 3 + err2prec ( 3.14 , eps / 0.0075 ) ) ) ; exp2p = exp ( exp2p . multiply ( new BigDecimal ( 2 ) ) ) ; BigDecimal c = exp2p . subtract ( BigDecimal . ONE ) ; exps = divideRound ( 1 , c ) ; for ( int npr = 2 ; npr <= kmax ; npr ++ ) { /* the error estimate above for npr=1 is the worst case of
                     * the absolute error created by an error in 2pi. So we can
                     * safely re-use the exp2p value computed above without
                     * reassessment of its error.
                     */ c = powRound ( exp2p , npr ) . subtract ( BigDecimal . ONE ) ; c = multiplyRound ( c , ( BigInteger . valueOf ( npr ) ) . pow ( n ) ) ; c = divideRound ( 1 , c ) ; exps = exps . add ( c ) ; } } else { /* since the argument n is at least 9 here, the drop
                 * of the terms is at rather constant pace at least 10^-3, for example
                 * 0.0096, 0.5e-7, 0.3e-11, 0.6e-15 etc. We want these terms
                 * fall below eps/10.
                 */ int kmax = ( 1 + mc . getPrecision ( ) ) / 3 ; eps /= kmax ; /* need an error of eps for 2/(exp(2pi)-1)*(1+4*Pi/8/(1-exp(-2pi)) = 0.0096
                 * at k=7 or = 0.00766 at k=13 for example.
                 * The absolute error is 0.017*err(pi) at k=9, 0.013*err(pi) at k=13, 0.012 at k=17
                 */ BigDecimal twop = pi ( new MathContext ( 3 + err2prec ( 3.14 , eps / 0.017 ) ) ) ; twop = twop . multiply ( new BigDecimal ( 2 ) ) ; BigDecimal exp2p = exp ( twop ) ; BigDecimal c = exp2p . subtract ( BigDecimal . ONE ) ; exps = divideRound ( 1 , c ) ; c = BigDecimal . ONE . subtract ( divideRound ( 1 , exp2p ) ) ; c = divideRound ( twop , c ) . multiply ( new BigDecimal ( 2 ) ) ; c = divideRound ( c , n - 1 ) . add ( BigDecimal . ONE ) ; exps = multiplyRound ( exps , c ) ; for ( int npr = 2 ; npr <= kmax ; npr ++ ) { c = powRound ( exp2p , npr ) . subtract ( BigDecimal . ONE ) ; c = multiplyRound ( c , ( BigInteger . valueOf ( npr ) ) . pow ( n ) ) ; BigDecimal d = divideRound ( 1 , exp2p . pow ( npr ) ) ; d = BigDecimal . ONE . subtract ( d ) ; d = divideRound ( twop , d ) . multiply ( new BigDecimal ( 2 * npr ) ) ; d = divideRound ( d , n - 1 ) . add ( BigDecimal . ONE ) ; d = divideRound ( d , c ) ; exps = exps . add ( d ) ; } } exps = exps . multiply ( new BigDecimal ( 2 ) ) ; return ftrm . subtract ( exps , mc ) ; } }
Riemann zeta function.

@param n  The positive integer argument.
32
@param mc Specification of the accuracy of the result.
@return zeta(n).
Riemann zeta function .