在线 O(1) 逆元(任意模数)

https://www.luogu.com/article/mgtvf9pq 的实现.

提交记录 https://qoj.ac/submission/418152.

代码
#include <bits/stdc++.h>
#include "inv.h"
using i64=long long;
using u32=unsigned;
using u64=unsigned long long;
using std::cin;
using std::cout;
//https://qoj.ac/submission/346829
inline u32 _udiv64(u64 u,u32 v){
	u32 u0=u,u1=u>>32,q,r;
	__asm__("divl %[v]" : "=a"(q), "=d"(r) : [v] "r"(v), "a"(u0), "d"(u1));
	return q;
}
int ttt=0;
struct ExgcdTable_i32{
    struct mat2x2{
        int a00,a01,a10,a11;
        mat2x2 operator*(const mat2x2&r)const{
            int b00=a00*r.a00+a01*r.a10;
            int b01=a00*r.a01+a01*r.a11;
            int b10=a10*r.a00+a11*r.a10;
            int b11=a10*r.a01+a11*r.a11;
            return {b00,b01,b10,b11};
        }
        std::tuple<int,int> apply(int b00,int b10)const{
            int c00=a00*b00+a01*b10;
            int c10=a10*b00+a11*b10;
            return {c00,c10};
        }
        void print(){
            cout<<a00<<","<<a01<<'\n';
            cout<<a10<<","<<a11<<'\n';
        }
    };
    int B,B2;
    std::vector<mat2x2> t1;
    std::vector<std::vector<mat2x2> > t2;
    std::vector<std::vector<std::pair<int,int> > > t3;
    mat2x2 bhgcd(int a,int b,int B)const{
        mat2x2 M={1,0,0,1}; 
        while(b>=B){
            int q=a/b,r=a%b;
            M=mat2x2{0,1,1,-q}*M,a=b,b=r;
        }
        return M;
    }
    ExgcdTable_i32(int _B):B(_B),B2(B*B),t1(B2),t2(3*B),t3(3*B){
        for(int i=B;i<B2;++i){t1[i]=bhgcd(B2-1,i,B);}
        for(int i=0;i<3*B;++i){t2[i].resize(i+1),t3[i].resize(i+1);}
        for(int i=1;i<3*B;++i){t2[i][0]={1,0,0,1},t3[i][0]={1,0};}
        for(int a=1;a<3*B;++a){
            for(int b=1;b<=a;++b){
                t2[a][b]=t2[b][a%b]*mat2x2{0,1,1,-a/b};
                t3[a][b]={t2[a][b].a00,t2[a][b].a01};
            }
        }
    }
    mat2x2 get(int a,int b)const{
        if(a<3*B){return t2[a][b];}
        if(a>i64(B)*b){return get(b,a%b)*mat2x2{0,1,1,-a/b};}
        auto M=t1[(b*i64(B2))/(a+1)];
        auto [x0,y0]=M.apply(a,b);
        if(x0<0){M.a00=-M.a00,M.a01=-M.a01,x0=-x0;}
        if(y0<0){M.a10=-M.a10,M.a11=-M.a11,y0=-y0;}
        if(x0<y0){std::swap(M.a00,M.a10),std::swap(M.a01,M.a11),std::swap(x0,y0);}
        return get(x0,y0)*M;
    }
    static int brute_calc(int a,int b){
        int uu=0,vv=1;
        while(b){
            std::tie(a,b,uu,vv)=std::tuple{b,u32(a)%u32(b),vv,uu-(u32(a)/u32(b))*vv};
        }
        return uu;
    }
    int calc(int a,int b)const{
        int uu=0,vv=1;
        while(a>=3*B){
            if(u32(a)>u64(B)*b){
                std::tie(a,b,uu,vv)=std::tuple{b,u32(a)%u32(b),vv,uu-(u32(a)/u32(b))*vv};
            }
            else{
                auto Mm=t1[_udiv64(b*u64(B2),a)];
                std::tie(a,b)=Mm.apply(a,b);
                std::tie(uu,vv)=Mm.apply(uu,vv);
                if(b<0){vv=-vv,b=-b;}
                if(a<b){std::swap(uu,vv),std::swap(a,b);}
            }
        }
        auto [kil,jok]=t3[a][b];
        return kil*uu+jok*vv;
    }
    int inv(int x,int mod)const{
        int res=calc(mod,x);
        return res<0?res+mod:res;
    }
    //x,y,r
    static std::tuple<int,int,int> exgcd(int a,int b){
        ++ttt;
        if(b==0){return {1,0,a};}
        auto [xx,yy,r]=exgcd(b,a%b);
        return {yy,xx-yy*(a/b),r};
    }
    static int inv_brute(int a,int mod){
        auto [x,y,r]=exgcd(mod,a);
        return y<0?y+mod:y;
    }
};
ExgcdTable_i32 t(150);
int mod;
void init(int p){
	mod=p;
}
int inv(int x){
	return t.inv(x,mod);
}

