Example usage for java.math BigDecimal ulp

List of usage examples for java.math BigDecimal ulp

Introduction

In this page you can find the example usage for java.math BigDecimal ulp.

Prototype

public BigDecimal ulp() 

Source Link

Document

Returns the size of an ulp, a unit in the last place, of this BigDecimal .

Usage

From source file:org.nd4j.linalg.util.BigDecimalMath.java

/**
 * The trigonometric tangent./*from  w  w w  . ja va 2  s  .c om*/
 *
 * @param x the argument in radians.
 * @return the tan(x)
 */
static public BigDecimal tan(final BigDecimal x) {
    if (x.compareTo(BigDecimal.ZERO) == 0) {
        return BigDecimal.ZERO;
    } else if (x.compareTo(BigDecimal.ZERO) < 0) {
        return tan(x.negate()).negate();
    } else {
        /* reduce modulo pi
         */
        BigDecimal res = modpi(x);
        /* absolute error in the result is err(x)/cos^2(x) to lowest order
         */
        final double xDbl = res.doubleValue();
        final double xUlpDbl = x.ulp().doubleValue() / 2.;
        final double eps = xUlpDbl / 2. / Math.pow(Math.cos(xDbl), 2.);
        if (xDbl > 0.8) {
            /* tan(x) = 1/cot(x) */
            BigDecimal co = cot(x);
            MathContext mc = new MathContext(err2prec(1. / co.doubleValue(), eps));
            return BigDecimal.ONE.divide(co, mc);
        } else {
            final BigDecimal xhighpr = scalePrec(res, 2);
            final BigDecimal xhighprSq = multiplyRound(xhighpr, xhighpr);
            BigDecimal resul = xhighpr.plus();
            /* x^(2i+1) */
            BigDecimal xpowi = xhighpr;
            Bernoulli b = new Bernoulli();
            /* 2^(2i) */
            BigInteger fourn = new BigInteger("4");
            /* (2i)! */
            BigInteger fac = new BigInteger("2");
            for (int i = 2;; i++) {
                Rational f = b.at(2 * i).abs();
                fourn = fourn.shiftLeft(2);
                fac = fac.multiply(new BigInteger("" + (2 * i))).multiply(new BigInteger("" + (2 * i - 1)));
                f = f.multiply(fourn).multiply(fourn.subtract(BigInteger.ONE)).divide(fac);
                xpowi = multiplyRound(xpowi, xhighprSq);
                BigDecimal c = multiplyRound(xpowi, f);
                resul = resul.add(c);
                if (Math.abs(c.doubleValue()) < 0.1 * eps) {
                    break;
                }
            }
            MathContext mc = new MathContext(err2prec(resul.doubleValue(), eps));
            return resul.round(mc);
        }
    }
}

From source file:org.nd4j.linalg.util.BigDecimalMath.java

/**
 * The trigonometric co-tangent./*from  w  w  w  . j a  v a2s  .  co  m*/
 *
 * @param x the argument in radians.
 * @return the cot(x)
 */
static public BigDecimal cot(final BigDecimal x) {
    if (x.compareTo(BigDecimal.ZERO) == 0) {
        throw new ArithmeticException("Cannot take cot of zero " + x.toString());
    } else if (x.compareTo(BigDecimal.ZERO) < 0) {
        return cot(x.negate()).negate();
    } else {
        /* reduce modulo pi
         */
        BigDecimal res = modpi(x);
        /* absolute error in the result is err(x)/sin^2(x) to lowest order
         */
        final double xDbl = res.doubleValue();
        final double xUlpDbl = x.ulp().doubleValue() / 2.;
        final double eps = xUlpDbl / 2. / Math.pow(Math.sin(xDbl), 2.);
        final BigDecimal xhighpr = scalePrec(res, 2);
        final BigDecimal xhighprSq = multiplyRound(xhighpr, xhighpr);
        MathContext mc = new MathContext(err2prec(xhighpr.doubleValue(), eps));
        BigDecimal resul = BigDecimal.ONE.divide(xhighpr, mc);
        /* x^(2i-1) */
        BigDecimal xpowi = xhighpr;
        Bernoulli b = new Bernoulli();
        /* 2^(2i) */
        BigInteger fourn = new BigInteger("4");
        /* (2i)! */
        BigInteger fac = BigInteger.ONE;
        for (int i = 1;; i++) {
            Rational f = b.at(2 * i);
            fac = fac.multiply(new BigInteger("" + (2 * i))).multiply(new BigInteger("" + (2 * i - 1)));
            f = f.multiply(fourn).divide(fac);
            BigDecimal c = multiplyRound(xpowi, f);
            if (i % 2 == 0) {
                resul = resul.add(c);
            } else {
                resul = resul.subtract(c);
            }
            if (Math.abs(c.doubleValue()) < 0.1 * eps) {
                break;
            }
            fourn = fourn.shiftLeft(2);
            xpowi = multiplyRound(xpowi, xhighprSq);
        }
        mc = new MathContext(err2prec(resul.doubleValue(), eps));
        return resul.round(mc);
    }
}

