高精度四则运算:从入门到入土

· · 算法·理论

Copyright (c) 2025 complete_binary_tree.

转载请标明出处。

0 前言

好像没人写高精度的各种运算的算法介绍,于是我来写下吧。(其实主要是四则运算介绍)

:::warning[本教程不是什么]{open} 本教程不是封装教程,提供的模板代码安全性和正确性不能保证。如果你 hack 了我的模板,你可以私信我,我会尝试 debug。

文中代码前后可能不一样(因为后期发现一些 bug 可能会更改前期的函数,有时候没有及时更新),请以文末的模板为准。

本教程仅介绍较为常用的高精度四则运算知识(以及前置知识),不常用的(比如 O(n^{\log_2^3}) 的 Karatsuba 分治乘法)不做介绍。

本教程仅作为算法介绍不保证代码实现的效率。如果你追求效率和实用性,你可以看看masonxiong.《目前最高效的高精度整数模板发布》(2025-09-01)。 :::

1 高精度介绍

C++ 对整数运算支持最多到 long long64 位),在部分电脑上可以做到 __int128128 位)。它们能表示的数太小了(<10^{40}),所以高精度应运而生。

高精度的思想就是用数组来存数:把数组的每一个元素当作一位,当这一位到达 base 后就进一位。通常这个 base10 的倍数(比较好读入输出),但是根据 C++ 数据类型的上限,理论每个 base 可以取 2^{32} 左右,但是由于实现上的原因基本取不到这个上限。

一般来说,高精度存储都是从低位往高位存。这样存运算较为方便。

当然,整数还有一个性质就是有符号,所以我们还需要一个来存符号的变量,通常是 bool 类型。

所以一个高精度结构体可以是:

struct BIGINT {
    int num[100005];        //存数字 
    bool f;                 //存符号 
    ll len;                 //存长度
};
//base = 10

2 数字转换:如何将常见的数据类型和高精数互相转换

整数-高精

现在我们有一个整数。我们要把它从低到高一位位拆解。

根据我们小学就学过的知识,如果我们要把一个整数 a 转化为 base 进制下的数,那么最低位就是 a \bmod \text{base},次低位就是 \left\lfloor\frac a {\text{base}}\right\rfloor \bmod \text{base},……,第 i 位就是 \left\lfloor\frac a {\text{base}^i}\right\rfloor \bmod \text{base}

我们可以用循环来模拟以上过程:

void trans( int a ){
    num[0] = 0;
    f = len = 0; //清空 
    if( a < 0 ) a = -a, f = 1; //负数 
    while( a ){
        num[len] = a % 10;  //当前位 
        a /= 10;            //这样子第 len 位就可以除 len 次 base 
        len ++;             //数字长度 +1 
    }
    if( !len ) len = 1;
}

反过来,我们把高精度转化为整数也一样,就是 \text{num}_0 \times \text{base}^0 + \text{num}_1 \times \text{base}^1+\dots+\text{num}_{\text{len-1}} \times \text{base}^{\text{len-1}}。这个也可以用循环模拟:

int getI() const{
    int ret = 0;
    for( int i = len - 1; i >= 0; -- i ){
        ret = ret * 10;     //这样子在 i 位一共会乘上 i 次 base 
        ret = ret + num[i]; //统计当前位 
    }
    if( f ) ret = -ret;
    return ret;
}

字符串-高精

我们假设字符串是 10 进制的。

如果 base = 10,那么字符串前后反转一下,每一位对应的数字就是 num 数组对应位上的数字了。高精转字符串也一样,反转一下就好了。

void trans( string a ){
    num[0] = 0;
    f = len = 0; //清空 
    reverse( a.begin(), a.end() );  //反转
    if( a.back() == '-' ){ //负数
        a.pop_back();   //删掉符号 
        f = 1; 
    }
    for( int i = 0; i < a.size(); ++ i )
        num[i] = a[i] - '0'; //注意 - '0' 后才是对应数字 
    len = a.size();
    while( len > 0 && num[len - 1] == 0 ) -- len;
    if( !len ) len = 1;
}
string getS() const{
    string ret;
    for( int i = 0; i < len; ++ i )
        ret += char( num[i] + '0' ); //string 可以用加法进行拼接
    if( f ) ret += '-';
    reverse( ret.begin(), ret.end() ); 
    if( ret.size() == 0 ) ret = "0"; //特判 0 
    return ret;
}

通过以上两种转换,我们就可以做到输入输出高精度数字了。

3 高精度数字的比较

我们小学就学过:

所以我们只要实现两个正数的比较即可。

两个正数的比较:

所以我们就可以写出以下代码:

bool abs_less( const BIGINT& b ) const{ //绝对值更小的
    if( len < b.len ) return 1;
    if( len > b.len ) return 0;
    for( int i = len - 1; i >= 0; -- i ){
        if( b.num[i] < num[i] ) return 0;
        if( num[i] < b.num[i] ) return 1;
    } 
    return 0;
}
bool operator==( const BIGINT& b ) const{ //相等的 
    if( f != b.f ) return 0;
    if( len != b.len ) return 0;
    for( int i = len - 1; i >= 0; -- i ){
        if( b.num[i] != num[i] ) return 0;
    } 
    return 1;
}
bool operator<( const BIGINT& b ) const{ //小于 
    if( f && !b.f ) return 1;
    if( !f && b.f ) return 0;
    if( f && b.f ) return !abs_less( b ) && !( *this == b );
    return abs_less( b );
}
bool operator>( const BIGINT& b ) const{ //大于 
    return !( *this < b ) && !( *this == b );
}

4 高精度的加减法运算

4.1 简单加法

我们考虑两个正数相加。

我们小学就学过竖式计算:

   1  1   4
+  9  1   9
------------
 (10)(2)(13)
------------ // 进行进位
 1 0  3   3

