浮点数转换

浮点数进制转换是一个很难的问题，几十年来，学术界和工业界都在不停地优化这个过程。rsc总结了这些算法的特点，提取出了非常简洁的浮点数打印和解析算法。相比于其它算法，该算法更容易理解，同时性能也非常出色。

浮点数

IEEE 754定义了浮点数的格式，它将一个二进制浮点数拆成三部分：符号位、阶码和尾数。

┌────────────┬──────────────────┬────────────────────────────────────────────┐
│ sign 1 bit │ exponent 11 bits │            fraction 52 bits                │
└────────────┴──────────────────┴────────────────────────────────────────────┘
63           62                 52                                           0

以 float64 为例，记符号位为 $s$，阶码字段为 $E$，fraction 字段为 $F$。$E$ 和 $F$ 都是编码中的无符号整数：

$$0\le E\le 2047,0\le F<2^{52}.$$

当 $1\le E\le 2046$ 时，编码表示规格化数。此时有效尾数有一个隐含前导位 $1$，整个64位表示的数值为

$$(-1{)}^s\cdot \left(1+F\cdot 2^{-52}\right)\cdot 2^{E-1023}$$

真实指数为 $e=E-1023$（使用bias避免了显式的符号位），规格化数的指数范围是 $e\in [-1022,1023]$。注意，$e$的下界不是 $-1023$，因为 $E=0$（$e=-1023$）表示零和非规格化数；上界不是 $1024$，因为 $E=2047$ 表示无穷大和 NaN。

为了让解析和打印尽量转化为整数运算，我们需要将有效尾数改写成整数风格。有效尾数 $1+F\cdot 2^{-52}$ 放大 $2^{52}$ 倍后得到整数尾数 $m=2^{52}+F$。为了保持数值不变，指数相应减去 52，因此同一个规格化数可以写成

$$(-1{)}^s\cdot m\cdot 2^{e-52}$$

令 $e^{\prime }=e-52$，规格化数便可写成：

$$(-1{)}^s\cdot m\cdot 2^{e^{\prime }}m\in [2^{52},2^{53}),e^{\prime }\in [-1074,971]$$

进制转换

十进制表示的数可记为 $d\cdot 10^p$，二进制表示的数可记为 $m\cdot 2^e$，其中 $m,d\in {\mathbb{Z}}_{>0}$，$e,p\in \mathbb{Z}$。解析与打印的核心，是在这两类有限表示之间建立对应关系。若二者数值相同，则有

$$m\cdot 2^e=d\cdot 10^p$$

将上述等式分别对 $m$ 和 $d$ 求解，得到两个互逆的公式：

$$\begin{array}{c}m=d\cdot 10^p\cdot 2^{-e}\end{array}\text{\hspace{1em}(1)}$$

$$\begin{array}{c}d=m\cdot 2^e\cdot 10^{-p}\end{array}\text{\hspace{1em}(2)}$$

将指数符号另行处理后，我们可以得到一个通用的计算形式：

$$y=x\cdot 2^e\cdot 10^p,x\in {\mathbb{Z}}_{>0},e,p\in \mathbb{Z}\text{\hspace{1em}(3)}$$

三元组 $(x,e,p)$ 可视为 $y$ 的精确表示。然而浮点数不需要绝对精确，在足够的有限精度约束下，这个式子可以高效地计算出来。

未舍入缩放

浮点数运算定义为先进行无限精度运算，再舍入到最接近的实际浮点数。实际的浮点数硬件不会使用无限精度进行运算，而是保留比有效位多几位的精度，保证后续能正确舍入即可。因此，我们可以模拟硬件来计算公式 $(\text{3})$ ：控制缩放因子$e$或者$p$，让结果保存在一个64位的整数中，其高位承载缩放后的整数部分，低两位分别保存 Guard 和 Sticky，从而后续基于该式子的计算在整数域中进行，并且使用方可以基于低两位信息准确判断边界进行舍入。

