List of usage examples for java.math BigDecimal ulp
public BigDecimal ulp()
From source file:Main.java
public static void main(String[] args) { BigDecimal bg1, bg2, bg3, bg4; bg1 = new BigDecimal("123"); bg2 = new BigDecimal("1.23"); // assign ulp value of bg1, bg2 to bg3, bg4 bg3 = bg1.ulp(); bg4 = bg2.ulp();/* ww w .j av a 2 s.c o m*/ System.out.println("ULP value of " + bg1 + " is " + bg3); System.out.println("ULP value of " + bg2 + " is " + bg4); }
From source file:com.github.jessemull.microflex.stat.statbigdecimal.GeometricMeanBigDecimalTest.java
/** * Corrects any rounding errors due to differences in the implementation of * the statistic between the Apache and MicroFlex libraries * @param BigDecimal the first result * @param BigDecimal the second result * @return corrected results */// w ww. j av a2 s . c o m private static BigDecimal[] correctRoundingErrors(BigDecimal bd1, BigDecimal bd2) { BigDecimal[] array = new BigDecimal[2]; int scale = mc.getPrecision(); while (!bd1.equals(bd2) && scale > mc.getPrecision() / 4) { bd1 = bd1.setScale(scale, RoundingMode.HALF_DOWN); bd2 = bd2.setScale(scale, RoundingMode.HALF_DOWN); if (bd1.subtract(bd1.ulp()).equals(bd2)) { bd1 = bd1.subtract(bd1.ulp()); } if (bd1.add(bd1.ulp()).equals(bd2)) { bd1 = bd1.add(bd1.ulp()); } scale--; } array[0] = bd1; array[1] = bd2; return array; }
From source file:org.nd4j.linalg.util.BigDecimalMath.java
/** * The integer root./*from w w w. j a va 2s . com*/ * * @param n the positive argument. * @param x the non-negative argument. * @return The n-th root of the BigDecimal rounded to the precision implied by x, x^(1/n). */ static public BigDecimal root(final int n, final BigDecimal x) { if (x.compareTo(BigDecimal.ZERO) < 0) { throw new ArithmeticException("negative argument " + x.toString() + " of root"); } if (n <= 0) { throw new ArithmeticException("negative power " + n + " of root"); } if (n == 1) { return x; } /* start the computation from a double precision estimate */ BigDecimal s = new BigDecimal(Math.pow(x.doubleValue(), 1.0 / n)); /* this creates nth with nominal precision of 1 digit */ final BigDecimal nth = new BigDecimal(n); /* Specify an internal accuracy within the loop which is * slightly larger than what is demanded by eps below. */ final BigDecimal xhighpr = scalePrec(x, 2); MathContext mc = new MathContext(2 + x.precision()); /* Relative accuracy of the result is eps. */ final double eps = x.ulp().doubleValue() / (2 * n * x.doubleValue()); for (;;) { /* s = s -(s/n-x/n/s^(n-1)) = s-(s-x/s^(n-1))/n; test correction s/n-x/s for being * smaller than the precision requested. The relative correction is (1-x/s^n)/n, */ BigDecimal c = xhighpr.divide(s.pow(n - 1), mc); c = s.subtract(c); MathContext locmc = new MathContext(c.precision()); c = c.divide(nth, locmc); s = s.subtract(c); if (Math.abs(c.doubleValue() / s.doubleValue()) < eps) { break; } } return s.round(new MathContext(err2prec(eps))); }
From source file:org.nd4j.linalg.util.BigDecimalMath.java
/** * The hypotenuse./*from w ww .j a v a 2 s. co m*/ * * @param x the first argument. * @param y the second argument. * @return the square root of the sum of the squares of the two arguments, sqrt(x^2+y^2). */ static public BigDecimal hypot(final BigDecimal x, final BigDecimal y) { /* compute x^2+y^2 */ BigDecimal z = x.pow(2).add(y.pow(2)); /* truncate to the precision set by x and y. Absolute error = 2*x*xerr+2*y*yerr, * where the two errors are 1/2 of the ulps. Two intermediate protectio digits. */ BigDecimal zerr = x.abs().multiply(x.ulp()).add(y.abs().multiply(y.ulp())); MathContext mc = new MathContext(2 + err2prec(z, zerr)); /* Pull square root */ z = sqrt(z.round(mc)); /* Final rounding. Absolute error in the square root is (y*yerr+x*xerr)/z, where zerr holds 2*(x*xerr+y*yerr). */ mc = new MathContext(err2prec(z.doubleValue(), 0.5 * zerr.doubleValue() / z.doubleValue())); return z.round(mc); }
From source file:org.nd4j.linalg.util.BigDecimalMath.java
/** * The hypotenuse.// ww w .j a va2s . c o m * * @param n the first argument. * @param x the second argument. * @return the square root of the sum of the squares of the two arguments, sqrt(n^2+x^2). */ static public BigDecimal hypot(final int n, final BigDecimal x) { /* compute n^2+x^2 in infinite precision */ BigDecimal z = (new BigDecimal(n)).pow(2).add(x.pow(2)); /* Truncate to the precision set by x. Absolute error = in z (square of the result) is |2*x*xerr|, * where the error is 1/2 of the ulp. Two intermediate protection digits. * zerr is a signed value, but used only in conjunction with err2prec(), so this feature does not harm. */ double zerr = x.doubleValue() * x.ulp().doubleValue(); MathContext mc = new MathContext(2 + err2prec(z.doubleValue(), zerr)); /* Pull square root */ z = sqrt(z.round(mc)); /* Final rounding. Absolute error in the square root is x*xerr/z, where zerr holds 2*x*xerr. */ mc = new MathContext(err2prec(z.doubleValue(), 0.5 * zerr / z.doubleValue())); return z.round(mc); }
From source file:org.nd4j.linalg.util.BigDecimalMath.java
/** * The exponential function.//w w w.j av a 2 s .co m * * @param x the argument. * @return exp(x). * The precision of the result is implicitly defined by the precision in the argument. * 16 * In particular this means that "Invalid Operation" errors are thrown if catastrophic * cancellation of digits causes the result to have no valid digits left. */ static public BigDecimal exp(BigDecimal x) { /* To calculate the value if x is negative, use exp(-x) = 1/exp(x) */ if (x.compareTo(BigDecimal.ZERO) < 0) { final BigDecimal invx = exp(x.negate()); /* Relative error in inverse of invx is the same as the relative errror in invx. * This is used to define the precision of the result. */ MathContext mc = new MathContext(invx.precision()); return BigDecimal.ONE.divide(invx, mc); } else if (x.compareTo(BigDecimal.ZERO) == 0) { /* recover the valid number of digits from x.ulp(), if x hits the * zero. The x.precision() is 1 then, and does not provide this information. */ return scalePrec(BigDecimal.ONE, -(int) (Math.log10(x.ulp().doubleValue()))); } else { /* Push the number in the Taylor expansion down to a small * value where TAYLOR_NTERM terms will do. If x<1, the n-th term is of the order * x^n/n!, and equal to both the absolute and relative error of the result * since the result is close to 1. The x.ulp() sets the relative and absolute error * of the result, as estimated from the first Taylor term. * We want x^TAYLOR_NTERM/TAYLOR_NTERM! < x.ulp, which is guaranteed if * x^TAYLOR_NTERM < TAYLOR_NTERM*(TAYLOR_NTERM-1)*...*x.ulp. */ final double xDbl = x.doubleValue(); final double xUlpDbl = x.ulp().doubleValue(); if (Math.pow(xDbl, TAYLOR_NTERM) < TAYLOR_NTERM * (TAYLOR_NTERM - 1.0) * (TAYLOR_NTERM - 2.0) * xUlpDbl) { /* Add TAYLOR_NTERM terms of the Taylor expansion (Eulers sum formula) */ BigDecimal resul = BigDecimal.ONE; /* x^i */ BigDecimal xpowi = BigDecimal.ONE; /* i factorial */ BigInteger ifac = BigInteger.ONE; /* TAYLOR_NTERM terms to be added means we move x.ulp() to the right * for each power of 10 in TAYLOR_NTERM, so the addition wont add noise beyond * whats already in x. */ MathContext mcTay = new MathContext(err2prec(1., xUlpDbl / TAYLOR_NTERM)); for (int i = 1; i <= TAYLOR_NTERM; i++) { ifac = ifac.multiply(new BigInteger("" + i)); xpowi = xpowi.multiply(x); final BigDecimal c = xpowi.divide(new BigDecimal(ifac), mcTay); resul = resul.add(c); if (Math.abs(xpowi.doubleValue()) < i && Math.abs(c.doubleValue()) < 0.5 * xUlpDbl) { break; } } /* exp(x+deltax) = exp(x)(1+deltax) if deltax is <<1. So the relative error * in the result equals the absolute error in the argument. */ MathContext mc = new MathContext(err2prec(xUlpDbl / 2.)); return resul.round(mc); } else { /* Compute exp(x) = (exp(0.1*x))^10. Division by 10 does not lead * to loss of accuracy. */ int exSc = (int) (1.0 - Math.log10(TAYLOR_NTERM * (TAYLOR_NTERM - 1.0) * (TAYLOR_NTERM - 2.0) * xUlpDbl / Math.pow(xDbl, TAYLOR_NTERM)) / (TAYLOR_NTERM - 1.0)); BigDecimal xby10 = x.scaleByPowerOfTen(-exSc); BigDecimal expxby10 = exp(xby10); /* Final powering by 10 means that the relative error of the result * is 10 times the relative error of the base (First order binomial expansion). * This looses one digit. */ MathContext mc = new MathContext(expxby10.precision() - exSc); /* Rescaling the powers of 10 is done in chunks of a maximum of 8 to avoid an invalid operation 17 * response by the BigDecimal.pow library or integer overflow. */ while (exSc > 0) { int exsub = Math.min(8, exSc); exSc -= exsub; MathContext mctmp = new MathContext(expxby10.precision() - exsub + 2); int pex = 1; while (exsub-- > 0) { pex *= 10; } expxby10 = expxby10.pow(pex, mctmp); } return expxby10.round(mc); } } }
From source file:org.nd4j.linalg.util.BigDecimalMath.java
/** * The natural logarithm./*from w w w . java2 s . c o m*/ * * @param x the argument. * @return ln(x). * The precision of the result is implicitly defined by the precision in the argument. */ static public BigDecimal log(BigDecimal x) { /* the value is undefined if x is negative. */ if (x.compareTo(BigDecimal.ZERO) < 0) { throw new ArithmeticException("Cannot take log of negative " + x.toString()); } else if (x.compareTo(BigDecimal.ONE) == 0) { /* log 1. = 0. */ return scalePrec(BigDecimal.ZERO, x.precision() - 1); } else if (Math.abs(x.doubleValue() - 1.0) <= 0.3) { /* The standard Taylor series around x=1, z=0, z=x-1. Abramowitz-Stegun 4.124. * The absolute error is err(z)/(1+z) = err(x)/x. */ BigDecimal z = scalePrec(x.subtract(BigDecimal.ONE), 2); BigDecimal zpown = z; double eps = 0.5 * x.ulp().doubleValue() / Math.abs(x.doubleValue()); BigDecimal resul = z; for (int k = 2;; k++) { zpown = multiplyRound(zpown, z); BigDecimal c = divideRound(zpown, k); if (k % 2 == 0) { resul = resul.subtract(c); } else { resul = resul.add(c); } if (Math.abs(c.doubleValue()) < eps) { break; } } MathContext mc = new MathContext(err2prec(resul.doubleValue(), eps)); return resul.round(mc); } else { final double xDbl = x.doubleValue(); final double xUlpDbl = x.ulp().doubleValue(); /* Map log(x) = log root[r](x)^r = r*log( root[r](x)) with the aim * to move roor[r](x) near to 1.2 (that is, below the 0.3 appearing above), where log(1.2) is roughly 0.2. */ int r = (int) (Math.log(xDbl) / 0.2); /* Since the actual requirement is a function of the value 0.3 appearing above, * we avoid the hypothetical case of endless recurrence by ensuring that r >= 2. */ r = Math.max(2, r); /* Compute r-th root with 2 additional digits of precision */ BigDecimal xhighpr = scalePrec(x, 2); BigDecimal resul = root(r, xhighpr); resul = log(resul).multiply(new BigDecimal(r)); /* error propagation: log(x+errx) = log(x)+errx/x, so the absolute error * in the result equals the relative error in the input, xUlpDbl/xDbl . */ MathContext mc = new MathContext(err2prec(resul.doubleValue(), xUlpDbl / xDbl)); return resul.round(mc); } }
From source file:org.nd4j.linalg.util.BigDecimalMath.java
/** * Power function./*w ww . j a v a2 s.c om*/ * * @param x Base of the power. * @param y Exponent of the power. * @return x^y. * The estimation of the relative error in the result is |log(x)*err(y)|+|y*err(x)/x| */ static public BigDecimal pow(final BigDecimal x, final BigDecimal y) { if (x.compareTo(BigDecimal.ZERO) < 0) { throw new ArithmeticException("Cannot power negative " + x.toString()); } else if (x.compareTo(BigDecimal.ZERO) == 0) { return BigDecimal.ZERO; } else { /* return x^y = exp(y*log(x)) ; */ BigDecimal logx = log(x); BigDecimal ylogx = y.multiply(logx); BigDecimal resul = exp(ylogx); /* The estimation of the relative error in the result is |log(x)*err(y)|+|y*err(x)/x| */ double errR = Math.abs(logx.doubleValue() * y.ulp().doubleValue() / 2.) + Math.abs(y.doubleValue() * x.ulp().doubleValue() / 2. / x.doubleValue()); MathContext mcR = new MathContext(err2prec(1.0, errR)); return resul.round(mcR); } }
From source file:org.nd4j.linalg.util.BigDecimalMath.java
/** * Trigonometric sine./* w w w. ja v a 2 s . c o m*/ * * @param x The argument in radians. * @return sin(x) in the range -1 to 1. */ static public BigDecimal sin(final BigDecimal x) { if (x.compareTo(BigDecimal.ZERO) < 0) { return sin(x.negate()).negate(); } else if (x.compareTo(BigDecimal.ZERO) == 0) { return BigDecimal.ZERO; } else { /* reduce modulo 2pi */ BigDecimal res = mod2pi(x); double errpi = 0.5 * Math.abs(x.ulp().doubleValue()); int val = 2 + err2prec(FastMath.PI, errpi); MathContext mc = new MathContext(val); BigDecimal p = pi(mc); mc = new MathContext(x.precision()); if (res.compareTo(p) > 0) { /* pi<x<=2pi: sin(x)= - sin(x-pi) */ return sin(subtractRound(res, p)).negate(); } else if (res.multiply(new BigDecimal("2")).compareTo(p) > 0) { /* pi/2<x<=pi: sin(x)= sin(pi-x) */ return sin(subtractRound(p, res)); } else { /* for the range 0<=x<Pi/2 one could use sin(2x)=2sin(x)cos(x) * to split this further. Here, use the sine up to pi/4 and the cosine higher up. */ if (res.multiply(new BigDecimal("4")).compareTo(p) > 0) { /* x>pi/4: sin(x) = cos(pi/2-x) */ return cos(subtractRound(p.divide(new BigDecimal("2")), res)); } else { /* Simple Taylor expansion, sum_{i=1..infinity} (-1)^(..)res^(2i+1)/(2i+1)! */ BigDecimal resul = res; /* x^i */ BigDecimal xpowi = res; /* 2i+1 factorial */ BigInteger ifac = BigInteger.ONE; /* The error in the result is set by the error in x itself. */ double xUlpDbl = res.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) (res.precision() / Math.log10(1.0 / res.doubleValue())) / 2; MathContext mcTay = new MathContext(err2prec(res.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(res).multiply(res).negate(); 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. */ mc = new MathContext(res.precision()); return resul.round(mc); } } } /* sin */ }
From source file:org.nd4j.linalg.util.BigDecimalMath.java
/** * Trigonometric cosine./*from w w w .j a v a2 s . c o m*/ * * @param x The argument in radians. * @return cos(x) in the range -1 to 1. */ static public BigDecimal cos(final BigDecimal x) { if (x.compareTo(BigDecimal.ZERO) < 0) { return cos(x.negate()); } else if (x.compareTo(BigDecimal.ZERO) == 0) { return BigDecimal.ONE; } else { /* reduce modulo 2pi */ BigDecimal res = mod2pi(x); double errpi = 0.5 * Math.abs(x.ulp().doubleValue()); int val = +err2prec(FastMath.PI, errpi); MathContext mc = new MathContext(val); BigDecimal p = pi(mc); mc = new MathContext(x.precision()); if (res.compareTo(p) > 0) { /* pi<x<=2pi: cos(x)= - cos(x-pi) */ return cos(subtractRound(res, p)).negate(); } else if (res.multiply(new BigDecimal("2")).compareTo(p) > 0) { /* pi/2<x<=pi: cos(x)= -cos(pi-x) */ return cos(subtractRound(p, res)).negate(); } else { /* for the range 0<=x<Pi/2 one could use cos(2x)= 1-2*sin^2(x) * to split this further, or use the cos up to pi/4 and the sine higher up. throw new ProviderException("Unimplemented cosine ") ; */ if (res.multiply(new BigDecimal("4")).compareTo(p) > 0) { /* x>pi/4: cos(x) = sin(pi/2-x) */ return sin(subtractRound(p.divide(new BigDecimal("2")), res)); } else { /* Simple Taylor expansion, sum_{i=0..infinity} (-1)^(..)res^(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 * res.ulp().doubleValue() * res.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+1)/(2k+1)! below this value. * x^(2k) < xUlpDbl; (2k)*log(x) < log(xUlpDbl); */ int k = (int) (Math.log(xUlpDbl) / Math.log(res.doubleValue())) / 2; 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(res).multiply(res).negate(); 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. */ mc = new MathContext(err2prec(resul.doubleValue(), xUlpDbl)); return resul.round(mc); } } } }