DZY Loves Math系列

Posted 2021-01-30 zhoushuyu

link

好久没写数学题了，再这样下去吃枣药丸啊。

找一套应该还比较有意思的数学题来做。

[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;
}