例如, $10/3$ 的二进制表示为 $11.010101\dots$。它的未舍入形式为$11.01$，通过缩放直接将其表示为整数$1101_2$。定义uscale是计算式$(\text{3})$的未舍入形式，$uscale(x,e,p)=⟨x\cdot e\cdot p⟩$，这可以通过大数运算来实现：

func uscale(x uint64, e, p int) unrounded {
    num :=   mul(big(4), big(x), pow(2, max(0, e)),  pow(10, max(0, p)))
    denom := mul(                pow(2, max(0, -e)), pow(10, max(0, -p)))
    div, mod := divmod(num, denom)
    return unrounded(div.uint64() | bool2[uint64](!mod.isZero()))
}

浮点数解析

要解析一个浮点数字符串，如 "3.14e10"，可以先做词法分析，得到

$$d=314,p=10-2=8.$$

该字符串表示的数值为 $d\cdot 10^p=314\cdot 10^8$。

由公式 $(\text{1})$ 可知，对 $d\cdot 10^p$ 乘以$2^e$ 做二进制伸缩($m=d\cdot 10^p\cdot 2^e$)，然后计算 $round(m)\cdot 2^{-e}$ 即可求得尾数，最后将$e$和$m$按照IEEE 754格式打包即可。然而，这个计算过程中有$m$和$e$两个未知数，怎么办呢？解决办法是先根据约束限定$e$的取值，然后再计算出$m$后，再进行修正即可得到最终的数值。伸缩比例只需使尾数 $m$ 占 53 个二进制位即可保证值的准确性，因此 $m\in [2^{52},2^{53})$。于是伸缩指数 $e$ 应满足：

$$d\cdot 10^p\cdot 2^e\in [2^{52},2^{53})=2^{52}\cdot [1,2)$$

取以 2 为底的对数：

$${\log }_2(d\cdot 10^p)+e\in 52+[0,1)$$

移项后，

$$\begin{array}{rl}e& \in 52-{\log }_2(d\cdot 10^p)+[0,1)\\ e& \in 52-({\log }_2(d\cdot 10^p)-[0,1))\end{array}$$

由于 $e$ 必须是整数，该区间中唯一的整数为

$$\begin{array}{rl}e& =52-\lfloor {\log }_2(d\cdot 10^p)\rfloor \\ e& =52-\lfloor {\log }_{2}d+p\cdot {\log }_{2}10\rfloor \end{array}\text{\hspace{1em}(4)}$$

因为$e$只影响缩放系数，不影响最终计算值的精度，因此我们并不需要一个准确的$e$，而是只需要得到e的近似值。${\log }_{2}d$ 的计算开销比较大，我们可以近似 ${\log }_{2}d\approx bits(d)-1$ 。因此公式 $(\text{4})$ 的计算可以进一步进行优化：

$$\begin{array}{rl}e& =52-\lfloor bits(d)-1+p\cdot {\log }_{2}10\rfloor \\ & =53-bits(d)-\lfloor p\cdot {\log }_{2}10\rfloor \end{array}\text{\hspace{1em}(5)}$$

确定 $e$ 后，代回 $m=d\cdot 10^p\cdot 2^e$求得$m$，即可构造出float64。

需要注意的是，$lo{g}_{2}d$ 和近似值 $bits-1$ 可能相差一倍，因此计算得到$m$如果超过范围$[2^{52},2^{53})$，则需要进行修正：缩小 $m$ 和 $e$ ，使得 $m=u/2，e=e-1$。（$m$ 的缩小是要在未舍入值上进行，以保证舍入信息准确）。

func Parse(d uint64, p int) float64 {
	b := bits.Len64(d)
	e := min(1074, 53-b-log2Pow10(p))
	u := uscale(d<<(64-b), e, p)
	m := u.round()
	if m >= 1<<53 {
		m, e = u.rsh(1).round(), e-1
	}
	return pack64(m, -e)
}

