浅谈 01Exp 的实现原理

引言

朴素的 Ln Exp 在模数非常小（小于长度）的情况下会比较拟人，于是我们需要扩域，于是经过我一晚上+一早上的研究现在至少得到了 \(O(n\log nF(\log^2n))\) 的保证正确的 Exp，\(F(x)\) 表示我们要保持 \(x\) 位精度所要付出的代价，这个是通过牛迭得到的精度上界。

如果使用朴素的递推，似乎需要 \(O(p^n)\) 的模数级别的精度，但是牛迭则说明我们至多需要 \(O(p^{\log^2 n})\) 的精度，因为 Exp牛迭 的式子里面有 \(\ln\) ，\(\ln\) 每次的精度指数会下降 \(O(\log n)\) ，我们需要 \(O(\log n)\) 轮，于是需要的初始精度指数为 \(O(\log^2 n)\) 的。

现在的问题来到关于这个精度，究竟最低限度需要多少？实践对拍表明，似乎只需要 \(O(\log n)\) 位！

猜想：

关于我们究竟需要多少精度，我们发现当

1）Ln/Exp 的实际值域太小？
这个猜想被否定了，因为经过最基本的观察我们发现，Ln/Inv/Exp 真值都是指数级的爆炸。

2）卷的多项式项数太少了？
这个猜想被否定了，因为项数多和项数少似乎没有本质上的区别，基本判定都是 需要值域*长度<保留值域。具体的，如果我们过程保留的是 \(p^k\) ，我们需要 \(p^t\) 的，次数是 \(p^l\) 的(不足做下取整，也就是最多的 \(p\) 因子的指数)，我们需要 \(t+l\le k\) ，这个显然是下界，但是很奇怪的是同时为上界，非常紧。

3）卷的多项式个数太少了？
这个猜想被否定了，我尝试将100个长度为2000的多项式的Ln形式加起来，然后Exp和直接乘积做比较，对拍得到没有差别。

4）卷的模数有特殊性？
这个猜想被否定了，我们模 \(p^k\) 意义下的 \(p\) 我们尝试了不同的，都过拍了。

关键代码保存一下：

模二意义下的对拍：

点击查看代码

#include <bits/stdc++.h>
#define m_p(a,b) make_pair(a,b)
#define pb push_back
#define ll long long
#define ull unsigned long long
#define ld long double
#define inf 0x7FFFFFFF
#define inff 9223372036854775807
#define rep(i,l,r) for(int i=l;i<r;++i)
#define repp(i,l,r) for(int i=l;i<=r;++i)
#define per(i,r,l) for(int i=r-1;i>=l;--i)
#define pper(i,r,l) for(int i=r;i>=l;--i)
#define pii pair<int,int>
#define fi first
#define se second
#define p_q priority_queue
#define all(x) x.begin(),x.end()
#define rall(x) x.rbegin(),x.rend()
#define ls(x) ((x)<<1)
#define rs(x) ((x)<<1|1)
#define lb(x) ((x)&-(x))
#define lg(x) (31^__builtin_clz(x))
#define vi vector<int>
#define vii vector<pii >
#define vT vector<T>
#define mm1 mint(1)
using i64=long long;
using u32=unsigned;
using u64=unsigned long long;
constexpr u32 p=2;
using poly=std::vector<u32>;
using namespace std;
const unsigned int ep = ((1<<18)-1);
poly _like_inv(const poly&f){
    assert(!f.empty()&&f.front()==1);
    int n=f.size();
    poly res(n);
    res[0]=1;
    for(int i=1;i<n;++i){
        u32 tmp=0;
        for(int j=0;j<i;++j){
            tmp+=res[j]*f[i-j];
        }
        res[i]=-tmp;
    }
    return res;
}
constexpr u32 inv(u32 x){
    u32 iv=2-x;
    for(int i=0;i<4;++i){
        iv*=2-x*iv;
    }
    return iv;
}
poly _like_mul_trunc(poly a,poly b){
    int n=a.size();
    poly c(n);
    for(int i=0;i<n;++i){
        for(int j=0;j<=i;++j){
            c[i]+=a[j]*b[i-j];
        }
    }
    return c;
}
poly mul_trunc(poly a,poly b){
    int n=a.size();
    poly c(n);
    for(int i=0;i<n;++i){
        for(int j=0;j<=i;++j){
            c[i]+=a[j]*b[i-j];
        }
        c[i]&=ep;
    }
    return c;
}
poly _like_ln(poly a){
    auto p=_like_inv(a);
    for(int i=0;i<int(a.size());++i){
        a[i]*=i;
    }
    return _like_mul_trunc(a,p);
}
poly _like_exp(poly a){
    int n=a.size();
    poly res(n);
    res[0]=1;
    for(int i=1;i<n;++i){
        u32 tmp=0;
        for(int j=0;j<i;++j){
            tmp+=res[j]*a[i-j];
        }
        int ctz=__builtin_ctz(i);
        if((tmp&((1<<ctz)-1))!=0){
            std::cerr<<tmp<<" "<<i<<" "<<ctz<<"bad bad00\n";
        }
        res[i]=(tmp>>ctz)*inv(i>>ctz);
    }
    for(auto&x:res){
        x&=ep;
    }
    return res;
}
poly _like_add(poly a,poly b){
    int n=a.size();
    for(int i=0;i<n;++i){
        a[i]+=b[i];
    }
    return a;
}
ostream& operator << (ostream & cout,poly p){
    cout<<"(";
    for(auto x:p)cout<<x<<',';
    cout<<")";
    return cout;
}

