好久没写数学题了，再这样下去吃枣药丸啊。
找一套应该还比较有意思的数学题来做。
[bzoj3309]DZY Loves Math
简单推一下。
[sum_{i=1}^nsum_{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}^nfrac nTfrac mTsum_{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
",ans);
}
return 0;
}
[bzoj3462]DZY Loves Math II
先把(S)质因数分解，如果有质因子出现超过一次则无解。
假设分解出来有(k)个质因子(p_1,p_2...p_k)，我们现在是要用这些质数拼凑出(n)，也就是求(prod_{i=1}^kfrac{1}{1-x^{p_i}})的(x^n)项系数。
每个质因子在最终方案里的出现次数一定是(xfrac{S}{p_i}+y)的形式，其中(0le 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
",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}dvarphi(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
",ans);return 0;
}
[bzoj3512]DZY Loves Math IV
我们令(n=xy)，其中(x)含有(n)中出现的每个质因子恰好一次。
[varphi(ni)=yvarphi(xi)\=yvarphi(i)varphi(frac{x}{gcd(i,x)})gcd(i,x)
]
利用(sum_{d|n}varphi(d)=n)
[=yvarphi(i)varphi(frac{x}{gcd(i,x)})sum_{d|i,d|x}varphi(d)\=yvarphi(i)sum_{d|i,d|x}varphi(frac xd)
]
令(S(n,m)=sum_{i=1}^mvarphi(ni))
[S(n,m)=sum_{i=1}^myvarphi(i)sum_{d|i,d|x}varphi(frac xd)\=ysum_{d|x}varphi(frac xd)sum_{i=1}^{m/d}varphi(di)\=ysum_{d|x}varphi(frac xd)S(d,lfloorfrac md
floor)
]
当(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
",ans);return 0;
}
[bzoj3560]DZY Loves Math V
考虑每个质数的贡献。对于一个质数(p)，记(X=prod_{i=1}^nsum_{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
",ans);return 0;
}
[bzoj3561]DZY Loves Math VI
[sum_{i=1}^nsum_{j=1}^m(frac{ij}{gcd(i,j)})^{gcd(i,j)}\=sum_{d=1}^nfrac{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^dsum_{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(nlog 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
",ans);return 0;
}