for abs(x) < eps,  -0.5 * x + 1.0 == 1.0 for single precision without "internal representation".
I cannot find John D. Cook's article, give me a hand.

[quote][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)吗？多点运算会增加精度吗？ 有点费解！[/quote]
数值计算和数学的化简是不一样的
这样会相对误差小的

[quote][quote][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)吗？多点运算会增加精度吗？ 有点费解！[/quote]
数值计算和数学的化简是不一样的
这样会相对误差小的[/quote]


能举个实例吗，想亲眼看看到底怎么减少相对误差的。

[quote]for abs(x) < eps,  -0.5 * x + 1.0 == 1.0 for single precision without "internal representation".
I cannot find John D. Cook's article, give me a hand.[/quote]

我也找了半天
应该是这篇吧
http://www.johndcook.com/blog/2010/06/07/math-library-functions-that-seem-unnecessary/

汗，你的编译器是gfortran吧，我刚用gfortran测试了下，根你结果一致。
开始的结果我用的是CVF

下面是gfortran版本，明天去测试下gfortran4.5.1：
liuying@liuying-laptop:~$ gfortran -v
Using built-in specs.
Target: i486-linux-gnu
Configured with: ../src/configure -v --with-pkgversion='Ubuntu 4.4.3-4ubuntu5' --with-bugurl=file:///usr/share/doc/gcc-4.4/README.Bugs --enable-languages=c,c++,fortran,objc,obj-c++ --prefix=/usr --enable-shared --enable-multiarch --enable-linker-build-id --with-system-zlib --libexecdir=/usr/lib --without-included-gettext --enable-threads=posix --with-gxx-include-dir=/usr/include/c++/4.4 --program-suffix=-4.4 --enable-nls --enable-clocale=gnu --enable-libstdcxx-debug --enable-plugin --enable-objc-gc --enable-targets=all --disable-werror --with-arch-32=i486 --with-tune=generic --enable-checking=release --build=i486-linux-gnu --host=i486-linux-gnu --target=i486-linux-gnu
Thread model: posix
gcc version 4.4.3 (Ubuntu 4.4.3-4ubuntu5)

[quote]
汗，你的编译器是gfortran吧，我刚用gfortran测试了下，根你结果一致。
开始的结果我用的是CVF

下面是gfortran版本，明天去测试下gfortran4.5.1：
liuying@liuying-laptop:~$ gfortran -v
Using built-in specs.
Target: i486-linux-gnu
Configured with: ../src/configure -v --with-pkgversion='Ubuntu 4.4.3-4ubuntu5' --with-bugurl=file:///usr/share/doc/gcc-4.4/README.Bugs --enable-languages=c,c++,fortran,objc,obj-c++ --prefix=/usr --enable-shared --enable-multiarch --enable-linker-build-id --with-system-zlib --libexecdir=/usr/lib --without-included-gettext --enable-threads=posix --with-gxx-include-dir=/usr/include/c++/4.4 --program-suffix=-4.4 --enable-nls --enable-clocale=gnu --enable-libstdcxx-debug --enable-plugin --enable-objc-gc --enable-targets=all --disable-werror --with-arch-32=i486 --with-tune=generic --enable-checking=release --build=i486-linux-gnu --host=i486-linux-gnu --target=i486-linux-gnu
Thread model: posix
gcc version 4.4.3 (Ubuntu 4.4.3-4ubuntu5)[/quote]

你说的是哪段程序？得到什么结果？

program main
 implicit none
 integer :: n=365,i,j
 real(kind=4) :: x=0.06,y,z

y=0.0 
do i=1,n,1
   z=1.0
   do j=1,i,1
      z=z*(n-j+1)*1.0/(n*(i-j+1))*x
   enddo
   y=y+z
! print*,z
enddo
y=(y/x)*n*100
print*,y
end program

g95 编译，结果是 37614.047
请检验

[quote]既然jstzhurj兄那么有热情, 改了那么多个程序, 我也来再献丑抛砖引玉吧.
我直接对 (1+x/nd)**nd 作多项式展开. 显然这样做也是基于x/nd是一个小值. 这样的话多项式只需要求前面几项就能达到楼主的要求.

program main
  implicit none
  integer, parameter :: SP=4
  integer(kind=4) :: nd=365, i
  real(kind=SP) :: x=0.06, y, z, temp

  y = 1.0_sp + x
  i = 1
  temp = real(nd, kind=sp)
  DO
    temp = temp*real((nd-i), kind=sp)/real(i+1, kind=sp)
    z = temp*((x/real(nd, kind=sp))**(i+1))
    y = y + z
    i = i +1
    IF((z)<1e-8) EXIT
  END DO
  write(*, *) i-1 
  write(*, *) y
  write(*, *) 100.*(y-1.0_sp)/(x/real(nd, kind=sp))

end program


运行结果
           4
   1.061831
   37614.00
Press any key to continue[/quote]

你这还是级数展开吧？多项式展开是有限项，不用设z小于多少才退出，见俺的程序

[quote]program main
 implicit none
 integer :: n=365,i,j
 real(kind=4) :: x=0.06,y,z

y=0.0 
do i=1,n,1
   z=1.0
   do j=1,i,1
      z=z*(n-j+1)*1.0/(n*(i-j+1))*x
   enddo
   y=y+z
! print*,z
enddo
y=(y/x)*n*100
print*,y
end program

g95 编译，结果是 37614.047
请检验[/quote]

程序需要优化！
1）有必要二重循环吗？做了很多重复计算。
2）z小于一定值退出是非常必要的，不设退出条件，花了很多时间做了很多无用功，精度却没有提高。

总之，不能满足30楼的要求：时间复杂度尽量小。

[quote][quote]program main
 implicit none
 integer :: n=365,i,j
 real(kind=4) :: x=0.06,y,z

y=0.0 
do i=1,n,1
   z=1.0
   do j=1,i,1
      z=z*(n-j+1)*1.0/(n*(i-j+1))*x
   enddo
   y=y+z
! print*,z
enddo
y=(y/x)*n*100
print*,y
end program

g95 编译，结果是 37614.047
请检验[/quote]

程序需要优化！
1）有必要二重循环吗？做了很多重复计算。
2）z小于一定值退出是非常必要的，不设退出条件，花了很多时间做了很多无用功，精度却没有提高。

总之，不能满足30楼的要求：时间复杂度尽量小。
[/quote]

思路是二项式定理展开，不知道怎么优化，貌似 x=9.06 也可以工作，求验证，先保证精度再考虑时间复杂度