poly A[1005];
poly LnA[1005],Sm,Exp,Mul;
void test(){
    std::mt19937_64 rng(random_device{}());
    for(int t=0;t<int(2000000);++t){
        int num = 100;
        int n=1100+2;
        std::cout<<n<<std::endl;
        rep(i,0,num){
            A[i].resize(n);
            A[i][0] = 1;
        }
        rep(i,0,num){
            rep(j,1,n){
                A[i][j] = (rng()&ep);
            }
        }
        // A0[0]=A1[0]=A2[0]=1;

        Sm.resize(n);
        Mul.resize(n);
        Exp.resize(n);
        rep(i,0,n)Sm[i] = 0;
        rep(i,0,n)Mul[i] = 0;
        rep(i,0,n)Exp[i] = 0;
        Mul[0] = 1;
        rep(i,0,num){
            LnA[i] = _like_ln(A[i]);
            Sm = _like_add(Sm,LnA[i]);
            Mul = mul_trunc(Mul,A[i]);
        }
        Exp = _like_exp(Sm);
        // poly lhs = mul_trunc(mul_trunc(A0,A1),A2);
        // poly rhs = _like_exp(_like_add(_like_add(_like_ln(A0),_like_ln(A1)),_like_ln(A2)));
        // poly tmp1 = _like_ln(A0);
        // poly tmp2 = _like_ln(A1);
        cout<<"tmp1 : "<<Exp<<endl;
        cout<<"tmp2 : "<<Mul<<endl;
        // cout<<lhs<<' '<<rhs<<endl;
        assert(Exp==Mul);
        cout<<"ok"<<endl;
        // if(t%int(1e6)==0){
        //     for(auto x:A0){
        //         cout<<x<<" ";
        //     }
        //     cout<<"\n";
        //     for(auto x:A1){
        //         cout<<x<<" ";
        //     }
        //     cout<<"\n";
        //     for(auto x:mul_trunc(A0,A1)){
        //         cout<<x<<" ";
        //     }
        //     cout<<"\n";
        //     cout<<std::endl;
        // }
    }
}
void solve(){
    poly A0{1,1,0,1,1},A1{1,0,0,1,1};
    //1 1 0 1 1
    //      1 1 0 1 1
    //        1 1 0 1 1
    //1 1 0 2 3 1 1 2 1
    //_like_add(_like_ln(A0),_like_ln(A1))
    //_like_ln(mul_trunc(A0,A1))
    for(auto x:_like_exp(_like_add(_like_ln(A0),_like_ln(A1)))){
        cout<<x<<" ";
    }
    //0 1 4294967295 4 4294967295
    //0 1 4294967295 7 3
    //0 1 4294967295 7 3
}
int main(){
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);
    test();
    return 0;
}

模 3 意义下的对拍：

点击查看代码

#include <bits/stdc++.h>
#define m_p(a,b) make_pair(a,b)
#define pb push_back
#define ll long long
#define ull unsigned long long
#define ld long double
#define inf 0x7FFFFFFF
#define inff 9223372036854775807
#define rep(i,l,r) for(int i=l;i<r;++i)
#define repp(i,l,r) for(int i=l;i<=r;++i)
#define per(i,r,l) for(int i=r-1;i>=l;--i)
#define pper(i,r,l) for(int i=r;i>=l;--i)
#define pii pair<int,int>
#define fi first
#define se second
#define p_q priority_queue
#define all(x) x.begin(),x.end()
#define rall(x) x.rbegin(),x.rend()
#define ls(x) ((x)<<1)
#define rs(x) ((x)<<1|1)
#define lb(x) ((x)&-(x))
#define lg(x) (31^__builtin_clz(x))
#define vi vector<int>
#define vii vector<pii >
#define vT vector<T>
#define mm1 mint(1)
#define int ll
// using i64=long long;
// using u32=unsigned;
// using u64=unsigned long long;
// constexpr u32 p=2;
using namespace std;

