主题:[讨论]x*log(1.0+x)/((1.0+x)-1.0) 的意义?
jstzhurj
[专家分:4680] 发布于 2010-10-14 22:08:00
有必要开一个新帖专门讨论一下。
asymptotic老师在他的帖子《基本功训练 04】浮点数计算过程中的误差》中给出了计算ln(1+x)代码,如下:
! compute ln(1 + x)
function Log_1Px(x)
implicit none
real, intent(in):: x
real:: Log_1Px
if ( 1.0 + x == 1.0 ) then
Log_1Px = x !* (1.0 - 0.5 * x)
else
! the parenthesis is necessary, must not be omitted.
Log_1Px = x * log(1.0 + x) / ((1.0 + x) - 1.0)
end if
return
end function
代码中等式 Log_1Px = x * log(1.0 + x) / ((1.0 + x) - 1.0)的意义何在? 与直接用 Log_1Px = log(1.0 + x) 有什么区别?在什么情况下必须用Log_1Px = x * log(1.0 + x) / ((1.0 + x) - 1.0) ?
回复列表 (共19个回复)
11 楼
jstzhurj [专家分:4680] 发布于 2010-10-15 22:44:00
[quote]在计算机计算时,1.0应该是用0.9999999999.....代替的,这样的话,机器算会发现0.99999997与1.0000001一样的接近0.99999999,但笔算,后者更接近1.0。[/quote]
能不能不要胡说八道?好好学习一下算术去!
12 楼
doctorlive [专家分:800] 发布于 2010-10-16 10:05:00
你这种人就不应该回复你,任何问题,无论回答的对错,都不应该侮辱回答者。回答你,至少说明别人是想了的,至于对错,是水平问题。
13 楼
yeg001 [专家分:14390] 发布于 2010-10-16 11:01:00
doctorlive兄 也是提供了一个想法而已, 大家来这里交流学习的. 昨晚吃西瓜吃撑了, 想了一下这个问题. 抱着抛砖引玉学习观摩的心态, 我就这个计算写写自己的想法吧.
楼主是想讨论 x*log(1.0+x)/((1.0+x)-1.0) 比 log(1.0+x) 精度有没有奉献和使用的必要性.(如果误会了请指正)
首先我们赋值一个实数x0(我们要计算的真值), 但由于实数二进制记数会发生一定截断. 但这部分截断引起的计算误差(与期望计算值之间的误差)与讨论的问题无关, 于是把x0被计算机截断后记录下来的值记为x, 此时真值不再是log(1.0+x0)而是log(1.0+x).
以下讨论使用rel_error代表相对误差, delta代表偏离值.
先考察 1.0+x 由于记录位数的限制. 由于x是小值, 加上1.0之后原本的有效二进制数位被属于1.0的大量占据, 使得实际加到1.0上的有效值偏离x.
即 x=x'+delta ..... (1)
其中x'表示加到1.0上的有效值, 显然[color=FF0000]x'/=0[/color], 否者(1+x)==1(这个等号判断的争论先放下)并非讨论范围; delta是由于进行加法而丢失的小值. 显然由于计算机处理机制delta可能为正(舍去delta)可能为负(向前进位后舍去). 而且当x>0, x'>=0; 当x<0, x'<0.
于是, ((1.0+x)-1.0)=x'
原来直接计算算法的相对误差
rel_error1 = ABS( (log(1+x')-log(1+x))/log(1+x) )
另一种算法的相对误差
rel_error2 = ABS( ([color=FF0000](x/x')[/color]*log(1+x')-log(1+x))/log(1+x) )
上面两式中的log(1+x)是理想真值, 非计算机结果.
把(1)代入第二种算法. 得到
rel_error2 = ABS( (log(1+x')-log(1+x)+[color=FF0000](delta/x')[/color]*log(1+x'))/log(1+x) )
可见两相对误差差别在于后者多了一项 [color=FF0000](delta/x')[/color]*log(1+x')
为了讨论方便, 假设x>0
此时 log(1+x')-log(1+x)=log( (1+x')/(1+x) )=log(1-delta/(1+x))
当delta>0时, log(1-delta/(1+x))<0, 而 [color=FF0000](delta/x')[/color]*log(1+x') > 0
可见, 这两项刚好异号, 也就是说分子绝对值是缩小的.
当delta<0时, log(1-delta/(1+x))>0, 而 [color=FF0000](delta/x')[/color]*log(1+x') < 0
可见, 同样是异号.
结论, 也就是第二种算法能够进一步减少绝对误差(都除以同一个真值故相对误差也就缩小了).
对于x<0的讨论就不进行了, 应该大同小异.
我把doctorlive兄的代码稍微改动一下, 加入一个log(1+x)以作比较.
program Log1
implicit none
integer::i
real:: x
real:: Log_1Px
do i=1,6
read(*,*) x
if ( 1.0 + x == 1.0 ) then
Log_1Px = x !* (1.0 - 0.5 * x)
else
! the parenthesis is necessary, must not be omitted.
Log_1Px = x * log(1.0 + x) / ((1.0 + x) - 1.0)
write(*,*) x,(1.0 + x) - 1.0
end if
write(*,*) Log_1Px, log(1.0+x)
enddo
end program
在cvf debug模式下计算
由于单精度浮点, 有效数字就6~7位, 所以如果输入值太小就直接被舍弃掉了, 于是关心1.e-5, 1.e-6, 1.e-7, 1.e-8 (其实不一定要1.0开头)其他更小的数被前面的if拦截了.
结果: [color=FF0000](不小心在write的地方在1.0后面加了D0, 更正一下结果)[/color]
1.e-5
9.9999997E-06 9.9999997E-06
9.9999497E-06 9.9999497E-06
1.e-6
1.0000000E-06 1.0000000E-06
9.9999954E-07 9.9999954E-07
1.e-7
1.0000000E-07 1.0000000E-07
9.9999994E-08 9.9999994E-08
1.e-8
9.9999999E-09 9.9999999E-09
9.9999999E-09 9.9999999E-09
这几个值没看出两种算法的差别. [color=FF0000]还有一个很有趣的事情是, 我用gfortran算了一下, 结果跟cvf的差别很大, 没有加任何其他开关, 直接gfortran编译. 这个问题再请教一下别人.[/color]
14 楼
yeg001 [专家分:14390] 发布于 2010-10-16 11:30:00
早上整理推导弄了大半个上午, 可能就是几位前辈没有时间慢慢分析的原因.
个人觉得只要不涉及到这么极端的情况, 提高精度求解就算了, 偶做计算真的很少考虑这些底层函数. 有点耿耿于怀的是, 三角函数的精度是不是不够??
15 楼
jstzhurj [专家分:4680] 发布于 2010-10-16 12:02:00
[quote]你这种人就不应该回复你,任何问题,无论回答的对错,都不应该侮辱回答者。回答你,至少说明别人是想了的,至于对错,是水平问题。[/quote]
这种低级错误,实在忍无可忍,你是有了先入为主的概念,硬往上面凑,有没有认认真真的算呢?有没有量级的概念呢?
16 楼
jstzhurj [专家分:4680] 发布于 2010-10-16 12:13:00
[quote]doctorlive兄 也是提供了一个想法而已, 大家来这里交流学习的. 昨晚吃西瓜吃撑了, 想了一下这个问题. 抱着抛砖引玉学习观摩的心态, 我就这个计算写写自己的想法吧.
楼主是想讨论 x*log(1.0+x)/((1.0+x)-1.0) 比 log(1.0+x) 精度有没有奉献和使用的必要性.(如果误会了请指正)
首先我们赋值一个实数x0(我们要计算的真值), 但由于实数二进制记数会发生一定截断. 但这部分截断引起的计算误差(与期望计算值之间的误差)与讨论的问题无关, 于是把x0被计算机截断后记录下来的值记为x, 此时真值不再是log(1.0+x0)而是log(1.0+x).
以下讨论使用rel_error代表相对误差, delta代表偏离值.
先考察 1.0+x 由于记录位数的限制. 由于x是小值, 加上1.0之后原本的有效二进制数位被属于1.0的大量占据, 使得实际加到1.0上的有效值偏离x.
即 x=x'+delta ..... (1)
其中x'表示加到1.0上的有效值, 显然[color=FF0000]x'/=0[/color], 否者(1+x)==1(这个等号判断的争论先放下)并非讨论范围; delta是由于进行加法而丢失的小值. 显然由于计算机处理机制delta可能为正(舍去delta)可能为负(向前进位后舍去). 而且当x>0, x'>=0; 当x<0, x'<0.
于是, ((1.0+x)-1.0)=x'
原来直接计算算法的相对误差
rel_error1 = ABS( (log(1+x')-log(1+x))/log(1+x) )
另一种算法的相对误差
rel_error2 = ABS( ([color=FF0000](x/x')[/color]*log(1+x')-log(1+x))/log(1+x) )
上面两式中的log(1+x)是理想真值, 非计算机结果.
把(1)代入第二种算法. 得到
rel_error2 = ABS( (log(1+x')-log(1+x)+[color=FF0000](delta/x')[/color]*log(1+x'))/log(1+x) )
可见两相对误差差别在于后者多了一项 [color=FF0000](delta/x')[/color]*log(1+x')
为了讨论方便, 假设x>0
此时 log(1+x')-log(1+x)=log( (1+x')/(1+x) )=log(1-delta/(1+x))
当delta>0时, log(1-delta/(1+x))<0, 而 [color=FF0000](delta/x')[/color]*log(1+x') > 0
可见, 这两项刚好异号, 也就是说分子绝对值是缩小的.
当delta<0时, log(1-delta/(1+x))>0, 而 [color=FF0000](delta/x')[/color]*log(1+x') < 0
可见, 同样是异号.
结论, 也就是第二种算法能够进一步减少绝对误差(都除以同一个真值故相对误差也就缩小了).
对于x<0的讨论就不进行了, 应该大同小异.
我把doctorlive兄的代码稍微改动一下, 加入一个log(1+x)以作比较.
program Log1
implicit none
integer::i
real:: x
real:: Log_1Px
do i=1,6
read(*,*) x
if ( 1.0 + x == 1.0 ) then
Log_1Px = x !* (1.0 - 0.5 * x)
else
! the parenthesis is necessary, must not be omitted.
Log_1Px = x * log(1.0 + x) / ((1.0 + x) - 1.0)
write(*,*) x,(1.0 + x) - 1.0
end if
write(*,*) Log_1Px, log(1.0+x)
enddo
end program
在cvf debug模式下计算
由于单精度浮点, 有效数字就6~7位, 所以如果输入值太小就直接被舍弃掉了, 于是关心1.e-5, 1.e-6, 1.e-7, 1.e-8 (其实不一定要1.0开头)其他更小的数被前面的if拦截了.
结果: [color=FF0000](不小心在write的地方在1.0后面加了D0, 更正一下结果)[/color]
1.e-5
9.9999997E-06 9.9999997E-06
9.9999497E-06 9.9999497E-06
1.e-6
1.0000000E-06 1.0000000E-06
9.9999954E-07 9.9999954E-07
1.e-7
1.0000000E-07 1.0000000E-07
9.9999994E-08 9.9999994E-08
1.e-8
9.9999999E-09 9.9999999E-09
9.9999999E-09 9.9999999E-09
这几个值没看出两种算法的差别. [color=FF0000]还有一个很有趣的事情是, 我用gfortran算了一下, 结果跟cvf的差别很大, 没有加任何其他开关, 直接gfortran编译. 这个问题再请教一下别人.[/color][/quote]
你说的完全对,我的本意就是讨论这个式子对精度有没有贡献?有没有必要使用?不妨用双进度,甚至四倍精度再来比较一下,再把输出结果放大10^N(若干)看看。
17 楼
yeg001 [专家分:14390] 发布于 2010-10-16 13:19:00
我觉得是否碰到这个用法的好处要看delta是否等于0, 可能选取的几个数比较特别, 截断的时候delta等于0所以么显现出来. 那就要另外花时间取检查了.
刚才跟dongyunxun网友聊了一下gfortran的事, 可能是编译器机制和函数算法的不同导致的出入.
18 楼
jstzhurj [专家分:4680] 发布于 2010-10-19 12:50:00
[em18]
19 楼
jstzhurj [专家分:4680] 发布于 2010-10-19 23:22:00
有一个问题,谁有能够精确计算log(1+x)的程序,至少保证小数点后面40位甚至更多位数精确?
我来回复