所以,我们可以先把两个数的对应数位加起来,最后再跑一遍进位。

显然每一位最多向前进 1,而两个数相加最多是 18,所以不会出现进位 > 1 的情况。所以两个数相加,得出的数长度最多是两个数长度之和 +1。所以时间复杂度是 O(|\text{len}|) 的。

参考代码:

BIGINT simple_add( const BIGINT& b ) const{
    BIGINT c( *this );
    c.len = max( len, b.len );
    for( int i = len; i <= c.len; ++ i ) c.num[i]=0; //防止 c.num 没清空,多清一位是要进位 
    for( int i = 0; i < b.len; ++ i ){
        c.num[i] += b.num[i]; //加法 
    }
    for( int i = 0; i < c.len; ++ i ){ //进位 
        c.num[i + 1] += c.num[i] / 10;
        c.num[i] %= 10;
    }
    if( c.num[c.len] ) ++ c.len; //判断最高位
    return c; 
}

4.2 简单减法

有两个正数 a,b,那么 a - b 有几种可能呢?

这样我们发现只需要实现 a>b 时的 a-b 就行了。

减法和加法一样,直接减就好了,最后像处理进位一样处理退位就行了。(考虑进位后每位最多 -10,所以最多只要借 1 位)

这样的时间复杂度也是 O(|\text{len}|) 的。

参考代码:

BIGINT operator-() const{ return BIGINT( *this, 1 ); }
BIGINT simple_minus( const BIGINT& b ) const{
    if( *this == b ) return BIGINT( 0 ); //相等 
    if( b > *this ) return -( b.simple_minus( *this ) ); //b > a
    //剩下 a > b 
    BIGINT c( *this );
    for( int i = 0; i < b.len; ++ i ){
        c.num[i] = c.num[i] - b.num[i]; //直接减 
    }
    for( int i = 0; i < len; ++ i ){
        if( c.num[i] < 0 ) c.num[i + 1] --, c.num[i] += 10; //退位 
    }
    while( c.len && !c.num[c.len - 1] ) -- c.len; //重新计算长度
    if( c.len == 0 ) ++ c.len;
    return c;
}

4.3 加法和减法的所有情况

先说加法 a+b

参考代码:

BIGINT operator+( const BIGINT& b ) const{cerr<<1;
    if( !f && !b.f ) return this->simple_add( b );
    if( f && b.f ) return -( ( -*this ).simple_add( -b ) );
    if( !f && b.f ) return simple_minus( -b );
    return b.simple_minus( -(*this) );
}

然后是减法 a-b,只要直接调用 a+(-b) 即可。

BIGINT operator-( const BIGINT& b )const{
    return *this + ( -b );
}

5 暴力高精乘

众所周知,小学的时候我们就学过竖式计算乘法:

       1  1   4 (a)
*      5  1   4 (b)
-----------------
       4  4 (16)
   1   1  4
 5 5 (20)
-----------------
 5 6 (25) 8 (16)
----------------- //进位
 5 8   5  9   6 (c)

观察这个式子,我们发现 a[i] * b[j] 最终会加到 c[i+j-1] 上,于是模拟即可。

注意到一位数乘一位数最多 81,在 |a| \leq 10^4 时不会超过 8.1 \times 10^5,加上进位不会超过 10^6,所以不会爆 int

乘法的符号问题比较简单,就留给读者思考了。

时间复杂度 O(\text{len}_a \cdot \text{len}_b)

参考实现:

BIGINT operator*( const BIGINT& b ) const{
    BIGINT c;
    c.len = 1;
    c.f = f ^ b.f; //符号 
    for( int i = 0; i < len + b.len + 10; ++ i ) c.num[i] = 0;
    for( int i = 0; i < len; ++ i ){
        for( int j = 0; j < b.len; ++ j ){
            c.num[i + j] += num[i] * b.num[j];
        }
    }
    c.len = len + b.len;
    for( int i = 0; i < c.len; ++ i ){ //进位 
        c.num[i + 1] += c.num[i] / 10;
        c.num[i] %= 10;
    }
    while( c.num[c.len - 1] >= 10 ){ //进位 
        c.num[c.len] += c.num[c.len - 1] / 10;
        c.num[c.len - 1] %= 10;
        ++ c.len;
    }
    while( c.len && !c.num[c.len - 1] ) -- c.len; //重新计算长度
    return c;
}

6 暴力高精除/取模

众所周知,小学的时候我们也学过竖式除法:

            4  5 //商
114 /------------
      5  1  4  0
      4  5  6
     ------------
         5  8  0
         5  7  0
     ------------
            1  0 //余数

实际上它是这样运作的:

 5140
-114  //省略 0
------
 4000
-114
------
 2860
 ...
------//一共四次,所以 c[1] = 4
 0580 //减不了 1140了
- 114 //那就减 114 吧
------
...
------
 0010 //减不了 114 了
      //实际上是余数

于是我们就可以得到结论:

由于每位最多减 9 次,所以时间复杂度 \Theta((\text{len}_a-\text{len}_b) \times \text{len}_b)=O(\text{len}_a^2)|a|\ge |b| 时,|a|<|b|O(\max\{\text{len}_a,\text{len}_b\}))。

所以我们就可以实现(细节有点多):