ll getmod(ll x,int d){
	ll rlt = 1;
	rep(i,0,d)rlt *= x;
	return rlt;
}
const int Base = 3,Power = 13;
const int mod = getmod(Base,Power);
const int phimod = mod/Base*(Base-1);
const int ep = getmod(Base,6);
ll pw(ll x,int d){
	ll t = 1;
	for(;d;d>>=1,x=x*x%mod)if(d&1)t = t*x%mod;
	return t;
}
struct m3{
	int v;
	m3(){v=0;}
	m3(int x){v=(x%mod+mod)%mod;}
	const m3 inv(){
		assert(v%3!=0);
		return pw(v,phimod-1);
	}
	const m3 operator + (const m3 y)const{return (v+y.v)%mod;}
	m3& operator += (const m3 y){return *this=*this+y;}
	const m3 operator - (const m3 y)const{return (v-y.v+mod)%mod;}
	m3& operator -= (const m3 y){return *this=*this-y;}
	const m3 operator * (const m3 y)const{return (v*y.v)%mod;}
	m3& operator *= (const m3 y){return *this=*this*y;}
	const m3 operator / (m3 y)const{
		assert(y.v!=0);
		m3 yy = y;
		m3 rlt = v;
		while(rlt.v%Base==0&&y.v%Base==0){
			rlt.v /= Base;
			y.v /= Base;
		}
		if(y.v%Base==0){
			cout<<"as :: "<<v<<" / "<<yy.v<<endl;
			assert(false);
		}
		return rlt*y.inv();
	}
	m3& operator /= (const m3 y){return *this=*this/y;}
	const bool operator == (const m3 rhs){return v==rhs.v;}
	const m3 operator - ()const {return m3(-v);}
};
ostream & operator << (ostream &cout,m3 x){
	cout<<x.v;
	return cout;
}
using poly=std::vector<m3>;
poly _like_inv(const poly&f){
    assert(!f.empty()&&f.front().v==1);
    int n=f.size();
    poly res(n);
    res[0]=1;
    for(int i=1;i<n;++i){
        m3 tmp=0;
        for(int j=0;j<i;++j){
            tmp+=res[j]*f[i-j];
        }
        res[i] = -tmp;
    }
    return res;
}
// constexpr m3 inv(m3 x){
//     m3 iv=2-x;
//     for(int i=0;i<4;++i){
//         iv*=2-x*iv;
//     }
//     return iv;
// }
poly _like_mul_trunc(poly a,poly b){
    int n=a.size();
    poly c(n);
    for(int i=0;i<n;++i){
        for(int j=0;j<=i;++j){
            c[i]+=a[j]*b[i-j];
        }
    }
    return c;
}
poly mul_trunc(poly a,poly b){
    int n=a.size();
    poly c(n);
    for(int i=0;i<n;++i){
        for(int j=0;j<=i;++j){
            c[i]+=a[j]*b[i-j];
        }
        c[i].v %= ep;
    }
    return c;
}
poly _like_ln(poly a){
    auto p=_like_inv(a);
    for(int i=0;i<(int)a.size();++i){
        a[i]*=i;
    }
    return _like_mul_trunc(a,p);
}
poly _like_exp(poly a){
    int n=a.size();
    poly res(n);
    res[0]=1;
    for(int i=1;i<n;++i){
        m3 tmp=0;
        for(int j=0;j<i;++j){
            tmp+=res[j]*a[i-j];
        }
        // int ctz=__builtin_ctz(i);
        // if((tmp&((1<<ctz)-1))!=0){
        //     std::cerr<<tmp<<" "<<i<<" "<<ctz<<"bad bad00\n";
        // }
        // res[i]=(tmp>>ctz)*inv(i>>ctz);
		tmp /= i;
		res[i] = tmp;
    }
    for(auto&x:res){
        x.v %= ep;
    }
    return res;
}
poly _like_add(poly a,poly b){
    int n=a.size();
    for(int i=0;i<n;++i){
        a[i]+=b[i];
    }
    return a;
}
ostream& operator << (ostream & cout,poly p){
    cout<<"(";
    for(auto x:p)cout<<x<<',';
    cout<<")";
    return cout;
}

void test(){
    std::mt19937_64 rng(random_device{}());
    // std::mt19937_64 rng(0);
    for(int t=0;t<100000;++t){
        // int n=rng()%5000+2;
		int n = 2187*3;
        std::cout<<n<<std::endl;
        poly A0(n),A1(n);
        A0[0]=A1[0]=1;
        for(int i=1;i<n;++i){
            A0[i]=((rng())%ep);
            A1[i]=((rng())%ep);
        }
        poly lhs = mul_trunc(A0,A1);
		// cout<<"lhs ok"<<endl;
		poly tmp = _like_add(_like_ln(A0),_like_ln(A1));
		// cout<<"tmp ok"<<endl;
        poly rhs = _like_exp(tmp);
		// cout<<"rhs ok"<<endl;
        rep(i,0,n)assert(lhs[i]==rhs[i]);
		// cout<<"lhs : "<<lhs<<endl;
		// cout<<"rhs : "<<rhs<<endl;
        cout<<"ok"<<endl;

		// poly C(n),D(n);
		// // cout<<"A0 : "<<A0<<endl;

		// C = _like_ln(A0);
		// // cout<<"C : "<<C<<endl;

		// D = _like_exp(C);
		// // cout<<"D : "<<D<<endl;
		// // assert(A0[0]==D[0]);
		// rep(i,0,n)assert(A0[i]==D[i]);
    }
}
void solve(){
    poly A0{1,1,0,1,1},A1{1,0,0,1,1};
    for(auto x:_like_exp(_like_add(_like_ln(A0),_like_ln(A1)))){
        cout<<x<<" ";
    }
}
signed main(){
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);
    test();
    return 0;
}