在fortran中快速向上/向下翻倍？

Question

在 Fortran 中是否有快速向上/向下舍入的方法？

由于正双数的位表示的线性顺序，可以如下实现舍入。

pinf 和 ninf 是全局常量，分别为 +/- 无穷大

function roundup(x)
    double precision ,intent(in) :: x
    double precision :: roundup

    if (isnan(x))then 
        roundup = pinf
        return
    end if 
    if (x==pinf)then
        roundup = pinf
        return 
    end if
    if (x==ninf)then
        roundup = ninf
        return 
    end if
    if (x>0)then
        roundup = transfer((transfer(x,1_8)+1_8),1d0)
    else if (x<0) then
        roundup = transfer((transfer(x,1_8)-1_8),1d0)
    else
        if (transfer(x,1_8)==Z'0000000000000000')then
            roundup = transfer((transfer(x,1_8)+1_8),1d0)
        else
            roundup = transfer((transfer(-x,1_8)+1_8),1d0)
        end if 
    end if 
end function roundup

我觉得这不是最好的方法，因为它很慢，但它几乎只使用位操作。

另一种方法是使用乘法和一些 epsilon eps = epsilon (1d0)

function roundup2(x)
    double precision ,intent(in) :: x
    double precision :: roundup2
    if (isnan(x)) then 
        roundup2 = pinf
        return 
    else if (x>=eps) then 
        roundup2 = x*(1d0+eps)
    else if (x<=-eps) then 
        roundup2 = x*(1d0-eps)
    else 
        roundup2 = eps
    end if 
end function roundup2

对于某些x，两个函数返回相同的结果 (1d0, 158d0)，对于某些不返回 (0.1d0, 15d0)。

第一个函数更准确，但比第二个慢 3.6 倍 （10^9 轮测试中 11.1 对 3.0 秒）

    print * ,x,y,abs(x-y)
    do i = 1, 1000000000
        x = roundup(x)
        !y = roundup2(y)
    end do 
    print * ,x,y,abs(x-y)

在不检查 NaN/Infinities 的情况下，第一个功能测试需要 8.5 秒 (-20%)。

我非常努力地使用圆形功能，并且在程序配置文件中需要很多时间。是否有跨平台的方法可以更快地进行舍入而不丢失精度？

更新

该问题怀疑当时的 roundup 和 rounddown 调用无法重新排序。我没有提到四舍五入以保持话题简短。

提示： 第一个函数使用两个transfer 函数和一个加法。在第二种情况下，它比一次乘法和一次加法要慢。当它对数字的位没有任何作用时，为什么传输成本如此之高？是否可以用更快的函数代替传输或完全避免添加调用？

Answer 1

我建议您查看 Fortran 标准 IEEE 浮点内在模块（IEEE_ARITHMETIC、IEEE_FEATURES、IEEE_EXCEPTIONS）。这些提供了 IEEE_SET_ROUNDING_MODE，您可以在其中设置后续操作的舍入模式。理想情况下，您会使用 IEEE_GET_ROUNDING_MODE 获取当前模式并保存，设置新模式，执行操作，然后恢复模式。

一些注意事项 - 更改处理器舍入模式本身是一项缓慢的操作，但如果您只做一次，然后再做很多轮，那将是一个胜利。并非所有当前的 Fortran 编译器都支持 IEEE 内在模块，但大多数合理的应该都支持。您可能需要告诉编译器您正在使用 IEEE 环境 - 对于 Intel Fortran，请使用“-fp-model strict”。

Answer 2

如果我正确理解您想要做什么，如果您将 +/- 无穷大作为参数，“最近的”内在函数不是做您想要的吗？

http://gcc.gnu.org/onlinedocs/gfortran/NEAREST.html#NEAREST

如果编译器以良好的性能实现它，这可能会起作用。如果您希望 NaN 舍入为 Inf，则必须将其添加到包装器中。

至于为什么 roundup2 更快，我无法确定你的机器上发生了什么，但我可以说两件事：

1. roundup2 中的加法可能已优化（如果 eps 是参数？），所以实际上只有乘法。
2. 如果传输确实有任何作用，那么很容易显着减慢函数的速度，因为函数本身很短。如果转移只是对 x 进行了多余的复制，这甚至可能是正确的。