回 帖 发 新 帖 刷新版面

主题:[讨论]x*log(1.0+x)/((1.0+x)-1.0) 的意义?

有必要开一个新帖专门讨论一下。

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 楼

[quote]在计算机计算时,1.0应该是用0.9999999999.....代替的,这样的话,机器算会发现0.99999997与1.0000001一样的接近0.99999999,但笔算,后者更接近1.0。[/quote]
能不能不要胡说八道?好好学习一下算术去!

12 楼

你这种人就不应该回复你,任何问题,无论回答的对错,都不应该侮辱回答者。回答你,至少说明别人是想了的,至于对错,是水平问题。

13 楼

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 楼

早上整理推导弄了大半个上午, 可能就是几位前辈没有时间慢慢分析的原因.
个人觉得只要不涉及到这么极端的情况, 提高精度求解就算了, 偶做计算真的很少考虑这些底层函数. 有点耿耿于怀的是, 三角函数的精度是不是不够??

15 楼

[quote]你这种人就不应该回复你,任何问题,无论回答的对错,都不应该侮辱回答者。回答你,至少说明别人是想了的,至于对错,是水平问题。[/quote]

这种低级错误,实在忍无可忍,你是有了先入为主的概念,硬往上面凑,有没有认认真真的算呢?有没有量级的概念呢?

16 楼

[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 楼

我觉得是否碰到这个用法的好处要看delta是否等于0, 可能选取的几个数比较特别, 截断的时候delta等于0所以么显现出来. 那就要另外花时间取检查了. 
刚才跟dongyunxun网友聊了一下gfortran的事, 可能是编译器机制和函数算法的不同导致的出入.

18 楼

[em18]

19 楼


有一个问题,谁有能够精确计算log(1+x)的程序,至少保证小数点后面40位甚至更多位数精确?

我来回复

您尚未登录,请登录后再回复。点此登录或注册