From source file:org.nd4j.linalg.util.BigDecimalMath.java

/**
 * The inverse trigonometric sine./*from   w  w  w .j a  v a2s  .  co m*/
 *
 * @param x the argument.
 * @return the arcsin(x) in radians.
 */
static public BigDecimal asin(final BigDecimal x) {
    if (x.compareTo(BigDecimal.ONE) > 0 || x.compareTo(BigDecimal.ONE.negate()) < 0) {
        throw new ArithmeticException("Out of range argument " + x.toString() + " of asin");

    } else if (x.compareTo(BigDecimal.ZERO) == 0) {
        return BigDecimal.ZERO;
    } else if (x.compareTo(BigDecimal.ONE) == 0) {
        /* arcsin(1) = pi/2
         */
        double errpi = Math.sqrt(x.ulp().doubleValue());
        MathContext mc = new MathContext(err2prec(3.14159, errpi));

        return pi(mc).divide(new BigDecimal(2));

    } else if (x.compareTo(BigDecimal.ZERO) < 0) {
        return asin(x.negate()).negate();

    } else if (x.doubleValue() > 0.7) {
        final BigDecimal xCompl = BigDecimal.ONE.subtract(x);
        final double xDbl = x.doubleValue();
        final double xUlpDbl = x.ulp().doubleValue() / 2.;
        final double eps = xUlpDbl / 2. / Math.sqrt(1. - Math.pow(xDbl, 2.));

        final BigDecimal xhighpr = scalePrec(xCompl, 3);
        final BigDecimal xhighprV = divideRound(xhighpr, 4);
        BigDecimal resul = BigDecimal.ONE;
        /* x^(2i+1) */
        BigDecimal xpowi = BigDecimal.ONE;
        /* i factorial */
        BigInteger ifacN = BigInteger.ONE;
        BigInteger ifacD = BigInteger.ONE;

        for (int i = 1;; i++) {
            ifacN = ifacN.multiply(new BigInteger("" + (2 * i - 1)));
            ifacD = ifacD.multiply(new BigInteger("" + i));

            if (i == 1) {
                xpowi = xhighprV;
            } else {
                xpowi = multiplyRound(xpowi, xhighprV);
            }
            BigDecimal c = divideRound(multiplyRound(xpowi, ifacN),
                    ifacD.multiply(new BigInteger("" + (2 * i + 1))));
            resul = resul.add(c);
            /* series started 1+x/12+... which yields an estimate of the sums error
             */

            if (Math.abs(c.doubleValue()) < xUlpDbl / 120.) {
                break;
            }

        }
        /* sqrt(2*z)*(1+...)
         */
        xpowi = sqrt(xhighpr.multiply(new BigDecimal(2)));
        resul = multiplyRound(xpowi, resul);
        MathContext mc = new MathContext(resul.precision());
        BigDecimal pihalf = pi(mc).divide(new BigDecimal(2));
        mc = new MathContext(err2prec(resul.doubleValue(), eps));

        return pihalf.subtract(resul, mc);

    } else {
        /* absolute error in the result is err(x)/sqrt(1-x^2) to lowest order
         */
        final double xDbl = x.doubleValue();
        final double xUlpDbl = x.ulp().doubleValue() / 2.;
        final double eps = xUlpDbl / 2. / Math.sqrt(1. - Math.pow(xDbl, 2.));
        final BigDecimal xhighpr = scalePrec(x, 2);
        final BigDecimal xhighprSq = multiplyRound(xhighpr, xhighpr);
        BigDecimal resul = xhighpr.plus();
        /* x^(2i+1) */
        BigDecimal xpowi = xhighpr;
        /* i factorial */
        BigInteger ifacN = BigInteger.ONE;
        BigInteger ifacD = BigInteger.ONE;

        for (int i = 1;; i++) {
            ifacN = ifacN.multiply(new BigInteger("" + (2 * i - 1)));
            ifacD = ifacD.multiply(new BigInteger("" + (2 * i)));
            xpowi = multiplyRound(xpowi, xhighprSq);
            BigDecimal c = divideRound(multiplyRound(xpowi, ifacN),
                    ifacD.multiply(new BigInteger("" + (2 * i + 1))));
            resul = resul.add(c);

            if (Math.abs(c.doubleValue()) < 0.1 * eps) {
                break;
            }

        }
        MathContext mc = new MathContext(err2prec(resul.doubleValue(), eps));

        return resul.round(mc);

    }
}

