DZY Loves Math系列
好久没写数学题了，再这样下去吃枣药丸啊。
找一套应该还比较有意思的数学题来做。
[bzoj3309]DZY Loves Math
简单推一下。
\[\sum_{i=1}^n\sum_{j=1}^mf(\gcd(i,j))\\=\sum_{d=1}^nf(d)\sum_{i=1}^{n/d}\mu(i)\frac n{id}\frac m{id}\\=\sum_{T=1}^n\frac nT\frac mT\sum_{d|T}f(d)\mu(\frac Td)
\]
令\(h(T)=\sum_{d|T}f(d)\mu(\frac Td)\)，我们现在需要\(h\)的前缀和。
随便分析一下可知\(h((p_1p_2...p_k)^t)=(-1)^{k+1}\)否则都是\(0\)。具体过程见这篇博客。
线性筛。求每个数的最小质因子的幂次，判断是否幂次相等即可。
#include<cstdio>
#include<algorithm>
using namespace std;
int gi(){
int x=0,w=1;char ch=getchar();
while ((ch<'0'||ch>'9')&&ch!='-') ch=getchar();
if (ch=='-') w=0,ch=getchar();
while (ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+ch-'0',ch=getchar();
return w?x:-x;
}
const int N = 1e7+5;
int zhi[N],pri[N],tot,a[N],low[N],h[N];
int main(){
for (int i=2;i<N;++i){
if (!zhi[i]) pri[++tot]=i,a[i]=h[i]=1,low[i]=i;
for (int j=1;i*pri[j]<N;++j){
zhi[i*pri[j]]=1;
if (i%pri[j]==0){
a[i*pri[j]]=a[i]+1;
low[i*pri[j]]=low[i]*pri[j];
if (i==low[i]) h[i*pri[j]]=1;
else h[i*pri[j]]=a[i/low[i]]==a[low[i]*pri[j]]?-h[i/low[i]]:0;
break;
}
a[i*pri[j]]=1;
low[i*pri[j]]=pri[j];
h[i*pri[j]]=a[i]==1?-h[i]:0;
}
}
for (int i=1;i<N;++i) h[i]+=h[i-1];
int T=gi();while(T--){
int n=gi(),m=gi();if(n>m)swap(n,m);
long long ans=0;
for (int i=1,j;i<=n;i=j+1)
j=min(n/(n/i),m/(m/i)),ans+=1ll*(h[j]-h[i-1])*(n/i)*(m/i);
printf("%lld\n",ans);
}
return 0;
}
[bzoj3462]DZY Loves Math II
先把\(S\)质因数分解，如果有质因子出现超过一次则无解。
假设分解出来有\(k\)个质因子\(p_1,p_2...p_k\)，我们现在是要用这些质数拼凑出\(n\)，也就是求\(\prod_{i=1}^k\frac{1}{1-x^{p_i}}\)的\(x^n\)项系数。
每个质因子在最终方案里的出现次数一定是\(x\frac{S}{p_i}+y\)的形式，其中\(0\le y<\frac{S}{p_i}\)。
把\(x\)跟\(y\)单独考虑。考虑\(x\)，相当于是往\(k\)个盒子里面放若干个球，可以隔板法计算方案。组合数的\(n\)可能会非常大，但考虑到\(k\)非常小所以可以\(O(k)\)地计算组合数。而考虑\(y\)，发现总值域不超过\(kS\)，所以做个多重背包统计一下方案数即可。
复杂度\(O(\sqrt S+Sk^2+qk^2)\)。
#include<cstdio>
#include<algorithm>
#include<cstring>
using namespace std;
#define ll long long
ll gi(){
ll x=0,w=1;char ch=getchar();
while ((ch<'0'||ch>'9')&&ch!='-') ch=getchar();
if (ch=='-') w=0,ch=getchar();
while (ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+ch-'0',ch=getchar();
return w?x:-x;
}
const int mod = 1e9+7;
int S,q,x,p[10],t,f[2][14000005],all,sum,inv[10],jcn;
void add(int &x,int y){x+=y;x>=mod?x-=mod:x;}
int cal(ll x){
x%=mod;int res=1;
for (int i=1;i<t;++i) res=1ll*res*(x+i)%mod;
return res;
}
int main(){
S=gi();q=gi();x=S;
for(int i=2;i*i<=x;++i)while(x%i==0)p[++t]=i,x/=i;
if(x>1)p[++t]=x;
for(int i=1;i<t;++i)if(p[i]==p[i+1]){while(q--)puts("0");return 0;}
f[0][0]=1;
for(int i=1,now=1,lst=0;i<=t;++i,lst=now,now^=1){
memcpy(f[now],f[lst],sizeof(f[now]));
for(int j=0;j<=i*S;++j){
if(j>=p[i])add(f[now][j],f[now][j-p[i]]);
if(j>=S)add(f[now][j],mod-f[lst][j-S]);
}
}
for(int i=1;i<=t;++i)sum+=p[i];
inv[1]=jcn=1;
for(int i=2;i<t;++i)inv[i]=1ll*inv[mod%i]*(mod-mod/i)%mod,jcn=1ll*jcn*inv[i]%mod;
while(q--){
ll n=gi()-sum;if(n<0){puts("0");continue;}int ans=0;
for(int i=n%S;i<=t*S&&i<=n;i+=S)
ans=(ans+1ll*f[t&1][i]*cal((n-i)/S))%mod;
printf("%lld\n",1ll*ans*jcn%mod);
}
return 0;
}
[bzoj3481]DZY Loves Math III
枚举\(x\)，统计有多少个\(y\)满足条件，显然就是\(\sum_{x=0}^{P-1}\gcd(x,P)[\gcd(x,P)|Q]\)。
令\(D=\gcd(P,Q)\)，上式改为枚举\(\gcd(x,P)\)，化成\(\sum_{d|D}d\varphi(\frac Pd)\)。
分别考虑每个质数的贡献。只要看\(D\)中该质数的幂次与\(P\)中是否相等就可以了。
#include<cstdio>
#include<algorithm>
#include<ctime>
#include<map>
using namespace std;
#define ll long long
ll gi(){
ll x=0,w=1;char ch=getchar();
while ((ch<'0'||ch>'9')&&ch!='-') ch=getchar();
if (ch=='-') w=0,ch=getchar();
while (ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+ch-'0',ch=getchar();
return w?x:-x;
}
ll mul(ll x,ll y,ll mod){
return (x*y-(ll)((long double)x/mod*y+0.5)*mod+mod)%mod;
}
ll fastpow(ll x,ll y,ll mod){
ll res=1;
while (y) {if (y&1) res=mul(res,x,mod);x=mul(x,x,mod);y>>=1;}
return res;
}
ll f[]={2,3,5,7,11,13,17,19,23,29};
bool Miller_Rabin(ll p){
for (int i=0;i<10;++i){
if (p<=f[i]) break;
if (fastpow(f[i],p-1,p)!=1) return false;
ll pp=p-1;
while (~pp&1){
pp>>=1;ll y=fastpow(f[i],pp,p);
if (mul(y,y,p)==1&&y!=1&&y!=p-1) return false;
}
}
return true;
}
ll Pollard_Rho(ll n){
ll x=0,y=0,t=1,q=1,a=1+rand()*rand()%(n-1);
for (int k=2;true;k<<=1,y=x,q=1){
for (int i=1;i<k;++i){
x=(mul(x,x,n)+a)%n;
q=mul(q,abs(x-y),n);
if (!(i&127)){
t=__gcd(q,n);
if (t>1) return t;
}
}
if (t>1||(t=__gcd(q,n))>1) break;
}
if (t==n) t=1;return t;
}
map<ll,ll>P,Q,D;map<ll,ll>::iterator it;
const ll mod = 1e9+7;ll fg,ans=1;
void fact(ll n,map<ll,ll>&M){
if (n==0) {fg=1;return;}if (n==1) return;
if (Miller_Rabin(n)) {++M[n];return;}
ll p=n;while (p==n) p=Pollard_Rho(p);
fact(p,M);fact(n/p,M);
}
int main(){
int n=gi();srand(20020415);
for (int i=1;i<=n;++i) fact(gi(),P);
for (int i=1;i<=n;++i) fact(gi(),Q);
if (fg) Q=P;
for (it=Q.begin();it!=Q.end();++it){
ll p=(*it).first;
if (P.find(p)!=P.end()) D[p]=min(P[p],Q[p]);
}
for (it=P.begin();it!=P.end();++it){
ll p=(*it).first,res;
if (D[p]<P[p]) res=(p-1)%mod*(D[p]+1)%mod*fastpow(p,P[p]-1,mod)%mod;
else res=((p-1)%mod*D[p]+p)%mod*fastpow(p,P[p]-1,mod)%mod;
ans=ans*res%mod;
}
printf("%lld\n",ans);return 0;
}
[bzoj3512]DZY Loves Math IV
我们令\(n=xy\)，其中\(x\)含有\(n\)中出现的每个质因子恰好一次。
\[\varphi(ni)=y\varphi(xi)\\=y\varphi(i)\varphi(\frac{x}{\gcd(i,x)})\gcd(i,x)
\]
利用\(\sum_{d|n}\varphi(d)=n\)
\[=y\varphi(i)\varphi(\frac{x}{\gcd(i,x)})\sum_{d|i,d|x}\varphi(d)\\=y\varphi(i)\sum_{d|i,d|x}\varphi(\frac xd)
\]
令\(S(n,m)=\sum_{i=1}^m\varphi(ni)\)
\[S(n,m)=\sum_{i=1}^my\varphi(i)\sum_{d|i,d|x}\varphi(\frac xd)\\=y\sum_{d|x}\varphi(\frac xd)\sum_{i=1}^{m/d}\varphi(di)\\=y\sum_{d|x}\varphi(\frac xd)S(d,\lfloor\frac md\rfloor)
\]
当\(m=1\)时就是求\(\varphi\)的前缀和，直接杜教筛，剩下的记忆化。
#include<cstdio>
#include<algorithm>
#include<map>
using namespace std;
int gi(){
int x=0,w=1;char ch=getchar();
while ((ch<'0'||ch>'9')&&ch!='-') ch=getchar();
if (ch=='-') w=0,ch=getchar();
while (ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+ch-'0',ch=getchar();
return w?x:-x;
}
const int N = 2e6+5;
const int mod = 1e9+7;
int zhi[N],pri[N],tot,phi[N],low[N],sum[N],ans;
map<int,int>M[N],P;
int S(int n){
if (n<N) return sum[n];
if (P[n]) return P[n];
int res=0;
for (int i=2,j;i<=n;i=j+1)
j=n/(n/i),res=(res+1ll*(j-i+1)*S(n/i))%mod;
return P[n]=((1ll*n*(n+1)>>1)-res+mod)%mod;
}
int S(int n,int m){
if (m==0) return 0;
if (n==1) return S(m);
if (M[n][m]) return M[n][m];
int res=0;
for (int i=1;i*i<=n;++i)
if (n%i==0){
res=(res+1ll*phi[n/i]*S(i,m/i))%mod;
if (i*i<n) res=(res+1ll*phi[i]*S(n/i,m/(n/i)))%mod;
}
return M[n][m]=res;
}
int main(){
int n=gi(),m=gi();phi[1]=low[1]=1;
for (int i=2;i<N;++i){
if (!zhi[i]) pri[++tot]=i,phi[i]=i-1,low[i]=i;
for (int j=1;i*pri[j]<N;++j){
zhi[i*pri[j]]=1;
if (i%pri[j]==0){
phi[i*pri[j]]=phi[i]*pri[j];
low[i*pri[j]]=low[i];
break;
}
phi[i*pri[j]]=phi[i]*(pri[j]-1);
low[i*pri[j]]=low[i]*pri[j];
}
}
for (int i=1;i<N;++i) sum[i]=(sum[i-1]+phi[i])%mod;
for (int i=1;i<=n;++i) ans=(ans+1ll*(i/low[i])*S(low[i],m))%mod;
printf("%d\n",ans);return 0;
}
[bzoj3560]DZY Loves Math V
考虑每个质数的贡献。对于一个质数\(p\)，记\(X=\prod_{i=1}^n\sum_{j=0}^{k_i}p^j\)，那么它对这个答案的贡献就是\(\frac{(X-1)(p-1)}{p}+1\)。
#include<cstdio>
#include<algorithm>
using namespace std;
int gi(){
int x=0,w=1;char ch=getchar();
while ((ch<'0'||ch>'9')&&ch!='-') ch=getchar();
if (ch=='-') w=0,ch=getchar();
while (ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+ch-'0',ch=getchar();
return w?x:-x;
}
const int N = 1e7+5;
const int mod = 1e9+7;
int zhi[N],pri[N],tot,d[N],inv[N],pro[N],n,ans=1;
int main(){
inv[1]=1;
for (int i=2;i<N;++i){
if (!zhi[i]) pri[++tot]=i,d[i]=i;
for (int j=1;i*pri[j]<N;++j){
zhi[i*pri[j]]=1;d[i*pri[j]]=pri[j];
if (i%pri[j]==0) break;
}
inv[i]=1ll*inv[mod%i]*(mod-mod/i)%mod;
pro[i]=1;
}
n=gi();
for (int i=1;i<=n;++i){
int x=gi();
while(x>1){
int p=d[x],s=1;
while (x%p==0) x/=p,s=(1ll*s*p+1)%mod;
pro[p]=1ll*pro[p]*s%mod;
}
}
for (int i=1;i<=tot;++i){
int p=pri[i];
ans=(1ll*(pro[p]-1)*(p-1)%mod*inv[p]+1)%mod*ans%mod;
}
printf("%d\n",ans);return 0;
}
[bzoj3561]DZY Loves Math VI
\[\sum_{i=1}^n\sum_{j=1}^m(\frac{ij}{\gcd(i,j)})^{\gcd(i,j)}\\=\sum_{d=1}^n\frac{1}{d^d}\sum_{i=1}^{n/d}\mu(i)(id)^{2d}S(\frac n{id},d)S(\frac m{id},d)\\=\sum_{d=1}^nd^d\sum_{i=1}^{n/d}\mu(i)i^{2d}S(\frac n{id},d)S(\frac m{id},d)
\]
其中\(S(n,m)=\sum_{i=1}^ni^m\)
使用一些高妙的方法可以把复杂度做到\(O(n\log n)\)。
#include<cstdio>
#include<algorithm>
using namespace std;
int gi(){
int x=0,w=1;char ch=getchar();
while ((ch<'0'||ch>'9')&&ch!='-') ch=getchar();
if (ch=='-') w=0,ch=getchar();
while (ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+ch-'0',ch=getchar();
return w?x:-x;
}
const int N = 5e5+5;
const int mod = 1e9+7;
int n,m,pw[N],zhi[N],pri[N],tot,mu[N],sum[N],ans;
int fastpow(int a,int b){
int res=1;
while (b) {if (b&1) res=1ll*res*a%mod;a=1ll*a*a%mod;b>>=1;}
return res;
}
int main(){
n=gi();m=gi();if(n<m)swap(n,m);
for (int i=1;i<=n;++i) pw[i]=1;
mu[1]=1;
for (int i=2;i<=n;++i){
if (!zhi[i]) pri[++tot]=i,mu[i]=mod-1;
for (int j=1;i*pri[j]<=n;++j){
zhi[i*pri[j]]=1;
if (i%pri[j]==0) break;
mu[i*pri[j]]=mod-mu[i];
}
}
for (int d=1;d<=n;++d){
int res=0;
for (int i=1;i<=n/d;++i) pw[i]=1ll*pw[i]*i%mod,sum[i]=(sum[i-1]+pw[i])%mod;
for (int i=1;i<=n/d;++i) res=(res+1ll*mu[i]*pw[i]%mod*pw[i]%mod*sum[n/d/i]%mod*sum[m/d/i])%mod;
ans=(ans+1ll*res*fastpow(d,d))%mod;
}
printf("%d\n",ans);return 0;
}
