请教一个ln(x)的收敛性较高的幂级数展开式

[ls251544415]

用泰勒级数展开式验证了结果
x不大的时候，结果准确性还行
x稍微大一点，结果的正确性就差远了，必须把累加的次数加倍，精度才会稍微提高一点，这样效率很低。

不知道有没有人晓得<math.h>中log()函数的源代码的原理是怎样的

[mLee79]

好像是 log( (1+x) / (1-x) ) 吧...
用 log( x ) 展开也容易的, 稍微动点手脚, 比如直接按浮点数IEEE格式, 提取尾数和阶码, 很容易的就使 .5 <= x < 1 , 收敛很快地...

[ls251544415]

[Quote=引用 1 楼 mlee79 的回复:]
好像是 log( (1+x) / (1-x) ) 吧...
用 log( x ) 展开也容易的, 稍微动点手脚, 比如直接按浮点数IEEE格式, 提取尾数和阶码, 很容易的就使 .5 <= x < 1 , 收敛很快地...
[/Quote]

先前试过了
ln( (1+x)/(1-x) ) = 2∑{ [x^(2n+1)] / (2n+1) } （ n∈[0,∞) ）

(1+x)/(1-x) == 100 且 n==100时
误差达到了0.1

[mLee79]

log( (1+x) / (1-x) ) 只是比 log( x ) 快一点点...

第一步当然是 提取尾数和阶码, 把 x 写成 a * 2^N ( .5 <= a < 1 ) , 这是很容易的事情, 然后结果就是
log( a ) + N * log( 2 ) , 对这个范围的 a 显然不管啥办法都很简单的说...

[Aniao]

如果你对这个要求高的话，建议去搜索专业的数学库

[mengqingsong_2009]

mathlab

[gbb21]

http://www.math.com/tables/expansion/log.htm

[ls251544415]

mLee79，有例子就通俗易懂了
这个转换巧妙的很呐。多谢

[mLee79]

这个精度和速度应该都够了:
#include <stdio.h> 
#include <stdlib.h>
#include <assert.h>
#include <math.h>

#define LOG2		(0.69314718055994529)

#ifdef _MSC_VER
typedef __int64			int64;
#define FRACMASK		(0xFFFFFFFFFFFFFI64)
#else
typedef long long int64;
#define FRACMASK		(0xFFFFFFFFFFFFFLL)
#endif

double mLog( double x )
{
	int    N , i;
	double r = 0. , x0;
	
	assert( x > 0 );
	
	N = (int)((*(int64*)&x) >> 52) - 1022;
	(*(int64*)&x) -= (int64)N << 52;

	if( ((*(int64*)&x) & FRACMASK ) == 0 )
		return (N-1) * LOG2;

	x0 = x = (x-1)/(x+1);	/* -.333... < x0 < 0 */
	for( i = 1; x0 < -1e-18; i+=2 , x0 *= x * x )
		r += x0 / i;
return 2 * r + N * LOG2;
}

int main()
{
	double x , r1 , r2;

	for( x = .000000001; x < 10000000; x *= 10 )
	{
		r1 = mLog( x );
		r2 = log ( x );
		
		printf( "%.8g : %.8g  %.8g  %.8g\n" , x , r1 , r2 , r1 - r2 );
	}

	for( x = .0001234567; x < 10000; x *= 1.234567 )
	{
		r1 = mLog( x );
		r2 = log ( x );

		printf( "%.8g : %.8g  %.8g  %.8g\n" , x , r1 , r2 , r1 - r2 );
	}

return 0;
}

[mLee79]

看起来好像很别扭, 还是这样吧:

double mLog( double x )
{
	int    N , i;
	double r = 0. , x0;
	
	assert( x > 0 );
	
	N = (int)((*(int64*)&x) >> 52) - 1022;
	(*(int64*)&x) -= (int64)N << 52;

	if( ((*(int64*)&x) & FRACMASK ) == 0 )
		return (N-1) * LOG2;

	x0 = x = (1-x)/(1+x);	/* 0 < x < .333... */
	for( i = 1; x0 > 1e-18; i+=2 , x0 *= x * x )
		r += x0 / i;
return 2 * -r + N * LOG2;
}