Java tutorial
//package com.java2s; /* * Copyright 2013 Valentyn Kolesnikov * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. */ import java.math.BigDecimal; import java.math.MathContext; import java.math.RoundingMode; public class Main { private static final int SCALE = 18; public static BigDecimal log10(BigDecimal b) { final int NUM_OF_DIGITS = SCALE + 2; // need to add one to get the right number of dp // and then add one again to get the next number // so I can round it correctly. MathContext mc = new MathContext(NUM_OF_DIGITS, RoundingMode.HALF_EVEN); // special conditions: // log(-x) -> exception // log(1) == 0 exactly; // log of a number lessthan one = -log(1/x) if (b.signum() <= 0) { throw new ArithmeticException("log of a negative number! (or zero)"); } else if (b.compareTo(BigDecimal.ONE) == 0) { return BigDecimal.ZERO; } else if (b.compareTo(BigDecimal.ONE) < 0) { return (log10((BigDecimal.ONE).divide(b, mc))).negate(); } StringBuilder sb = new StringBuilder(); // number of digits on the left of the decimal point int leftDigits = b.precision() - b.scale(); // so, the first digits of the log10 are: sb.append(leftDigits - 1).append("."); // this is the algorithm outlined in the webpage int n = 0; while (n < NUM_OF_DIGITS) { b = (b.movePointLeft(leftDigits - 1)).pow(10, mc); leftDigits = b.precision() - b.scale(); sb.append(leftDigits - 1); n++; } BigDecimal ans = new BigDecimal(sb.toString()); // Round the number to the correct number of decimal places. ans = ans.round(new MathContext(ans.precision() - ans.scale() + SCALE, RoundingMode.HALF_EVEN)); return ans; } public static BigDecimal pow(BigDecimal savedValue, BigDecimal value) { BigDecimal result = null; result = exp(ln(savedValue, 32).multiply(value), 32); return result; } /** * Compute e^x to a given scale. Break x into its whole and fraction parts * and compute (e^(1 + fraction/whole))^whole using Taylor's formula. * * @param x * the value of x * @param scale * the desired scale of the result * @return the result value */ public static BigDecimal exp(BigDecimal x, int scale) { // e^0 = 1 if (x.signum() == 0) { return BigDecimal.valueOf(1); } // If x is negative, return 1/(e^-x). else if (x.signum() == -1) { return BigDecimal.valueOf(1).divide(exp(x.negate(), scale), scale, BigDecimal.ROUND_HALF_EVEN); } // Compute the whole part of x. BigDecimal xWhole = x.setScale(0, BigDecimal.ROUND_DOWN); // If there isn't a whole part, compute and return e^x. if (xWhole.signum() == 0) { return expTaylor(x, scale); } // Compute the fraction part of x. BigDecimal xFraction = x.subtract(xWhole); // z = 1 + fraction/whole BigDecimal z = BigDecimal.valueOf(1).add(xFraction.divide(xWhole, scale, BigDecimal.ROUND_HALF_EVEN)); // t = e^z BigDecimal t = expTaylor(z, scale); BigDecimal maxLong = BigDecimal.valueOf(Long.MAX_VALUE); BigDecimal result = BigDecimal.valueOf(1); // Compute and return t^whole using intPower(). // If whole > Long.MAX_VALUE, then first compute products // of e^Long.MAX_VALUE. while (xWhole.compareTo(maxLong) >= 0) { result = result.multiply(intPower(t, Long.MAX_VALUE, scale)).setScale(scale, BigDecimal.ROUND_HALF_EVEN); xWhole = xWhole.subtract(maxLong); Thread.yield(); } return result.multiply(intPower(t, xWhole.longValue(), scale)).setScale(scale, BigDecimal.ROUND_HALF_EVEN); } /** * Compute the natural logarithm of x to a given scale, x > 0. */ public static BigDecimal ln(BigDecimal x, int scale) { // Check that x > 0. if (x.signum() <= 0) { throw new IllegalArgumentException("x <= 0"); } // The number of digits to the left of the decimal point. int magnitude = x.toString().length() - x.scale() - 1; if (magnitude < 3) { return lnNewton(x, scale); } // Compute magnitude*ln(x^(1/magnitude)). else { // x^(1/magnitude) BigDecimal root = intRoot(x, magnitude, scale); // ln(x^(1/magnitude)) BigDecimal lnRoot = lnNewton(root, scale); // magnitude*ln(x^(1/magnitude)) return BigDecimal.valueOf(magnitude).multiply(lnRoot).setScale(scale, BigDecimal.ROUND_HALF_EVEN); } } /** * Compute e^x to a given scale by the Taylor series. * * @param x * the value of x * @param scale * the desired scale of the result * @return the result value */ private static BigDecimal expTaylor(BigDecimal x, int scale) { BigDecimal factorial = BigDecimal.valueOf(1); BigDecimal xPower = x; BigDecimal sumPrev; // 1 + x BigDecimal sum = x.add(BigDecimal.valueOf(1)); // Loop until the sums converge // (two successive sums are equal after rounding). int i = 2; do { // x^i xPower = xPower.multiply(x).setScale(scale, BigDecimal.ROUND_HALF_EVEN); // i! factorial = factorial.multiply(BigDecimal.valueOf(i)); // x^i/i! BigDecimal term = xPower.divide(factorial, scale, BigDecimal.ROUND_HALF_EVEN); // sum = sum + x^i/i! sumPrev = sum; sum = sum.add(term); ++i; Thread.yield(); } while (sum.compareTo(sumPrev) != 0); return sum; } /** * Compute x^exponent to a given scale. Uses the same algorithm as class * numbercruncher.mathutils.IntPower. * * @param x * the value x * @param exponent * the exponent value * @param scale * the desired scale of the result * @return the result value */ public static BigDecimal intPower(BigDecimal x, long exponent, int scale) { // If the exponent is negative, compute 1/(x^-exponent). if (exponent < 0) { return BigDecimal.valueOf(1).divide(intPower(x, -exponent, scale), scale, BigDecimal.ROUND_HALF_EVEN); } BigDecimal power = BigDecimal.valueOf(1); // Loop to compute value^exponent. while (exponent > 0) { // Is the rightmost bit a 1? if ((exponent & 1) == 1) { power = power.multiply(x).setScale(scale, BigDecimal.ROUND_HALF_EVEN); } // Square x and shift exponent 1 bit to the right. x = x.multiply(x).setScale(scale, BigDecimal.ROUND_HALF_EVEN); exponent >>= 1; Thread.yield(); } return power; } /** * Compute the natural logarithm of x to a given scale, x > 0. Use Newton's * algorithm. */ private static BigDecimal lnNewton(BigDecimal x, int scale) { int sp1 = scale + 1; BigDecimal n = x; BigDecimal term; // Convergence tolerance = 5*(10^-(scale+1)) BigDecimal tolerance = BigDecimal.valueOf(5).movePointLeft(sp1); // Loop until the approximations converge // (two successive approximations are within the tolerance). do { // e^x BigDecimal eToX = exp(x, sp1); // (e^x - n)/e^x term = eToX.subtract(n).divide(eToX, sp1, BigDecimal.ROUND_DOWN); // x - (e^x - n)/e^x x = x.subtract(term); Thread.yield(); } while (term.compareTo(tolerance) > 0); return x.setScale(scale, BigDecimal.ROUND_HALF_EVEN); } /** * Compute the integral root of x to a given scale, x >= 0. Use Newton's * algorithm. * * @param x * the value of x * @param index * the integral root value * @param scale * the desired scale of the result * @return the result value */ public static BigDecimal intRoot(BigDecimal x, long index, int scale) { // Check that x >= 0. if (x.signum() < 0) { throw new IllegalArgumentException("x < 0"); } int sp1 = scale + 1; BigDecimal n = x; BigDecimal i = BigDecimal.valueOf(index); BigDecimal im1 = BigDecimal.valueOf(index - 1); BigDecimal tolerance = BigDecimal.valueOf(5).movePointLeft(sp1); BigDecimal xPrev; // The initial approximation is x/index. x = x.divide(i, scale, BigDecimal.ROUND_HALF_EVEN); // Loop until the approximations converge // (two successive approximations are equal after rounding). do { // x^(index-1) BigDecimal xToIm1 = intPower(x, index - 1, sp1); // x^index BigDecimal xToI = x.multiply(xToIm1).setScale(sp1, BigDecimal.ROUND_HALF_EVEN); // n + (index-1)*(x^index) BigDecimal numerator = n.add(im1.multiply(xToI)).setScale(sp1, BigDecimal.ROUND_HALF_EVEN); // (index*(x^(index-1)) BigDecimal denominator = i.multiply(xToIm1).setScale(sp1, BigDecimal.ROUND_HALF_EVEN); // x = (n + (index-1)*(x^index)) / (index*(x^(index-1))) xPrev = x; x = numerator.divide(denominator, sp1, BigDecimal.ROUND_DOWN); Thread.yield(); } while (x.subtract(xPrev).abs().compareTo(tolerance) > 0); return x; } }