3 回答

TA貢獻(xiàn)1826條經(jīng)驗(yàn) 獲得超6個(gè)贊
Java Number Cruncher:《 Java數(shù)值計(jì)算程序員指南》提供了使用牛頓方法的解決方案。這本書的源代碼在這里。以下內(nèi)容摘自第12.5章大十進(jìn)制函數(shù)(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);
}

TA貢獻(xiàn)1830條經(jīng)驗(yàn) 獲得超3個(gè)贊
一種關(guān)系不好的小算法,適用于大量數(shù)字log(AB) = log(A) + log(B)
。以下是以10為基數(shù)的方法(您可以將其轉(zhuǎn)換為任何其他對(duì)數(shù)基數(shù)):
計(jì)算答案中的小數(shù)位數(shù)。那是對(duì)數(shù)的必不可少的部分,再加上一個(gè)。例如:
floor(log10(123456)) + 1
是6,因?yàn)?23456有6位數(shù)字。如果您需要的只是對(duì)數(shù)的整數(shù)部分,可以在這里停止:只需從步驟1的結(jié)果中減去1。
要獲得對(duì)數(shù)的小數(shù)部分,請(qǐng)將數(shù)字除以
10^(number of digits)
,然后使用math.log10()
(或其他方法;計(jì)算對(duì)數(shù)的對(duì)數(shù))(如果沒有其他可用方法,則使用簡(jiǎn)單的序列近似),然后將其加到整數(shù)部分。示例:要獲取的小數(shù)部分log10(123456)
,請(qǐng)計(jì)算math.log10(0.123456) = -0.908...
并將其添加到步驟1的結(jié)果中6 + -0.908 = 5.092
,即log10(123456)
。請(qǐng)注意,您基本上只是在大數(shù)前面加上小數(shù)點(diǎn);在您的用例中可能有一種優(yōu)化此方法的好方法,對(duì)于很大的數(shù)字,您甚至不必費(fèi)心去抓住所有數(shù)字-log10(0.123)
近似于log10(0.123456789)
。

TA貢獻(xiàn)1876條經(jīng)驗(yàn) 獲得超7個(gè)贊
這是超級(jí)快的,因?yàn)椋?/p>
沒有 toString()
無BigInteger數(shù)學(xué)(牛頓/續(xù)分?jǐn)?shù))
甚至沒有實(shí)例化一個(gè)新的 BigInteger
僅使用固定數(shù)量的非??焖俚牟僮?/p>
一次通話大約需要20微秒(每秒大約5萬(wàn)次通話)
但:
僅適用于 BigInteger
解決方法BigDecimal(未測(cè)試速度):
移動(dòng)小數(shù)點(diǎn),直到值> 2 ^ 53
使用toBigInteger()(div內(nèi)部使用)
該算法利用了可以將對(duì)數(shù)計(jì)算為指數(shù)和尾數(shù)對(duì)數(shù)之和的事實(shí)。例如:
12345有5位數(shù)字,因此以10為底的對(duì)數(shù)介于4和5之間。log(12345)= 4 + log(1.2345)= 4.09149 ...(以10為底的對(duì)數(shù))
該函數(shù)計(jì)算以2為底的對(duì)數(shù),因?yàn)檎业秸加玫奈粩?shù)很簡(jiǎn)單。
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!
}
添加回答
舉報(bào)