\(k\ge\varphi(m)\)時,\(b^k\equiv b^{k\%\varphi(m)+\varphi(m)}(\bmod m\))html
推式子合併同餘方程。
http://www.javashuo.com/article/p-xqxomktl-dw.htmlc++
根號大暴力,就是細節不少。
http://www.javashuo.com/article/p-zlwoggrh-ea.html算法
對於模數不是質數的,惟一分解爲\(\prod p_i^{k_i}\),算出組合數對每個\(p_i^{k_i}\)取模的結果,用CRT合併。
問題在於求\(n!\%p^k\)。
將\(n!\)中全部\(p\)因子提出來,因子個數\(\sum\limits_{i=1}\lfloor\frac{n}{p^i}\rfloor\)。
剩下的數的貢獻單獨考慮,記爲\(F_n\),發現它在模意義下以\(p^k\)爲週期,\(fac_n=F_{p^k}^{\lfloor\frac{n}{p^k}\rfloor}F_{n\%p^k}fac_{\lfloor\frac n p\rfloor}\),遞歸處理。
時間複雜度\(\sum p_i^{k_i}\)。
洛谷P4720 【模板】擴展盧卡斯數組
#include<bits/stdc++.h> #define LL long long #define R register LL using namespace std; LL n,m,YL,F[1000009]; void exgcd(R a,R b,R&x,R&y){ if(!b){x=1;y=0;return;} exgcd(b,a%b,y,x);y-=a/b*x; } inline LL qpow(R b,R k,R p){ R a=1; for(;k;k>>=1,b=b*b%p) if(k&1ll)a=a*b%p; return a; } LL fac(R n,R p,R k){ return n?qpow(F[k],n/k,k)*F[n%k]%k*fac(n/p,p,k)%k:1; } inline LL inv(R n,R p){ R x,y;exgcd(n,p,x,y); return x<=0?x+p:x; } inline LL cnt(R n,R p){ R k=0; for(n/=p;n;n/=p)k+=n; return k; } inline LL C(R n,R m,R p,R k){ for(R i=F[0]=1;i<=k;++i) F[i]=i%p?F[i-1]*i%YL:F[i-1]; return fac(n,p,k)*inv(fac(m,p,k),k)%k*inv(fac(n-m,p,k),k)%k*qpow(p,cnt(n,p)-cnt(m,p)-cnt(n-m,p),k)%k; } inline LL CRT(R a,R p){ return YL/p*inv(YL/p,p)%YL*a%YL; } int main(){ cin>>n>>m>>YL; R now=YL,lim=sqrt(YL),ans=0; for(R p=2;p<=lim;++p){ if(now%p)continue; R k=1;while(now%p==0)now/=p,k*=p; ans=(ans+CRT(C(n,m,p,k),k))%YL; lim=sqrt(now); } if(now>1)ans=(ans+CRT(C(n,m,now,now),now))%YL; cout<<ans<<endl; return 0; }
背板子
洛谷P4718 【模板】Pollard-Rho算法
辣雞卡常題,參考了一下討論版裏的卡常技巧
upd:模數爲long long
級別的乘法取模原理
若是NOI能用什麼__int128
的話倒還好說
可是別忘了long double
和__int128
的位長是同樣的啊!所以用它作long long
級別的運算的精度損失是徹底能夠接受的。
咱們知道\(a*b\bmod p=a*b-\lfloor\frac{a*b}p\rfloor*p\)
而後咱們用long double
先除後乘,考慮浮點數除法會丟掉小數點後面一個很小的值,而以後強制轉換整形是要強制捨去的,因此要加一個eps。函數
#include<bits/stdc++.h> #define SL __int128 #define LL long long #define RG register #define R RG int using namespace std; LL t,n,ans; inline LL Mul(LL a,LL b,LL p){//擓的玄學乘法取模 LL d=((long double)a/p*b+1e-8); LL r=a*b-d*p; return r<0?r+p:r; } inline LL Gcd(LL a,LL b){//擓的builtin展轉相減 if(!a||!b)return a|b; int t=__builtin_ctzll(a|b); a>>=__builtin_ctzll(a); do{ b>>=__builtin_ctzll(b); if(a>b)swap(a,b); b-=a; }while(b); return a<<t; } inline LL Qpow(LL b,LL k,LL p,LL a=1){ for(;k;k>>=1,b=(SL)b*b%p) if(k&1ll)a=(SL)a*b%p; return a; } inline LL MR(LL n){ if(n==2)return 1; if(n==1||(1&n)==0)return 0; static int a[]={2,3,7,13,61,24251}; LL u=n-1,x,y;R k=0; while((1&u)==0)u>>=1,++k; for(R i=0;i<6&&a[i]<n;++i){ x=Qpow(a[i],u,n); for(R j=0;j<k;++j,x=y){ y=(SL)x*x%n; if(y==1&&x!=1&&x!=n-1)return 0; } if(x!=1)return 0; } return 1; } inline LL F(LL a,LL c,LL n){ LL t=Mul(a,a,n)+c; return t<n?t:t-n; } inline LL PR(LL n){ if((1&n)==0)return 2; LL c=rand(),a=rand(),b=a,g; do{ a=F(a,c,n);b=F(F(b,c,n),c,n); g=Gcd(n,abs(a-b)); if(g!=1&&g!=n)return g; }while(a!=b); return g; } void find(LL n){ if(n==1||n<=ans)return; if(MR(n)){ans=max(ans,n);return;} LL d=n; while(d==n)d=PR(n); while(n%d==0)n/=d; find(n);find(d); } int main(){ srand(time(NULL)); cin>>t; while(t--){ cin>>n;ans=1;find(n); if(ans==n)puts("Prime"); else cout<<ans<<'\n'; } return 0; }
積性函數的性質:兩個積性函數的狄利克雷卷積仍是積性函數。
有一些推式子題裏須要經過構造積性函數來加快函數求值。
經常使用的卷積式
\(\mu*\textbf1=\epsilon\Leftrightarrow\sum\limits_{d|n}\mu(d)=[n==1]\)
\(\varphi*\textbf1=\textbf{id}\Leftrightarrow\sum\limits_{d|n}\varphi(d)=n\)
\(\mu*\textbf{id}=\varphi\Leftrightarrow\sum\limits_{d|n}\frac{n}{d}\mu(d)=\varphi(n)\)
\(\textbf1*\textbf1=\textbf{d};\textbf{d}*\textbf1=\sigma\)(感受還有不少都是能夠相互推出來的)ui
杜教篩:求積性函數\(f\)的前綴和(記爲\(S\))。
構造積性函數\(g\),前提是\(g\)和\(F=f*g\)的前綴和均可以\(O(1)\)算。
\(\sum\limits_{i=1}^nF(i)=\sum\limits_{i=1}^n\sum\limits_{xy=i}f(x)g(y)=\sum\limits_{y=1}^ng(y)\sum\limits_{x=1}^{\lfloor\frac{n}{y}\rfloor}f(x)=g(1)S(x)+\sum\limits_{y=2}^ng(y)S(\lfloor\frac n y\rfloor)\)
\(g(1)S(x)=\sum\limits_{i=1}^nF(i)-\sum\limits_{y=2}^ng(y)S(\lfloor\frac n y\rfloor)\)
看到了能夠整除分塊+遞歸的地方。
預處理\(F,g\)到\(n^\frac 2 3\)的前綴和,時間複雜度\(O(n^\frac 2 3)\)。
洛谷P4213 【模板】杜教篩(Sum)spa
#include<bits/stdc++.h> #define LL long long #define UI unsigned int #define RG register #define R RG int using namespace std; const LL N=2147483648,M=1.7e6,L=2000; int pr[M],cnt;bool vis[M]; struct Dat{LL p;int u;}s[M],t[L]; Dat&S(UI n){ if(n<M)return s[n]; UI x=N/n;if(vis[x])return t[x];vis[x]=1; Dat&ans=t[x]; ans.p=(LL)n*(n+1)>>1;ans.u=1; for(UI l=2,r;l<=n;l=r+1){ r=n/(n/l);Dat&ret=S(n/l); ans.p-=(r-l+1)*ret.p; ans.u-=(r-l+1)*ret.u; } return ans; } int main(){ s[1]=(Dat){1,1};vis[1]=1; for(R i=2;i<M;++i){ if(!vis[i])s[pr[++cnt]=i]=(Dat){i-1,-1}; for(R j=1,x;j<=cnt&&(x=i*pr[j])<M;++j){ vis[x]=1; if(i%pr[j]==0){s[x].p=s[i].p*pr[j];break;} s[x].p=s[i].p*(pr[j]-1);s[x].u=-s[i].u; } s[i].p+=s[i-1].p;s[i].u+=s[i-1].u; } memset(vis,0,M); R t,n;cin>>t; while(t--){ cin>>n; Dat ret=S(n); cout<<ret.p<<' '<<ret.u<<endl; memset(vis,0,L); } return 0; }
一樣是求積性函數前綴和,不過擴展性很強。
yyb巨佬告訴我這是今年新鮮出爐的篩法,吊打了某洲閣篩(話說寫這句話的時候是今年的最後一天了)
核心思想是體如今Step2中的將最小質因子相同的數的貢獻一塊兒算。
應用前提:對於\(p\in P\),\(f(p)\)是一個簡單多項式,\(f(p^e)\)能夠快速計算。比某杜教篩不知道高到哪裏去了。
時間複雜度玄學,可是跑得快
http://www.javashuo.com/article/p-qfqnlnsq-cn.html
http://www.javashuo.com/article/p-ocbisnol-do.html
http://www.javashuo.com/article/p-pghfthdh-da.html
http://www.javashuo.com/article/p-utazeyfm-cq.html
http://www.javashuo.com/article/p-kgecsqlp-cb.html (都是高級操做)code
把多項式拆開,對每一個\(f(x)=x^k\)的項分開算貢獻。
把\(\sqrt N\)範圍的質數篩出來,記爲\(P\)。預處理質數的函數前綴和\(fs\)。
先對於每一個\(i\le \sqrt N\)求\(\sum\limits_{x=1}^{\lfloor\frac N i\rfloor}[x\in P]f(x)\)。
設\(g(n,j)\)爲\(x\in[1,n]\)中知足\(x\in P\)或\(x\)的最小質因子小於\(P_j\)的\(f(x)\)之和
因而顯然有\(g(n,|P|)=\sum\limits_{x=1}^{n}[x\in P]f(x)\)。
實現過程:滾動數組,外層枚舉\(j\),直接用一維數組保存每個\(n=\lfloor\frac N i\rfloor\)當前的\(g(n,j)\)。
內層枚舉\(i\),有\(g(n,j)=g(n,j-1)-f(P_j)\cdot(g(\lfloor\frac{n}{P_j}\rfloor,j-1)-fs(j-1))\),結合定義用容斥思想理解。htm
開始算答案了。
設\(S(n,j)\)爲\(x\in[1,n]\)中知足\(x\)的最小質因子大於等於\(P_j\)的\(f(x)\)之和(注意和\(g\)的區別)
因而咱們要求的答案便是\(S(N,1)+f(1)\)(\(S\)裏面不包含\(f(1)\))
求\(S(n,j)\)時,分質數和合數算貢獻
質數:\(g(n,|P|)-fs(j-1)\)
合數:枚舉最小質因子\(P_j\)以及它的次數\(e\),相同的一塊兒算(注意不包含\(f(1)\),所以還有一個尾項)
\(\sum\limits_{k=j}^{P_k^2\le n}\sum\limits_{e=1}^{P_k^{e+1}\le n}S(\lfloor\frac{n}{P_k^e}\rfloor,k+1)f(P_k^e)+f(P_k^{e+1})\)
LOJ6053 簡單的函數blog
#include<bits/stdc++.h> #define LL long long #define RG register #define R RG LL using namespace std; const LL N=2e5+9,YL=1e9+7; LL n,Sq,p,m,pr[N],id1[N],id2[N],w[N],fx[N],gx[N],g1[N]; bool np[N]; inline LL Getid(R x){ return x<=Sq?id1[x]:id2[n/x]; } void Sieve(){ for(R i=2;i<=Sq;++i){ if(np[i])continue; pr[++p]=i;fx[p]=(fx[p-1]+i)%YL; for(R j=i*i;j<=Sq;j+=i)np[j]=1; } } LL S(R n,R j){ if(n<=1||pr[j]>n)return 0; R x=Getid(n); R ret=(gx[x]-fx[j-1]-g1[x]+j-1+2*YL)%YL; for(R k=j;k<=p&&pr[k]*pr[k]<=n;++k){ R p1=pr[k],p2=p1*p1; for(R e=1;p2<=n;++e,p1=p2,p2*=pr[k]) ret=(ret+S(n/p1,k+1)*(pr[k]^e)+(pr[k]^(e+1)))%YL; } return ret; } int main(){ cin>>n; if(n==1)return puts("1"),0; Sq=sqrt(n);Sieve(); for(R i=1,j;i<=n;i=j+1){ R x=w[++m]=n/i;j=n/x; (x<=Sq?id1[x]:id2[n/x])=m; gx[m]=(x+2)%YL*((x-1)%YL)%YL*((YL+1)>>1)%YL; g1[m]=(x-1)%YL; } for(R j=1;j<=p;++j) for(R i=1;i<=m&&pr[j]*pr[j]<=w[i];++i){ R x=Getid(w[i]/pr[j]); gx[i]=((gx[i]-pr[j]*(gx[x]-fx[j-1]))%YL+YL)%YL; g1[i]=(g1[i]-g1[x]+j-1+YL)%YL; } return cout<<S(n,1)+3<<endl,0; }