基本三角方程的浮点异常

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

我知道 sindcosd 是非标准函数,但它们通过适当的转换链接到目标文件。我之前测试过这个,它适用于其他使用度数而不是弧度的代码:

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。