基本三角方程的浮点异常

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 = 1
(gdb) print rr
$2 = 6367443.5
(gdb) print geo(1)
$3 = 50.350000000000001
(gdb) print geo(2)
$4 = -127.55800000000001
(gdb)

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

Answer 1

sind 和 cosd 函数确实返回 REAL（单精度）浮点数。

但是subroutine geo2xyz对此一无所知。
对于IMPLICIT REAL*8，它将假定这些函数返回一个双精度浮点数。

因此，您必须在使用前将 sind 和 cosd 正确声明为 real，和/或将它们的返回类型更改为 real*8。

有趣的是，这个程序可能在旧的 x86 架构上运行，因为返回的单精度值在 ST0 寄存器中被提升为双精度值。它不适用于使用 SSE2 并注册 XMM0 的 64 位英特尔：在这种架构中没有提升双倍。