Java Number Cruncher：《 Java数值计算程序员指南》提供了使用牛顿方法的解决方案。这本书的源代码在这里。以下内容摘自第12.5章大十进制函数（p330＆p331）：
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); }}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);}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;}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);}
