我的方法yeg001思想是一样的，我只不过把100.*nd/x并到计算中了！我的程序在nd=365不变的情况下，x约小于89.4运行良好，直至x约大于89.4结果超界。

[quote]yeg001  网友的思想我是看懂了，能很好的解决某在原帖中提出的问题，但若是将 x = 9.06 ，则我电脑 CPU 是无法完成这个任务了。
jstzhurj 网友在 23 楼提供的程序，看来还不错，但我没有看明白算法，请明示。

我更倾向于这样的代码：不仅 x 小的时候可以，比如 x = 0.05， 同时 x 较大的时候亦可。同样，对 ND 也是如此。 另外，代码尽量优化，时间复杂度尽量小。

这样一来，类似 Taylor 展开就不好使了，因为 Taylor 收敛的快慢是对某点而言的，要类似于“一致收敛”的算法来计算，则适用面更广。

感谢 以上两位网友参与讨论，也希望别的网友能提出你们的看法。

[/quote]


找到yeg001程序的问题了，yeg001把系数单独出来，即temp是不收敛的，在x比较小的时候temp，temp还没超界，z就收敛到要求了！！！而x比较大的时候，z还没达到要求，temp已经超界了！！！！因此，系数不能单独出来，虽然程序看起来易懂一些！！！系数要和对应的级数项要合并在一起考虑！！

为了避免出现浮点异常的问题，再对23楼的程序作一点改进。

program main

 implicit none
 integer :: nd=365,i
 real(kind=4) :: x=0.06,y,z,term,esp=1.e-8
 
 y=x/nd
 term=100./y
 z=0.
 i=1
 do  while(.true.)
  term=term*nd*y/i
  z=z+term
  if (term/z.lt.esp) exit
  i=i+1
  nd=nd-1
 enddo
 write(*,*) z
 
end program

本人最近在看浮点数的相关知识、Elementary Function （比如 sin, cos, exp, log, atan）的实现（即：计算机算出数值所用的算法），可能的问题。

本人的原意是看看 log 函数，到底如何？现在看来，log(x) 在 x = 1 附近的相对误差值得人们注意。同样的，三角函数，比如 sin(x) 在 x 很大时，系统内置的函数也不管用了；这个已经有人专门成文论述如何处理了。我们进行数值计算，知道这些总是有好处的。

