50

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.

5

11 Answers 11

27

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);
}
7
  • 5
    Why not use Math.log() as the first approximation? Jan 15, 2010 at 23:54
  • 14
    The call to Thread.yield() should not be there. If your aim is to make a computationally intensive thread a "good citizen", then you could replace it with some code to test the Thread's "interrupted" flag and bail out. But a call to Thread.yield() interferes with normal thread scheduling and could make the method run very slowly ... depending on what else is going on.
    – Stephen C
    Dec 19, 2011 at 1:17
  • 9
    Note that this answer is not complete, the code for exp() and intRoot() is missing. Jan 31, 2013 at 0:06
  • 1
    You can used .precision() instead of toString().length()
    – user502187
    Feb 15, 2017 at 1:56
  • 1
    @MaartenBodewes exp() and intRoot() github.com/javadev/calc/blob/master/src/main/java/com/github/…
    – kane
    Oct 23, 2019 at 21:51
9

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).

2
  • 2
    How does this approach not work for arbitrary precision? You give me a number and a tolerance, and I can use that algorithm to calculate its logarithm, with absolute error guaranteed to be less than your tolerance. I'd say that means it works for arbitrary precision.
    – kquinn
    Apr 14, 2009 at 2:02
  • My simple non-optimized implementation for BigInteger, in tune with this answer, and generalizable to BigDecimal, here stackoverflow.com/questions/6827516/logarithm-for-biginteger/…
    – leonbloy
    Nov 2, 2011 at 15:08
8

This one is super fast, because:

  • No toString()
  • No BigInteger math (Newton's/Continued fraction)
  • Not even instantiating a new BigInteger
  • Only uses a fixed number of very fast operations

One call takes about 20 microseconds (about 50k calls per second)

But:

  • Only works for BigInteger

Workaround for BigDecimal (not tested for speed):

  • Shift the decimal point until the value is > 2^53
  • Use toBigInteger() (uses one div internally)

This algorithm makes use of the fact that the log can be calculated as the sum of the exponent and the log of the mantissa. eg:

12345 has 5 digits, so the base 10 log is between 4 and 5. log(12345) = 4 + log(1.2345) = 4.09149... (base 10 log)


This function calculates base 2 log because finding the number of occupied bits is trivial.

public double log(BigInteger val)
{
    // Get the minimum number of bits necessary to hold this value.
    int n = val.bitLength();

    // Calculate the double-precision fraction of this number; as if the
    // binary point was left of the most significant '1' bit.
    // (Get the most significant 53 bits and divide by 2^53)
    long mask = 1L << 52; // mantissa is 53 bits (including hidden bit)
    long mantissa = 0;
    int j = 0;
    for (int i = 1; i < 54; i++)
    {
        j = n - i;
        if (j < 0) break;

        if (val.testBit(j)) mantissa |= mask;
        mask >>>= 1;
    }
    // Round up if next bit is 1.
    if (j > 0 && val.testBit(j - 1)) mantissa++;

    double f = mantissa / (double)(1L << 52);

    // Add the logarithm to the number of bits, and subtract 1 because the
    // number of bits is always higher than necessary for a number
    // (ie. log2(val)<n for every val).
    return (n - 1 + Math.log(f) * 1.44269504088896340735992468100189213742664595415298D);
    // Magic number converts from base e to base 2 before adding. For other
    // bases, correct the result, NOT this number!
}
3
  • 1
    Out of curiosity, why is 1.44269504088896340735992468100189213742664595415298D so long? Java's floating points only have a precision of 16, so 1.44269504088896340735992468100189213742664595415298D == 1.4426950408889634 in Java (and most languages with floating point precision). Can still confirm it's working pretty well though, so +1 from me. Apr 24, 2018 at 11:53
  • 1
    It's what windows calculator gave me, and I'm lazy Apr 25, 2018 at 14:11
  • @KevinCruijssen That's integer precision in decimal digits. When talking about fractional precision it's a completely different ballgame, because of using base-2 fractions, some of which convert to repeating. There is no single figure for fractional decimal precision, because basically there is no fractional decimal.
    – user207421
    Oct 23, 2019 at 23:07
5

You could decompose it using

log(a * 10^b) = log(a) + b * log(10)

Basically b+1 is going to be the number of digits in the number, and a will be a value between 0 and 1 which you could compute the logarithm of by using regular double arithmetic.

Or there are mathematical tricks you can use - for instance, logarithms of numbers close to 1 can be computed by a series expansion

ln(x + 1) = x - x^2/2 + x^3/3 - x^4/4 + ...

Depending on what kind of number you're trying to take the logarithm of, there may be something like this you can use.

EDIT: To get the logarithm in base 10, you can divide the natural logarithm by ln(10), or similarly for any other base.

2
  • I found an algorithm that works on the first equn you give, but the second gives the natural log.
    – masher
    Apr 14, 2009 at 0:04
  • oops, yeah, I should have mentioned that - the series is for the natural log. I'll make an edit.
    – David Z
    Apr 14, 2009 at 0:37
4

This is what I've come up with:

//http://everything2.com/index.pl?node_id=946812        
public BigDecimal log10(BigDecimal b, int dp)
{
    final int NUM_OF_DIGITS = dp+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),dp)).negate();

    StringBuffer sb = new StringBuffer();
    //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() + dp, RoundingMode.HALF_EVEN));
    return ans;
}
4

If all you need is to find the powers of 10 in the number you can use:

public int calculatePowersOf10(BigDecimal value)
{
    return value.round(new MathContext(1)).scale() * -1;
}
3

A Java implementation of Meower68 pseudcode which I tested with a few numbers:

