一个BigDecimal的对数

50

我如何计算一个BigDecimal的对数?有人知道可以使用哪些算法吗?

到目前为止,我的谷歌搜索只提供了(无用的)将BigDecimal转换为double并使用Math.log的想法。

我会提供所需答案的精度。

编辑:任何底数都可以。如果在x进制中更容易,我会采用x进制。


对数的基是什么?2、10、还是e? - paxdiablo
2
任何进制。一旦我有一个实现,进制之间的转换就变得微不足道了。 - masher
我已经在这里提供了解决方案 https://dev59.com/wmfWa4cB1Zd3GeqPcgHw#22556217 - softawareblog.com
我需要这个。有人测试过给出的答案的性能吗? - AD - Stop Putin -
请参见 https://dev59.com/Umw15IYBdhLWcg3wGn9c#7982137。 - leonbloy
11个回答

27

Java数字计算专家:Java程序员的数值计算指南 提供了使用牛顿法的解决方案,书中源代码在此处可获得。以下内容摘自第12.5大十进制函数章节(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);
}

5
为什么不把Math.log()作为第一个近似值? - Peter Lawrey
14
调用 Thread.yield() 应该不在这里。如果你的目的是让一个计算密集型线程成为一个“好公民”,那么你可以将其替换为一些代码来测试线程的“中断”标志并退出。但是调用 Thread.yield() 会干扰正常的线程调度,并且可能会使方法运行非常缓慢...这取决于其他事情发生了什么。 - Stephen C
9
请注意,这个答案不完整,缺少exp()intRoot()的代码。 - Maarten Bodewes
1
你可以使用.precision()代替toString().length()。 - user502187
1
@MaartenBodewes exp()intRoot() https://github.com/javadev/calc/blob/master/src/main/java/com/github/calc/BigDecimalUtil.java - kane
显示剩余2条评论

9
一个非常适用于大数的hacky算法使用关系式log(AB) = log(A) + log(B)。以下是如何在以10为底的情况下进行操作(你可以轻松地将其转换为任何其他对数基数):
  1. 统计答案中的十进制数字数量,这就是对数的整数部分,再加上1。例如:floor(log10(123456)) + 1得到6,因为123456有6个数字。

  2. 如果只需要对数的整数部分,则可以在此处停止:只需从步骤1的结果中减去1即可。

  3. 要获取对数的小数部分,请将数字除以10^(数字的位数),然后使用math.log10()(或其他任何方法;如果没有其他方法可用,则使用简单的级数近似)计算该数字的对数,并将其添加到整数部分。例如:要获取log10(123456)的小数部分,请计算math.log10(0.123456) = -0.908...,并将其加到步骤1的结果上:6 + -0.908 = 5.092,这就是log10(123456)。请注意,您基本上只是在大数前面加上一个小数点;在您的用例中可能有一种很好的优化方法,对于非常大的数字,您甚至不需要担心获取所有数字--log10(0.123)log10(0.123456789)的一个很好的近似值。


2
这种方法为什么不能用于任意精度?你给我一个数字和一个公差,我可以使用该算法计算其对数,并保证绝对误差小于您的公差。我认为这意味着它适用于任意精度。 - kquinn
我的BigInteger简单非优化实现,与这个答案相一致,并且可以推广到BigDecimal,链接在这里https://dev59.com/Umw15IYBdhLWcg3wGn9c#7982137。 - leonbloy

8
这个算法非常快,因为:
  • 没有toString()
  • 没有使用BigInteger函数(牛顿/连分数)
  • 甚至没有实例化一个新的BigInteger
  • 只使用了一定数量的非常快速的运算
一个调用大约需要20微秒(每秒约50,000次调用)
但是:
  • 仅适用于BigInteger
对于BigDecimal的解决办法(未经过速度测试):
  • 将小数点移位,直到值 > 2^53
  • 使用toBigInteger()(内部使用一个div

这个算法利用了对数可以被计算为指数和尾数的对数之和这一事实。比如:

12345 有5位数字,所以以10为底的对数在4和5之间。 log(12345) = 4 + log(1.2345) = 4.09149... (以10为底的对数)


这个函数计算的是以2为底的对数,因为找到占用位数的数量是微不足道的。
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!
}

1
出于好奇,为什么 1.44269504088896340735992468100189213742664595415298D 这么长?Java 的浮点数只有 16 位精度,所以在 Java 中 1.44269504088896340735992468100189213742664595415298D == 1.4426950408889634(大多数具有浮点精度的语言也是如此)。不过我还是可以确认它工作得非常好,所以给个赞。 - Kevin Cruijssen
1
这是Windows计算器给我的结果,而我很懒。 - Mark Jeronimus
@KevinCruijssen 这是十进制数字的整数精度。当谈论小数精度时,情况完全不同,因为使用基于2的分数,其中一些转换为重复。没有小数十进制精度的单个数字,因为基本上没有小数十进制。 - user207421

6
您可以使用分解来实现它。
log(a * 10^b) = log(a) + b * log(10)

基本上,b+1 是数字的位数,a 是 0 到 1 之间的一个值,你可以通过使用常规的 double 算术计算其对数。或者有一些数学技巧可以使用 - 例如,接近 1 的数的对数可以通过级数展开计算。
ln(x + 1) = x - x^2/2 + x^3/3 - x^4/4 + ...

根据你要取对数的数字类型,可能会有类似这样的方法可用。

编辑: 要获取以10为底的对数,可以将自然对数除以ln(10),或者对于任何其他底数也可以使用类似的方法。


我发现一个算法适用于你提供的第一个等式,但第二个等式给出了自然对数。 - masher
哦,是的,我应该提到这一点——这个系列是关于自然对数的。我会进行编辑。 - David Z

4
如果您只需要找出数字中的10的幂次,可以使用以下方法:

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

4
这是我想到的内容:
//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;
}

3

这是一个Meower68伪代码的Java实现,我用几个数字进行了测试:

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

我一直在寻找这个东西,最终采用了连分数的方法。 连分数可以在这里这里找到。

代码:

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

进行对数的伪代码算法。

假设我们想要计算 x 的 log_n。

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;
 }

大的while循环可能看起来有些困惑。
每次循环中,您可以将输入平方,也可以取基数的平方根;无论哪种方式,都必须将分数除以2。我发现将输入平方,并保持基数不变更准确。
如果输入为1,则完成。任何基数的对数为1,这意味着我们不需要再添加任何内容。
如果(result + fraction)不大于result,则我们已经达到了数字系统精度的极限。我们可以停止。
显然,如果您正在使用具有任意多位精度的系统,则需要放入其他内容以限制循环。

1

虽然这是一个旧问题,但我认为这个答案更可取。它具有很高的精度,并支持几乎任何大小的参数。

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;
}

核心逻辑(logBigInteger 方法)是从我另一个答案中复制的。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接