BZOJ 3512 DZY Loves Math IV

题解网上都有。。。那我写下这个式子的证明吧。。。

phi(n*k)=Σ(d|i,d|k)phi(n/d)*phi(k) （|miu[n]|==1）

既然miu[n]==1，那么n的每个质因子的次数都是1，d也是这样。

令n*k=∏pi，gcd(n,k)=∏qi，那么phi(n*k)=n*k*∏(1-1/pi)。

而原式=phi(n/d)*phi(k)=(n/d)*∏(1-1/pi)（pi为只在n中出现的质因数）*∏(1-1/pi)（属于(n,k)的n/d的质因数）*k*∏(1-1/pi)（k中的质因数）。

        =n*k*∏pi（n*k中所有的质因数）/d*∏(1-1/pi)（属于(n,k)的n/d的质因数）。

前面那一坨在结论中出现，我们现在只需要证Σ后面那一坨=1。

重新叙述一下，即是证对于一个集合s，里面的每一个元素都可以表示成1/pi或（1-1/pi），求所有集合的∏的Σ。

这么想，有n个事件，每个事件有1/pi的概率发生，求所有情况的概率的和。

那不就是1吗。

就好了。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<map>
#include<cmath>
#define maxn 5000050
#define mod 1000000007
#define inv 500000004
using namespace std;
int n,m,prime[maxn],tot=0,phi[maxn],phd[maxn],ans=0,q[maxn];
bool vis[maxn];
map <long long,int> mp1,mp2;
void get_table()
{
    phi[1]=1;q[1]=1;
    for (int i=2;i<=maxn-50;i++)
    {
        if (!vis[i])
        {
            vis[i]=true;q[i]=i;
            prime[++tot]=i;phi[i]=i-1;
        }
        for (int j=1;j<=tot && i*prime[j]<=maxn-50;j++)
        {
            vis[i*prime[j]]=true;
            if (i%prime[j]) {phi[i*prime[j]]=phi[i]*(prime[j]-1);q[i*prime[j]]=q[i]*prime[j];}
            else {phi[i*prime[j]]=phi[i]*prime[j];q[i*prime[j]]=q[i];break;}
        }
    }
    for (int i=1;i<=maxn-50;i++) phd[i]=(phi[i]+phd[i-1])%mod;
}
int ask2(int m)
{
    if (m<=maxn-50) return phd[m];
    if (mp2.find(m)!=mp2.end()) return mp2[m];
    int ret=(long long)m*(m+1)%mod*inv%mod,l=2,r;
    while (l<=m)
    {
        r=m/(m/l);
        ret=(ret-(long long)ask2(m/l)*(r-l+1)%mod+mod)%mod;
        l=r+1;
    }
    mp2[m]=ret;return ret;
}
int ask1(int m,int n)
{
      if (!m) return 0;
      long long now=(long long)n*mod+m;
    if (mp1.find(now)!=mp1.end()) return mp1[now];
    int ret=0;
    if (n==1) ret=ask2(m);
    else
    {
        int top=sqrt(n);
        for (int i=1;i<=top;i++)
        {
             if (n%i) continue;
             ret=(ret+(long long)phi[n/i]*ask1(m/i,i)%mod)%mod;
            if ((i!=top) || (top*top!=n)) ret=(ret+(long long)phi[i]*ask1(m/(n/i),n/i)%mod)%mod;
        }
    }
    mp1[now]=ret;return ret;
}
int main()
{
    scanf("%d%d",&n,&m);
    get_table();
    for (int i=1;i<=n;i++)
        ans=(ans+(long long)ask1(m,q[i])*i/q[i]%mod)%mod;
    printf("%d\n",ans);
    return 0;
}