浮点数打印

如果不考虑性能和最优可读性，二进制格式的浮点数打印成十进制格式是一个比较简单的问题。然而，实现正确高效的浮点数打印则非常困难。根据”四舍六入五成双”的舍入规则，浮点数实际上表示一段区间，这个区间内的十进制数可能有多个，无论选择哪一个，它们都可以表示这个浮点数。为了让转换稳定易读，我们必须选择满足这三个约束的十进制输出：

• 正确：转成十进制表示的数字可以再完美转回二进制浮点数；

• 易读：用最短的十进制表示这个浮点数

  e.g. 如果 xxxxxx.19e10 和 xxxxxx.2e11 都能表示这个浮点数，那就选择后者；

• 稳定：同一个浮点数转换成的十进制表示要唯一，不能在不同时候输出不同的表示；

那如何实现这个转换呢？

由公式 $(\text{2})$ 可知，$m$ 和 $e$ 是输入，$d$ 和 $p$ 是输出。一个等式存在两个未知数，需要先把其中一个确定，才能求解出另一个。如果我们能确定$p$，使用uscale 和 $m\cdot 2^e$ 运算就可以得到$d$。因此，问题归结为：怎么选择正确的 $p$？浮点数在数轴上是离散的，每个 $f$ 都”管辖”一段连续区间：左右邻居的中点之间的所有实数都会舍入到它。我们把这段区间的宽度称为浮点数的 footprint，记作 $footprint(f)$。对绝大多数浮点数，左右邻居距离都是 $2^e$，$footprint(f)$ 是宽度为 $2^e$ 的对称区间。但是有一个特殊例外：如果$f$ 恰为 2 的幂，左侧邻居的指数比 $f$ 低一阶，此时左右两个区间中点构成的区间宽度为：

$$\underset{\text{left}}{\underbrace{\frac{1}{2}\cdot 2^{e-1}}}+\underset{\text{right}}{\underbrace{\frac{1}{2}\cdot 2^e}}=\frac{3}{4}\cdot 2^e$$

因此$footprint(f)$ 存在 $2^e$ 和 $\frac{3}{4}\cdot 2^e$ 两种情况，我们将其称作对称形和倾斜形。

最短十进制

同一个数值 $d\cdot 10^p$ 可以有多种写法：

我们可以通过伸缩$p$来控制$d$ 的长度，而”最短打印”等价于在所有合法$p$ 中取最大值。$p$的合法性可由$footprint$ 给出：选定 $p$ 后，十进制刻度 $d\cdot 10^p$只要落在 $footprint(f)$内，输出的十进制就能正确解析回 $f$。

那么，如何得到 $p$ ？枚举所有 $p$ 不现实，需要一个能直接给出 $p$ 的方法。考察

$$\frac{\text{footprint}(f)}{10^p}\approx \frac{2^e}{10^p}$$

该式子是 $f$ 在二进制下的区间宽度与候选十进制步长的比值，它表示这个二进数区间内能容纳多少个十进制数（偏斜形下 footprint 为 $\frac{3}{4}\cdot 2^e$，与 $2^e$ 仅差一个常数因子，故用 $\approx$），这些十进制数都等价于该二进制数。我们需要找到一个明确的最短的十进制数，因此要对伸缩条件进行限制，使得这个区间中的数字个数：

• 不小于 1，使区间至少容下一个刻度，候选不落空；
• 不超过 10，使区间至多包含 10 个连续整数；其中至多一个0结尾的数，防止多个等长最短数引起歧义。

这两点约束要求 $\frac{\text{footprint}(f)}{10^p}\in [1,10)$

现在，$p$ 可以按 $footprint$ 的两种形态分别求解。

