【洛谷3768】简单的数学题【莫比乌斯反演】【杜教筛】【小学奥数】
傳送門
題意:給定p,Np,Np,N,求
∑i=1N∑j=1Nijgcd(i,j)modp\sum_{i=1}^{N}\sum_{j=1}^{N}ijgcd(i,j)\text{ }mod \text{ }pi=1∑N?j=1∑N?ijgcd(i,j)?mod?p
ppp為質數,在1e91e91e9左右
N≤1e10N \leq 1e10N≤1e10
神仙題
前置芝士:杜教篩
懶得重新寫,所以順便講一下
杜教篩可以在非線性時間求積性函數的前綴和
設所求函數為f(n)f(n)f(n),S(n)為前綴和S(n)為前綴和S(n)為前綴和
S(n)=∑i=1nf(i)S(n)=\sum_{i=1}^{n}f(i)S(n)=i=1∑n?f(i)
構造一個可以快速計算前綴和的函數g(n)g(n)g(n),且和f(n)f(n)f(n)的狄利克雷卷積(f?g)(n)(f*g)(n)(f?g)(n)的前綴和可以快速求出
考慮
∑i=1n(f?g)(i)\sum_{i=1}^{n}(f*g)(i)i=1∑n?(f?g)(i)
暴力拆開
∑i=1n∑d∣if(d)g(id)\sum_{i=1}^{n}\sum_{d \mid i}f(d)g(\frac{i}ze8trgl8bvbq)i=1∑n?d∣i∑?f(d)g(di?)
枚舉ddd
∑d=1ng(d)∑d∣if(id)\sum_{d=1}^{n}g(d)\sum_{d \mid i}f(\frac{i}ze8trgl8bvbq)d=1∑n?g(d)d∣i∑?f(di?)
∑d=1ng(d)S(?nd?)\sum_{d=1}^{n}g(d)S(\lfloor\frac{n}ze8trgl8bvbq\rfloor)d=1∑n?g(d)S(?dn??)
我們要求的是S(n)S(n)S(n),所以可以求出g(1)S(n)g(1)S(n)g(1)S(n)
即
∑d=1ng(d)S(?nd?)?∑d=2ng(d)S(?nd?)\sum_{d=1}^{n}g(d)S(\lfloor\frac{n}ze8trgl8bvbq\rfloor)-\sum_{d=2}^{n}g(d)S(\lfloor\frac{n}ze8trgl8bvbq\rfloor)d=1∑n?g(d)S(?dn??)?d=2∑n?g(d)S(?dn??)
換成上面
∑i=1n(f?g)(i)?∑d=2ng(d)S(?nd?)\sum_{i=1}^{n}(f*g)(i)-\sum_{d=2}^{n}g(d)S(\lfloor\frac{n}ze8trgl8bvbq\rfloor)i=1∑n?(f?g)(i)?d=2∑n?g(d)S(?dn??)
左邊直接求,右邊整除分塊套遞歸
復雜度是O(N34)O(N^\frac{3}{4})O(N43?)
如果線性篩出N23N^\frac{2}{3}N32?以內的并用mapmapmap記憶化,可以達到O(N23)O(N^\frac{2}{3})O(N32?)
正文
∑i=1N∑j=1Nijgcd(i,j)\sum_{i=1}^{N}\sum_{j=1}^{N}ijgcd(i,j)i=1∑N?j=1∑N?ijgcd(i,j)
套路
∑d=1N∑i=1N∑j=1N[gcd(i,j)=d]ijd\sum_{d=1}^{N}\sum_{i=1}^{N}\sum_{j=1}^{N}[gcd(i,j)=d]ijdd=1∑N?i=1∑N?j=1∑N?[gcd(i,j)=d]ijd
∑d=1Nd∑i=1N∑j=1N[gcd(i,j)=d]ij\sum_{d=1}^{N}d\sum_{i=1}^{N}\sum_{j=1}^{N}[gcd(i,j)=d]ijd=1∑N?di=1∑N?j=1∑N?[gcd(i,j)=d]ij
右邊可以反演
設
f(n)=∑i=1N∑j=1N[gcd(i,j)=n]ijf(n)=\sum_{i=1}^{N}\sum_{j=1}^{N}[gcd(i,j)=n]ijf(n)=i=1∑N?j=1∑N?[gcd(i,j)=n]ij
F(d)=∑d∣nf(n)=∑i=1N∑j=1N[d∣i][d∣j]ijF(d)=\sum_{d\mid n}f(n)=\sum_{i=1}^{N}\sum_{j=1}^{N}[d\mid i][d\mid j]ijF(d)=d∣n∑?f(n)=i=1∑N?j=1∑N?[d∣i][d∣j]ij
記sum(n)=n(n+1)2sum(n)=\frac{n(n+1)}{2}sum(n)=2n(n+1)?
F(d)=d2(sum?Nd?)2F(d)=d^2(sum\lfloor\frac{N}ze8trgl8bvbq\rfloor)^2F(d)=d2(sum?dN??)2
f(d)=∑d∣nF(n)μ(nd)=∑d∣nn2(sum?Nn?)2μ(nd)f(d)=\sum_{d\mid n}F(n)\mu(\frac{n}ze8trgl8bvbq)=\sum_{d\mid n}n^2(sum\lfloor\frac{N}{n}\rfloor)^2\mu(\frac{n}ze8trgl8bvbq)f(d)=d∣n∑?F(n)μ(dn?)=d∣n∑?n2(sum?nN??)2μ(dn?)
原來求的是
∑d=1Ndf(d)\sum_{d=1}^{N}df(d)d=1∑N?df(d)
∑d=1Nd∑d∣nn2(sum?Nn?)2μ(nd)\sum_{d=1}^{N}d\sum_{d\mid n}n^2(sum\lfloor\frac{N}{n}\rfloor)^2\mu(\frac{n}ze8trgl8bvbq)d=1∑N?dd∣n∑?n2(sum?nN??)2μ(dn?)
枚舉nnn
∑n=1Nn2(sum?Nn?)2∑d∣ndμ(nd)\sum_{n=1}^{N}n^2(sum\lfloor\frac{N}{n}\rfloor)^2\sum_{d \mid n}d\mu(\frac{n}ze8trgl8bvbq)n=1∑N?n2(sum?nN??)2d∣n∑?dμ(dn?)
右邊就是φ\varphiφ
∑n=1N(sum?Nn?)2n2φ(n)\sum_{n=1}^{N}(sum\lfloor\frac{N}{n}\rfloor)^2n^2\varphi(n)n=1∑N?(sum?nN??)2n2φ(n)
左邊是個整除分塊,所以求出右邊就可以了
#undef f
設f(n)=n2φ(n)f(n)=n^2\varphi(n)f(n)=n2φ(n)
顯然這是個積性函數,考慮杜教篩
(f?g)(n)=∑d∣ng(d)f(nd)=∑d∣ng(d)n2d2φ(nd)(f*g)(n)=\sum_{d \mid n}g(d)f(\frac{n}ze8trgl8bvbq)=\sum_{d \mid n}g(d)\frac{n^2}{d^2}\varphi(\frac nd)(f?g)(n)=d∣n∑?g(d)f(dn?)=d∣n∑?g(d)d2n2?φ(dn?)
令g(n)=n2g(n)=n^2g(n)=n2就可以消掉
則
(f?g)(n)=∑d∣nn2φ(nd)=n2∑d∣nφ(d)=n3(f*g)(n)=\sum_{d \mid n}n^2\varphi(\frac nd)=n^2\sum_{d \mid n}\varphi(d)=n^3(f?g)(n)=d∣n∑?n2φ(dn?)=n2d∣n∑?φ(d)=n3
眾所周知
12+22+32+...+n2=(n)(n+1)(2n+1)61^2+2^2+3^2+...+n^2=\frac{(n)(n+1)(2n+1)}{6}12+22+32+...+n2=6(n)(n+1)(2n+1)?
13+23+33+...+n3=(1+2+3+...+n)21^3+2^3+3^3+...+n^3=(1+2+3+...+n)^213+23+33+...+n3=(1+2+3+...+n)2
兩個前綴和都可以O(1)O(1)O(1)算出
這個問題就解決了
#include <iostream> #include <cstdio> #include <cstring> #include <cctype> #include <map> #define MAXN 10000005 using namespace std; typedef long long ll; int MOD,inv2,inv6; inline int add(const int& x,const int& y){return x+y>=MOD? x+y-MOD:x+y;} inline int dec(const int& x,const int& y){return x-y<0? x-y+MOD:x-y;} inline int mul(const int& x,const int& y){return (ll)x*y%MOD;} inline int qpow(int a,int p) {int ans=1;while (p){if (p&1) ans=mul(ans,a);a=mul(a,a);p>>=1;}return ans; } const int N=10000000; int np[MAXN],pl[1000005],cnt; int phi[MAXN],f_sum[MAXN]; void init() {np[1]=phi[1]=1;for (int i=2;i<=N;++i){if (!np[i]) phi[pl[++cnt]=i]=i-1;int x;for (int j=1;(x=i*pl[j])<=N;++j){np[x]=1;if (i%pl[j]==0){phi[x]=mul(phi[i],pl[j]);break;}phi[x]=mul(phi[i],pl[j]-1);}}for (int i=1;i<=N;i++) f_sum[i]=add(f_sum[i-1],mul(phi[i],mul(i,i))); } inline int calc(const int& n){return mul(mul(n,n+1),inv2);} inline int calc2(const int& n){return mul(mul(n,n+1),mul(2*n+1,inv6));} inline int calc3(const int& n){int t=calc(n);return mul(t,t);} map<int,int> m; int getsum(ll n) {if (n<=N) return f_sum[n];if (m[n]) return m[n];int ans=calc3(n%MOD);for (ll l=2,r;l<=n;l=r+1){r=n/(n/l);int tmp=getsum(n/l);ans=dec(ans,mul(dec(calc2(r%MOD),calc2((l-1)%MOD)),tmp));}return m[n]=ans; } //int gcd(int a,int b) //{ // return b? gcd(b,a%b):a; //} int main() {ll n;scanf("%d%lld",&MOD,&n);inv2=(MOD+1)>>1;inv6=qpow(6,MOD-2);init(); int ans=0; // for (int i=1;i<=n;i++) // for (int j=1;j<=n;j++) // ans=(ans+(ll)i*j%MOD*gcd(i,j))%MOD; // cerr<<ans<<'\n'; // ans=0;for (ll l=1,r;l<=n;l=r+1){r=n/(n/l);ans=add(ans,mul(calc3((n/l)%MOD),dec(getsum(r),getsum(l-1))));}printf("%d\n",ans);return 0; }總結
以上是生活随笔為你收集整理的【洛谷3768】简单的数学题【莫比乌斯反演】【杜教筛】【小学奥数】的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 更换笔记本电脑内置电池笔记本电脑如何换电
- 下一篇: 电脑中常用的软件_开发效率软件开发电脑软