From source file:org.nd4j.linalg.util.BigDecimalMath.java

/**
 * The inverse trigonometric tangent./*w  w  w. java 2s.c  om*/
 *
 * @param x the argument.
 * @return the principal value of arctan(x) in radians in the range -pi/2 to +pi/2.
 */
static public BigDecimal atan(final BigDecimal x) {
    if (x.compareTo(BigDecimal.ZERO) < 0) {
        return atan(x.negate()).negate();

    } else if (x.compareTo(BigDecimal.ZERO) == 0) {
        return BigDecimal.ZERO;
    } else if (x.doubleValue() > 0.7 && x.doubleValue() < 3.0) {
        /* Abramowitz-Stegun 4.4.34 convergence acceleration
         * 2*arctan(x) = arctan(2x/(1-x^2)) = arctan(y). x=(sqrt(1+y^2)-1)/y
         * This maps 0<=y<=3 to 0<=x<=0.73 roughly. Temporarily with 2 protectionist digits.
         */
        BigDecimal y = scalePrec(x, 2);
        BigDecimal newx = divideRound(hypot(1, y).subtract(BigDecimal.ONE), y);
        /* intermediate result with too optimistic error estimate*/
        BigDecimal resul = multiplyRound(atan(newx), 2);
        /* absolute error in the result is errx/(1+x^2), where errx = half of the ulp. */

        double eps = x.ulp().doubleValue() / (2.0 * Math.hypot(1.0, x.doubleValue()));
        MathContext mc = new MathContext(err2prec(resul.doubleValue(), eps));

        return resul.round(mc);

    } else if (x.doubleValue() < 0.71) {
        /* Taylor expansion around x=0; Abramowitz-Stegun 4.4.42 */
        final BigDecimal xhighpr = scalePrec(x, 2);
        final BigDecimal xhighprSq = multiplyRound(xhighpr, xhighpr).negate();
        BigDecimal resul = xhighpr.plus();
        /* signed x^(2i+1) */
        BigDecimal xpowi = xhighpr;
        /* absolute error in the result is errx/(1+x^2), where errx = half of the ulp.
         */

        double eps = x.ulp().doubleValue() / (2.0 * Math.hypot(1.0, x.doubleValue()));

        for (int i = 1;; i++) {
            xpowi = multiplyRound(xpowi, xhighprSq);
            BigDecimal c = divideRound(xpowi, 2 * i + 1);
            resul = resul.add(c);

            if (Math.abs(c.doubleValue()) < 0.1 * eps) {
                break;
            }

        }
        MathContext mc = new MathContext(err2prec(resul.doubleValue(), eps));

        return resul.round(mc);

    } else {
        /* Taylor expansion around x=infinity; Abramowitz-Stegun 4.4.42 */
        /* absolute error in the result is errx/(1+x^2), where errx = half of the ulp.
         */
        double eps = x.ulp().doubleValue() / (2.0 * Math.hypot(1.0, x.doubleValue()));
        /* start with the term pi/2; gather its precision relative to the expected result
         */
        MathContext mc = new MathContext(2 + err2prec(3.1416, eps));
        BigDecimal onepi = pi(mc);
        BigDecimal resul = onepi.divide(new BigDecimal(2));
        final BigDecimal xhighpr = divideRound(-1, scalePrec(x, 2));
        final BigDecimal xhighprSq = multiplyRound(xhighpr, xhighpr).negate();
        /* signed x^(2i+1) */
        BigDecimal xpowi = xhighpr;

        for (int i = 0;; i++) {
            BigDecimal c = divideRound(xpowi, 2 * i + 1);
            resul = resul.add(c);

            if (Math.abs(c.doubleValue()) < 0.1 * eps) {
                break;
            }
            xpowi = multiplyRound(xpowi, xhighprSq);

        }
        mc = new MathContext(err2prec(resul.doubleValue(), eps));

        return resul.round(mc);

    }
}

