在 C 中，如何将浮点数或双精度数除以 2 的 i 次幂？

Posted 2023-04-14

【中文标题】在 C 中，如何将浮点数或双精度数除以 2 的 i 次幂？

【英文标题】：In C, how to divide a float or double number by 2 raised to the power i?

【发布时间】：2020-08-16 19:10:00

【问题描述】：

我用 C 语言为 Atmel 的微控制器 SAM E70 编写了代码，它处理 32 位宽的整数值。为了进一步计算，我将整数值标准化为 0...1.0，如下所示：

#define DIV4294967296 ((double) 1.0) / ((double) 4294967296.0)
.
.
double doubleValue;
doubleValue = ((double) intValue) * DIV4294967296;

我知道我可以从doubleValue 的指数中减去 32，从而避免更昂贵的乘法。我知道ldexp() 允许将指数乘以 2 的幂 i 但我找不到任何可以让我显式读取、操作和写回 a 的指数的东西双倍的。执行所有这些步骤实际上可能并不比执行乘法更快，因此从指数中直接减去 32 是理想的。这通常如何在 C 中完成？更重要的是，ARM 的 Cortex V7 指令集如何做到最好？

附录：回答 Eric 的问题，这是 Atmel Studio 7 向我展示的反汇编代码，用于使用 ldexp、scalbn 以及与 0x1p-32 的乘法：

uint32_t intV = 123456;
 ldr    r3, [pc, #424]
 str    r3, [r7, #28]
double doubleV0 = ((double) intV) * DIV4096;
 ldr    r3, [r7, #36]        
 vmov   s15, r3      
 vcvt.f64.u32   d7, s15      
 vldr   d6, [pc, #272]       
 vmul.f64   d7, d7, d6       
 vstr   d7, [r7, #24]       
double doubleV1 = ldexp(intV, -32);
 ldr    r3, [r7, #28]
 vmov   s15, r3
 vcvt.f64.u32   d7, s15
 mvn    r0, #31
 vmov.f64   d0, d7 
 ldr    r3, [pc, #408]
 blx    r3
 vstr   d0, [r7, #16]
double doubleV2 = scalbn(intV, -32);
 ldr    r3, [r7, #28]
 vmov   s15, r3
 vcvt.f64.u32   d7, s15
 mvn    r0, #31 
 vmov.f64   d0, d7
 ldr    r3, [pc, #384]
 blx    r3
 vstr   d0, [r7, #8]
double doubleV3 = intV * 0x1p-32;
 ldr    r3, [r7, #28]
 vmov   s15, r3
 vcvt.f64.u32   d7, s15
 vldr   d6, [pc, #164]
 vmul.f64   d7, d7, d6
 vstr   d7, [r7]

看起来这些都不匹配任何 ARM 指令（就像 C 函数 fabs() 直接编译为汇编指令 vabs）。 ldexp 和 scalbn 的编码方式相同。与0x1p-32 的乘法的编码方式与我最初提出的问题的乘法相同。

附录 2： 根据 chqrlie 的建议显示它编译成的代码：

double doubleV4 = ((double) intV);
 vstr   d7, [r7]    
*(uint64_t *)&doubleV4 -= 32ULL << 52;
 mov    r3, r7       
 ldrd   r2, r3, [r3]         
 mov    r1, r7       
 adds   r4, r2, #0       
 adc    r5, r3, #4261412864      
 strd   r4, r5, [r1]

在我看来这是最便宜的实现。

最终判决：我喜欢 chqrlie 的回答，因为它可能对我们当中乘法太慢的人有用。不过，在我的情况下，我运行了一个基于中断的例程并测量了我的初始代码和 chqrlie 的替代代码的执行时间，如果最佳优化 (-O3) 与 GCC 9.3.1 一起使用，它们的运行时间完全相同。

【问题讨论】：

您的编译器为ldexp(intValue, −32) 生成了什么？它为scalbn(intValue, −32) 生成了什么？还是intValue * 0x1p-32？

问题是“如何划分float 或double”但示例代码是uint32_t intV ... uint32_t intV。 dividend的真实类型是什么？

您是否尝试过将您的数字除以 1u

如果您可以断言字节顺序，您可以直接对 32 位半部分甚至包含指数的 16 位部分进行操作。这可能会生成更少的代码。

【参考方案1】：

如果您可以断言 double 是使用 IEEE 754 double-precision binary floating-point format: binary64 存储的，具有与 64 位整数相同的字节顺序和对齐要求，并且它的值足够大以使结果仍然是正常值，你可以直接用这个表达式破解表示，它应该编译成 2 或 3 条指令：

*(uint64_t *)&doubleValue -= 32ULL << 52;

然而，这种形式的类型双关语可能会给积极的优化器带来麻烦，因为它违反了 C 别名规则，因为类型 double 的值是通过指向不同类型的指针而不是字符指针来访问的。可以通过union 使用更好的类型双关语形式，它可以在大多数编译器中正常工作：

union  double d; uint64_t u;  u = doubleValue;
u.u -= 32ULL << 52;
doubleValue = u.d;

要完全避免 C 别名问题，您可以使用memcpy：

uint64_t u;
memcpy(&u, &doubleValue, sizeof u);
u -= 32ULL << 52;
memcpy(&doubleValue, &u, sizeof u);

一个好的优化编译器应该将这些 memcpy 调用转换为单个指令。

【讨论】：

...如果double 对齐要求不超过uint64_t...（也许在ARM 中这不是问题。）

不过，基于类型的别名可能会咬你一口。最好正确地执行双关语。

你不能像这样访问错误类型的对象。您可以将其存储到联合中并在联合中进行操作，或者 memcpy 将其存储到 `uint64_t` 对象并返回

要使代码符合C的别名规则，可以使用double foo(int x) union double d; uint64_t u; t = x ; t.u -= (uint64_t) 32 << 52; return t.d; 。

@chux-ReinstateMonica：你的意思是如果uint64_t的对齐要求不超过double的对齐要求？