【学习笔记】使用魔数快速求立方根

简介

介绍使用魔数0x2a517d47快速求立方根 x 3 {\sqrt[3]{x}} 3x 的C语言实现和公式的推导。

代码

float MagicCubeRoot(float x)
{float xthird = 0.333f * x;int i = *(int*)&x;i = (0x2a517d47 + (0.333f * i));x = *(float*)&i;x = 0.667f * x + xthird / (x * x);return x;
}

代码用于快速计算立方根 x 3 {\sqrt[3]{x}} 3x 。

代码中核心部分是

i = (0x2a517d47 + (0.333f * i));

该行代码就完成了计算立方根。

另外再使用一次牛顿迭代法提高下精度

x = 0.667f * x + xthird / (x * x);

所以整个计算过程就是

1.i = *(int*)&x;

将输入的数转换成整数

2.i = (0x2a517d47 + (0.333f * i));

通过魔数完成立方根的计算。

3.x = *(float*)&i;

转换回浮点数。

4.x = 0.667f * x + xthird / (x * x);

使用一次牛顿迭代法提高下精度。

完成快速计算立方根 x 3 {\sqrt[3]{x}} 3x 。

如果需要提高精度，可以多进行一次牛顿迭代。

float MagicCubeRoot(float x)
{float xthird = 0.333f * x;int i = *(int*)&x;i = (0x2a517d47 + (0.333f * i));x = *(float*)&i;x = 0.667f * x + xthird / (x * x);x = 0.667f * x + xthird / (x * x);return x;
}

注意，上面的函数只支持正数，传进来的值x需要大于或等于0。

如果要支持负数，可以增加如下判断处理。

float MagicCubeRoot(float x)
{if (x < 0){x = -x;float xthird = 0.333f * x;int i = *(int*)&x;i = (0x2a517d47 + (0.333f * i));//i = (int) (0x2a517d3c + (0.333f * i));x = *(float*)&i;x = 0.667f * x + xthird / (x * x);//x = 0.667f * x + xthird / (x * x);x = -x;}else{float xthird = 0.333f * x;int i = *(int*)&x;i = (0x2a517d47 + (0.333f * i));//i = (int) (0x2a517d3c + (0.333f * i));x = *(float*)&i;x = 0.667f * x + xthird / (x * x);//x = 0.667f * x + xthird / (x * x);}return x;
}

牛顿迭代法计算方式

一般计算立方根可以使用牛顿迭代法。

x n + 1 = x n − f ( x n ) f ′ ( x n ) x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} xn+1=xn−f′(xn)f(xn)

我们计算立方根的公式：

y = x 3 = x 1 3 y = {\sqrt[3]{x}} = x^{\frac 13} y=3x =x31

所以

y 3 = x {{y}^3} = x y3=x

构建以y为自变量的函数方程为

f ( y ) = y 3 − x = 0 f(y) = {{y}^3} - x = 0 f(y)=y3−x=0

f ′ ( y ) = 3 y 2 = 0 f'(y) = {3}{y^2} = 0 f′(y)=3y2=0

将 f ( y ) f(y) f(y)和 f ′ ( y ) f'(y) f′(y)带入