From source file:org.nd4j.linalg.util.BigDecimalMath.java

/**
 * The hyperbolic cosine./*from   ww  w .  j  a  v a2  s . co  m*/
 *
 * @param x The argument.
 * @return The cosh(x) = (exp(x)+exp(-x))/2 .
 */
static public BigDecimal cosh(final BigDecimal x) {
    if (x.compareTo(BigDecimal.ZERO) < 0) {
        return cos(x.negate());
    } else if (x.compareTo(BigDecimal.ZERO) == 0) {
        return BigDecimal.ONE;
    } else {
        if (x.doubleValue() > 1.5) {
            /* cosh^2(x) = 1+ sinh^2(x).
             */
            return hypot(1, sinh(x));

        } else {
            BigDecimal xhighpr = scalePrec(x, 2);
            /* Simple Taylor expansion, sum_{0=1..infinity} x^(2i)/(2i)! */
            BigDecimal resul = BigDecimal.ONE;
            /* x^i */
            BigDecimal xpowi = BigDecimal.ONE;
            /* 2i factorial */
            BigInteger ifac = BigInteger.ONE;
            /* The absolute error in the result is the error in x^2/2 which is x times the error in x.
             */

            double xUlpDbl = 0.5 * x.ulp().doubleValue() * x.doubleValue();
            /* The error in the result is set by the error in x^2/2 itself, xUlpDbl.
             * We need at most k terms to push x^(2k)/(2k)! below this value.
             * x^(2k) < xUlpDbl; (2k)*log(x) < log(xUlpDbl);
             */

            int k = (int) (Math.log(xUlpDbl) / Math.log(x.doubleValue())) / 2;
            /* The individual terms are all smaller than 1, so an estimate of 1.0 for
             * the absolute value will give a safe relative error estimate for the indivdual terms
             */
            MathContext mcTay = new MathContext(err2prec(1., xUlpDbl / k));

            for (int i = 1;; i++) {
                /* TBD: at which precision will 2*i-1 or 2*i overflow?
                 */
                ifac = ifac.multiply(new BigInteger("" + (2 * i - 1)));
                ifac = ifac.multiply(new BigInteger("" + (2 * i)));
                xpowi = xpowi.multiply(xhighpr).multiply(xhighpr);
                BigDecimal corr = xpowi.divide(new BigDecimal(ifac), mcTay);
                resul = resul.add(corr);

                if (corr.abs().doubleValue() < 0.5 * xUlpDbl) {
                    break;
                }

            } /* The error in the result is governed by the error in x itself.
              */
            MathContext mc = new MathContext(err2prec(resul.doubleValue(), xUlpDbl));

            return resul.round(mc);

        }
    }
}

From source file:org.nd4j.linalg.util.BigDecimalMath.java

/**
 * The hyperbolic sine./* ww  w  .  j  av  a2 s.  co m*/
 *
 * @param x the argument.
 * @return the sinh(x) = (exp(x)-exp(-x))/2 .
 */
