在线O(1)求逆元

2022-01-11 15:04:04 阅读：425 来源： 互联网

~~怎么还有厉害的在线O(1)求逆元，不过常数确实有点儿太大了~~

~~本文大部分搬运于这里~~

相信大家都做过 POJ2478 这道题吧，这道题的 Farey 序列 $F_n$ 包含了分子分母不大于 $n$ 且互质的数。该分数可以为 $0$ 和 $1$。

嗯我们现在要把 $F_{\sqrt [3]p}$求出来，然后有一个妙妙定理，就可以在线地 $O(1)$ 求逆元。

定理：给定一个分数 $\frac a b(0 \leq \frac a b \leq 1) $和一个 $ n(n > 1)$，你都可以在 $F_{n-1}$ 中找到一个分数使得 $|\frac a b-\frac x y| \leq \frac 1 {yn}$，并且这个分数一定是 $\frac a b$ 在序列中使用 useful binary search 后往前或往后的第一个数。

证明？注意到 $F_{n-1}$ 实际上包含了分子分母都不大于 $n$ 的所有分数，只不过将其去重了。所以我们对每一个分母都讨论一遍。

差不多就是要证明对于所有不大于 $n$ 的数对 $(x,y)$，$[\frac x y-\frac 1 {ny},\frac x y+\frac 1 {ny}]$ 的并包含了 $[0,1]$。也就是对于任意一个区间，一定有另一个在其右侧区间与其相交。

通分：$[\frac {nx-1} {ny},\frac {nx+1} {ny}]$。

那么我们接下来需要证明对于 $x_1<y_1$ 一定有 $x_2<y_2$ 满足 $\frac {x_1}{y_1} \ne \frac {x_2}{y_2}$ 且 $\frac {nx_1+1}{ny_1}>\frac {nx_2-1}{ny_2}$。

还是通分：$nx_1y_2+y_2>nx_2y_1-y_1$，也就是 $y_1+y_2>n(x_2y_1-x_1y_2)$。很明显一定有 $x_2=x_1-1,y_2=y_1-1$ 满足条件。

啥？$x_1=0,y_1=1$？右端点为 $\frac 1 n$，那么存在 $\frac 1 {n-1}$ 满足 $\frac 1 {n-1} - \frac 1 {n(n-1)}=\frac 1 n$。

所以自然就是二分查找后向前或向后的第一个数。

~~以上内容均为口胡~~

我们发现这个结论 $|\frac a b - \frac x y| \leq \frac 1 {yn}$，通分之后就是 $| ay-bx | \leq \lfloor \frac b n \rfloor$。这里把 $b$ 用 $p$ 替换一下。

上面的结论相当于 $ay \equiv t \pmod p$，且 $t \leq \lfloor \frac p n \rfloor$，相当于 $a^{-1} \equiv yt^{-1} \pmod p$。

所以只需要找到这个分数之后，处理 $t$ 的逆元就好了啊。

因为 $t \leq \lfloor \frac p n \rfloor$，我们预处理一个长度为 $n$ 的 Farey 肯定需要至少 $O(n^2)$ 的时间，$n^2=\frac p n$，得到 $n=\sqrt [3]p$。所以需要预处理 $F_{\sqrt [3]p}$ 和 $p^{\frac 2 3}$ 以内的数的逆元。

那么如何 $O(n^2)$ 预处理 Farey 序列和 $O(1)$ 二分查找？

~~我们可以请教神圣二分帝国皇帝 Um_nik~~

对于 Farey 序列中的任意两个数 $\frac {x_1}{y_1}$ 和 $\frac {x_2}{y_2}$，有 $\frac {x_1}{y_1}-\frac {x_2}{y_2}=\frac {x_1y_2-x_2y_1}{y_1y_2}$。在分母最大的情况下不大于 $p^{\frac 2 3}$，分子最小不小于 $1$，所以只需要将原分数乘上 $p^{\frac 2 3}$ 后向下取整，一定互不相同。

所以只需要把 $\frac x y$ 映射到 $\lfloor \frac {x \times p^{\frac 2 3}} y \rfloor$ 后桶排序就可以啦。

代码压过行，还卡过常，实际上写起来也不是很长。

目测常数是正常离线求逆的 $2 \sim 3$ 倍。

#include<iostream>
#include<cctype>
#include<cmath>
typedef unsigned ui;
typedef __uint128_t L;
typedef unsigned long long ull;
const ui M1=1000,M2=M1*M1;
ui T,m1,m2,MOD,len,fra[M2+5],inv[M2+5],sum[M2+5];
ui q[M2+5];ull p[M2+5];
double invm1,INV;
char buf[1<<22|1],*p1=buf;
inline char Getchar(){
	return*p1=='\0'&&(std::cin.read(p1=buf,1<<22)),*p1++;
}
struct FastMod{
	ull b,m;
	FastMod(ull b):b(b),m(ull((L(1)<<64)/b)){}
	friend inline ull operator%(const ull&a,const FastMod&mod){
		ull r=a-(L(mod.m)*a>>64)*mod.b;
		return r>=mod.b?r-mod.b:r;
	}
}mod(2);
inline ull abs(const ull&a){
	return a>>63?-a:a;
}
inline ui read(){
	ui n(0);char s;
	while(!isdigit(s=Getchar()));
	while(n=n*10+(s&15),isdigit(s=Getchar()));
	return n;
}
inline void init(){
	ui i,j,x;
	for(i=1;i^m1;++i){
		const double&INV=1./i+1e-15;
		for(j=0;j^i;++j)if(!sum[x=1ull*j*m2*INV])sum[x]=1,fra[x]=m1*i+j;
	}
	for(i=0;i<=m2;++i){
		if(sum[i])++len,q[len]=fra[i]*invm1,p[len]=1ull*(fra[i]-q[len]*m1)*MOD;
		if(i)sum[i]+=sum[i-1];inv[i]=i>1?1ull*(MOD-MOD/i)*inv[MOD%i]%mod:i;
	}
}
inline ui Inv(const ui&a){
	static ui q,k;static ull t;
	if(a<=m2)return inv[a];if(MOD-a<=m2)return MOD-inv[MOD-a];k=sum[ui(a*INV)];
	if(k<=len){
		q=::q[k];t=1ull*a*q-p[k];
		if(abs(t)<=m2)return 1ull*q*(t>>63?MOD-inv[-t]:inv[t])%mod;
	}
	if(++k<=len){
		q=::q[k];t=1ull*a*q-p[k];
		if(abs(t)<=m2)return 1ull*q*(t>>63?MOD-inv[-t]:inv[t])%mod;
	}
	return-1;
}
signed main(){
	std::ios::sync_with_stdio(false);std::cin.tie(0);std::cout.tie(0);
	ui k,x,sum(0);T=read();MOD=read();mod=FastMod(MOD);x=k=read();
	m1=pow(MOD,1./3)+1;m2=m1*m1;invm1=1./m1+1e-15;INV=1.*m2/MOD+1e-15;init();
	while(T--)sum=(sum+1ull*x*Inv(read()))%mod,x=1ull*x*k%mod;std::cout<<sum;
}

标签：frac,在线,leq,逆元,ui,m1,mod,MOD
来源： https://www.cnblogs.com/lmpp/p/15788297.html

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

ICode9

在线O(1)求逆元