List of usage examples for java.math MathContext MathContext
public MathContext(String val)
From source file:org.nd4j.linalg.util.BigDecimalMath.java
/** * The trigonometric tangent.//w w w . j a v a2s .com * * @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 ww. j av a2 s. c o 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.// w w w .j av a2 s.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.//from ww w. j av a 2 s .com * * @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 va2s . c o 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.// w ww .j a va 2 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 .ja v a 2 s . com * * @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./* ww w . jav a 2 s . c o 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.//from ww w . j a v a 2 s. c om * * @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 w w. j a va2 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); } }