pair< BIGINT, BIGINT > divide( const BIGINT& b ) const{
    //余数返回和 C++ 取模一致,符号与 this 一致 
    if( b == 0 ) { //除 0 
        cerr << "BIGINT:Devide by zero!" << endl;
        exit( 32 );
    }
    if( abs_less( b ) ) return make_pair( BIGINT( 0 ), *this ); 
    BIGINT c( *this ), d; //c模 d商
    d.len = c.len - b.len + 1;
    //if( d.len <= 0 ) return make_pair( BIGINT( 0 ), *this );
    for( int i = 0; i <= d.len; ++ i ) d.num[i] = 0;
    for( int i = c.len - b.len; i >= 0; -- i ){
        while( 1 ){ //循环减直到 c < b * 10^i
            bool _great = 1; //c >= b * 10^i 判断 
            if( c.len <= i + b.len ) //注意多一位也是大 
                for( int j = b.len - 1; j >= 0; -- j ){ //c >= b * 10^i 判断
                    if( c.num[i + j] > b.num[j] ) break;
                    if( c.num[i + j] < b.num[j] ){
                        _great = 0;
                        break;
                    }
                }
            if( !_great ) break; //小于,退
            for( int j = b.len - 1; j >= 0; -- j ){ //减法 
                c.num[i + j] -= b.num[j];
            }
            for( int j = 0; j < b.len; ++ j ){ //退位 
                if( c.num[i + j] < 0 ){
                    c.num[i + j] += 10;
                    c.num[i + j + 1] --;
                }
            }
            while( c.len && !c.num[c.len - 1] ) -- c.len; //更新长度 
            d.num[i] ++;
        }
    }
    while( d.len && !d.num[d.len - 1] ) -- d.len; //更新长度 
    d.f = f ^ b.f;
    if( c.len == 0 && c.f == 1 ) c.f = 0; //注意 -0 
    if( c.len == 0 ) c.len = 1; //len=1 
    if( d.len == 0 && d.f == 1 ) d.f = 0;
    if( d.len == 0 ) d.len = 1;
    return make_pair( d, c );
}
BIGINT operator/( const BIGINT& b ) const{
    return divide( b ).first;
}
BIGINT operator%( const BIGINT& b ) const{
    BIGINT ret = divide( b ).second;
    if( b > BIGINT( 0 ) && ret < BIGINT( 0 ) ) ret = ret + b; //好心的变成正数 
    if( b < BIGINT( 0 ) && ret > BIGINT( 0 ) ) ret = ret + b;
    return ret;
}

:::success[你可以在这里检验目前为止的学习成果]{open}

P1932 A+B A-B A*B A/B A%B Problem

:::

:::info[这是到此为止的封装模板长度]

6.3 \text K
#include<iostream>
#include<string>
#include<algorithm>
#include<cstring>

using std::ostream;
using std::istream;
using std::string;
using std::pair;
using std::reverse;
using std::max;
using std::cerr;
using std::endl;
using std::make_pair;

struct BIGINT {
    int num[100005];        //存数字 
    bool f;                 //存符号 
    size_t len;             //存长度

    BIGINT() : f( 0 ), len( 1 ) {num[0] = 0;}
    BIGINT( const int& a ){ trans( a ); }
    BIGINT( const BIGINT& b, const bool changef = 0 ) : f( b.f ^ changef ), len( b.len ) { memcpy( num, b.num, sizeof( unsigned ) * b.len ); }

    inline size_t size(){ return len; }

    friend ostream& operator<<( ostream& out, const BIGINT& a ){
        return out << a.getS();
    }
    friend istream& operator>>( istream& in, BIGINT& a ){
        string s;
        in >> s;
        a.trans( s );
        return in;
    }

    void trans( int a ){
        f = 0, len = 0; //清空 
        if( a < 0 ) a = -a, f = 1; //负数 
        while( a ){
            num[len] = a % 10;  //当前位 
            a /= 10;            //这样子第 len 位就可以除 len 次 base 
            len ++;             //数字长度 +1 
        }
        if( !len ) len = 1;
    }
    int getI() const{
        int ret = 0;
        for( int i = len - 1; i >= 0; -- i ){
            ret = ret * 10;     //这样子在 i 位一共会乘上 i 次 base 
            ret = ret + num[i]; //统计当前位 
        }
        if( f ) ret = -ret;
        return ret;
    }
    void trans( string a ){
        f = len = 0; //清空 
        reverse( a.begin(), a.end() );  //反转
        if( a.back() == '-' ){ //负数
            a.pop_back();   //删掉符号 
            f = 1; 
        }
        for( int i = 0; i < a.size(); ++ i )
            num[i] = a[i] - '0'; //注意 - '0' 后才是对应数字 
        len = a.size();
        while( len > 0 && num[len - 1] == 0 ) -- len;
        if( !len ) len = 1;
    }
    string getS() const{
        string ret;
        for( int i = 0; i < len; ++ i )
            ret += char( num[i] + '0' ); //string 可以用加法进行拼接
        if( f ) ret += '-';
        reverse( ret.begin(), ret.end() ); 
        if( ret.size() == 0 ) ret = "0"; //特判 0 
        return ret;
    }

    bool abs_less( const BIGINT& b ) const{ //绝对值更小的
        if( len < b.len ) return 1;
        if( len > b.len ) return 0;
        for( int i = len - 1; i >= 0; -- i ){
            if( b.num[i] < num[i] ) return 0;
            if( num[i] < b.num[i] ) return 1;
        } 
        return 0;
    }
    bool operator==( const BIGINT& b ) const{ //相等的 
        if( f != b.f ) return 0;
        if( len != b.len ) return 0;
        for( int i = len - 1; i >= 0; -- i ){
            if( b.num[i] != num[i] ) return 0;
        } 
        return 1;
    }
    bool operator<( const BIGINT& b ) const{ //小于 
        if( f && !b.f ) return 1;
        if( !f && b.f ) return 0;
        if( f && b.f ) return !abs_less( b ) && !( *this == b );
        return abs_less( b );
    }
    bool operator>( const BIGINT& b ) const{ //大于 
        return !( *this < b ) && !( *this == b );
    }