这是本人写出的程序：  具体的实现是： 先化成 幂指函数 -- 然后调用能保证 x = 1 附近 log(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


subroutine PowerOperator
  include "link_fnl_shared.h"
  use ALNREL_INT
  implicit none
  integer, save:: ND = 365
  real:: x, y
  real:: u
  procedure( real ):: Log_1Px
  
  x = 0.06
  u = x / real(ND, kind(x))

  ! method 0 -- estimate 
  ! cannot get 6-7 decimal significant bit
  y = 100.0 * (exp(x) - 1.0) / u
  write(unit = *, fmt = *) "method 0: ", y 
    
  ! method 1 -- direct comput, 
  ! cannot get 6-7 decimal significant bit
  y = 100.0 * ((1.0 + u) ** ND - 1.0) / u
  write(unit = *, fmt = *) "method 1: ", y 
  
  ! method 2
  ! log -- cann't control the relative error near (u) zero
  y = ND * log(1.0 + u)
  y = exp(y)
  y = 100.0 * (y - 1.0) / u
  write(unit = *, fmt = *) "method 2: ", y

  ! method 3 -- Taylor Series
  ! good for small x, but when x larger,
  ! it convergence slowly
  y = ND * (u - 0.5 * u * u + u ** 3 / 3.0 - 0.25 * u ** 4)
  y = exp(y)
  y = 100.0 * (y - 1.0) / u
  write(unit = *, fmt = *) "method 3: ", y

  ! method 4
  y = ND * Log_1Px(u)
  y = exp(y)
  y = 100.0 * (y - 1.0) / u
  write(unit = *, fmt = *) "method 4: ", y
        
  ! method 5
  y = ND * ALNREL(u)
  y = exp(y)
  y = 100.0 * (y - 1.0) / u
  write(unit = *, fmt = *) "method 0: ", y        
         
  return
end subroutine

[quote]本人最近在看浮点数的相关知识、Elementary Function （比如 sin, cos, exp, log, atan）的实现（即：计算机算出数值所用的算法），可能的问题。

本人的原意是看看 log 函数，到底如何？现在看来，log(x) 在 x = 1 附近的相对误差值得人们注意。同样的，三角函数，比如 sin(x) 在 x 很大时，系统内置的函数也不管用了；这个已经有人专门成文论述如何处理了。我们进行数值计算，知道这些总是有好处的。

这是本人写出的程序：  具体的实现是： 先化成 幂指函数 -- 然后调用能保证 x = 1 附近 log(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


subroutine PowerOperator
  include "link_fnl_shared.h"
  use ALNREL_INT
  implicit none
  integer, save:: ND = 365
  real:: x, y
  real:: u
  procedure( real ):: Log_1Px
  
  x = 0.06
  u = x / real(ND, kind(x))

  ! method 0 -- estimate 
  ! cannot get 6-7 decimal significant bit
  y = 100.0 * (exp(x) - 1.0) / u
  write(unit = *, fmt = *) "method 0: ", y 
    
  ! method 1 -- direct comput, 
  ! cannot get 6-7 decimal significant bit
  y = 100.0 * ((1.0 + u) ** ND - 1.0) / u
  write(unit = *, fmt = *) "method 1: ", y 
  
  ! method 2
  ! log -- cann't control the relative error near (u) zero
  y = ND * log(1.0 + u)
  y = exp(y)
  y = 100.0 * (y - 1.0) / u
  write(unit = *, fmt = *) "method 2: ", y

  ! method 3 -- Taylor Series
  ! good for small x, but when x larger,
  ! it convergence slowly
  y = ND * (u - 0.5 * u * u + u ** 3 / 3.0 - 0.25 * u ** 4)
  y = exp(y)
  y = 100.0 * (y - 1.0) / u
  write(unit = *, fmt = *) "method 3: ", y

  ! method 4
  y = ND * Log_1Px(u)
  y = exp(y)
  y = 100.0 * (y - 1.0) / u
  write(unit = *, fmt = *) "method 4: ", y
        
  ! method 5
  y = ND * ALNREL(u)
  y = exp(y)
  y = 100.0 * (y - 1.0) / u
  write(unit = *, fmt = *) "method 0: ", y        
         
  return
end subroutine[/quote]


楼主的方法里面似乎没有采用yeg001和我使用的方法，这个方法完全避开了exp, log的运算，因而结果似乎更精确。

1.0 + x == 1.0 
这个条件废了
浮点数不能用相等比较吧

楼主再解释一下这个语句：

 Log_1Px = x * log(1.0 + x) / ((1.0 + x) - 1.0)

这对提高精度有什么好处？

(1) for single precision, epsilon = 2 ** (-23)
(2) when abs(x) < epsilon, then 1.0 + x == 1.0, so the relative error for Log_1Px is
0.5 * epsilon
(3) when 1.0 + x /= x, then the relative error for expression 
 Log_1Px = x * log(1.0 + x) / ((1.0 + x) - 1.0)
is 5.0 * epsilon < 6.0E-7

so the relative error for function Log_1Px is 6.0E-7.
When we use fortran programming, the "internal precision" will help you to decrease rounding error. 

[quote](1) for single precision, epsilon = 2 ** (-23)
(2) when abs(x) < epsilon, then 1.0 + x == 1.0, so the relative error for Log_1Px is
0.5 * epsilon
(3) when 1.0 + x /= x, then the relative error for expression 
 Log_1Px = x * log(1.0 + x) / ((1.0 + x) - 1.0)
is 5.0 * epsilon < 6.0E-7

so the relative error for function Log_1Px is 6.0E-7.
When we use fortran programming, the "internal precision" will help you to decrease rounding error. 
[/quote]
其实这个问题John D. Cook的文章里提到过
lz的这句
 if ( 1.0 + x == 1.0 ) then
    Log_1Px = x 

根据他的说法最好改为
 if ( abs(x)<eps ) then
    Log_1Px = (-0.5*x + 1.0)*x; 

印象中是这样的，文章名好像翻译过来叫看似毫无必要的数学函数

[quote](1) for single precision, epsilon = 2 ** (-23)
(2) when abs(x) < epsilon, then 1.0 + x == 1.0, so the relative error for Log_1Px is
0.5 * epsilon
(3) when 1.0 + x /= x, then the relative error for expression 
 Log_1Px = x * log(1.0 + x) / ((1.0 + x) - 1.0)
is 5.0 * epsilon < 6.0E-7

so the relative error for function Log_1Px is 6.0E-7.
When we use fortran programming, the "internal precision" will help you to decrease rounding error. 
[/quote]

Log_1Px = x * log(1.0 + x) / ((1.0 + x) - 1.0) 不是 Log_1Px = log(1.0 + x)吗？多点运算会增加精度吗？ 有点费解！