• 对称形：$\text{footprint}(f)=2^e$。

  $$\frac{2^e}{10^p}\in [1,10)⟺p=\lfloor e{\log }_{10}2\rfloor .$$

  无理数 ${\log }_{10}2$ 用一个分数近似：

  // log10Pow2(e) returns ⌊e * log₁₀2⌋.
  func log10Pow2(e int) int {
      // log₁₀ 2 ≈ 0.30102999566 ≈ 78913 / 2^18
      return (e * 78913) >> 18
  }

• 偏斜形（$f$ 是 2 的幂）：$\text{footprint}(f)=\frac{3}{4}\cdot 2^e$。

  $$\frac{(3/4)\cdot 2^e}{10^p}\in [1,10)⟺p=\lfloor e{\log }_{10}2-{\log }_{10}(4/3)\rfloor .$$

  近似计算：

  // ⌊log₁₀(3/4) + log₁₀(2**e)⌋ = ⌊e*(log₁₀2)-(log₁₀(4/3))⌋.
  func skewed(e int) int {
      return (e*631305 - 261663) >> 21
  }

可靠性

现在，我们可以分情况把 $p$ 算出来了，这还不够，$p$ 对应的 $footprint(f)$ 可能包含多个连续的十进制整数。我们要得到最终的十进制数字，一方面要在这个范围内选择合适的数字，另一方面还要考虑边界情况：如果 $p$ 对应的 $footprint(f)$ 范围内只有一个十进制数且这个数还位于中点上，就需要要考虑这个数是否会舍入到 $f$。如果不能舍入到 $f$ ，那要考虑”这个区间可能内没有合适的数” （出现这种情况从根本上就否定了这个算法）。记 $r=\text{footprint}/10^p\in [1,10)$。$r>1$ 时十进制在区间内严格大于一个刻度。因此即使两个端点都恰为整数(大部分情况不可能)，奇尾数排除两个端点候选后，区间内仍至少含一个十进制数可供选择；

然而， $r=1$ 时候选区间宽度恰好一个刻度时，则需要考虑这个唯一的整数会不会存在端点上，否则奇尾数把端点候选排除后区间内候选的数就不存在了。那 $r$ 究竟会不会真的等于 1？

• 偏斜形下 $r=\frac{3}{4}\cdot 2^e\cdot 10^p$。要它等于 1，分母里的 3 得被 $10^p$ 约掉，但 $10^p=2^p\cdot 5^p$ 只有因子 2 和 5，3 永远无法被约掉。

• 对称形下 $r=2^e\cdot 10^p$。要它等于 1，只有 $e=p=0$，此时缩放是恒等变换，十进制刻度直接就是二进制整数，本来就不存在歧义。

所以在所有有意义的情形下 $r$ 都严格大于 1：如果尾数是奇数，即使把端点（中点）排除掉，footprint 区间内仍然存在可选的十进制数。

实现

把这些逻辑片段拼到一起，最短打印就是下面这个函数：

func (u unrounded) nudge(δ int) unrounded { return u + unrounded(δ) }
func (u unrounded) ceil() uint64    { return uint64((u + 3) >> 2) }
func (u unrounded) floor() uint64     { return uint64((u + 0) >> 2) }
func Short(f float64) (d uint64, p int) {
    m, e := unpack64(f)

    // 按 footprint 形态算 p，使比值 footprint/10^p 落到 [1, 10)
    var left uint64
    z := 11
    if m == 1<<63 && e > minExp {
        p = -skewed(e + z)
        left = m - 1<<(z-2) // left = m - 1/4 * 2**(e+z)
    } else {
        if e < minExp {
            z = 11 + (minExp - e)
        }
        p = -log10Pow2(e + z)
        left = m - 1<<(z-1) // left = m - 1/2 * 2**(e+z)
    }
    right := m + 1<<(z-1)

    // 端点开闭按尾数奇偶判定
    odd := int(m>>z) & 1
    pre := prescale(e, p, log2Pow10(p))
    dmin := uscale(left, pre).nudge(+odd).ceil()
    dmax := uscale(right, pre).nudge(-odd).floor()

    // 优先取 0 结尾的候选：dmax/10 是 ≤ dmax 的最大十倍刻度，
    // 若它仍 ≥ dmin，说明区间内有 0 结尾的整数，trim 末尾零升一阶即最短
    if d = dmax / 10; d*10 >= dmin {
        return trimZeros(d, -(p - 1))
    }
    // 没有 0 结尾的候选；区间内若仍有多个整数，按 m·2^e 的实际值正确舍入；
    // 否则 dmin == dmax，区间内仅一个整数，直接返回
    if d = dmin; d < dmax {
        d = uscale(m, pre).round()
    }
    return d, -p
}

