在浮子上施加（某些情况下）周期性边界条件的有效方法？

Posted 2021-05-04

通过简单地执行以下操作，可以非常有效地对整数强加一些periodic boundary conditions (PBC)的情况：

myWrappedWithinPeriodicBoundary = myUIntValue & mask

当边界是半开放范围[0，upperBound]时，这是有效的，其中（不包括）upperBound是2 ^ exp，这样

mask = (1 << exp) - 1

例如：

let pbcUpperBoundExp = 2 // so the periodic boundary will be [0, 4)
let mask = (1 << pbcUpperBoundExp) - 1
for x in -7 ... 7 { print(x & mask, terminator: " ") }

（在Swift中）将打印：

1 2 3 0 1 2 3 0 1 2 3 0 1 2 3

问题：是否有任何（大致相似的）有效方法在浮点数（32或64位IEEE-754）上施加（某些情况下）PBC？

答案

有几种合理的方法：

1. fmod(x,1)
2. modf(x,&dummy) - 具有静态知道其除数的优势，但在我的测试中来自libc.so.6甚至与-ffast-math
3. x-floor(x)（Jens在评论中建议） - 直接支持负面输入
4. 手动bit-twiddling直接实现
5. 手动bit-twiddling实现floor

前两个保留了他们输入的标志;如果它是负数，你可以加1。

两位操作非常相似：您确定哪些有效位对应于整数部分，并将它们（用于直接实现）或其余（用于实现floor）屏蔽掉。直接实现可以通过浮点除法完成，也可以通过手动重新组装double来完成;前者即使硬件CLZ快28％。 floor实现可以立即重建double：floor永远不会更改其参数的指数，除非它返回0.大约需要20行C。

以下时间是使用double和gcc -O3，时序循环超过代表性输入，其中内联操作代码。

fmod: 41.8 ns
modf: 19.6 ns
floor: 10.6 ns

使用-ffast-math：

fmod: 26.2 ns
modf: 30.0 ns
floor: 21.9 ns

位操作：

direct: 18.0 ns
floor: 20.6 ns

手动实施具有竞争力，但floor技术是最好的。奇怪的是，三个库函数中的两个在没有-ffast-math的情况下表现更好：即作为PLT函数调用而不是内联内置函数。

另一答案

我将这个答案添加到我自己的问题中，因为它在撰写本文时描述了我找到的最佳解决方案。它在Swift 4.1中（应该直接转换为C）并且已经在各种用例中进行了测试：

extension BinaryFloatingPoint {
    /// Returns the value after restricting it to the periodic boundary
    /// condition [0, 1).
    /// See https://forums.swift.org/t/why-no-fraction-in-floatingpoint/10337
    @_transparent
    func wrappedToUnitRange() -> Self {
        let fract = self - self.rounded(.down)
        // Have to clamp to just below 1 because very small negative values
        // will otherwise return an out of range result of 1.0.
        // Turns out this:
        if fract >= 1.0 { return Self(1).nextDown } else { return fract }
        // is faster than this:
        //return min(fract, Self(1).nextDown)
    }
    @_transparent
    func wrapped(to range: Range<Self>) -> Self {
        let measure = range.upperBound - range.lowerBound
        let recipMeasure = Self(1) / measure
        let scaled = (self - range.lowerBound) * recipMeasure
        return scaled.wrappedToUnitRange() * measure + range.lowerBound
    }
    @_transparent
    func wrappedIteratively(to range: Range<Self>) -> Self {
        var v = self
        let measure = range.upperBound - range.lowerBound
        while v >= range.upperBound { v = v - measure }
        while v < range.lowerBound { v = v + measure }
        return v
    }
}

在我的带有2 GHz Intel Core i7的MacBook Pro上，随机（有限）Double值上的一千万（可能是内联）调用wrapped(to range:)需要0.6秒，这大约是每秒1.66亿次调用（不是多线程）。静态知道或不知道的范围，或具有2的幂等的界限或量度，可以产生一些差异，但不会像人们想象的那样多。

wrappedToUnitRange()需要大约0.2秒，这意味着我的系统每秒有5亿次呼叫。

鉴于正确的情况，wrappedIteratively(to range:)和wrappedToUnitRange()一样快。

通过比较基线测试（没有包装一些值，但仍然使用它来计算例如一个简单的xor校验和）到一个值被包装的相同测试来进行计时。这些之间的时间差是我给包装调用的时间。

我使用了Swift开发工具链2018-02-21，使用-O -whole-module-optimization -static-stdlib -gnone进行编译。并且已经注意使测试相关，即使用不同分布的真实随机输入来防止死代码移除等。一般地编写包装函数，如BinaryFloatingPoint上的这个扩展，结果被优化为等效代码，就好像我有为浮动和双重编写单独的专用版本。

有兴趣看到比我更熟练的人进一步调查（C或Swift或任何其他语言并不重要）。

编辑：对于任何感兴趣的人，这里是simd float2的一些版本：

extension float2 {
    @_transparent
    func wrappedInUnitRange() -> float2 {
        return simd.fract(self)
    }
    @_transparent
    func wrappedToMinusOneToOne() -> float2 {
        let scaled = (self + float2(1, 1)) * float2(0.5, 0.5)
        let scaledFract = scaled - floor(scaled)
        let wrapped = simd_muladd(scaledFract, float2(2, 2), float2(-1, -1))
        // Note that we have to make sure the result is not out of bounds, like
        // simd fract does:
        let oneNextDown = Float(bitPattern:
            0b0_01111110_11111111111111111111111)
        let oneNextDownFloat2 = float2(oneNextDown, oneNextDown)
        return simd.min(wrapped, oneNextDownFloat2)
    }
    @_transparent
    func wrapped(toLowerBound lowerBound: float2,
                 upperBound: float2) -> float2
    {
        let measure = upperBound - lowerBound
        let recipMeasure = simd_precise_recip(measure)
        let scaled = (self - lowerBound) * recipMeasure
        let scaledFract = scaled - floor(scaled)
        // Note that we have to make sure the result is not out of bounds, like
        // simd fract does:
        let wrapped = simd_muladd(scaledFract, measure, lowerBound)
        let maxX = upperBound.x.nextDown // For some reason, this won't be
        let maxY = upperBound.y.nextDown // optimized even when upperBound is
        // statically known, and there is no similar simd function available.
        let maxValue = float2(maxX, maxY)
        return simd.min(wrapped, maxValue)
    }
}

我问了一些与simd相关的问题here可能会引起关注。

Aaditi：

从上面的Swift论坛主题可以看出：

// Note that tiny negative values like:
let x: Float = -1e-08
// May produce results outside the [0, 1) range:
let wrapped = x - floor(x)
print(wrapped < 1.0) // false
// which may result in out-of-bounds table accesses
// in common usage, so it's probably better to use:
let correctlyWrapped = simd_fract(x)
print(correctlyWrapped < 1.0) // true

我已经更新了代码来解释这个问题。