    BIGINT simple_add( const BIGINT& b ) const{
        BIGINT c( *this );
        c.len = max( len, b.len );
        for( int i = len; i <= c.len; ++ i ) c.num[i]=0; //防止 c.num 没清空,多清一位是要进位 
        for( int i = 0; i < b.len; ++ i ){
            c.num[i] += b.num[i]; //加法 
        }
        for( int i = 0; i < c.len; ++ i ){ //进位 
            c.num[i + 1] += c.num[i] / 10;
            c.num[i] %= 10;
        }
        if( c.num[c.len] ) ++ c.len; //判断最高位
        return c; 
    }
    BIGINT operator-() const{ return BIGINT( *this, 1 ); }
    BIGINT simple_minus( const BIGINT& b ) const{
        if( *this == b ) return BIGINT( 0 ); //相等 
        if( b > *this ) return -( b.simple_minus( *this ) ); //b > a
        //剩下 a > b 
        BIGINT c( *this );
        for( int i = 0; i < b.len; ++ i ){
            c.num[i] = c.num[i] - b.num[i]; //直接减 
        }
        for( int i = 0; i < c.len; ++ i ){
            if( c.num[i] < 0 ) c.num[i + 1] --, c.num[i] += 10; //退位 
        }
        while( c.len && !c.num[c.len - 1] ) -- c.len; //重新计算长度 
        if( c.len == 0 ) ++ c.len;
        return c;
    }
    BIGINT operator+( const BIGINT& b ) const{
        if( !f && !b.f ) return this->simple_add( b );
        if( f && b.f ) return -( ( -*this ).simple_add( -b ) );
        if( !f && b.f ) return simple_minus( -b );
        return b.simple_minus( -(*this) );
    }
    BIGINT operator-( const BIGINT& b ) const{
        return *this + ( -b );
    }

    BIGINT operator*( const BIGINT& b ) const{
        BIGINT c;
        c.len = 1;
        c.f = f ^ b.f; //符号 
        for( int i = 0; i < len + b.len + 10; ++ i ) c.num[i] = 0;
        for( int i = 0; i < len; ++ i ){
            for( int j = 0; j < b.len; ++ j ){
                c.num[i + j] += num[i] * b.num[j];
            }
        }
        c.len = len + b.len;
        for( int i = 0; i < c.len; ++ i ){ //进位 
            c.num[i + 1] += c.num[i] / 10;
            c.num[i] %= 10;
        }
        while( c.num[c.len - 1] >= 10 ){ //进位 
            c.num[c.len] += c.num[c.len - 1] / 10;
            c.num[c.len - 1] %= 10;
            ++ c.len;
        }
        while( c.len && !c.num[c.len - 1] ) -- c.len; //重新计算长度
        return c;
    }

    pair< BIGINT, BIGINT > divide( const BIGINT& b ) const{
        //余数返回和 C++ 取模一致,符号与 this 一致 
        if( b == 0 ) { //除 0 
            cerr << "BIGINT:Devide by zero!" << endl;
            exit( 32 );
        }
        if( abs_less( b ) ) return make_pair( BIGINT( 0 ), *this ); 
        BIGINT c( *this ), d; //c模 d商
        d.len = c.len - b.len + 1;
        //if( d.len <= 0 ) return make_pair( BIGINT( 0 ), *this );
        for( int i = 0; i <= d.len; ++ i ) d.num[i] = 0;
        for( int i = c.len - b.len; i >= 0; -- i ){
            while( 1 ){ //循环减直到 c < b * 10^i
                bool _great = 1; //c >= b * 10^i 判断 
                if( c.len <= i + b.len ) //注意多一位也是大 
                    for( int j = b.len - 1; j >= 0; -- j ){ //c >= b * 10^i 判断
                        if( c.num[i + j] > b.num[j] ) break;
                        if( c.num[i + j] < b.num[j] ){
                            _great = 0;
                            break;
                        }
                    }
                if( !_great ) break; //小于,退
                for( int j = b.len - 1; j >= 0; -- j ){ //减法 
                    c.num[i + j] -= b.num[j];
                }
                for( int j = 0; j < b.len; ++ j ){ //退位 
                    if( c.num[i + j] < 0 ){
                        c.num[i + j] += 10;
                        c.num[i + j + 1] --;
                    }
                }
                while( c.len && !c.num[c.len - 1] ) -- c.len; //更新长度 
                d.num[i] ++;
            }
        }
        while( d.len && !d.num[d.len - 1] ) -- d.len; //更新长度 
        d.f = f ^ b.f;
        if( c.len == 0 && c.f == 1 ) c.f = 0; //注意 -0 
        if( c.len == 0 ) c.len = 1; //len=1 
        if( d.len == 0 && d.f == 1 ) d.f = 0;
        if( d.len == 0 ) d.len = 1;
        return make_pair( d, c );
    }
    BIGINT operator/( const BIGINT& b ) const{
        return divide( b ).first;
    }
    BIGINT operator%( const BIGINT& b ) const{
        BIGINT ret = divide( b ).second;
        if( b > BIGINT( 0 ) && ret < BIGINT( 0 ) ) ret = ret + b; //好心的变成正数 
        if( b < BIGINT( 0 ) && ret > BIGINT( 0 ) ) ret = ret + b;
        return ret;
    }
};

:::

:::warning[好东西就要来了!]{open}

接下来就要上强度喽! :::

7 FFT/NTT 优化高精乘