这个做法常数有点大,用点黑科技才卡过去.

卡常技巧 (thanks killer_joke,FR)

无需完整处理第一部分,预处理一个B^2级别的数的全部半gcd.

表太大会因为缓存原因变慢.

用无符号除法,最后一段只记录两个系数而非 2x2 矩阵.

从 FR 提交里搬了个 u64 除 u32.

大体耗时对比:

方法 耗时(\(10^7\)次,on atc)
\(exgcd\) \(1030ms\)
\(exgcd\) (仅计算逆元版) \(900ms\)
快速幂(固定模数) \(1060ms\)
快速幂(蒙哥马利约减) \(990ms\)
\(this\) (计算完整矩阵) \(800ms\)
\(this\) (仅计算逆元) \(590ms\)
经典在线 \(O(1)\) 逆元(取qoj最优解) \(200ms\)

此外,这样求解 exgcd 的时候不能保证不互质时的值域最小性,但非平凡时满足 \(|y| < a\), \(|x| < b\).

测试代码

\(exgcd\)

std::tuple<int,int,int> exgcd(int a,int b){
    if(b==0){return {1,0,a};}
    auto [xx,yy,r]=exgcd(b,a%b);
    return {yy,xx-yy*(a/b),r};
}
int inv(int a,int mod){
    auto [x,y,r]=exgcd(mod,a);
    return y<0?y+mod:y;
}

\(exgcd\) (仅计算逆元)

int calc(int a,int b){
    int uu=0,vv=1;
    while(b){
        std::tie(a,b,uu,vv)=std::tuple{b,u32(a)%u32(b),vv,uu-(u32(a)/u32(b))*vv};
    }
    return uu;
}
int inv(int x,int mod){
    int res=calc(mod,x);
    return res<0?res+mod:res;
}

快速幂(固定模数)

constexpr u32 mul(u32 a,u32 b){
    return u64(a)*b%998244353;
}
constexpr u32 qpw(u32 a,u32 b,u32 r=1){
    for(;b;b>>=1,a=mul(a,a)){
        if(b&1){
            r=mul(r,a);
        }
    }
    return r;
}
constexpr u32 inv(u32 x){
    return qpw(a,998244351);
}

快速幂(蒙哥马利约减)

constexpr u32 get_nr(u32 m){
    u32 niv=2+m;
    for(int i=0;i<4;++i){niv*=2+m*niv;}
    return niv;
}
constexpr u32 mod=998244353,niv=get_nr(mod);
constexpr u32 reduce(u64 x){
    return (x+u64(u32(x)*niv)*mod)>>32;
}
constexpr u32 mul(u32 x,u32 y){
    return reduce(u64(x)*y);
}
constexpr u32 qpw(u32 a,u32 b,u32 r){
    for(;b;b>>=1,a=mul(a,a)){
        if(b&1){
            r=mul(r,a);
        }
    }
    return std::min(r,r-mod);
}
constexpr u32 inv(u32 x){
    constexpr u32 ou1=reduce(1);
    return qpw(x,mod-2,ou1);
}
EI曾说过

【保留】【不保留】 2024/2/11 18:18:46
给定 n 个数,问是否存在三个数的异或和为 0,能做到低于 O(n³) 或规约矩阵乘法吗/kel(视 logV 为 64)

JOHNKRAM 2024/2/11 18:19:30
这不枚举两个数就行了

【保留】【不保留】 2024/2/11 18:19:55
*优于 O(n²)

旧日余孽 2024/2/11 18:41:07
感觉低于n^2,很简单

旧日余孽 2024/2/11 18:41:31
我来给大家云一个至少能除以(logn)^0.1的办法

旧日余孽 2024/2/11 18:43:33
首先把64个bit,随机映射到 k = log (n(logn)^0.1) 个 bit

旧日余孽 2024/2/11 18:44:02
然后问题转化为判断是否有一组 a(i) ^ b(j) ^ c(i^j)=0

旧日余孽 2024/2/11 18:44:56
当然首先这个有哈希冲突,要把出现超过一次的位置都图图,但这个期望只消耗 n2/(log)0.1 的时间,因为冲突不多

旧日余孽 2024/2/11 18:45:46
然后解决任何一个 b * b 的小块,先把他们的线性关系消元出来,这只需要 b 的时间,然后查表

旧日余孽 2024/2/11 18:46:22
打表需要 2^ b^2 时间,所以可开到b = (logn)^.49

旧日余孽 2024/2/11 18:47:56
呜呜呜,没有人看我魔怔

posted @ 2024-05-23 16:38  QedDust  阅读(483)  评论(0)    收藏  举报