定点乘法算法

Question

我正在尝试将时间戳（仅秒的小数部分）从纳秒（单位为 10^-9 秒）重新调整为 NTP 时间戳的下半部分（单位为 2^-32 秒）。实际上，这意味着乘以 4.2949673。但是我需要在没有浮点数学的情况下进行，并且不使用大于 32 位的整数（事实上，我实际上是为 8 位微控制器编写的，所以即使是 32 位数学也很昂贵，尤其是除法）。

我提出了一些运行良好的算法，但我对数值方法没有任何真正的基础，所以我很感激任何关于如何改进它们的建议，或者任何其他可能的算法更准确和/或更快。

算法 1

uint32_t intts = (ns >> 16) * 281474 + (ns << 16) / 15259 + ns / 67078;

选择前两个常数是为了稍微低于而不是超过正确的数字，并且根据经验确定最终因子 67078 来纠正这个问题。在正确值的 +/- 4 NTP 单位内产生结果，即 +/- 1 ns - 可以接受，但残差会随 ns 变化。我想我可以添加另一个术语。

算法 2

uint32_t ns2 = (2 * ns) + 1;
uint32_t intts = (ns2 << 1)
  + (ns2 >> 3) + (ns2 >> 6) + (ns2 >> 8) + (ns2 >> 9) + (ns2 >> 10)
  + (ns2 >> 16) + (ns2 >> 18) + (ns2 >> 19) + (ns2 >> 20) + (ns2 >> 21)
  + (ns2 >> 22) + (ns2 >> 24) + (ns2 >> 30) + 3;

基于 4.2949673 的二进制展开（实际上是基于 2.14748365 的二进制展开，因为我是先加一加一来完成四舍五入的）。可能比算法 1 更快（我还没有完成基准测试）。 +3 是根据经验确定的，以消除截断所有这些低位的下冲，但它并没有做最好的工作。

Answer 1

我可能会说显而易见的……但是您是否在 interwebz 上搜索过定点数学库？他们有很多。这是 Flipcode 档案中的一个很好的 C++ 和 x86 实现：

http://www.flipcode.com/archives/Fixed_Point_Routines.shtml

Answer 2

uint32_t convert(uint32_t x) {
    const uint32_t chi = 0x4b82;
    const uint32_t clo = 0xfa09;
    const uint32_t round = 0x9525;
    const uint32_t xhi = x >> 16;
    const uint32_t xlo = x & 0xffff;
    uint32_t lowTerm = xlo*clo;
    uint32_t crossTerms = xhi*clo + xlo*chi;
    uint32_t rounded = crossTerms + (lowTerm >> 16) + round >> 16;
    uint32_t highTerm = xhi*chi;
    return (x << 2) + highTerm + rounded;
}

基本定点乘法，使用四个 16x16 -> 32 乘积模拟 32x32 -> 64 乘积。选择常量round 是为了使用简单的二分搜索来最小化错误。此表达式适用于整个有效范围内的 +/-0.6 NTP。

比例因子中的前导 4 在班次中处理。编译器通常可以为这类事情生成相当不错的代码，但如果需要，通常可以通过手写汇编对其进行简化。

如果你不需要这么高的精度，你可以去掉lowTerm和round，得到一个对 +/-1.15 NTP 有好处的答案：

uint32_t convert(uint32_t x) {
    const uint32_t chi = 0x4b82;
    const uint32_t clo = 0xfa09;
    const uint32_t xhi = x >> 16;
    const uint32_t xlo = x & 0xffff;
    uint32_t crossTerms = xhi*clo + xlo*chi;
    uint32_t highTerm = xhi*chi;
    return (x << 2) + highTerm + (crossTerms >> 16) + 1;
}