基本三角方程的浮点异常
Floating point exceptions for basic trigonometric equations
我在从地理坐标转换为地心坐标的子例程中遇到浮点异常问题。变量 geo(x)
作为纬度和经度对输入到子例程中。变量 xyz(x)
作为分量的三元组输出(x
-- 格林威治子午线和赤道;y
-- 90 度经度和赤道;z
-- 北极) .
subroutine geo2xyz(geo,xyz)
IMPLICIT REAL*8 (A-H,O-Z)
dimension geo(2),xyz(3)
rr=6367443.5
xyz(1)=rr*sind(90.-geo(1))*cosd(geo(2))
xyz(2)=rr*sind(90.-geo(1))*sind(geo(2))
xyz(3)=rr*cosd(90.-geo(1))
return
end
我知道 sind
和 cosd
是非标准函数,但它们通过适当的转换链接到目标文件。我之前测试过这个,它适用于其他使用度数而不是弧度的代码:
real function sind(x)
IMPLICIT REAL*8 (A-H,O-Z)
sind=sin(x*3.141592653589793d0/180.0d0)
return
end
real function cosd(x)
IMPLICIT REAL*8 (A-H,O-Z)
cosd=cos(x*3.141592653589793d0/180.0d0)
return
end
我试着用 GDB 遍历程序,但没能找出问题所在。好像方程里所有的变量都可以,但是xyz(1) = 1
,如果我用计算器计算就不是这样了。
Program received signal SIGFPE, Arithmetic exception.
0x0000000100001cb8 in geo2xyz (geo=..., xyz=...) at disslip.f:758
758 xyz(1)=rr*sind(90.-geo(1))*cosd(geo(2))
(gdb) print xyz(1)
= 1
(gdb) print rr
= 6367443.5
(gdb) print geo(1)
= 50.350000000000001
(gdb) print geo(2)
= -127.55800000000001
(gdb)
我在这里遗漏了什么导致浮点异常?我敢肯定这很简单,我对此很陌生。
sind 和 cosd 函数执行 return REAL
(单精度)浮点。
但是subroutine geo2xyz
对此一无所知。
使用 IMPLICIT REAL*8
它将假定这些函数 return 是双精度浮点数。
因此,在使用sind和cosd之前必须正确声明为real,and/or将它们的return类型改为real*8。
有趣的是,该程序可能适用于旧的 x86 架构,因为单精度值 returned 在 ST0 寄存器中被提升为双精度。它不适用于使用 SSE2 和寄存器 XMM0 的 64 位英特尔:在这种架构中没有提升到 double。
我在从地理坐标转换为地心坐标的子例程中遇到浮点异常问题。变量 geo(x)
作为纬度和经度对输入到子例程中。变量 xyz(x)
作为分量的三元组输出(x
-- 格林威治子午线和赤道;y
-- 90 度经度和赤道;z
-- 北极) .
subroutine geo2xyz(geo,xyz)
IMPLICIT REAL*8 (A-H,O-Z)
dimension geo(2),xyz(3)
rr=6367443.5
xyz(1)=rr*sind(90.-geo(1))*cosd(geo(2))
xyz(2)=rr*sind(90.-geo(1))*sind(geo(2))
xyz(3)=rr*cosd(90.-geo(1))
return
end
我知道 sind
和 cosd
是非标准函数,但它们通过适当的转换链接到目标文件。我之前测试过这个,它适用于其他使用度数而不是弧度的代码:
real function sind(x)
IMPLICIT REAL*8 (A-H,O-Z)
sind=sin(x*3.141592653589793d0/180.0d0)
return
end
real function cosd(x)
IMPLICIT REAL*8 (A-H,O-Z)
cosd=cos(x*3.141592653589793d0/180.0d0)
return
end
我试着用 GDB 遍历程序,但没能找出问题所在。好像方程里所有的变量都可以,但是xyz(1) = 1
,如果我用计算器计算就不是这样了。
Program received signal SIGFPE, Arithmetic exception.
0x0000000100001cb8 in geo2xyz (geo=..., xyz=...) at disslip.f:758
758 xyz(1)=rr*sind(90.-geo(1))*cosd(geo(2))
(gdb) print xyz(1)
= 1
(gdb) print rr
= 6367443.5
(gdb) print geo(1)
= 50.350000000000001
(gdb) print geo(2)
= -127.55800000000001
(gdb)
我在这里遗漏了什么导致浮点异常?我敢肯定这很简单,我对此很陌生。
sind 和 cosd 函数执行 return REAL
(单精度)浮点。
但是subroutine geo2xyz
对此一无所知。
使用 IMPLICIT REAL*8
它将假定这些函数 return 是双精度浮点数。
因此,在使用sind和cosd之前必须正确声明为real,and/or将它们的return类型改为real*8。
有趣的是,该程序可能适用于旧的 x86 架构,因为单精度值 returned 在 ST0 寄存器中被提升为双精度。它不适用于使用 SSE2 和寄存器 XMM0 的 64 位英特尔:在这种架构中没有提升到 double。