几种快速平方根算法及其拓展应用

Play_with · 2023-08-23 15:31:07 · 个人记录

前置芝士

全文均为实数，不考虑复数。

快速幂

微积分入门：1、2、3、4、5。（其实只影响一些细节性的理解）

单精度浮点数在计算机中的储存原理

如图

其中：

符号位（sign）用来表示正负，接下来简写为 S。若为正数，S=0，否则 S=1。由于后面的算法都是间接或直接计算平方根的，用计算机对负数开平方根失去了意义，所以之后所有的数都默认为正数；

指数位（exponent）用来表示该浮点数在二进制下的科学计数法的指数部分。为了能够表示出小数部分，它必须要区分正数和负数，所以会被偏置 127 以适应该性质，接下来简写为 E；

尾数位（mantissa）用来存放该数的尾数部分。由于尾数必以 1 开头，故忽略这个 1 不计。所以 23 位的位置可以存 24 位数，接下来简写为 M。

我们举个栗子 9.5，它将先被表示为二进制数 1001.1，即 2^3+2^0+2^{-1}，接着被分整数部分 2^3+2^0 和小数部分 2^{-1} 分别储存。
用公式写出来，就是类似这样：

x=(-1)^S\times(1.0+\sum_{i=1}^{23} b_{23-i}2^{-i})\times2^{E-127}=(-1)^S\times(1.0+\frac{M}{2^{23}})\times2^{E-127}

的式子，想感受具体转换就点这个或这个。

从上面也能看出，并非所有的浮点数都能毫无精度损失地储存在计算机中，例如 1.14 在单精度下就是 1.13999998569488525390625。

而整数可以视为没有小数部分的浮点数，于是就有了如下表示方法。

x=(-1)^S\times (E\times 2^{23}+M)

正文

还是先贴个图放着：

二分法

我们将开平方根的结果设为 y，则有 y=\sqrt x。整个函数如你所见，单调透顶，可以二分，时间复杂度为 \log_2{n}。

牛顿迭代法

在这里我们用 y=f(x) 来描述该函数，用 x_i 来代表我们所找到的近似根。这里只讲它最本来的应用：求根。
随便找一个函数，比如这个：

那么，会有

(x_0-x_1)f(x_0)'=f(x_0)

稍加变形，就有

