[学习笔记] 高精度除法（牛顿迭代法）

没想到竟然又在 blog 里面发了 oi 相关

由于大作业涉及到了这个知识而网上并没有看到什么好的教程决定自己写一篇。

有参考这篇博客和这个ppt，但是思路还是我自己整理的。

牛顿迭代法

这个是用来求方程 $f(x)=0$ 的根，这里提到的是用切线迭代，即 $x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}$。

这个是我们做高精度除法的核心步骤，但其方法并不是需要重点讲解的地方，所以这部分请移步wikipedia。

但知道这个方法后，我一时间也没反应过来这个和高精度除法有什么关联。

但实际上，我们要计算的是 $\lfloor\frac{a}{b}\rfloor$ ，那么如果我们能算出 $\frac{1}{b}$，再乘上 $a$ 稍作调整，即可得到答案。

而求 $\frac{1}{b}$ 则可用牛顿迭代法实现。

具体地，考虑 $f(x)=\frac{1}{x}-b$，则其零点就是 $\frac{1}{b}$。那么用牛顿迭代的式子：

$\begin{aligned} x_{n+1}&=x_n-\frac{\frac{1}{x_n}-b}{-\frac{1}{x_n^2}} \\ &=x_n+x_n-b\times x_n^2 \\ &=2x_n-b\times x_n^2 \end{aligned}$

则我们能实现加法和乘法即可。

当然，这只是初步的思路，实际上有很多具体的实现细节，比如说迭代的次数，起始和停止条件，误差和调整，等等。这也是我在大多教程中所没有读到或者读懂的。那么接下来就来讲讲这些。

问题的进一步转化

回到问题 $\lfloor\frac{a}{b}\rfloor$，其实会发现一个事情，就是求 $\frac{1}{b}$ 的话，如果我们只写了大整数类，那么是无法表示这个小数的。怎么办呢？我们考虑 $b$ 有 $m$ 位（这里的 $m$ 位是在压位后的 $P$ 进制下的，后同），求出 $\lfloor\frac{P^{2m}}{b}\rfloor$，再算 $\lfloor\frac{a\times \lfloor\frac{P^{2m}}{b}\rfloor}{P^{2m}}\rfloor$ 即可。

当然，这里涉及到两个问题：为什么选 $P^{2m}$ 以及算出来的结果和 $\lfloor\frac{a}{b}\rfloor$ 的关系。

前面已经提到，我们需要仅用大整数类来实现对 $b$ 求逆，考虑这样一件事情：

$\lfloor\frac{a}{b}\rfloor-\lceil\frac{a}{t}\rceil \leq\lfloor\frac{a}{b}-\frac{a}{t}\rfloor=\lfloor\frac{a(t-b)}{bt}\rfloor=\lfloor\frac{a(\frac{t}{b}-1)}{t}\rfloor \leq\lfloor\frac{a\lfloor\frac{t}{b}\rfloor}{t}\rfloor\leq\lfloor\frac{a\frac{t}{b}}{t}\rfloor=\lfloor\frac{a}{b}\rfloor$

其中最关键的是

$\lfloor\frac{a}{b}\rfloor-\lceil\frac{a}{t}\rceil \leq\lfloor\frac{a\lfloor\frac{t}{b}\rfloor}{t}\rfloor\leq\lfloor\frac{a}{b}\rfloor$

这意味着当我们选取一个 $t>a$ 时，求得的误差便 $\leq1$，利用余数微调即可。

而这里我们选取了 $t=P^{2m}$，所以我们需要 $a$ 的位数 $n \leq 2m$。

如果不满足的话也很简单，让 $a,b$ 左移 $n-2m$ 位即可。

而关于为什么 $P^{2m}$ ，我们暂且放在这里，待之后再谈。下面要做的就是求 $\lfloor\frac{P^{2m}}{b}\rfloor$ 了。

迭代的应用

设 $c=\lfloor\frac{P^{2m}}{b}\rfloor$，考虑如何利用迭代求得 $c$。

