List of usage examples for java.math MathContext getPrecision
public int getPrecision()
From source file:com.qcadoo.model.internal.units.UnitConversionImpl.java
private UnitConversionImpl(final String unitFrom, final String unitTo, final BigDecimal ratio, final MathContext mathContext) { Preconditions.checkNotNull(unitFrom); Preconditions.checkNotNull(ratio);/*from w ww . jav a2s . c o m*/ Preconditions.checkNotNull(mathContext); this.unitFrom = unitFrom; this.unitTo = unitTo; this.ratio = ratio.setScale(mathContext.getPrecision(), mathContext.getRoundingMode()); this.mathContext = mathContext; }
From source file:org.nd4j.linalg.util.BigDecimalMath.java
/** * Eulers constant./*from w w w .java 2 s. com*/ * * @param mc The required precision of the result. * @return 3.14159... */ static public BigDecimal pi(final MathContext mc) { /* look it up if possible */ if (mc.getPrecision() < PI.precision()) { return PI.round(mc); } else { /* Broadhurst \protect\vrule width0pt\protect\href{http://arxiv.org/abs/math/9803067}{arXiv:math/9803067} */ int[] a = { 1, 0, 0, -1, -1, -1, 0, 0 }; BigDecimal S = broadhurstBBP(1, 1, a, mc); return multiplyRound(S, 8); } }
From source file:org.nd4j.linalg.util.BigDecimalMath.java
/** * Euler-Mascheroni constant./* www .j a v a 2 s. co m*/ * * @param mc The required precision of the result. * @return 0.577... */ static public BigDecimal gamma(MathContext mc) { /* look it up if possible */ if (mc.getPrecision() < GAMMA.precision()) { return GAMMA.round(mc); } else { double eps = prec2err(0.577, mc.getPrecision()); /* Euler-Stieltjes as shown in Dilcher, Aequat Math 48 (1) (1994) 55-85 14 */ MathContext mcloc = new MathContext(2 + mc.getPrecision()); BigDecimal resul = BigDecimal.ONE; resul = resul.add(log(2, mcloc)); resul = resul.subtract(log(3, mcloc)); /* how many terms: zeta-1 falls as 1/2^(2n+1), so the * terms drop faster than 1/2^(4n+2). Set 1/2^(4kmax+2) < eps. * Leading term zeta(3)/(4^1*3) is 0.017. Leading zeta(3) is 1.2. Log(2) is 0.7 */ int kmax = (int) ((Math.log(eps / 0.7) - 2.) / 4.); mcloc = new MathContext(1 + err2prec(1.2, eps / kmax)); for (int n = 1;; n++) { /* zeta is close to 1. Division of zeta-1 through * 4^n*(2n+1) means divion through roughly 2^(2n+1) */ BigDecimal c = zeta(2 * n + 1, mcloc).subtract(BigDecimal.ONE); BigInteger fourn = new BigInteger("" + (2 * n + 1)); fourn = fourn.shiftLeft(2 * n); c = divideRound(c, fourn); resul = resul.subtract(c); if (c.doubleValue() < 0.1 * eps) { break; } } return resul.round(mc); } }
From source file:org.nd4j.linalg.util.BigDecimalMath.java
/** * The base of the natural logarithm.//from ww w .j a v a 2s. c o m * * @param mc the required precision of the result * @return exp(1) = 2.71828.... */ static public BigDecimal exp(final MathContext mc) { /* look it up if possible */ if (mc.getPrecision() < E.precision()) { return E.round(mc); } else { /* Instantiate a 1.0 with the requested pseudo-accuracy * and delegate the computation to the public method above. */ BigDecimal uni = scalePrec(BigDecimal.ONE, mc.getPrecision()); return exp(uni); } }
From source file:org.nd4j.linalg.util.BigDecimalMath.java
/** * The natural logarithm./*from www . j a v a 2s . c om*/ * * @param n The main argument, a strictly positive integer. * @param mc The requirements on the precision. * @return ln(n). */ static public BigDecimal log(int n, final MathContext mc) { /* the value is undefined if x is negative. */ if (n <= 0) { throw new ArithmeticException("Cannot take log of negative " + n); } else if (n == 1) { return BigDecimal.ZERO; } else if (n == 2) { if (mc.getPrecision() < LOG2.precision()) { return LOG2.round(mc); } else { /* Broadhurst \protect\vrule width0pt\protect\href{http://arxiv.org/abs/math/9803067}{arXiv:math/9803067} * Error propagation: the error in log(2) is twice the error in S(2,-5,...). */ int[] a = { 2, -5, -2, -7, -2, -5, 2, -3 }; BigDecimal S = broadhurstBBP(2, 1, a, new MathContext(1 + mc.getPrecision())); S = S.multiply(new BigDecimal(8)); S = sqrt(divideRound(S, 3)); return S.round(mc); } } else if (n == 3) { /* summation of a series roughly proportional to (7/500)^k. Estimate count * of terms to estimate the precision (drop the favorable additional * 1/k here): 0.013^k <= 10^(-precision), so k*log10(0.013) <= -precision * so k>= precision/1.87. */ int kmax = (int) (mc.getPrecision() / 1.87); MathContext mcloc = new MathContext(mc.getPrecision() + 1 + (int) (Math.log10(kmax * 0.693 / 1.098))); BigDecimal log3 = multiplyRound(log(2, mcloc), 19); /* log3 is roughly 1, so absolute and relative error are the same. The * result will be divided by 12, so a conservative error is the one * already found in mc */ double eps = prec2err(1.098, mc.getPrecision()) / kmax; Rational r = new Rational(7153, 524288); Rational pk = new Rational(7153, 524288); for (int k = 1;; k++) { Rational tmp = pk.divide(k); if (tmp.doubleValue() < eps) { break; } /* how many digits of tmp do we need in the sum? */ mcloc = new MathContext(err2prec(tmp.doubleValue(), eps)); BigDecimal c = pk.divide(k).BigDecimalValue(mcloc); if (k % 2 != 0) { log3 = log3.add(c); } else { log3 = log3.subtract(c); } pk = pk.multiply(r); } log3 = divideRound(log3, 12); return log3.round(mc); } else if (n == 5) { /* summation of a series roughly proportional to (7/160)^k. Estimate count * of terms to estimate the precision (drop the favorable additional * 1/k here): 0.046^k <= 10^(-precision), so k*log10(0.046) <= -precision * so k>= precision/1.33. */ int kmax = (int) (mc.getPrecision() / 1.33); MathContext mcloc = new MathContext(mc.getPrecision() + 1 + (int) (Math.log10(kmax * 0.693 / 1.609))); BigDecimal log5 = multiplyRound(log(2, mcloc), 14); /* log5 is roughly 1.6, so absolute and relative error are the same. The * result will be divided by 6, so a conservative error is the one * already found in mc */ double eps = prec2err(1.6, mc.getPrecision()) / kmax; Rational r = new Rational(759, 16384); Rational pk = new Rational(759, 16384); for (int k = 1;; k++) { Rational tmp = pk.divide(k); if (tmp.doubleValue() < eps) { break; } /* how many digits of tmp do we need in the sum? */ mcloc = new MathContext(err2prec(tmp.doubleValue(), eps)); BigDecimal c = pk.divide(k).BigDecimalValue(mcloc); log5 = log5.subtract(c); pk = pk.multiply(r); } log5 = divideRound(log5, 6); return log5.round(mc); } else if (n == 7) { /* summation of a series roughly proportional to (1/8)^k. Estimate count * of terms to estimate the precision (drop the favorable additional * 1/k here): 0.125^k <= 10^(-precision), so k*log10(0.125) <= -precision * so k>= precision/0.903. */ int kmax = (int) (mc.getPrecision() / 0.903); MathContext mcloc = new MathContext( mc.getPrecision() + 1 + (int) (Math.log10(kmax * 3 * 0.693 / 1.098))); BigDecimal log7 = multiplyRound(log(2, mcloc), 3); /* log7 is roughly 1.9, so absolute and relative error are the same. */ double eps = prec2err(1.9, mc.getPrecision()) / kmax; Rational r = new Rational(1, 8); Rational pk = new Rational(1, 8); for (int k = 1;; k++) { Rational tmp = pk.divide(k); if (tmp.doubleValue() < eps) { break; } /* how many digits of tmp do we need in the sum? */ mcloc = new MathContext(err2prec(tmp.doubleValue(), eps)); BigDecimal c = pk.divide(k).BigDecimalValue(mcloc); log7 = log7.subtract(c); pk = pk.multiply(r); } return log7.round(mc); } else { /* At this point one could either forward to the log(BigDecimal) signature (implemented) * or decompose n into Ifactors and use an implemenation of all the prime bases. * Estimate of the result; convert the mc argument to an absolute error eps * log(n+errn) = log(n)+errn/n = log(n)+eps */ double res = Math.log((double) n); double eps = prec2err(res, mc.getPrecision()); /* errn = eps*n, convert absolute error in result to requirement on absolute error in input */ eps *= n; /* Convert this absolute requirement of error in n to a relative error in n */ final MathContext mcloc = new MathContext(1 + err2prec((double) n, eps)); /* Padd n with a number of zeros to trigger the required accuracy in * the standard signature method */ BigDecimal nb = scalePrec(new BigDecimal(n), mcloc); return log(nb); } }
From source file:org.nd4j.linalg.util.BigDecimalMath.java
/** * The natural logarithm.// w w w. j ava2 s . c om * * @param r The main argument, a strictly positive value. * @param mc The requirements on the precision. * @return ln(r). */ static public BigDecimal log(final Rational r, final MathContext mc) { /* the value is undefined if x is negative. */ if (r.compareTo(Rational.ZERO) <= 0) { throw new ArithmeticException("Cannot take log of negative " + r.toString()); } else if (r.compareTo(Rational.ONE) == 0) { return BigDecimal.ZERO; } else { /* log(r+epsr) = log(r)+epsr/r. Convert the precision to an absolute error in the result. * eps contains the required absolute error of the result, epsr/r. */ double eps = prec2err(Math.log(r.doubleValue()), mc.getPrecision()); /* Convert this further into a requirement of the relative precision in r, given that * epsr/r is also the relative precision of r. Add one safety digit. */ MathContext mcloc = new MathContext(1 + err2prec(eps)); final BigDecimal resul = log(r.BigDecimalValue(mcloc)); return resul.round(mc); } }
From source file:org.nd4j.linalg.util.BigDecimalMath.java
/** * Riemann zeta function./* ww w .j a v a2 s. c o m*/ * * @param n The positive integer argument. * 32 * @param mc Specification of the accuracy of the result. * @return zeta(n). */ 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, (new BigInteger("" + 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, (new BigInteger("" + 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); } }
From source file:org.nd4j.linalg.util.BigDecimalMath.java
/** * Broadhurst ladder sequence.//w w w.j a va 2 s . c om * * @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} */ 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(new BigInteger("" + a[k]), (new BigInteger("" + (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); }
From source file:org.nd4j.linalg.util.BigDecimalMath.java
/** * Boost the precision by appending decimal zeros to the value. This returns a value which appears to have * a higher precision than the input.//w ww. j av a2 s .c o m * * @param x The input value * @param mc The requirement on the minimum precision on return. * @return The same value as the input but with increased (pseudo) precision. */ static public BigDecimal scalePrec(final BigDecimal x, final MathContext mc) { final int diffPr = mc.getPrecision() - x.precision(); if (diffPr > 0) { return scalePrec(x, diffPr); } else { return x; } }