本帖最后由 ocw^。 于 2020-11-14 21:30 编辑
编写一个进行50位有效数字进行计算的程序,出来很多bug,但不知道这些bug的原因和如何修改,求助
[Fortran] 纯文本查看 复制代码 module MyPrec
implicit none
integer, parameter :: rp = selected_real_kind(p = 50)!保存50位浮点数
end module
module MyHmlt
use MyPrec
implicit none
real(rp) :: t1, t2, p0, er1, er2, er3, er4!两个精度为50的浮点数
integer :: n0 !有效数字
end module
program play
use MyPrec
use MyHmlt
implicit none
character :: a
write(*,*)"输入运算符号"
read(*,*) a
t1 = 1.454578787_rp
t2 = 5.9897854685_rp
if(a == '+') then
p0 = 0._rp
p0 = t1 + t2
call jiaxiangduiwucha(t1, t2)
call youxiaoshuzi(p0, er1)
call jiafa(p0)
end if
end program
!subroutine guiling
! use MyPrec
! use MyHmlt
! implicit none
! t1 = 0._rp
! t2 = 0._rp
! er1 = 0._rp
! er2 = 0._rp
! er3 = 0._rp
! er4 = 0._rp
! p0 = 0._rp
! rp2 = 0
! n0 = 0
!end subroutine
subroutine jiafa
use MyPrec
use MyHmlt
implicit none
real(rp) :: x1
integer :: i, m, n
character(0:99) :: string, jieguo
write(x1,"(E100.50)") p0
write(string,"(I100)") x1
m = 0
n = 0
do i = 0, 99, 1
if(n == n0) then
exit
end if
if(string(i) /= ' ') then
jieguo(m) = string(i)
m = m + 1
if(string(i) == '.') then
cycle
end if
n = n + 1
end if
end do
open(30, file = "jiafa.txt", access = "append")
write(30, *) jieguo
close(30)
end subroutine
subroutine jiaxiangduiwucha
use MyPrec
use MyHmlt
implicit none
integer :: m1, m2
m1 = floor(log10(t1))
m2 = floor(log10(t2))
er1 = (0.5*10**(m1+1-50) + 0.5*10**(m2+1-50))/(abs(t1+t2))
write(*,*) "t1+t2的相对误差为 ", er1
end subroutine
subroutine jianxiangduiwucha
use MyPrec
use MyHmlt
implicit none
integer :: m1, m2
m1 = floor(log10(t1))
m2 = floor(log10(t2))
er2 = (0.5*10**(m1+1-50) + 0.5*10**(m2+1-50))/(abs(t1-t2))
write(*,*) "t1-t2的相对误差为 ", er2
end subroutine
subroutine chengxiangduiwucha
use MyPrec
use MyHmlt
implicit none
integer :: a1, a2
a1 = FLOOR(abs(t1*10**(-LOG10(t1))))
a2 = FLOOR(abs(t1*10**(-LOG10(t2))))
er3 = 10**(-49)/(2*a1) + 10**(-49)/(2*a2)
write(*,*) "t1*t2的相对误差为 ", er3
end subroutine
subroutine chuxiangduiwucha
use MyPrec
use MyHmlt
implicit none
integer :: a1, a2
a1 = FLOOR(abs(t1*10**(-LOG10(t1))))
a2 = FLOOR(abs(t2*10**(-LOG10(t2))))
er4 = abs(10**(-49)/(2*a1) - 10**(-49)/(2*a2))
write(*,*) "t1/t2的相对误差为 ", er4
end subroutine
subroutine youxiaoshuzi(er0)
use MyPrec
use MyHmlt
implicit none
real(rp) :: er0
integer :: a1
a1 = FLOOR(LOG10(p0))
a1 = FLOOR(abs(p0)*10**(-a1))
n0 = -FLOOR(LOG10(2*er0*(a1+1)))
write(*,*) "有效数字为", n0
end subroutine
|