假设我们已经求出了 $b$ 的最高 $k$ 位，即 $b'=\lfloor\frac{b}{P^{m-k}}\rfloor$ 的答案，其为 $c'=\lfloor\frac{P^{2k}}{b'}\rfloor$。

考虑函数 $f(x)=\frac{1}{x}-\frac{b}{P^{2m}}$。

则之前提到的迭代公式可化为 $x_{n+1}=2x_n-\frac{bx_n^2}{P^{2m}}$。（中间步骤就不详细写出了）

但我们前面已经算出的 $c'$ 是 $b$ 最高 $k$ 位，即 $b'$ 的答案，那么我们要将它左移 $P$ 进制下的 $m-k$ 位再代入，即 $x_n=c'\times P^{m-k}$。

然后得到 $c^*=x_{n+1}=2c'\times P^{m-k}-b\times c'^2\times P^{-2k}$。（化简过程略）

然后则可将 $c^*$ 作为 $c$ 的一个近似值来使用，而 $c'$ 则可通过递归实现。

那么关于选取 $P^{2m}$，实际上我也并未完全弄清楚，我猜测，可能是作为一个定义，以方便地进行迭代的计算以及下文的误差分析。

接下来我们来看看当前的误差。

误差分析

设 $c'=\frac{P^{m+k}}{b}+\xi$，则

$\begin{aligned} \frac{P^{2m}}{b}-c^*&=\frac{P^{2m}}{b}-2\times (\frac{P^{m+k}}{b}+\xi)\times P^{m-k}+b\times (\frac{P^{m+k}}{b}+\xi)^2\times P^{-2k} \\ &=\frac{P^{2m}}{b}-2\frac{P^{2m}}{b}-2\xi\times P^{m-k}+\frac{P^{2m}}{b}+2b\frac{\xi}{b}\times P^{m-k}+b\xi^2\times P^{-2k} \\ &=\frac{b\xi^2}{P^{2k}} \end{aligned}$

考虑到 $c'$ 的取值范围

$\frac{P^{m+k}}{b}-1\leq \frac{P^{2k}}{b'}-1\leq c' \leq \frac{P^{2k}}{b'}\leq\frac{P^{2k}}{\frac{b}{P^{m-k}}-1}=\frac{P^{m+k}}{b-P^{m-k}}$

因而 $\xi$ 的取值范围

$-1\leq \xi \leq \frac{P^{m+k}}{b-P^{m-k}}-\frac{P^{m+k}}{b}=\frac{P^{2m}}{b(b-P^{m-k})}$

而我们有 $b \geq P^{m-1}$ （$b$ 有 $m$ 位）

所以

$\frac{b\xi^2}{P^{2k}} \leq \frac{P^{4m-2k}}{b(b-P^{m-k})^2} \leq\frac{P^{m+1}}{(P^{k-1}-1)^2}= \frac{\frac{P^{m+1}}{P^{k-1}}}{\frac{(P^{k-1}-1)^2}{P^{k-1}}} \leq\frac{P^{m-k+2}}{P^{k-1}-2}$

且 $\frac{b\xi^2}{P^{2k}} \geq 0$

由于我们递归计算需要将 $k$ 的规模缩小，继续放缩

$\frac{P^{m-k+2}}{P^{k-1}-2}=P^{m-2k+3}\frac{P^{k-1}}{P^{k-1}-2}$

一般的高精度进制 $P\geq10$，且递归时保证 $k>1$，则 $\frac{P^{k-1}}{P^{k-1}-2} \leq2$。

则误差 $\leq 2\times P^{m-2k+3}$。

取 $k=\lceil\frac{m}{2}\rceil+2$，则 $m+4 \leq 2k \leq m+5$，误差 $\leq \frac{2}{P} <1$，至多调整一次。

而当 $k$ 足够小的时候，我们就可以进行暴力计算了。

关于暴力

一个简单的想法是当 $k$ 小到 $\text{long long}$ 范围下时，便直接计算。

然而也是 $naive$ 的想法，因为从上文可以看出，$m=5$ 时 $k$ 便不会缩小了。

所以此时暴力也是 $P^{10}$，即便压 $2$ 位也会炸。

而我写的时候参考的一种写法是：

当 $k\leq10$ 时，考虑倍增的思路，算出 $b,2b,2^2b,2^3b,\dots,2^nb$。

然后再从大到小考虑计算商。

时间复杂度

设计算 $\frac{P^{2m}}{b}$ 的时间为 $T(n)$，$n$ 为位数。

则每次递归时，位数的规模减半，化为 $T(\frac{n}{2})$。

当然，这里的减半也不完全是减半，当 $n$ 较小时减小的量也较少，所以这也是 $\leq 10$ 时选择暴力的原因。

而在一次计算的范围内，计算加法和乘法的复杂度为 $O(n\log n)$。

而最终的暴力（我的写法）时 ，到 $2^kb>P^{2m}$ 时停止，则计算的复杂度为 $O(n^2\log P)$以及较大的常数，但由于此处 $n$ 规模较小所以问题不大。

最终可以得到 $T(n)=T(\frac{n}{2})+O(n\log n)$，

计算可得最终复杂度为 $O(n\log n)$。（这里涉及到主定理），也不展开讲述。

code

建议不要参考代码

由于懒直接把几个函数蒯出来了

还是等大作业cr之后再发吧

虽然发了但是懒得写注释凑合着看吧

int2048 div(const int2048 &a, const int2048 &b) { //violence
    if (b == int2048(1)) return a;
    std::vector<int2048> sta;
    std::vector<int2048> s;
    int2048 ret = b, tem = int2048(1);
    while (ret <= a) {
        sta.push_back(ret);
        s.push_back(tem);
        ret += ret;
        tem += tem;
    }
    ret = a;tem = int2048(0);
    int len = sta.size() - 1;
    while(len >= int2048(0)) {
        int2048 tmp = sta[len];
        if (tmp <= ret) ret -= tmp, tem += s[len];
        len--;
    }
    return std::move(tem);
}

int2048 inverse(const int2048 &a) {
    int m = a.num.size();
    if (m <= 9) {
        return div(int2048(1) <<= (m << 1), a); 
    }
    int k = (m + 5) >> 1;

    //used to wrongly moved by base^(k+1)
    int2048 b = a;
    b >>= (m - k); //get b prime
    int2048 c = inverse(b); //get c prime 
    int2048 c1 = c;
    c1 += c;
    c1 <<= (m - k);
    b = a * c * c;
    b >>= (k << 1);
    int2048 d = c1 - b; //get c star
    d -= int2048(1); //quotient
    c = (int2048(1) <<= (m << 1)) - d * a; //remain
    //while (c >= a) d += int2048(1), c -= a;
    //while (c < int2048(0)) d -= int2048(1), c += a;
    if (c >= a) d += int2048(1);
    //if (c < 0) d -= int2048(1);
    return std::move(d);
}

int2048 &int2048::operator/=(int2048 b){
    int l1 = num.size(), l2 = b.num.size();
    if (l1 < l2) {
        (*this) = int2048(0 - (sign*b.sign < 0));
        return (*this);
    }
    int o = sign * b.sign;
    sign = sign > 0? sign : -sign;
    b.sign = b.sign > 0? b.sign : -b.sign;
    int2048 x = (*this), y = b;
    if (o < 0) x += y - int2048(1);
    if (l1 > (l2 << 1)) {
        int tmp = l1 - (l2 << 1) ;
        x <<= tmp;
        y <<= tmp;
        l2 = l1 - l2;
        l1 = l2 << 1;
    }
    int2048 z = inverse(y);
    y = x * z; //quotient
    y >>= (l2 << 1);
    x = b * y; //the quotient times b
    int2048 ret = (*this) - x;
    //while (ret >= b) y += int2048(1), ret -= b;
    //while (ret < int2048(0)) y -= int2048(1), ret += b;
    if (ret >= b) y += int2048(1);
    //if (ret < 0) y -= int2048(1);
    (*this) = y;
    sign = o;
    return (*this);
}