R调用Fortran sinc子程序出现精度异常与传参问题求助

【已编辑】

问题背景 #

我正在尝试实现R与Fortran的集成，在R环境中调用Fortran子程序。我是Fortran初学者，但已熟练掌握R语言及通用编程技能。

我编写了一个用于计算sinc函数的Fortran子程序，计划在R中调用，该子程序代码如下：

SUBROUTINE SINC(X,Y)
    IMPLICIT NONE
    DOUBLE PRECISION , INTENT(IN) :: X
    DOUBLE PRECISION , PARAMETER :: PI = 4.D0*ATAN(1.0)
    DOUBLE PRECISION , INTENT(OUT) :: Y
    
    IF(X == 0) THEN
        Y = 1
    ELSE
        Y = (SIN(PI*X))/(PI*X)
    END IF
END SUBROUTINE SINC

我编写了配套的Fortran测试代码，运行状态正常，测试代码如下：

PROGRAM DEMO
DOUBLE PRECISION A, B
PRINT *, "WHAT IS THE NUMBER YOU WANT TO SINC?"
READ *, A
CALL SINC(A, B)
PRINT *, "THE SINC IS", B
END PROGRAM DEMO

SUBROUTINE SINC(X,Y)
    IMPLICIT NONE
    DOUBLE PRECISION , INTENT(IN) :: X
    DOUBLE PRECISION , PARAMETER :: PI = 4.D0*ATAN(1.0)
    DOUBLE PRECISION , INTENT(OUT) :: Y
    
    IF(X == 0) THEN
        PRINT *, "X IS NULL"
        Y = 1
    ELSE
        PRINT *, "X IS NON NULL"
        Y = (SIN(PI*X))/(PI*X)
    END IF
END SUBROUTINE SINC

我使用如下命令将上述子程序（不含主程序部分）编译为.dll文件：

R CMD SHLIB sinc.f95

此前我若使用REAL类型而非DOUBLE PRECISION定义变量，运行时传入的Y参数会被原样返回，故障原因不明，对应的R端调用代码及返回结果如下：

dyn.load("D:\\Fortran\\sinc.dll")
.Fortran("sinc", X=as.double(0), Y=as.double(1))

$X
[1] 0

$Y
[1] 1

当前使用DOUBLE PRECISION定义变量时，子程序对1.1、2.5这类非整数值计算完全正常，但计算1、2、3、4等整数输入值时，返回结果为-2.7827534378485793E-008，而R原生计算同场景结果为3.898043e-17，二者虽理论上均趋近于0，但精度存在明显差异。

我使用gfortran作为编译器，需要定位上述问题的故障原因，找到能在R中正常调用该Fortran子程序的修正方案。

故障根因 #

两个异常的本质都是Fortran侧的类型/精度定义和R接口的传递规则不匹配：

• 单精度REAL类型返回值错乱：R的.Fortran()接口默认传递双精度（8字节）浮点数值，Fortran侧使用单精度REAL（通常为4字节）定义参数时，内存读写长度不匹配，会直接导致返回值异常，即观察到的Y参数被原样返回的现象。
• 整数输入精度偏差：PI的定义存在精度损失。4.D0*ATAN(1.0)中传入ATAN的参数1.0是单精度浮点数，函数返回值也为单精度，即便乘以双精度常量4.D0，最终得到的PI仅保留约7位有效数字，远低于双精度15-17位有效数字的要求。当输入为整数时，理论上sin(π*X)=0，但因为PI精度不足，计算得到的π*X与真实的π整数倍存在明显偏差，最终得到的sin值误差在1e-8量级，和R原生双精度PI计算出的1e-17量级误差形成明显差距。

修正方案 #

1. 所有和R交互的Fortran数值变量，统一使用DOUBLE PRECISION双精度类型定义，不要使用默认单精度的REAL类型。
2. 修正PI的定义，保证计算PI时所有常量均为双精度，将PI定义行修改为DOUBLE PRECISION, PARAMETER :: PI = 4.D0*ATAN(1.D0)，把ATAN的参数从单精度的1.0换成双精度的1.D0，得到完整双精度精度的PI值。
3. 优化边界判断逻辑：浮点数计算存在固有精度误差，不要直接用X == 0做精确相等判断，改为判断X的绝对值是否小于极小阈值，避免边界值判断失效。

修正后的Fortran子程序参考：

SUBROUTINE SINC(X,Y)
    IMPLICIT NONE
    DOUBLE PRECISION , INTENT(IN) :: X
    DOUBLE PRECISION , PARAMETER :: PI = 4.D0*ATAN(1.D0)
    DOUBLE PRECISION , INTENT(OUT) :: Y
    
    IF(ABS(X) < 1D-15) THEN
        Y = 1.0D0
    ELSE
        Y = SIN(PI*X)/(PI*X)
    END IF
END SUBROUTINE SINC

修改完成后重新执行R CMD SHLIB sinc.f95编译，再在R中调用即可得到和原生计算精度一致的结果。

内容的提问来源于stack exchange，提问作者YetAnotherUsr