同样的,我们小学的时候就学过 FFT 和 NTT。(???

如果你不会,那么看看我的学习笔记。

所以说我们就可以打出一个 NTT 板子(这里是递归实现,效率很低;循环实现见我的学习笔记,那个效率比较高):

typedef long long ll;
const int mod=998244353,N=(1<<21+5),maxn=1<<21;

inline ll fstp(ll a,ll b){
    ll ans=1,x=a;
    while(b){
        if(b&1)ans=ans*x%mod;
        x=x*x%mod,b>>=1;
    }
    return ans;
}

void Nearly_Totally_TLE(int n,ll* a,int type){
    if(n==1)return;
    ll*a1=new ll[n>>1],*a2=new ll[n>>1];
    for(int i=0;i<n;++i)if(i&1)a2[i/2]=a[i];else a1[i/2]=a[i];
    Nearly_Totally_TLE(n>>1,a1,type),Nearly_Totally_TLE(n>>1,a2,type);
    ll wmul=((type+1)?fstp(3,(mod-1)/n):fstp(fstp(3,mod-2),(mod-1)/n)),w=1;
    for(int i=0;i<(n>>1);++i)a[i]=(a1[i]+w*a2[i])%mod,a[i+(n>>1)]=((a1[i]-w*a2[i])%mod+mod)%mod,w=w*wmul%mod;
    delete[] a1;
    delete[] a2;
}

然后你就需要修改一下你的 operator*

BIGINT operator*( const BIGINT& b ) const{
    //计算长度 
    int n = 1;
    while( n < len + b.len ) n *= 2;

    //新建数组 
    ll* a = new ll[n], *c = new ll[n];
    for( int i = 0; i < n; ++ i ) a[i] = c[i] = 0; //初始化 

    //乘法 3 次 NTT 
    for( int i = 0; i < len; ++ i ) a[i] = num[i];
    for( int i = 0; i < b.len; ++ i ) c[i] = b.num[i];
    Nearly_Totally_TLE( n, a, 1 );
    Nearly_Totally_TLE( n, c, 1 );
    for( int i = 0; i < n; ++ i ) a[i] = a[i] * c[i] % mod;
    Nearly_Totally_TLE( n, a, -1 );
    for( int i = 0; i < n; ++ i ) a[i] = a[i] * fstp( n, mod - 2 ) % mod; //记得要除掉 n!

    //处理进位 
    for( int i = 0; i < n; ++ i ) a[i + 1] += a[i] / 10, a[i] %= 10;

    //处理答案 
    BIGINT d;
    d.len = n;
    d.f = f ^ b.f;
    for(int i = 0; i < n; ++ i ) d.num[i] = a[i];
    while( d.len > 1 && !d.num[d.len - 1] ) -- d.len; //重新计算长度 

    //删除数组(不然可能内存泄漏) 
    delete[] a;
    delete[] c;
    return d;
}

看过我那篇文章的人应该都看得懂吧

它的时间复杂度是 O(n\log n),分析在我的学习笔记中。

8 牛顿迭代法优化高精度除法

众所周知……好了这真不能是小学知识了。

我们求 \lfloor\frac ab\rfloor 可以使用牛顿迭代法求出 \frac 1b 的近似值,只要精度足够那么乘上 a 取整就是 \lfloor\frac ab\rfloor

牛顿迭代的式子是 x_1 = 2x_0-bx_0^2。关于这个式子的推导详见我的题解。这里主要讲如何实现。

我们要实现一个简单的高精度小数,并实现它的乘法。

8.1 如何用高精度整数实现一个简略的高精度小数

我们可以直接存储个位数的位置。具体地,可以定义一个 BIGFLOAT 结构体:

struct BIGFLOAT{
    BIGINT a;    //整数当小数存 
    ll zero_pos; //0 的位置 
};

定义一些构造函数:

    BIGFLOAT() : a( BIGINT( 0 ) ), zero_pos( 0 ){}
    BIGFLOAT( const BIGINT& b, const ll zero_p = 0 ) : a( b ), zero_pos( zero_p ){}
    BIGFLOAT( const BIGFLOAT& b ) : a( b.a ), zero_pos( b.zero_pos ){}

接下来,我们要定义左移右移函数(\times 10^x\times \frac1{10^x}):

    BIGFLOAT operator<<( const ll x ) const{ // *10^x
        if( a == 0 ) return BIGFLOAT( BIGINT( 0 ) );
        if( zero_pos >= x ) return BIGFLOAT( a, zero_pos - x ); //可以直接右移小数点 
        ll mov = x - zero_pos;
        BIGFLOAT b( BIGINT( 0 ) );
        b.a.len = a.len + mov;
        for( int i = b.a.len - 1; i >= 0; -- i ) b.a.num[i] = 0; //注意清空 
        for( int i = a.len - 1; i >= 0; -- i ) b.a.num[i + mov] = a.num[i]; //移位 
        b.a.f = a.f;
        return b;
    }
    BIGFLOAT operator>>( const ll x ) const{ // /10^x
        if( a == 0 ) return BIGFLOAT( BIGINT( 0 ) );
        ll mov = 0;
        while( mov < a.len && !a.num[mov] ) ++mov; //后缀 0 
        if( mov == a.len ) return BIGFLOAT( BIGINT( 0 ) ); //全 0 
        BIGFLOAT b( 0 );
        b.a.len = a.len - mov; //长度 
        for( int i = b.a.len - 1; i >= 0; -- i ) 
            b.a.num[i] = a.num[i + mov]; //右移 
        b.zero_pos = zero_pos + ( x - mov ); //计算个位数位置 
        return b;
    }

接下来是下取整函数:

    BIGINT floor() const{
        if( zero_pos >= a.len ) return BIGINT( 0 ); //小数点比整数还前面 
        BIGINT b( 0 );
        b.len = a.len - zero_pos;
        for( int i = b.len - 1; i >= 0; -- i ) b.num[i] = 0;
        for( int i = b.len - 1; i >= 0; -- i ) b.num[i] = a.num[i + zero_pos]; //同移位 
        b.f = a.f;
        return b;
    }

接下来我们就可以实现加减乘法了(注意加法还需要对齐):

    BIGFLOAT lft_mov( const ll x ) const{ //数不变,精度增加 x 位,用于对齐 
        if( a == 0 ) return BIGFLOAT( BIGINT( 0 ) ); //特判 0 
        BIGFLOAT b( BIGINT( 0 ) );
        b.a.len = a.len + x;
        for( int i = b.a.len - 1; i >= 0; -- i ) b.a.num[i] = 0; //清空 
        for( int i = a.len - 1; i >= 0; -- i ) b.a.num[i + x] = a.num[i]; //左移 
        b.zero_pos = x + zero_pos; //小数点一起移动 
        b.a.f = a.f;
        return b;
    }
    BIGFLOAT operator+( const BIGFLOAT& b ) const{ //加法
        BIGFLOAT c( BIGINT( 0 ) );
        c.zero_pos = max( zero_pos, b.zero_pos ); //小数点位置取最大,再对齐 
        ll m1 = c.zero_pos - zero_pos, m2 = c.zero_pos - b.zero_pos; //不定义这俩直接传参会 UB 
        c.a = lft_mov( m1 ).a + b.lft_mov( m2 ).a; //位对齐 
        return c;
    }
    BIGFLOAT operator-() const{ return BIGFLOAT( -a, zero_pos ); } //取反 
    BIGFLOAT operator-( const BIGFLOAT& b ) const{ //减法 
        return ( (*this) + ( -b ) ); //偷懒 
    }
    BIGFLOAT operator*( const BIGFLOAT& b ) const{ //乘法 
        BIGFLOAT c;
        c.a = a * b.a; //直接乘 
        c.zero_pos = zero_pos + b.zero_pos; //小数点移动 
        return c;
    }

8.2 牛顿迭代,启动

接下来就可以在结构体外写我们的牛顿迭代了。

注意在牛顿迭代的时候要限制精度,不然会爆内存。

牛顿迭代初始值需要是一个极小值。

最后记得也加一个极小值,不然有误差。

BIGFLOAT calc1_b( int n, const BIGFLOAT& b ){
    auto x0 = BIGFLOAT( BIGINT( int( 1 ) ) ) >> b.a.len; //初始值不能取太大,10^(-len) 
    while( n-- ){
        //迭代 
        x0 = BIGFLOAT( BIGINT( 2 ) ) * x0 - b * x0 * x0;

        //限制精度,精度太高会爆内存 
        x0 = x0 << b.a.len * 2;
        x0 = BIGFLOAT( x0.floor() );
        x0 = x0 >> b.a.len * 2;
    }
    return x0;
}
BIGINT devide_helper( const BIGINT& a, const BIGINT& b ){
    auto x = calc1_b( 20 , BIGFLOAT( b ) ); //得到迭代结果 
    x = x + ( BIGFLOAT( BIGINT( 1 ) ) >> ( b.len * 2 ) ); //加上极小值避免精度误差 
    return ( BIGFLOAT( a ) * x ).floor(); //返回相乘后下取整 
}

最后再把我们的除法和取模更新一下。

    BIGINT operator/( const BIGINT& b ) const{
        BIGINT c( 0 );
        c = devide_helper( abs(), b.abs() );
        c.f = f ^ b.f;
        return c;
    }
    BIGINT operator%( const BIGINT& b ) const{
        return *this - *this / b * b; //朴素取模 
    }

我才不会告诉你我上面写的高精度除法模板过不了模板题呢如果你想过高精度除法模板题,你需要巨量的常数优化。

:::info[高精度四则运算模板]

```cpp #include<bits/stdc++.h> using namespace std; typedef long long ll; typedef long long ll; const int mod=998244353,maxn=1<<21; inline ll fstp(ll a,ll b){ ll ans=1,x=a; while(b){ if(b&1)ans=ans*x%mod; x=x*x%mod,b>>=1; } return ans; } void Nearly_Totally_TLE(int n,ll* a,int type){ if(n==1)return; ll*a1=new ll[n>>1],*a2=new ll[n>>1]; for(int i=0;i<n;++i)if(i&1)a2[i/2]=a[i];else a1[i/2]=a[i]; Nearly_Totally_TLE(n>>1,a1,type),Nearly_Totally_TLE(n>>1,a2,type); ll wmul=((type+1)?fstp(3,(mod-1)/n):fstp(fstp(3,mod-2),(mod-1)/n)),w=1; for(int i=0;i<(n>>1);++i)a[i]=(a1[i]+w*a2[i])%mod,a[i+(n>>1)]=((a1[i]-w*a2[i])%mod+mod)%mod,w=w*wmul%mod; delete[] a1; delete[] a2; } struct BIGFLOAT; struct BIGINT; BIGFLOAT calc1_b( int n, const BIGFLOAT& b ); BIGINT devide_helper( const BIGINT& a, const BIGINT& b ); struct BIGINT { int num[100005]; //存数字 bool f; //存符号 ll len; //存长度 BIGINT() : f( 0 ), len( 1 ) { num[0] = 0; } BIGINT( const int& a ){ trans( a ); } BIGINT( const BIGINT& b, const bool changef = 0 ) : f( b.f ^ changef ), len( b.len ) { memcpy( num, b.num, sizeof( unsigned ) * b.len ); } inline size_t size(){ return len; } friend ostream& operator<<( ostream& out, const BIGINT& a ){ return out << a.getS(); } friend istream& operator>>( istream& in, BIGINT& a ){ string s; in >> s; a.trans( s ); return in; } void trans( int a ){ num[0] = 0; f = 0, len = 0; //清空 if( a < 0 ) a = -a, f = 1; //负数 while( a ){ num[len] = a % 10; //当前位 a /= 10; //这样子第 len 位就可以除 len 次 base len ++; //数字长度 +1 } if( !len ) len = 1; } int getI() const{ int ret = 0; for( int i = len - 1; i >= 0; -- i ){ ret = ret * 10; //这样子在 i 位一共会乘上 i 次 base ret = ret + num[i]; //统计当前位 } if( f ) ret = -ret; return ret; } void trans( string a ){ num[0] = 0; f = len = 0; //清空 reverse( a.begin(), a.end() ); //反转 if( a.back() == '-' ){ //负数 a.pop_back(); //删掉符号 f = 1; } for( int i = 0; i < a.size(); ++ i ) num[i] = a[i] - '0'; //注意 - '0' 后才是对应数字 len = a.size(); while( len > 0 && num[len - 1] == 0 ) -- len; if( !len ) len = 1; } string getS() const{ string ret; for( int i = 0; i < len; ++ i ) ret += char( num[i] + '0' ); //string 可以用加法进行拼接 if( f ) ret += '-'; reverse( ret.begin(), ret.end() ); if( ret.size() == 0 ) ret = "0"; //特判 0 return ret; } bool abs_less( const BIGINT& b ) const{ //绝对值更小的 if( len < b.len ) return 1; if( len > b.len ) return 0; for( int i = len - 1; i >= 0; -- i ){ if( b.num[i] < num[i] ) return 0; if( num[i] < b.num[i] ) return 1; } return 0; } bool operator==( const BIGINT& b ) const{ //相等的 if( f != b.f ) return 0; if( len != b.len ) return 0; for( int i = len - 1; i >= 0; -- i ){ if( b.num[i] != num[i] ) return 0; } return 1; } bool operator<( const BIGINT& b ) const{ //小于 if( f && !b.f ) return 1; if( !f && b.f ) return 0; if( f && b.f ) return !abs_less( b ) && !( *this == b ); return abs_less( b ); } bool operator>( const BIGINT& b ) const{ //大于 return !( *this < b ) && !( *this == b ); } BIGINT simple_add( const BIGINT& b ) const{ BIGINT c( *this ); c.len = max( len, b.len ); for( int i = len; i <= c.len; ++ i ) c.num[i]=0; //防止 c.num 没清空,多清一位是要进位 for( int i = 0; i < b.len; ++ i ){ c.num[i] += b.num[i]; //加法 } for( int i = 0; i < c.len; ++ i ){ //进位 c.num[i + 1] += c.num[i] / 10; c.num[i] %= 10; } if( c.num[c.len] ) ++ c.len; //判断最高位 return c; } BIGINT operator-() const{ return BIGINT( *this, 1 ); } BIGINT simple_minus( const BIGINT& b ) const{ if( *this == b ) return BIGINT( 0 ); //相等 if( b > *this ) return -( b.simple_minus( *this ) ); //b > a //剩下 a > b BIGINT c( *this ); for( int i = 0; i < b.len; ++ i ){ c.num[i] = c.num[i] - b.num[i]; //直接减 } for( int i = 0; i < c.len; ++ i ){ if( c.num[i] < 0 ) c.num[i + 1] --, c.num[i] += 10; //退位 } while( c.len && !c.num[c.len - 1] ) -- c.len; //重新计算长度 if( c.len == 0 ) ++ c.len; return c; } BIGINT operator+( const BIGINT& b ) const{ if( !f && !b.f ) return this->simple_add( b ); if( f && b.f ) return -( ( -*this ).simple_add( -b ) ); if( !f && b.f ) return simple_minus( -b ); return b.simple_minus( -(*this) ); } BIGINT operator-( const BIGINT& b ) const{ return *this + ( -b ); } BIGINT operator*( const BIGINT& b ) const{ //计算长度 int n = 1; while( n < len + b.len ) n *= 2; //新建数组 ll* a = new ll[n], *c = new ll[n]; for( int i = 0; i < n; ++ i ) a[i] = c[i] = 0; //初始化 //乘法 3 次 NTT for( int i = 0; i < len; ++ i ) a[i] = num[i]; for( int i = 0; i < b.len; ++ i ) c[i] = b.num[i]; Nearly_Totally_TLE( n, a, 1 ); Nearly_Totally_TLE( n, c, 1 ); for( int i = 0; i < n; ++ i ) a[i] = a[i] * c[i] % mod; Nearly_Totally_TLE( n, a, -1 ); for( int i = 0; i < n; ++ i ) a[i] = a[i] * fstp( n, mod - 2 ) % mod; //记得要除掉 n! //处理进位 for( int i = 0; i < n; ++ i ) a[i + 1] += a[i] / 10, a[i] %= 10; //处理答案 BIGINT d; d.len = n; d.f = f ^ b.f; for(int i = 0; i < n; ++ i ) d.num[i] = a[i]; while( d.len > 1 && !d.num[d.len - 1] ) -- d.len; //重新计算长度 //删除数组(不然可能内存泄漏) delete[] a; delete[] c; return d; } BIGINT abs() const{ if( f ) return -*this; return *this; } BIGINT operator/( const BIGINT& b ) const{ BIGINT c( 0 ); c = devide_helper( abs(), b.abs() ); c.f = f ^ b.f; return c; } BIGINT operator%( const BIGINT& b ) const{ return *this - *this / b * b; //朴素取模 } }; struct BIGFLOAT{ BIGINT a; ll zero_pos; //0 的位置 BIGFLOAT() : a( BIGINT( 0 ) ), zero_pos( 0 ){} BIGFLOAT( const BIGINT& b, const ll zero_p = 0 ) : a( b ), zero_pos( zero_p ){} BIGFLOAT( const BIGFLOAT& b ) : a( b.a ), zero_pos( b.zero_pos ){} BIGFLOAT operator<<( const ll x ) const{ // *10^x if( a == 0 ) return BIGFLOAT( BIGINT( 0 ) ); if( zero_pos >= x ) return BIGFLOAT( a, zero_pos - x ); //可以直接右移小数点 ll mov = x - zero_pos; BIGFLOAT b( BIGINT( 0 ) ); b.a.len = a.len + mov; for( int i = b.a.len - 1; i >= 0; -- i ) b.a.num[i] = 0; //注意清空 for( int i = a.len - 1; i >= 0; -- i ) b.a.num[i + mov] = a.num[i]; //移位 b.a.f = a.f; return b; } BIGFLOAT operator>>( const ll x ) const{ // /10^x if( a == 0 ) return BIGFLOAT( BIGINT( 0 ) ); ll mov = 0; while( mov < a.len && !a.num[mov] ) ++mov; //后缀 0 if( mov == a.len ) return BIGFLOAT( BIGINT( 0 ) ); //全 0 BIGFLOAT b( 0 ); b.a.len = a.len - mov; //长度 for( int i = b.a.len - 1; i >= 0; -- i ) b.a.num[i] = a.num[i + mov]; //右移 b.zero_pos = zero_pos + ( x - mov ); //计算个位数位置 return b; } BIGINT floor() const{ if( zero_pos >= a.len ) return BIGINT( 0 ); //小数点比整数还前面 BIGINT b( 0 ); b.len = a.len - zero_pos; for( int i = b.len - 1; i >= 0; -- i ) b.num[i] = 0; for( int i = b.len - 1; i >= 0; -- i ) b.num[i] = a.num[i + zero_pos]; //同移位 b.f = a.f; return b; } BIGFLOAT lft_mov( const ll x ) const{ //数不变,精度增加 x 位,用于对齐 if( a == BIGINT( 0 ) ) return BIGFLOAT( BIGINT( 0 ) ); //特判 0 BIGFLOAT b( BIGINT( 0 ) ); b.a.len = a.len + x; for( int i = b.a.len - 1; i >= 0; -- i ) b.a.num[i] = 0; //清空 for( int i = a.len - 1; i >= 0; -- i ) b.a.num[i + x] = a.num[i]; //左移 b.zero_pos = x + zero_pos; //小数点一起移动 b.a.f = a.f; return b; } BIGFLOAT operator+( const BIGFLOAT& b ) const{ //加法 BIGFLOAT c( BIGINT( 0 ) ); c.zero_pos = max( zero_pos, b.zero_pos ); //小数点位置取最大,再对齐 ll m1 = c.zero_pos - zero_pos, m2 = c.zero_pos - b.zero_pos; //不定义这俩直接传参会 UB c.a = lft_mov( m1 ).a + (b.lft_mov( m2 )).a; //位对齐 return c; } BIGFLOAT operator-() const{ return BIGFLOAT( -a, zero_pos ); } //取反 BIGFLOAT operator-( const BIGFLOAT& b ) const{ //减法 return ( (*this) + ( -b ) ); //偷懒 } BIGFLOAT operator*( const BIGFLOAT& b ) const{ //乘法 BIGFLOAT c; c.a = a * b.a; //直接乘 c.zero_pos = zero_pos + b.zero_pos; //小数点移动 auto zp = c.zero_pos; return c; } friend ostream& operator<<( ostream& out, const BIGFLOAT& b ){ if( b.a.f ) out << '-'; if( b.zero_pos >= b.a.len ) out << '.'; for( ll i = b.zero_pos - 1; i >= b.a.len && i >= 0; -- i ) out <<'0'; for( ll i = b.a.len - 1; i >= 0; -- i ){ out << b.a.num[i]; if( i == b.zero_pos ) out << '.'; } return out; } }; BIGFLOAT calc1_b( int n, const BIGFLOAT& b ){ auto x0 = BIGFLOAT( BIGINT( int( 1 ) ) ) >> b.a.len; //初始值不能取太大,10^(-len) while( n-- ){ //迭代 x0 = BIGFLOAT( BIGINT( 2 ) ) * x0 - b * x0 * x0; //限制精度,精度太高会爆内存 x0 = x0 << b.a.len * 2; x0 = BIGFLOAT( x0.floor() ); x0 = x0 >> b.a.len * 2; } return x0; } BIGINT devide_helper( const BIGINT& a, const BIGINT& b ){ auto x = calc1_b( 20 , BIGFLOAT( b ) ); //得到迭代结果 x = x + ( BIGFLOAT( BIGINT( 1 ) ) >> ( b.len * 2 ) ); //加上极小值避免精度误差 return ( BIGFLOAT( a ) * x ).floor(); //返回相乘后下取整 } ``` ::: ## 9 尾声 现在你已经学会了高精度四则运算。你可以轻松地通过需要它们的题目了! 你也可以自行完成下面的挑战: 1. 实现高精开根。 :::info[思路] 还是牛顿迭代。 已知 $f(x) = x^2 - a$,那么 $f(x)=0$ 时 $x = \sqrt a$。 用牛顿迭代法推方程可得: $$k=f'(x_0)=2x_0\\ x_0^2-a=2x_0^2+b \Rightarrow b = -(x_0^2+a)\\ y=2x_0x-(x_0^2+a)\\ y=0\Rightarrow x_1=\frac {x_0}{2} + \frac{a}{2x_0}$$ 然后直接迭代即可。 代码请读者自行实现。 ::: 2. 优化高精度模板常数。 :::info[思路] - 压位高精度; - 使用动态数组; - 使用更大的 NTT 模数/使用 FFT,使用循环形式; - 根据数据大小动态调整迭代次数; - 优化精度限制的实现; - 在数据较小的时候使用暴力; - $\dots
:::
  1. 实现高精度快速幂、\exp\ln\gcd 等常用函数。

    注意到高精度除法和取模非常慢,所以你的 \gcd 可能得用 更相减损术。

  2. 实现能表示近乎无限大数(但是精度有限)的高精度实数模板。

    我的思路是实现一个类似浮点数的,前面存真实数据、后面存指数的浮点数模板,这两个的数据类型都可以自定义(也就是说假设是 Decimal<class T,class T2> 那么可以 Decimal<BIGINT,Decimal<BIGINT,...>>)。

  3. 封装高精度模板以随时方便地使用。

感谢你的阅读。如果这篇文章对你有帮助,那么你可以点个赞。如果有什么疑问可以在评论区提问。

不要脸地求关注.jpg