基本三角方程的浮点异常

Question

我在从地理坐标转换为地心坐标的子例程中遇到浮点异常问题。变量 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)

我在这里遗漏了什么导致浮点异常？我敢肯定这很简单，我对此很陌生。

Answer 1

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。

基本三角方程的浮点异常

Floating point exceptions for basic trigonometric equations

floating-point

fortran

gdb

sigfpe