1.YY的GCD
[sum_{x=1}^nsum_{y=1}^m[gcd(x,y)==P]
]
[sum_{p}sum_{x=1}^{n/p}sum_{y=1}^{m/p}[gcd(x,y)==1]
]
[sum_psum_{d=1}^{n/p}mu(d)frac{n}{pd}frac{m}{pd}
]
[sum_{T=1}^nfrac{n}{T}frac{m}{T}sum_{p|T}mu(T/p)
]
设(f[n]=sum_{p|n}mu(n/p)),从(f[i])转移到(f[i*pr[j]]).
(i\% pr[j])
(pr[j])作为p,贡献为(mu[i]).
(pr[j])作为n/p,贡献为(-f[i]).
(i\% pr[j]==0)
(pr[j])作为p,贡献为(mu[i]).
(pr[j])作为n/p,贡献为0.
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=1e7+50;
int n,m;
ll pr[N],v[N],f[N],mu[N];
inline void Pre(){
const int M=1e7;
mu[1]=1;
for(int i=2;i<=M;++i){
if(!v[i]) v[i]=1,pr[++pr[0]]=i,mu[i]=-1;
for(int j=1;j<=pr[0];++j){
if(i*pr[j]>M) break;
v[i*pr[j]]=1;
if(i%pr[j]!=0) mu[i*pr[j]]=-mu[i];
else{mu[i*pr[j]]=0;break;}
}
}
for(int i=1;i<=pr[0];++i) for(int j=1;j*pr[i]<=M;++j) f[j*pr[i]]+=mu[j];
for(int i=1;i<=M;++i) f[i]+=f[i-1];
}
void solve(){
scanf("%d%d",&n,&m);
if(n>m) swap(n,m);
ll ans=0;
for(int l=1,r;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
ans+=(f[r]-f[l-1])*(n/l)*(m/l);
}
printf("%lld
",ans);
}
int main(){
Pre();
int T;scanf("%d",&T);
while(T--) solve();
return 0;
}
2.数表
[sum_{i=1}^nsum_{j=1}^msigma(gcd(i,j))
]
[sum_{g=1}^nsigma(g)sum_{i=1}^{n/g}sum_{g=1}^{m/g}[gcd(i,j)==1]
]
[sum_{d=1}^nmu(d)sum_{g=1}^{n/d}sigma(g)frac{n}{gd}frac{m}{gd}
]
[sum_{T=1}^nfrac{n}{T}frac{m}{T}sum_{d|T}mu(d)sigma(T/d)
]
每次要求(sigma(x)<=a)
设(f[n]=sum_{d|n}mu(d)sigma(n/d)).
把所有二元组((sigma(n),n))升序,把询问升序.
增量构造,每次更新d的倍数.
所以总更新次数(nln)
因为需要查询区间和,单点修改,BIT维护.
#include<bits/stdc++.h>
using namespace std;
const int Query=2e4+50;
const int N=1e5+50;
const int M=1e5;
const int mod=(1ll<<31)-1;
int n,m;
int pri[N],v[N],mu[N],c[N],sum[N],ret[N];
inline void add(int p,int v){
for(;p<N;p+=p&-p) c[p]+=v;
}
inline int get(int p,int ans=0){
for(;p;p-=p&-p) ans+=c[p];
return ans;
}
struct Node{int n,m,a,id;}q[N];
bool cmp(Node a,Node b){return a.a<b.a;}
#define pr pair<int,int>
#define fir first
#define sec second
pr ord[N];
inline void pre(){
mu[1]=1;sum[1]=1;
for(int i=2;i<=M;++i){
if(!v[i]) v[i]=1,pri[++pri[0]]=i,mu[i]=-1,sum[i]=1+i;
for(int j=1;j<=pri[0];++j){
if(i*pri[j]>M) break;
v[i*pri[j]]=1;
if(i%pri[j]!=0) mu[i*pri[j]]=-mu[i],sum[i*pri[j]]=sum[i]*sum[pri[j]];
else {mu[i*pri[j]]=0;sum[i*pri[j]]=sum[i]*sum[pri[j]]-sum[i/pri[j]]*pri[j];break;}
}
}
for(int i=1;i<=M;++i) ord[i]=make_pair(sum[i],i);
sort(ord+1,ord+M+1);
}
int main(){
pre();
int T;scanf("%d",&T);
for(int i=1;i<=T;++i) q[i].id=i,scanf("%d%d%d",&q[i].n,&q[i].m,&q[i].a);
sort(q+1,q+T+1,cmp);
int it=1;
for(int i=1;i<=T;++i){
n=q[i].n,m=q[i].m;
while(it<=M&&ord[it].fir<=q[i].a){
for(int j=1;ord[it].sec*j<=M;++j) add(ord[it].sec*j,mu[j]*ord[it].fir);
it++;
}
if(n>m) swap(n,m);
int ans=0;
for(int l=1,r;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
ans+=(get(r)-get(l-1))*(n/l)*(m/l);
}
ret[q[i].id]=ans&mod;
}
for(int i=1;i<=T;++i) printf("%d
",ret[i]);
return 0;
}
3.DZY loves Math
由于一个神奇的性质,导致本来无法线性筛的本题异常简单.
套路不用再讲,不套路的地方又不会讲.
所以鸽掉了.
4.约数个数和
都做过一次还不会,太丢人了...
[sum_{i=1}^nsum_{j=1}^md(ij)
]
[d(ij)=sum_{x|i}sum_{y|j}[gcd(x,y)==1]
]
解释的话,因为i,j互质所以可以看作质因子的自由组合.
[sum_{i=1}^nsum_{j=1}^msum_{x|i}sum_{y|j}[gcd(x,y)==1]
]
[sum_{i=1}^nsum_{j=1}^m[gcd(i,j)==1]frac{n}{i}frac{m}{j}
]
[sum_{d=1}^nmu(d)sum_{i=1}^{n/d}frac{n}{id}sum_{j=1}^{m/d}frac{m}{jd}
]
到这里后发现按照原来的套路搞不动了.
[sum_{d=1}^nmu(d)sum_{i=1}^{n/d}frac{frac{n}{d}}{i}sum_{j=1}^{m/d}frac{frac{m}{d}}{j}
]
设(f[n]=sum_{i=1}^nfrac{n}{i})
[sum_{d=1}^nmu(d)f(n/d)f(m/d)
]
问题在求(f[n]).
[f[n]=sum_{i=1}^nd(i)
]
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=5e4+50;
ll n,m;
ll pr[N],v[N],mu[N],F[N];
inline void pre(){
const int M=5e4;
mu[1]=1;
for(int i=2;i<=M;++i){
if(!v[i]) pr[++pr[0]]=i,v[i]=1,mu[i]=-1;
for(int j=1;j<=pr[0];++j){
if(i*pr[j]>M) break;
v[i*pr[j]]=1;
if(i%pr[j]!=0) mu[i*pr[j]]=-mu[i];
else {mu[i*pr[j]]=0;break;}
}
}
for(int i=1;i<=M;++i) mu[i]+=mu[i-1];
for(int D=1;D<=M;++D)
for(int l=1,r;l<=D;l=r+1)
r=D/(D/l),F[D]+=(r-l+1)*(D/l);
}
inline void solve(){
scanf("%lld%lld",&n,&m);
if(n>m) swap(n,m);
ll ans=0;
for(ll l=1,r;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
ans+=1ll*(mu[r]-mu[l-1])*F[n/l]*F[m/l];
}
printf("%lld
",ans);
return ;
}
int main(){
pre();
int T;scanf("%d",&T);
while(T--) solve();
return 0;
}
5.数字表格
[prod_{i=1}^nprod_{j=1}^mf[gcd(i,j)]
]
[prod_{d=1}^nf[d]^{sum_{g=1}^{n/d}mu(g)frac{n}{gd}frac{m}{gd}}
]
[prod_{g=1}^nprod_{d=1}^{n/d}(f[d]^{mu(g)})^{frac{n}{gd}frac{m}{gd}}
]
[prod_{T=1}^nprod_{d|T}(f[d]^{mu(T/d)})^{frac{n}{T}frac{m}{T}}
]
处理出前缀积就行了.
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=1e6+50;
const int mod=1e9+7;
ll n,m;
ll fib[N],mu[N],F[N],pr[N],v[N];
ll mgml(ll a,ll b,ll ans=1){
b=(b%(mod-1)+mod-1)%(mod-1);
for(;b;b>>=1,a=a*a%mod) if(b&1) ans=ans*a%mod;
return ans;
}
inline void pre(){
const int M=1e6;
fib[0]=0,fib[1]=1,mu[1]=1;
for(int i=0;i<=M;++i) F[i]=1;
for(ll i=2;i<=M;++i){
fib[i]=(fib[i-1]+fib[i-2])%mod;
if(!v[i]) v[i]=1,pr[++pr[0]]=i,mu[i]=-1;
for(int j=1;j<=pr[0];++j){
if(i*pr[j]>M) break;
v[i*pr[j]]=1;
if(i%pr[j]!=0) mu[i*pr[j]]=-mu[i];
else {mu[i*pr[j]]=0;break;}
}
}
for(ll i=1;i<=M;++i){
ll st[3]={mgml(fib[i],mod-2),1,fib[i]};
for(ll j=1;i*j<=M;++j) F[i*j]=F[i*j]*st[mu[j]+1]%mod;
}
for(ll i=1;i<=M;++i) F[i]=F[i-1]*F[i]%mod;
}
inline void solve(){
scanf("%lld%lld",&n,&m);
if(n>m) swap(n,m);
ll ans=1;
for(ll l=1,r;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
ans=ans*mgml(F[r]*mgml(F[l-1],mod-2)%mod,(n/l)*(m/l))%mod;
}
printf("%lld
",ans);
return;
}
int main(){
pre();
int T;scanf("%d",&T);
while(T--) solve();
return 0;
}
6.于神之怒加强版
[sum_{i=1}^nsum_{j=1}^mgcd(i,j)^k
]
[sum_{g=1}^nmu(g)sum_{d=1}^{n/g}d^kfrac{n}{gd}frac{m}{gd}
]
[sum_{T=1}^nfrac{n}{T}frac{m}{T}sum_{g|T}mu(g)(frac{T}{g})^k
]
设(f[n]=sum_{g|n}mu(g)(frac{n}{g})^k).
由于是积性函数卷积的形式,所以依然是积性的.
(i\% pr[j])
[f[i*pr[j]]=f[i]*f[pr[j]]
]
(i\% pr[j]==0)
(pr[j])只能呆在(n/g)里,放在g里就会使(mu(g)=0).
[f[i*pr[j]]=f[i]*E[pr[j]]
]
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=5e6+50;
const int mod=1e9+7;
ll n,m,k;
ll v[N],pr[N],mu[N],F[N],E[N];
ll mgml(ll a,ll b,ll ans=1){
for(;b;b>>=1,a=a*a%mod) if(b&1) ans=ans*a%mod;
return ans;
}
inline void pre(){
const int M=5e6;
mu[1]=1;E[1]=1; F[1]=1;
for(int i=2;i<=M;++i){
if(!v[i]) v[i]=1,pr[++pr[0]]=i,mu[i]=-1,E[i]=mgml(i,k),F[i]=E[i]-1;
for(int j=1;j<=pr[0];++j){
if(i*pr[j]>M) break;
v[i*pr[j]]=1;
E[i*pr[j]]=E[i]*E[pr[j]]%mod;
if(i%pr[j]!=0) mu[i*pr[j]]=-mu[i],F[i*pr[j]]=F[i]*F[pr[j]]%mod;
else {mu[i*pr[j]]=0;F[i*pr[j]]=F[i]*E[pr[j]]%mod;break;}
}
}
// for(int i=1;i<=M;++i) printf("F[%d]=%lld
",i,F[i]);
for(int i=1;i<=M;++i) F[i] =(F[i]+F[i-1])%mod;
}
inline void solve(){
scanf("%lld%lld",&n,&m);
if(n>m) swap(n,m);
ll ans=0;
for(ll l=1,r;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
ans=(ans+1ll*(F[r]-F[l-1]+mod)*(n/l)%mod*(m/l)%mod)%mod;
}
printf("%lld
",ans);
return;
}
int main(){
int T;scanf("%d%lld",&T,&k);
pre();
while(T--) solve();
return 0;
}
7.jzptab
[sum_{i=1}^nsum_{j=1}^mlcm(i,j)
]
[sum_{i=1}^nsum_{j=1}^mfrac{ij}{gcd(i,j)}
]
[sum_{d=1}^ndsum_{i=1}^{n/d}isum_{j=1}^{m/d}j[gcd(i,j)==1]
]
[sum_{T=1}^nsum_{i=1}^{n/T}isum_{j=1}^{m/T}j Tsum_{g|T}mu(g)g
]
设(f[n]=sum_{g|n}mu(g)g)
(f=mu * Id)
(i\% pr[j])
(f[i*pr[j]]=f[i]*f[pr[j]])
(i\% pr[j]==0)
(pr[j])只能在(n/g)里.
(f[i*pr[j]]=f[i])
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=1e7+50;
const int mod=1e8+9;
int n,m;
int pr[N],v[N],mu[N],f[N];
inline void pre(){
mu[1]=1;
f[1]=1;
const int M=1e7;
for(int i=2;i<=M;++i){
if(!v[i]) v[i]=1,pr[++pr[0]]=i,mu[i]=-1,f[i]=1-i;
for(int j=1;j<=pr[0];++j){
if(i*pr[j]>M) break;
v[i*pr[j]]=1;
if(i%pr[j]!=0) mu[i*pr[j]]=-mu[i],f[i*pr[j]]=f[pr[j]]*f[i];
else {mu[i*pr[j]]=0;f[i*pr[j]]=f[i];break;}
}
}
for(int i=1;i<=M;++i) f[i]=1ll*i*f[i]%mod;
for(int i=1;i<=M;++i) f[i]=(f[i]+f[i-1])%mod;
}
inline int calc(int D){
return 1ll*((1ll+(n/D))*(n/D)/2%mod)*((1ll+(m/D))*(m/D)/2%mod)%mod;
}
inline void solve(){
scanf("%lld%lld",&n,&m);
if(n>m) swap(n,m);
int ans=0;
for(int l=1,r;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
ans=(ans+1ll*calc(l)*(f[r]-f[l-1])%mod+mod)%mod;
}
printf("%lld
",(ans+mod)%mod);
}
signed main(){
pre();
int T;scanf("%lld",&T);
while(T--) solve();
return 0;
}
8.一个人的数论
再看这道题就感觉简单多了.
这不是你颓题解的理由
[sum_{i=1}^n[gcd(i,n)==1]i^k
]
[sum_{i=1}^ni^ksum_{d|i&&d|n}mu(d)
]
[sum_{d|n}mu(d)d^ksum_{i=1}^{n/d}i^k
]
[frac{1}{k}sum_{i=0}^kC(k+1,i)B_isum_{d|n}mu(d)d^k(frac{n}{d})^{k+1-i}
]
[frac{n^k}{k}sum_{i=0}^kC(k+1,i)B_isum_{d|n}mu(d)(frac{n}{d})^{1-i}
]
(mu*id^{1-i})是积性的.
[ans=frac{n^k}{k}sum_{i=0}^kC(k+1,i)B_iprod_{j=1}^{pr[0]}((pr[j]^{c[j]})^{1-i}-(pr[j]^{c[j]-1})^{1-i})
]
#include<bits/stdc++.h>
#define ll long long
#define int long long
using namespace std;
const int mod=1e9+7;
int k,w,p[1500],c[1500];
ll a[105][105];
ll mgml(ll a,ll b,ll ans=1){
b=(b%(mod-1)+mod-1)%(mod-1);
for(;b;b>>=1,a=a*a%mod) if(b&1) ans=ans*a%mod;
return ans;
}
int L;
inline void gauss(){
for(int i=1;i<=L;++i){
int maxi=i;
for(int j=i+1;j<=L;++j) if(a[j][i]>a[maxi][i]) maxi=j;
swap(a[maxi],a[i]);
if(!a[i]) continue;
ll rate=mgml(a[i][i],mod-2);
for(int j=1;j<=L+1;++j) a[i][j]=a[i][j]*rate%mod;
for(int j=1;j<=L;++j){
if(j==i) continue;
ll rate=a[j][i];
for(int d=1;d<=L+1;++d) a[j][d]=(a[j][d]-rate*a[i][d]%mod+mod)%mod;
}
}
}
inline int rd(register int x=0,register char ch=getchar(),register int f=1){
while(!isdigit(ch)) f=ch=='-'?-1:1,ch=getchar();
while(isdigit(ch)) x=(x<<1)+(x<<3)+ch-48,ch=getchar();
return x*f;
}
signed main(){
k=rd(); w=rd(); L=k+1;
for(int i=1;i<=w;++i) p[i]=rd(),c[i]=rd();
ll pre=0;
for(ll i=1;i<=k+1;++i){
pre=(pre+mgml(i,k))%mod;
a[i][L+1]= pre;
for(ll j=1,powi=i;j<=L;++j,powi=1ll*powi*i%mod) a[i][j]=powi;
}
gauss();
ll ans=0;
for(int i=1;i<=L;++i){
ll tmp=1;
for(int j=1;j<=w;++j){
tmp=tmp*mgml(p[j],i*c[j])%mod*(1-mgml(p[j],k-i)+mod)%mod;
}
ans=(ans+a[i][L+1]*tmp%mod+mod)%mod;
}
printf("%lld
",ans);
return 0;
}