函数体总体上分为三步：

• 求 $p$：按 footprint 是对称形还是偏斜形，调用 log10Pow2 或 skewed 算出使比值落到 $[1,10)$ 的那个 $p$。
• 取候选区间：把 footprint 的边界 $left,right$ 换算到 $10^p$ 的刻度上，整数化为 $d_{left},d_{right}$，得到 $f$ 的候选整数区间。
• 确定 $d$：候选区间 $[d_{left},d_{right}]$ 里有 0 结尾的整数就取它，trim 末尾零等价于把 $p$ 升到最大；否则按 $m\cdot 2^e$ 的实际值挑离它最近的值。

区间 $d_{left},d_{right}$ 是否是开区间由两个端点 $d_{left},d_{right}$ 的计算来确定。

高效计算uscale

$y=x\cdot 2^e\cdot 10^p$ 的计算难点在于 $x\cdot 10^p$，有限位数下的二进制下 $10^p$ 可能没有精确表示。理论上，我们可以将$10^p$转成一个无穷精度的二进制表示：$\mu \cdot 2^k$，其中$\mu$的整数部分为128位，k为二进制指数。

$$x\cdot 10^p=x\cdot \mu \cdot 2^k,\mu \in [2^{127},2^{128}),k\in \mathbb{Z}.$$

$\mu$的范围固定后(${\log }_2\mu \in [127,128)$)，也就唯一确定了指数k。结合公式 $(\text{7})$ 可求得$k$：

$$\begin{array}{rl}{\log }_2(10^p)& ={\log }_2(\mu \cdot 2^k)\\ p\cdot {\log }_{2}10& ={\log }_2\mu +k\\ k& =p\cdot {\log }_{2}10-{\log }_2\mu \\ & \stackrel{(\text{7})}{=}\lfloor p\cdot {\log }_{2}10\rfloor -127.\end{array}$$

记 $lp=\lfloor p\cdot {\log }_{2}10\rfloor$，则我们确定了指数 $k=lp-127$，原式：

$$y=x\cdot 2^e\cdot 10^p=x\cdot \mu \cdot 2^{e+k}=x\cdot \mu \cdot 2^{e+lp-127}\text{\hspace{1em}(6)}$$

现在我们将 $x\cdot 10^p$ 转成了计算 $x\cdot \mu$ 的问题。然而$\mu$ 是理想的精确实数，$x\cdot \mu$ 无法计算出来。为了高效地在整数域中进行计算，我们可以将 $\mu$ 近似为 128 位整数 $M$，并向上取整，引入误差 $\epsilon$：

$$\begin{array}{rl}M=\lceil \mu \rceil & =\mu +\epsilon ,0\le \epsilon <1\end{array}$$

精确值 $x\cdot \mu$ 与近似值 $x\cdot M$ 的误差为 $x\cdot \epsilon$。限定输入值 $x$ 为左对齐的64位整数，那么 $0\le x\cdot \epsilon <x<2^{64}$。这个误差不足以影响我们最终构造的待舍入数：因为 $x\cdot M$ 是一个很大的整数乘积(192位)，而我们只需要64位以内的待舍入数，右移取高位后，$x\cdot \epsilon$ 直接被丢弃，这样就高效地计算出了结果。