y − f ( y ) f ′ ( y ) = y − y 3 − x 3 y 2 y - \frac{f(y)}{f'(y)} = y - \frac{{{y}^{3}} - x}{{3}{y^2}} y−f′(y)f(y)=y−3y2y3−x

= 1 3 ( 2 y + x y 2 ) = \frac{1}{3}{(2y + \frac{x}{y^2})} =31(2y+y2x)

所以

y n + 1 = 1 3 ( 2 y n + x y n 2 ) y_{n+1} = \frac{1}{3}{(2y_n + \frac{x}{y_n^2})} yn+1=31(2yn+yn2x) = ( 2 y n + x ÷ y n 2 ) ÷ 3 (2{y_n} + {x}÷{y_n^2})÷3 (2yn+x÷yn2)÷3

实际计算，设 x = 2 x = 2 x=2， y 1 = 2 y_1 = 2 y1=2

y 2 = ( 2 × 2 + 2 ÷ 2 2 ) ÷ 3 = 1.5 y_2 = (2 × 2 + 2 ÷ 2^2) ÷ 3 = 1.5 y2=(2×2+2÷22)÷3=1.5

y 3 = ( 2 × 1.5 + 2 ÷ 1. 5 2 ) ÷ 3 = 1.296296 y_3 = (2 × 1.5 + 2 ÷ 1.5^2) ÷ 3 = 1.296296 y3=(2×1.5+2÷1.52)÷3=1.296296

y 4 = ( 2 × 1.296296 + 2 ÷ 1.29629 6 2 ) ÷ 3 = 1.260932 y_4 = (2 × 1.296296 + 2 ÷ 1.296296^2) ÷ 3 = 1.260932 y4=(2×1.296296+2÷1.2962962)÷3=1.260932

y 5 = ( 2 × 1.260932 + 2 ÷ 1.26093 2 2 ) ÷ 3 = 1.259922 y_5 = (2 × 1.260932 + 2 ÷ 1.260932^2) ÷ 3 = 1.259922 y5=(2×1.260932+2÷1.2609322)÷3=1.259922

经过多次迭代，计算出来的结果 2 3 ≈ 1.259922 {\sqrt[3]{2}} \approx 1.259922 32 ≈1.259922。

使用牛顿迭代法计算立方根需要大量的浮点数运算，运算速度会远远弱于整数运算。

算法讲解

接下来解释该算法的核心i = (0x2a517d47 + (0.333f * i));。

首先需要知道浮点数的基本概念。

32位浮点数基本概念

浮点数由三部分组成：符号，指数和尾数。

32位浮点数，用二进制表示

用公式表示就是

( − 1 ) s ( 1 + m ) 2 e (-1)^s(1+m)2^e (−1)s(1+m)2e

这里s是符号(sign)，e是指数(exponent)，m是尾数(mantissa)。

因为这里计算立方根只计算正数，负数可以先将负数转换成正数计算，之后再转换成负数。

所以从假设符号位是 0 ，公式简化为

( 1 + m ) 2 e (1+m)2^e (1+m)2e

指数部分 e e e的数值范围， − 127 ≤ e ≤ 128 -127 \leq e \leq 128 −127≤e≤128。

尾数部分 m m m的数值范围， 0 ≤ m < 1 0 \leq m < 1 0≤m<1。

如果将浮点数转换成整数时，整数的数值就是

M + L E M + LE M+LE

这里 E E E表示指数(exponent)， M M M表示尾数(mantissa)， L L L是 2 23 2^{23} 223。

将指数部分看做是整数，用 E E E来表示，那么范围是 0 ≤ E ≤ 255 0 \leq E \leq 255 0≤E≤255。 E E E如果减去127，范围变成-127 - 128，变成一个有符号数。

所以 E E E和 e e e的转换关系是 e = E − B e = E - B e=E−B， B B B是127。

如果将尾数部分看做是整数，用 M M M来表示。数值范围是 0 ≤ M < 2 23 0 \leq M < 2^{23} 0≤M<223。

所以 M M M和 m m m的转换关系是 m = M L m = \frac{M}{L} m=LM， L L L是 2 23 2^{23} 223。

推导过程

接下来进入推导过程。

给定一个 x x x，计算立方根 y y y：

y = x 3 = x 1 3 y = {\sqrt[3]{x}} = x^{\frac 13} y=3x =x31

首先对等式两边取以 2 为底的对数

log ⁡ 2 y = 1 3 log ⁡ 2 x \log_2 y = {\frac 13}\log_2 x log2y=31log2x

将x和y用浮点数替换：

log ⁡ 2 ( ( 1 + m y ) 2 e y ) = 1 3 ( log ⁡ 2 ( ( 1 + m x ) 2 e x ) ) \log_2 ((1+m_y)2^{e_y}) = {\frac 13}(\log_2 ((1+m_x)2^{e_x})) log2((1+my)2ey)=31(log2((1+mx)2ex))

log ⁡ 2 ( 1 + m y ) + e y = 1 3 ( log ⁡ 2 ( 1 + m x ) + e x ) \log_2 (1+m_y) + e_y = {\frac 13}(\log_2 (1+m_x) + e_x) log2(1+my)+ey=31(log2(1+mx)+ex)

算式两边都有这样的项

log ⁡ 2 ( 1 + v ) \log_2(1 + v) log2(1+v)

其中v的值范围0 < v v v < 1。

当v的取值在0到1之间时，这个函数和一条直线很接近：

方程式：

log ⁡ 2 ( 1 + v ) ≈ v + σ \log_2(1 + v) \approx v + \sigma log2(1+v)≈v+σ

其中 σ \sigma σ是一个常数，可以通过调整这个常数让两个曲线更加近似。

上面的方程

log ⁡ 2 ( 1 + m y ) + e y = 1 3 ( log ⁡ 2 ( 1 + m x ) + e x ) \log_2 (1+m_y) + e_y = {\frac 13}(\log_2 (1+m_x) + e_x) log2(1+my)+ey=31(log2(1+mx)+ex)

进行化简

m y + σ + e y ≈ 1 3 ( m x + σ + e x ) m_y + \sigma + e_y \approx {\frac 13}(m_x + \sigma + e_x) my+σ+ey≈31(mx+σ+ex)

接下来用整形解释下的指数和尾数来替代浮点数解释：

M y L + σ + E y − B ≈ 1 3 ( M x L + σ + E x − B ) \frac{M_y}{L} + \sigma + E_y - B \approx {\frac 13}(\frac{M_x}{L} + \sigma + E_x - B) LMy+σ+Ey−B≈31(LMx+σ+Ex−B)

进行化简

M y L + E y ≈ 1 3 ( M x L + σ + E x − B ) − σ + B \frac{M_y}{L} + E_y \approx {\frac 13}(\frac{M_x}{L} + \sigma + E_x - B) - \sigma + B LMy+Ey≈31(LMx+σ+Ex−B)−σ+B

M y L + E y ≈ 1 3 ( M x L + E x ) − 2 3 ( σ − B ) \frac{M_y}{L} + E_y \approx {\frac 13}(\frac{M_x}{L} + E_x) - \frac{2}{3}(\sigma - B) LMy+Ey≈31(LMx+Ex)−32(σ−B)

M y + L E y ≈ 2 3 L ( B − σ ) + 1 3 ( M x + L E x ) M_y + LE_y \approx {\frac 23}L(B - \sigma) + {\frac 13}(M_x + LE_x) My+LEy≈32L(B−σ)+31(Mx+LEx)

两边都得到了整形解释的值：

I y ≈ 2 3 L ( B − σ ) + 1 3 I x \mathbf{I_y} \approx {\frac 23}L(B - \sigma) + {\frac 13}\mathbf{I_x} Iy≈32L(B−σ)+31Ix

其中 2 3 L ( B − σ ) {\frac 23}L(B - \sigma) 32L(B−σ)是一个常数。

用代码表示就是：

i = K + i / 3;

K就是算法中的常数。通过选取合适的 σ \sigma σ的值，就能得到K值。

魔数的选定

计算魔数所用的 σ \sigma σ可以采用穷举搜索或者二分法寻找最优值。

关于 σ \sigma σ值的选定的详细内容，可以参考以下论文：

Lomont, Chris. Fast Inverse Square Root. February 2003.

McEniry, Charles. The Mathematics Behind the Fast Inverse Square Root Function Code . August 2007 .

这里选定的值为 σ \sigma σ = 0.0450465。

带入计算出K值：

K = 2 3 L ( B − σ ) = 2 3 2 23 ( 127 − 0.0450465 ) ≈ 709983559 = 0 x 2 a 517 d 47 K = {\frac 23}L(B - \sigma) = {\frac 23}{2^{23}}(127 - 0.0450465) \approx 709983559 = 0x2a517d47 K=32L(B−σ)=32223(127−0.0450465)≈709983559=0x2a517d47

这里的K值就是公式中的0x2a517d47。

参考资料

平方根倒数速算法
0x5f3759df

本文链接：https://blog.csdn.net/u012028275/article/details/113822421