Logarithm of a BigDecimal

Question

How can I calculate the logarithm of a BigDecimal? Does anyone know of any algorithms I can use?

My googling so far has come up with the (useless) idea of just converting to a double and using Math.log.

I will provide the precision of the answer required.

edit: any base will do. If it's easier in base x, I'll do that.

Answer 1

Java Number Cruncher: The Java Programmer's Guide to Numerical Computing provides a solution using Newton's Method. Source code from the book is available here. The following has been taken from chapter 12.5 Big Decimal Functions (p330 & p331):

/**  * 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 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; }  /**  * 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); } 

Answer 2

A hacky little algorithm that works great for large numbers uses the relation log(AB) = log(A) + log(B). Here's how to do it in base 10 (which you can trivially convert to any other logarithm base):

1. Count the number of decimal digits in the answer. That's the integral part of your logarithm, plus one. Example: floor(log10(123456)) + 1 is 6, since 123456 has 6 digits.

2. You can stop here if all you need is the integer part of the logarithm: just subtract 1 from the result of step 1.

3. To get the fractional part of the logarithm, divide the number by 10^(number of digits), then compute the log of that using math.log10() (or whatever; use a simple series approximation if nothing else is available), and add it to the integer part. Example: to get the fractional part of log10(123456), compute math.log10(0.123456) = -0.908..., and add it to the result of step 1: 6 + -0.908 = 5.092, which is log10(123456). Note that you're basically just tacking on a decimal point to the front of the large number; there is probably a nice way to optimize this in your use case, and for really big numbers you don't even need to bother with grabbing all of the digits -- log10(0.123) is a great approximation to log10(0.123456789).