static public BigDecimal sinh(final BigDecimal x) {
    if (x.compareTo(BigDecimal.ZERO) < 0) {
        return sinh(x.negate()).negate();
    } else if (x.compareTo(BigDecimal.ZERO) == 0) {
        return BigDecimal.ZERO;
    } else {
        if (x.doubleValue() > 2.4) {
            /* Move closer to zero with sinh(2x)= 2*sinh(x)*cosh(x).
             */
            BigDecimal two = new BigDecimal(2);
            BigDecimal xhalf = x.divide(two);

            BigDecimal resul = sinh(xhalf).multiply(cosh(xhalf)).multiply(two);
            /* The error in the result is set by the error in x itself.
             * The first derivative of sinh(x) is cosh(x), so the absolute error
             * in the result is cosh(x)*errx, and the relative error is coth(x)*errx = errx/tanh(x)
             */

            double eps = Math.tanh(x.doubleValue());
            MathContext mc = new MathContext(err2prec(0.5 * x.ulp().doubleValue() / eps));

            return resul.round(mc);

        } else {
            BigDecimal xhighpr = scalePrec(x, 2);
            /* Simple Taylor expansion, sum_{i=0..infinity} x^(2i+1)/(2i+1)! */
            BigDecimal resul = xhighpr;
            /* x^i */
            BigDecimal xpowi = xhighpr;
            /* 2i+1 factorial */
            BigInteger ifac = BigInteger.ONE;
            /* The error in the result is set by the error in x itself.
             */

            double xUlpDbl = x.ulp().doubleValue();
            /* The error in the result is set by the error in x itself.
             * We need at most k terms to squeeze x^(2k+1)/(2k+1)! below this value.
             * x^(2k+1) < x.ulp; (2k+1)*log10(x) < -x.precision; 2k*log10(x)< -x.precision;
             * 2k*(-log10(x)) > x.precision; 2k*log10(1/x) > x.precision
             */

            int k = (int) (x.precision() / Math.log10(1.0 / xhighpr.doubleValue())) / 2;
            MathContext mcTay = new MathContext(err2prec(x.doubleValue(), xUlpDbl / k));

            for (int i = 1;; i++) {
                /* TBD: at which precision will 2*i or 2*i+1 overflow?
                 */
                ifac = ifac.multiply(new BigInteger("" + (2 * i)));
                ifac = ifac.multiply(new BigInteger("" + (2 * i + 1)));
                xpowi = xpowi.multiply(xhighpr).multiply(xhighpr);
                BigDecimal corr = xpowi.divide(new BigDecimal(ifac), mcTay);
                resul = resul.add(corr);

                if (corr.abs().doubleValue() < 0.5 * xUlpDbl) {
                    break;
                }

            } /* The error in the result is set by the error in x itself.
              */
            MathContext mc = new MathContext(x.precision());

            return resul.round(mc);

        }
    }
}

From source file:org.nd4j.linalg.util.BigDecimalMath.java

/**
 * The hyperbolic tangent./* w  ww.j  a va 2s . c o m*/
 *
 * @param x The argument.
 * @return The tanh(x) = sinh(x)/cosh(x).
 */
static public BigDecimal tanh(final BigDecimal x) {
    if (x.compareTo(BigDecimal.ZERO) < 0) {
        return tanh(x.negate()).negate();
    } else if (x.compareTo(BigDecimal.ZERO) == 0) {
        return BigDecimal.ZERO;
    } else {
        BigDecimal xhighpr = scalePrec(x, 2);
        /* tanh(x) = (1-e^(-2x))/(1+e^(-2x)) .
         */
        BigDecimal exp2x = exp(xhighpr.multiply(new BigDecimal(-2)));
        /* The error in tanh x is err(x)/cosh^2(x).
         */

        double eps = 0.5 * x.ulp().doubleValue() / Math.pow(Math.cosh(x.doubleValue()), 2.0);
        MathContext mc = new MathContext(err2prec(Math.tanh(x.doubleValue()), eps));

        return BigDecimal.ONE.subtract(exp2x).divide(BigDecimal.ONE.add(exp2x), mc);

    }
}

From source file:org.nd4j.linalg.util.BigDecimalMath.java

/**
 * The inverse hyperbolic sine./*from  ww w  .  j  av  a  2s  .  co m*/
 *
 * @param x The argument.
 * @return The arcsinh(x) .
 */
static public BigDecimal asinh(final BigDecimal x) {
    if (x.compareTo(BigDecimal.ZERO) == 0) {
        return BigDecimal.ZERO;
    } else {
        BigDecimal xhighpr = scalePrec(x, 2);
        /* arcsinh(x) = log(x+hypot(1,x))
         */
        BigDecimal logx = log(hypot(1, xhighpr).add(xhighpr));
        /* The absolute error in arcsinh x is err(x)/sqrt(1+x^2)
         */

        double xDbl = x.doubleValue();

        double eps = 0.5 * x.ulp().doubleValue() / Math.hypot(1., xDbl);
        MathContext mc = new MathContext(err2prec(logx.doubleValue(), eps));

        return logx.round(mc);

    }
}

From source file:org.nd4j.linalg.util.BigDecimalMath.java

/**
 * The inverse hyperbolic cosine.//  w  w  w  .  j a  va  2  s . c o  m
 *
 * @param x The argument.
 * @return The arccosh(x) .
 */
