c - 定点乘法算法

Question

我正在尝试将时间戳(仅秒的小数部分)从纳秒(单位为 10^-9 秒)重新调整为 NTP 时间戳的下半部分(单位为 2^-32 秒)。实际上，这意味着乘以 4.2949673。但我需要在没有 float 学的情况下完成它，并且不使用大于 32 位的整数(事实上，我实际上是为 8 位微 Controller 编写这个，所以即使是 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

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;
}