public static BigDecimal log(int base_int, BigDecimal x) {
        BigDecimal result = BigDecimal.ZERO;

        BigDecimal input = new BigDecimal(x.toString());
        int decimalPlaces = 100;
        int scale = input.precision() + decimalPlaces;

        int maxite = 10000;
        int ite = 0;
        BigDecimal maxError_BigDecimal = new BigDecimal(BigInteger.ONE,decimalPlaces + 1);
        System.out.println("maxError_BigDecimal " + maxError_BigDecimal);
        System.out.println("scale " + scale);

        RoundingMode a_RoundingMode = RoundingMode.UP;

        BigDecimal two_BigDecimal = new BigDecimal("2");
        BigDecimal base_BigDecimal = new BigDecimal(base_int);

        while (input.compareTo(base_BigDecimal) == 1) {
            result = result.add(BigDecimal.ONE);
            input = input.divide(base_BigDecimal, scale, a_RoundingMode);
        }

        BigDecimal fraction = new BigDecimal("0.5");
        input = input.multiply(input);
        BigDecimal resultplusfraction = result.add(fraction);
        while (((resultplusfraction).compareTo(result) == 1)
                && (input.compareTo(BigDecimal.ONE) == 1)) {
            if (input.compareTo(base_BigDecimal) == 1) {
                input = input.divide(base_BigDecimal, scale, a_RoundingMode);
                result = result.add(fraction);
            }
            input = input.multiply(input);
            fraction = fraction.divide(two_BigDecimal, scale, a_RoundingMode);
            resultplusfraction = result.add(fraction);
            if (fraction.abs().compareTo(maxError_BigDecimal) == -1){
                break;
            }
            if (maxite == ite){
                break;
            }
            ite ++;
        }

        MathContext a_MathContext = new MathContext(((decimalPlaces - 1) + (result.precision() - result.scale())),RoundingMode.HALF_UP);
        BigDecimal roundedResult = result.round(a_MathContext);
        BigDecimal strippedRoundedResult = roundedResult.stripTrailingZeros();
        //return result;
        //return result.round(a_MathContext);
        return strippedRoundedResult;
    }
2

I was searching for this exact thing and eventually went with a continued fraction approach. The continued fraction can be found at here or here

Code:

import java.math.BigDecimal;
import java.math.MathContext;

public static long ITER = 1000;
public static MathContext context = new MathContext( 100 );
public static BigDecimal ln(BigDecimal x) {
    if (x.equals(BigDecimal.ONE)) {
        return BigDecimal.ZERO;
    }

    x = x.subtract(BigDecimal.ONE);
    BigDecimal ret = new BigDecimal(ITER + 1);
    for (long i = ITER; i >= 0; i--) {
    BigDecimal N = new BigDecimal(i / 2 + 1).pow(2);
        N = N.multiply(x, context);
        ret = N.divide(ret, context);

        N = new BigDecimal(i + 1);
        ret = ret.add(N, context);

    }

    ret = x.divide(ret, context);
    return ret;
}
2

Pseudocode algorithm for doing a logarithm.

Assuming we want log_n of x

double fraction, input;
int base;
double result;

result = 0;
base = n;
input = x;

while (input > base){
  result++;
  input /= base;
}
fraction = 1/2;
input *= input;   

while (((result + fraction) > result) && (input > 1)){
  if (input > base){
    input /= base;
    result += fraction;
  }
  input *= input;
  fraction /= 2.0;
 }

The big while loop may seem a bit confusing.

On each pass, you can either square your input or you can take the square root of your base; either way, you must divide your fraction by 2. I find squaring the input, and leaving the base alone, to be more accurate.

If the input goes to 1, we're through. The log of 1, for any base, is 0, which means we don't need to add any more.

if (result + fraction) is not greater than result, then we've hit the limits of precision for our numbering system. We can stop.

Obviously, if you're working with a system which has arbitrarily many digits of precision, you will want to put something else in there to limit the loop.

1

Old question, but I actually think this answer is preferable. It has good precision and supports arguments of practically any size.

private static final double LOG10 = Math.log(10.0);

/**
 * Computes the natural logarithm of a BigDecimal 
 * 
 * @param val Argument: a positive BigDecimal
 * @return Natural logarithm, as in Math.log()
 */
public static double logBigDecimal(BigDecimal val) {
    return logBigInteger(val.unscaledValue()) + val.scale() * Math.log(10.0);
}

private static final double LOG2 = Math.log(2.0);

/**
 * Computes the natural logarithm of a BigInteger. Works for really big
 * integers (practically unlimited)
 * 
 * @param val Argument, positive integer
 * @return Natural logarithm, as in <tt>Math.log()</tt>
 */
public static double logBigInteger(BigInteger val) {
    int blex = val.bitLength() - 1022; // any value in 60..1023 is ok
    if (blex > 0)
        val = val.shiftRight(blex);
    double res = Math.log(val.doubleValue());
    return blex > 0 ? res + blex * LOG2 : res;
}

The core logic (logBigInteger method) is copied from this other answer of mine.

0

I created a function for BigInteger but it can be easily modified for BigDecimal. Decomposing the log and using some properties of the log is what I do but I get only double precision. But it works for any base. :)

public double BigIntLog(BigInteger bi, double base) {
    // Convert the BigInteger to BigDecimal
    BigDecimal bd = new BigDecimal(bi);
    // Calculate the exponent 10^exp
    BigDecimal diviser = new BigDecimal(10);
    diviser = diviser.pow(bi.toString().length()-1);
    // Convert the BigDecimal from Integer to a decimal value
    bd = bd.divide(diviser);
    // Convert the BigDecimal to double
    double bd_dbl = bd.doubleValue();
    // return the log value
    return (Math.log10(bd_dbl)+bi.toString().length()-1)/Math.log10(base);
}

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.