static public BigDecimal acosh(final BigDecimal x) {
    if (x.compareTo(BigDecimal.ONE) < 0) {
        throw new ArithmeticException("Out of range argument cosh " + x.toString());
    } else if (x.compareTo(BigDecimal.ONE) == 0) {
        return BigDecimal.ZERO;
    } else {
        BigDecimal xhighpr = scalePrec(x, 2);
        /* arccosh(x) = log(x+sqrt(x^2-1))
         */
        BigDecimal logx = log(sqrt(xhighpr.pow(2).subtract(BigDecimal.ONE)).add(xhighpr));
        /* The absolute error in arcsinh x is err(x)/sqrt(x^2-1)
         */

        double xDbl = x.doubleValue();

        double eps = 0.5 * x.ulp().doubleValue() / Math.sqrt(xDbl * xDbl - 1.);
        MathContext mc = new MathContext(err2prec(logx.doubleValue(), eps));

        return logx.round(mc);

    }
}

From source file:org.nd4j.linalg.util.BigDecimalMath.java

/**
 * The Gamma function.//from w  ww  . j a va 2 s .c  om
 *
 * @param x The argument.
 * @return Gamma(x).
 */
static public BigDecimal Gamma(final BigDecimal x) {
    /* reduce to interval near 1.0 with the functional relation, Abramowitz-Stegun 6.1.33
     */
    if (x.compareTo(BigDecimal.ZERO) < 0) {
        return divideRound(Gamma(x.add(BigDecimal.ONE)), x);
    } else if (x.doubleValue() > 1.5) {
        /* Gamma(x) = Gamma(xmin+n) = Gamma(xmin)*Pochhammer(xmin,n).
         */
        int n = (int) (x.doubleValue() - 0.5);
        BigDecimal xmin1 = x.subtract(new BigDecimal(n));

        return multiplyRound(Gamma(xmin1), pochhammer(xmin1, n));

    } else {
        /* apply Abramowitz-Stegun 6.1.33
         */
        BigDecimal z = x.subtract(BigDecimal.ONE);
        /* add intermediately 2 digits to the partial sum accumulation
         */
        z = scalePrec(z, 2);
        MathContext mcloc = new MathContext(z.precision());
        /* measure of the absolute error is the relative error in the first, logarithmic term
         */

        double eps = x.ulp().doubleValue() / x.doubleValue();
        BigDecimal resul = log(scalePrec(x, 2)).negate();

        if (x.compareTo(BigDecimal.ONE) != 0) {
            BigDecimal gammCompl = BigDecimal.ONE.subtract(gamma(mcloc));
            resul = resul.add(multiplyRound(z, gammCompl));

            for (int n = 2;; n++) {
                /* multiplying z^n/n by zeta(n-1) means that the two relative errors add.
                 * so the requirement in the relative error of zeta(n)-1 is that this is somewhat
                 * smaller than the relative error in z^n/n (the absolute error of thelatter is the
                 * absolute error in z)
                 */
                BigDecimal c = divideRound(z.pow(n, mcloc), n);
                MathContext m = new MathContext(err2prec(n * z.ulp().doubleValue() / 2. / z.doubleValue()));
                c = c.round(m);
                /* At larger n, zeta(n)-1 is roughly 1/2^n. The product is c/2^n.
                 * The relative error in c is c.ulp/2/c . The error in the product should be small versus eps/10.
                 * Error from 1/2^n is c*err(sigma-1).
                 * We need a relative error of zeta-1 of the order of c.ulp/50/c. This is an absolute
                 * error in zeta-1 of c.ulp/50/c/2^n, and also the absolute error in zeta, because zeta is
                 * of the order of 1.
                 */

                if (eps / 100. / c.doubleValue() < 0.01) {
                    m = new MathContext(err2prec(eps / 100. / c.doubleValue()));
                } else {
                    m = new MathContext(2);
                }
                /* zeta(n) -1 */
                BigDecimal zetm1 = zeta(n, m).subtract(BigDecimal.ONE);
                c = multiplyRound(c, zetm1);

                if (n % 2 == 0) {
                    resul = resul.add(c);
                } else {
                    resul = resul.subtract(c);
                }
                /* alternating sum, so truncating as eps is reached suffices
                 */

                if (Math.abs(c.doubleValue()) < eps) {
                    break;
                }

            }
        }
        /* The relative error in the result is the absolute error in the
         * input variable times the digamma (psi) value at that point.
         */
        double psi = 0.5772156649;

        double zdbl = z.doubleValue();

        for (int n = 1; n < 5; n++) {
            psi += zdbl / n / (n + zdbl);
        }
        eps = psi * x.ulp().doubleValue() / 2.;
        mcloc = new MathContext(err2prec(eps));

        return exp(resul).round(mcloc);

    }
}