怎么還有厲害的在線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;
}