x_1=x_0-\frac{f(x_0)}{f(x_0)'}

而且接下来还有 x_2，x_3 等等，它们都满足这个式子，即

x_{i+1}=x_i-\frac{f(x_i)}{f(x_i)'}

于是你就可以无穷逼近这个根，就像这样，

需要注意的是，假如是求一般的函数，一定要在单调有根的区间内使用该方法，否则会有漏根的情况。同时还要注意极值点，因为该点的 f(x)'=0，此时作牛顿迭代就会失去意义。

在求平方根时则要变形一下函数，假设你要求 n 的平方根，则将 f(x)=x^2-n 带入求根。显然它在第一、四象限是单调的，所以最终结果当然可以使 f(x)逼近 0 ，所以此时的 x 就近似是 n 的平方根，时间复杂度近似是常数复杂度，但仍然没有接下来的算法快。

卡马克快速倒数平方根算法

重头戏来了，不过这次先贴上那段来自 Quake3 的著名代码:

float Q_rsqrt( float number )
{
    long i;
    float x2, y;
    const float threehalfs = 1.5F;
    x2 = number * 0.5F;
    y  = number;
    i  = * ( long * ) &y; // evil floating point bit level hacking
    i  = 0x5f3759df - ( i >> 1 ); // what the fuck?
    y  = * ( float * ) &i;
    y  = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
//  y  = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed
    return y;
}

它是用来求 \frac{1}{\sqrt x} 的，比直接调用库函数并手动倒数快了近 4 倍。

先分析一下为什么这段代码这么高效吧。它几乎只进行了加减和乘，少数用除法的时候都是靠位运算实现的；而在 CPU 中，加减的速度跟位运算相当，乘法的速度比加减慢近 10 倍，除法的速度比加减慢至少 20 倍，且数据类型精数越高，效率越低。

进入正题

假设你要求 y 的倒数平方根，则有 y=\frac{1}{\sqrt x}。

看这三行，先忽略中间一行~~和它的注释~~。

    i  = * ( long * ) &y; // evil floating point bit level hacking
    i  = 0x5f3759df - ( i >> 1 ); // what the fuck?
    y  = * ( float * ) &i;

它将输入的数的整数单独部分取了出来。

然后再分析一下中间那行哔格爆表的代码。

想一想，函数 y=\frac{1}{\sqrt x} 相当于 y=x^{-2}。但是有烦人的指数一直在 x 头上，这让人很不好处理。再联想一下，在求极限或者别的什么奇怪的问题的时候，会遇到指数难以处理的情况。我们一般会取对数，把乘方变成乘法，把乘法变成加法，把高维降成低维。

所以我们对 y=x^{-\frac{1}{2}} 取对数，有

\log_2 y=-\frac{1}{2}\log_2 x

我们将它们都视为单精度浮点数，把一开始储存浮点数的式子套进去，就会有

\log_2((1.0+\frac{M_y}{2^{23}})\times2^{E_y-127})=-\frac{1}{2}\log_2((1.0+\frac{M_x}{2^{23}})\times2^{E_x-127})

化简一下，则会有

(E_y-127)+\log_2(1.0+\frac{M_y}{2^{23}})=-\frac{1}{2}(E_y-127)--\frac{1}{2}\log_2(1.0+\frac{M_x}{2^{23}})

为了进一步推导，我们先来看看近似函数的一些知识。根据泰勒展开，当 x 很小时，\ln(x+1)\approx x；（如下图）

所以我们推测 \log_2(x+1)\approx x+k。具体的 k 值对于这篇博客来说过于复杂了，细节可以参考这篇宝藏般的博客。原码所取的 k=0.0450465。（如下图）

我们可以将 \log_2(x+1)\approx x+k 代入之前所推导的那个式子 (E_y-127)+\log_2(1.0+\frac{M_y}{2^{23}})=-\frac{1}{2}(E_y-127)-\log_2(1.0+\frac{M_x}{2^{23}})，化简得

E_y\times 2^{23}+M_y=\frac{3}{2}(127-k)\times2^{23}-\frac{1}{2}(E_x\times2^{23}+M_x)

注意看，这个等式的左边正好就是 y 的整数部分，右边括号内正好就是 x 的整数部分，\frac{3}{2}(127-k)\times2^{23} 正好也等于 0x5f3759df，所以那句

    i  = 0x5f3759df - ( i >> 1 );               // what the fuck?

其实就是在求 y 的值！那后面几行

    y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
//  y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

是在作甚？回忆一下，浮点运算会丢失部分精度；用近似函数计算也会有误差。再回忆一下上一节的牛顿迭代法，它可以做什么？无限逼近真实值啊！所以这两行就是在作牛顿迭代降低误差！这也是为什么之前可以用近似函数来计算。事实上这个算法的优越性就在于它可以更快地接近函数近似值。于是整个算法是完全的常数复杂度（因为不超过 4 次就已经远远超过精度需求了），而且常数还小得离谱。

一开始也说了，这个算法是求 \frac{1}{\sqrt x} 的。若要用来求 \sqrt x，为了答案能更精确且过程更迅速，我们选择在最后返回 \frac{1}{y}，而非一开始就输入 \frac{1}{x}（因为 \frac{1}{\sqrt{\frac{1}{x}}}=\sqrt{x}）。

拓展

我们其实可以选择把上面的公式推广，不单单只是求 y=x^{-\frac{1}{2}} 的函数值。现在，我们假设要求 y=x^a 的函数值。还是先取对数，则

\log_2 y=a\log_2 x

再套入单精度浮点数的模板，就有

\log_2((1.0+\frac{M_y}{2^{23}})\times2^{E_y-127})=a\log_2((1.0+\frac{M_x}{2^{23}})\times2^{E_x-127})

再化简，并带入 \log_2(x+1)\approx x+k，可得

E_y\times 2^{23}+M_y=(1-a)(127-k)\times2^{23}+a(E_x\times2^{23}+M_x)

又是这个熟悉的形式。但这只是一个近似解，还要做牛顿迭代。

因为你要求 y=x^a 的函数值，但你又不能用这个直接做迭代法，所以你想变形这个函数。设最终结果为 n，你把它变形成了 y=x^{\frac{1}{a}}-n。于是

y'=\frac{1}{a} x^{\frac{1}{a}-1}

x_{i+1}=x_i-\frac{y_i}{y_i'}=x_i(1-a+anx_i^{-\frac{1}{a}})

从公式也能看出来，\frac{1}{a} 得是个整数才能使用这个算法，而且计算 x^{\frac{1}{a}} 多数时候还得用快速幂。所以这个算法的长处是：快速计算 a 次方根或它的倒数，即形如 x^{\pm\frac{1}{b}} 的式子，其中 b 为正整数。不过，最好是计算形如 x^{-\frac{1}{b}} 的式子，稍后会解释为什么。

为了之后称呼方便，我们称这个算法为拓展卡马克算法。

不过这个算法还有拓展的余地：求任何浮点数次幂的值。还是假设我们要求 y=x^a，那么浮点数 a 可以写成 b+c 的形式，其中左边是整数部分，右边为小数部分。现在有两种处理方法，先讲简单的那种。

再将 c 视为最简真分数 \frac{d}{e}，于是我们要求 y=x^{b+\frac{d}{e}}。整数部分直接交给快速幂，小数部分可以视为 y=(x^{\frac{1}{e}})^d；x^\frac{1}{e} 可以用拓展卡马克解决，接着也可以交给快速幂了。

接下来是更优的那种。先联想一下普通快速幂，它只能支持整数幂，将这个幂转换成二进制来处理。但浮点数也是可以用二进制表示的。例如 x^{0.75} 就可以利用二进制表示为 x^{\frac{1}{2}+\frac{1}{4}}，接着又可以用拓展卡马克解决。

现在，拓展卡马克算法就可以解决浮点数次幂了。

双精度

显然单精度不够用，所以带来双精度版本。

我们先来了解一下双精度浮点数在计算机中的储存方法：

具体规则与双精度相同，只有内存上的差异。其表示方法如下：

x=(-1)^S\times(1.0+\sum_{i=1}^{52} b_{52-i}2^{-i})\times2^{E-1023}=(-1)^S\times(1.0+\frac{M}{2^{52}})\times2^{E-1023}

同样地，超长整形也可以视为没有小数部分的双精度浮点数，所以可以如下表示：

x=(-1)^S\times(E\times 2^{52}+M)

接着再在草稿本上演算一下之前的步骤，就能发现：

E_y\times 2^{52}+M_y=(1-a)(1023-k)\times2^{52}+a(E_x\times2^{52}+M_x)

其余部分一致。也就是说，要想由单精度进化为双精度，只需要在原来的基础上改掉那个奇怪的数字和数据类型！

但是不要忘了精度问题！改成双精度后需要多做几次牛顿迭代。个人推荐 4 次，~~米四达狂喜~~。

至此，我们终于理解了这数串强大，神秘而又美丽的代码。现在来练练手吧！

我们拿求立方根举例子，把 a=\frac{1}{3} 带入可得如下代码。

#include<bits/stdc++.h>
using namespace std;
long double f(double number){//不开long double会#1 WA
    long long i;
    double y,x2=1.0*number/3;
    const double k = 2.0/3;
    y = number;
    i = *(long long*) &y;
    i = 0x2A9F8507753418BC + 1.0*i/3;
    y =  *(double*) &i;
    y = k*y+1.0*x2/(y*y);
    y = k*y+1.0*x2/(y*y);
    y = k*y+1.0*x2/(y*y);
    y = k*y+1.0*x2/(y*y);
    return y;
}
long long n;
int main(){
    cin>>n;
    cout<<(long long)f(n);
    return 0;
}

不过你肯定发现了，在牛顿迭代时出现了万恶的除法运算。而且回顾之前的迭代式子

x_{i+1}=x_i-\frac{y_i}{y_i'}=x_i(1-a+anx_i^{-\frac{1}{a}})

会发现：只要计算的不是形如 x^{-\frac{1}{b}} 的式子，就一定会出现除法运算。所以我们之前发现的拓展卡马克算法为了优化，可以在计算 x^{\frac{1}{b}} 时先算出 x^{-\frac{1}{b}}，最后再输出倒数，尽可能地减少除法运算。

后记

笔者第一次写这么复杂的东西，可能有很多细节或错误需要完善，敬请私信。如果以后有时间会补代码上去。

写完翻日报的时候才发现日报里有类似的，幸好他后面没讲原理和拓展算法。不过他前面的牛顿迭代法讲得比我好，看不懂我的请看他的。