现在，问题的核心变成了计算 $x\cdot M$ 的高位。由于 $M$ 只依赖 $10^p$，可以把各个 $10^p$ 对应的 $M$ 预先存入表中。CPU 通常只提供 64x64 乘法，因此要把 $M$ 拆成高 64 位的 $M_{hi}$ 和低 64 位的 $M_{lo}$：

$$M=M_{hi}\cdot 2^{64}+M_{lo},0\le M_{lo}<2^{64}.$$

分别计算

$$\begin{array}{rl}x\cdot M_{hi}& =hi\cdot 2^{64}+mid\\ x\cdot M_{lo}& =mi{d}_2\cdot 2^{64}+lo\end{array}$$

代回精确值：

$$\begin{array}{rl}\mu & =M-\epsilon =M_{hi}\cdot 2^{64}+M_{lo}-\epsilon \\ x\cdot \mu & =x\cdot M_{hi}\cdot 2^{64}+x\cdot M_{lo}-x\cdot \epsilon \\ & =(hi\cdot 2^{64}+mid)\cdot 2^{64}+mi{d}_2\cdot 2^{64}+lo-x\cdot \epsilon \\ & =hi\cdot 2^{128}+(mid+mi{d}_2)\cdot 2^{64}+lo-x\cdot \epsilon .\end{array}$$

未舍入数要在低位保留两个 bits，用于记录 Guard 和 Sticky 信息，因此要把 $y$ 放到 $y\cdot 2^2$ 的整数刻度。若暂时忽略低位修正，只由最高段给出候选值，则

$$\begin{array}{rl}{y}_0& =hi\cdot 2^{128}\cdot 2^{e+k}\cdot 2^2\\ & =hi\cdot 2^{e+k+130}\\ & =hi\cdot 2^{e+lp-127+130}\\ & =hi\cdot 2^{e+lp+3}.\end{array}$$

这里的 $+3$ 来自三处位数偏移：$\mu$ 被规格化到 128 位，所以 $k=lp-127$；未舍入数要保留两位，所以指数多出 $+2$；取 $hi$ 等价于已经越过完整乘积的低 128 位。

经过这样处理，可以实现最终的快速uscale版本了：

// prescale returns the scaling constants for e, p. lp must be log2Pow10(p).
func prescale(e, p, lp int) scaler {
	return scaler{pm: pow10Tab[p-pow10Min], s: -(e + lp + 3)}
}
// uscale returns unround(x * 2**e * 10**p).
// The caller should pass c = prescale(e, p, log2Pow10(p))
// and should have left-justified x so its high bit is set.
func uscale(x uint64, c scaler) unrounded {
	hi, mid := bits.Mul64(x, c.pm.hi)
	mid2, _ := bits.Mul64(x, c.pm.lo)
	mid, carry := bits.Add64(mid, mid2, 0)
	hi += carry
	sticky := bool2[unrounded](mid != 0 || hi&((1<<c.s)-1) != 0)
	return unrounded(hi>>c.s) | sticky
}

Math

• 整数范围
• $lo{g}_d$

向下取整

设 $x\in \mathbb{R}$，$m,k\in \mathbb{Z}$。若存在 $y\in [m,m+1)$，满足 $k=x-y$，则

$$\begin{array}{c}k=\lfloor x\rfloor -m\end{array}\text{\hspace{1em}(7)}$$

证明：

$$\begin{array}{rl}k=x-y,y\in [m,m+1)& ⟺y=x-k,m\le y<m+1\\ & ⟺m\le x-k<m+1\\ & ⟺m+k\le x<m+k+1\\ & ⟺\lfloor x\rfloor =m+k\\ & ⟺k=\lfloor x\rfloor -m.\end{array}$$

references

• https://hachyderm.io/@rsc