TA貢獻1934條經(jīng)驗 獲得超2個贊
Java Number Cruncher：《 Java數(shù)值計算程序員指南》提供了使用牛頓方法的解決方案。這本書的源代碼在這里。以下內(nèi)容摘自第12.5章大十進制函數(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貢獻1794條經(jīng)驗 獲得超7個贊
一種關(guān)系不好的小算法，適用于大量數(shù)字log(AB) = log(A) + log(B)。以下是以10為基數(shù)的方法（您可以將其轉(zhuǎn)換為任何其他對數(shù)基數(shù)）：
計算答案中的小數(shù)位數(shù)。那是對數(shù)的必不可少的部分，再加上一個。例如：floor(log10(123456)) + 1是6，因為123456有6位數(shù)字。
如果您需要的只是對數(shù)的整數(shù)部分，可以在這里停止：只需從步驟1的結(jié)果中減去1。
要獲得對數(shù)的小數(shù)部分，請將數(shù)字除以10^(number of digits)，然后使用math.log10()（或其他方法；計算對數(shù)的對數(shù)）（如果沒有其他可用方法，則使用簡單的序列近似），然后將其加到整數(shù)部分。示例：要獲取的小數(shù)部分log10(123456)，請計算math.log10(0.123456) = -0.908...并將其添加到步驟1的結(jié)果中6 + -0.908 = 5.092，即log10(123456)。請注意，您基本上只是在大數(shù)前面加上小數(shù)點；在您的用例中可能有一種優(yōu)化此方法的好方法，對于很大的數(shù)字，您甚至不必費心去抓住所有數(shù)字- log10(0.123)近似于log10(0.123456789)。
TA貢獻1874條經(jīng)驗 獲得超12個贊
這是超級快的，因為：
沒有 toString()
無BigInteger數(shù)學(xué)（牛頓/續(xù)分?jǐn)?shù)）
甚至沒有實例化一個新的 BigInteger
僅使用固定數(shù)量的非?？焖俚牟僮?/p>
一次通話大約需要20微秒（每秒大約5萬次通話）
但：
僅適用于 BigInteger
解決方法BigDecimal（未測試速度）：
移動小數(shù)點，直到值> 2 ^ 53
使用toBigInteger()（div內(nèi)部使用）
該算法利用了可以將對數(shù)計算為指數(shù)和尾數(shù)對數(shù)之和的事實。例如：
12345有5位數(shù)字，因此以10為底的對數(shù)介于4和5之間。log（12345）= 4 + log（1.2345）= 4.09149 ...（以10為底的對數(shù)）
該函數(shù)計算以2為底的對數(shù)，因為找到占用的位數(shù)很簡單。
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!
}
舉報