主题:[讨论]基本功训练 04】浮点数计算过程中的误差
asymptotic
[专家分:16630] 发布于 2010-10-13 09:10:00
用 kind = 4 的 real、integer 计算如下表达式:
100 * ( (1 + x / ND) ** ND - 1 ) / ( x / ND )
其中,ND = 365, x = 0.06; 精确答案 -- 37614.05。
编写程序,精确到有效数字第六位,第七位可以有区别。
提示:这个问题主要是设计算法,不能考虑更高精度的浮点数、整数,
以减少浮点数计算过程中的误差(当然,若是这样,则要求至少 14 位
有效数字正确)。
回复列表 (共57个回复)
51 楼
yeg001 [专家分:14390] 发布于 2010-10-15 10:08:00
[quote][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]
我是做二项式展开, 不是级数展开. 项是有限但没必要一直算下去, 达到精度就退出. 结果输出的4表示循环进行4次已经达到精度.
52 楼
jstzhurj [专家分:4680] 发布于 2010-10-15 10:20:00
[quote][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 也可以工作,求验证,先保证精度再考虑时间复杂度[/quote]
是,能工作,x=0.000001,n=365000000还能正常工作吗?
53 楼
yrliu [专家分:750] 发布于 2010-10-15 10:41:00
[quote][quote][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 也可以工作,求验证,先保证精度再考虑时间复杂度[/quote]
是,能工作,x=0.000001,n=365000000还能正常工作吗?[/quote]
那是不行了。那又要精确又要时间复杂度小,只好不在一个算法上吊死,程序长点,多用几种算法呗
54 楼
jstzhurj [专家分:4680] 发布于 2010-10-15 11:26:00
[quote][quote][quote][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 也可以工作,求验证,先保证精度再考虑时间复杂度[/quote]
是,能工作,x=0.000001,n=365000000还能正常工作吗?[/quote]
那是不行了。那又要精确又要时间复杂度小,只好不在一个算法上吊死,程序长点,多用几种算法呗[/quote]
1)有必要二重循环吗?做了很多重复计算。
2)z小于一定值退出是非常必要的,不设退出条件,花了很多时间做了很多无用功,精度却没有提高。
从这两点考虑优化即可。
55 楼
yeg001 [专家分:14390] 发布于 2010-10-15 12:38:00
拜读了asymptotic的回帖和dongyuanxun贴的链接的文章. 觉得有时间应该看看一些基本函数在极端情况下的处理. 不知道两位有没有这类基本函数的精度分析的文章呢? 最好是pdf的.
前段时间写了一个在二维空间中选取多边形区域的点的程序. 因为知道夹角当时我就用了一个旋转矩阵的方法得到另一条限制直线. 但后续的测试中发现这个直线方程并不准确. 于是用别的方法构造直线就能很好满足精度要求了. 所以我的理解是误差是由旋转矩阵中cos, sin函数结果精度不够引起的. 时间比较紧, 实现了功能之后就没再细究这个问题了. 这两天偶尔上来看到这个问题, 经大家讨论, 我觉得或许应该对基本函数的实现和精度有一个更好的把握.
56 楼
yrliu [专家分:750] 发布于 2010-10-19 15:20:00
[quote
1)有必要二重循环吗?做了很多重复计算。
2)z小于一定值退出是非常必要的,不设退出条件,花了很多时间做了很多无用功,精度却没有提高。
从这两点考虑优化即可。[/quote]
嗯,确实做了很多无用功,现在优化成下面的程序可以了吧?
program main
implicit none
real(kind=4):: a=100,b=0.06,n=3650000
real(kind=4):: nn,i,x,tmp,res=0.0,saveRes=0.0
tmp=n; nn=n+1.0
x=b/n
do i=2,n
tmp=tmp*(nn/i-1.0)*x
res=res+tmp
if (saveRes .eq. res) then
print*,i,tmp,res
exit
endif
saveRes=res
enddo
print *,a*(res+n)
end
57 楼
jstzhurj [专家分:4680] 发布于 2010-10-19 16:38:00
这回相当不错。
我来回复