在线 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
呜呜呜，没有人看我魔怔
浙公网安备 33010602011771号