您的位置:首页 > 其它

BZOJ 3512 DZY Loves Math IV

2017-03-06 10:01 417 查看
题解网上都有。。。那我写下这个式子的证明吧。。。

phi(n*k)=Σ(d|i,d|k)phi(n/d)*phi(k) (|miu
|==1)

既然miu
==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;
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: