针对 P1001 A+B Problem 的大炮使用手册

· · 题解

大炮维护日志:修改了大炮的运行逻辑,特别优化了其在对蚊子实施打击时的使用体验。同时,同步更新了《大炮使用手册》,润色了一些句子,并增添了对大炮新功能的说明。

思路:

观察数据范围,发现答案竟然在 [-2^{31},2^{31}-1] 范围内,为了防止溢出,我们考虑使用高精度算法,然而朴素的高精度算法时间复杂度为 O(n),也许需要一定程度的优化。
我们知道,FFT(快速傅里叶变换)可以高效地计算多项式卷积,因此被用来优化乘法的计算。借用这种算法,我们尝试对加法进行优化。具体思路借用了这道题。
FFT 具有线性性质。这也是我们应用其优化运算的基础。具体地,对于两个序列 ab,有 \operatorname{FFT}(a\times b)=\operatorname{FFT}(a)\times \operatorname{FFT}(b),同理,有 \operatorname{FFT}(a+ b)=\operatorname{FFT}(a)+\operatorname{FFT}(b),考虑利用这个性质,对 ab 进行 FFT,然后将结果相加,最后再对相加后的结果逆 FFT 即得答案。

常数优化:

高精度运算的一种简单且有效的优化方式就是压位,然而 FFT 在计算过程中涉及浮点运算,这会带来严重的精度问题,尤其是在使用压位高精的时候。所以我们只压两位。对于本题,\dfrac{1}{2} 的常数已经足够了。

使用体验优化:

为了方便大炮的使用,我对其进行了封装,同时内附一份标准的 FFT 优化高精度乘法运算模板。

#include<cmath>
#include<iostream>
#include<cstring>
//********一些辅助函数与FFT实现
const double PI=4*atan(1);
template<typename T>
void Swap(T &a,T &b){
    T c=a;
    a=b;
    b=c;
    return;
}
template<typename T>
T Max(const T &a,const T &b){
    return a<b?b:a;
}
typedef long long ll;
struct comp{
    double real,imag;
    comp operator+(const comp &x)const{
        return {real+x.real,imag+x.imag};
    }
    comp operator-(const comp &x)const{
        return {real-x.real,imag-x.imag};
    }
    comp operator*(const comp &x)const{
        return {real*x.real-imag*x.imag,real*x.imag+x.real*imag};
    }
    comp operator/(const unsigned &x)const{
        return {real/(double)x,imag/(double)x};
    }
};
void FFT(comp *f,unsigned n,int rev){
    for(unsigned i=1,j=n>>1,k;i<n-1;i++){//位逆序置换
        if(i<j)
            Swap(f[i],f[j]);
        k=n>>1;
        while(j>=k){
            j-=k;
            k>>=1;
        }
        j+=k;
    }
    for(unsigned l=2;l<=n;l<<=1){//蝶形运算
        double arg=2*PI*rev/l;
        comp wn={cos(arg),sin(arg)};
        for(unsigned i=0;i<n;i+=l){
            comp w={1,0};
            for(unsigned j=0;j<(l>>1);j++){
                comp f1=f[i+j];
                comp f2=f[i+j+(l>>1)];
                f[i+j]=f1+w*f2;
                f[i+j+(l>>1)]=f1-w*f2;
                w=w*wn;
            }
        }
    }
    if(!~rev)
        for(unsigned i=0;i<n;i++)
            f[i]=f[i]/n;
}
//********高精度整数类
#define BASE 100//因为FFT的精度问题严重,我们只压2位
template<const unsigned Size>
class bigint{
private:
    unsigned len;
    int num[Size];
    void init(){
        memset(num,0,sizeof(num));
        len=1;
    }
    bool abs_greater_equal(const bigint &a)const{
        if(len!=a.len)
            return len>a.len;
        for(int i=len;i;i--)
            if(num[i]!=a.num[i])
                return num[i]>a.num[i];
        return 1;
    }
public:
    bigint(){
        init();
    }
    void get_num(std::string s){
        init();
        int f=0;
        unsigned slen=s.length();
        if(s[0]=='-')
            num[0]=f=1;
        len=0;
        unsigned temp=0,w=1;
        for(int i=slen-1;i>=f;i--){
            temp+=(s[i]^48)*w;
            w=(w<<1)+(w<<3);
            if(w==BASE||i==f){
                num[++len]=(int)temp;
                temp=0;
                w=1;
            }
        }
        if(temp||len==0)
            num[++len]=temp;
    }
    bool operator<(const bigint &a)const{
        if(num[0]&&!a.num[0])
            return 1;
        if(!num[0]&&a.num[0])
            return 0;
        if(num[0]){
            if(len!=a.len)
                return len>a.len;
            for(int i=len;i;i--)
                if(num[i]!=a.num[i])
                    return num[i]>a.num[i];
        }
        else{
            if(len!=a.len)
                return len<a.len;
            for(int i=len;i;i--)
                if(num[i]!=a.num[i])
                    return num[i]<a.num[i];
        }
        return 0;
    }
    bigint operator+(const bigint &a)const{
        bigint res;
        if(len==1&&num[1]==0){
            res=a;
            return res;
        }
        if(a.len==1&&a.num[1]==0){
            res=*this;
            return res;
        }
        if(num[0]==a.num[0]){
            res.num[0]=num[0];
            unsigned len_sum=1;
            while(len_sum<len+a.len)
                len_sum<<=1;
            comp *fa=new comp[len_sum]();
            comp *fb=new comp[len_sum]();
            for(unsigned i=0;i<len;i++)
                fa[i]={(double)num[i+1],0};
            for(unsigned i=0;i<a.len;i++)
                fb[i]={(double)a.num[i+1],0};
            FFT(fa,len_sum,1);
            FFT(fb,len_sum,1);
            for(unsigned i=0;i<len_sum;i++)
                fa[i]=fa[i]+fb[i];
            FFT(fa,len_sum,-1);
            res.len=Max(len,a.len);
            ll temp=0;
            for(unsigned i=0;i<res.len;i++){
                ll val=(ll)round(fa[i].real)+temp;
                res.num[i+1]=(int)(val%BASE);
                temp=val/BASE;
            }
            if(temp)
                res.num[++res.len]=temp;
            while(res.len>1&&res.num[res.len]==0)
                res.len--;
            delete[] fa;
            delete[] fb;
        }
        else{
            if(abs_greater_equal(a)){
                res.num[0]=num[0];
                unsigned len_sum=1;
                while(len_sum<len+a.len)
                    len_sum<<=1;
                comp *fa=new comp[len_sum]();
                comp *fb=new comp[len_sum]();
                for(unsigned i=0;i<len;i++)
                    fa[i]={(double)num[i+1],0};
                for(unsigned i=0;i<a.len;i++)
                    fb[i]={(double)a.num[i+1],0};
                FFT(fa,len_sum,1);
                FFT(fb,len_sum,1);
                for(unsigned i=0;i<len_sum;i++)
                    fa[i]=fa[i]-fb[i];
                FFT(fa,len_sum,-1);
                res.len=Max(len,a.len);
                ll temp=0;
                for(unsigned i=0;i<res.len;i++){
                    ll val=(ll)round(fa[i].real)+temp;
                    if(val<0){
                        val+=BASE;
                        temp=-1;
                    }
                    else
                        temp=0;//借位
                    res.num[i+1]=(int)(val%BASE);
                }
                if(temp)
                    res.num[++res.len]=temp;
                while(res.len>1&&res.num[res.len]==0)
                    res.len--;
                delete[] fa;
                delete[] fb;
            }
            else{
                res.num[0]=a.num[0];
                unsigned len_sum=1;
                while(len_sum<len+a.len)
                    len_sum<<=1;
                comp *fa=new comp[len_sum]();
                comp *fb=new comp[len_sum]();
                for(unsigned i=0;i<len;i++)
                    fa[i]={(double)num[i+1],0};
                for(unsigned i=0;i<a.len;i++)
                    fb[i]={(double)a.num[i+1],0};
                FFT(fa,len_sum,1);
                FFT(fb,len_sum,1);
                for(unsigned i=0;i<len_sum;i++)
                    fa[i]=fb[i]-fa[i];
                FFT(fa,len_sum,-1);
                res.len=Max(len,a.len);
                ll temp=0;
                for(unsigned i=0;i<res.len;i++){
                    ll val=(ll)round(fa[i].real)+temp;
                    if(val<0){
                        val+=BASE;
                        temp=-1;
                    }
                    else
                        temp=0;
                    res.num[i+1]=(int)(val%BASE);
                }
                if(temp)
                    res.num[++res.len]=temp;
                while(res.len>1&&res.num[res.len]==0)
                    res.len--;
                delete[] fa;
                delete[] fb;
            }
            if(res.len==1&&res.num[1]==0)
                res.num[0]=0;
        }
        return res;
    }
    bigint operator*(const bigint &a)const{
        bigint res;
        if((len==1&&num[1]==0)||(a.len==1&&a.num[1]==0))
            return res;
        res.num[0]=num[0]^a.num[0];
        unsigned len_sum=1;
        while(len_sum<len+a.len)
            len_sum<<=1;
        comp *fa=new comp[len_sum]();
        comp *fb=new comp[len_sum]();
        for(unsigned i=0;i<len;i++)
            fa[i]={(double)num[i+1],0};
        for(unsigned i=0;i<a.len;i++)
            fb[i]={(double)a.num[i+1],0};
        FFT(fa,len_sum,1);
        FFT(fb,len_sum,1);
        for(unsigned i=0;i<len_sum;i++)
            fa[i]=fa[i]*fb[i];
        FFT(fa,len_sum,-1);
        res.len=len+a.len;
        ll temp=0;
        for(unsigned i=0;i<res.len;i++){
            ll val=(ll)(fa[i].real+0.5)+temp;
            res.num[i+1]=(int)(val%BASE);
            temp=val/BASE;
        }
        if(temp)
            res.num[++res.len]=temp;
        while(res.len>1&&res.num[res.len]==0)
            res.len--;
        delete[] fa;
        delete[] fb;
        return res;
    }
    void read(){
        init();
        std::string s;
        char ch=getchar();
        while(ch<'0'||ch>'9'){
            if(ch=='-')
                s.push_back('-');
            ch=getchar();
        }
        while(ch>='0'&&ch<='9'){
            s.push_back(ch);
            ch=getchar();
        }
        get_num(s);
    }
    void print(){
        if(num[0])
            putchar('-');
        bool leading_zero=1;
        for(int i=len;i;i--){
            if(leading_zero)
                printf("%d",num[i]);
            else
                printf("%02d",num[i]);
            leading_zero=0;
        }
        putchar('\n');
        return;
    }
};
//********程序主体
const int N=1<<10,M=100;
int n,m;
bigint<114514> a,b,c;
int main(){
    a.read();
    b.read();
    c=a+b;
    c.print();
    return 0;
}

写在结尾:

我们成功把 O(1) 的朴素算法优化到了 O(n),接着进一步将 O(n) 的高精度算法优化到了 FFT 的 O(n \log n),这体现了我们对 OI 知识的灵活迁移运用。希望大家能(可千万别)将这种思想运用到学习生活中。