Bzoj3512 DZY Loves Math IV

Time Limit: 15 Sec  Memory Limit: 128 MB
Submit: 447  Solved: 223

Description

给定n,m,求 模10^9+7的值。

Input

仅一行,两个整数n,m。

Output

仅一行答案。

Sample Input

100000 1000000000

Sample Output

857275582
数据规模:
1<=n<=10^5,1<=m<=10^9,本题共4组数据。

HINT

 

Source

By Jc 鸣谢杜教

数学问题 莫比乌斯反演 杜教筛

数学题的麻烦在于,如果你拖久了不更博,就不得不再推一遍式子了。


注意到n比较小,只有1e5,这是一个可以枚举的范围,我们来看看如果枚举n会怎么样:
$S(n,m)=sum_{i=1}^{m} varphi(i*n)$
有一个结论:
若$|mu(n)|=1$,即n是无平方数,则:
$$S(n,m)=sum_{i=1}^{m} varphi(i) sum_{d|gcd(i,n)} varphi(frac{n}{d})$$
然而博主并不知道如何能想到这个结论。手推一下发现它是对的。
http://www.cnblogs.com/ziliuziliu/p/6508562.html
↑这里有一个关于它的神奇证明

于是当n是无平方数的时候可以这么求:
$$S(n,m)=sum_{d|n} varphi(frac{n}{d}) sum_{i=1}^{lfloor frac{m}{d} floor} varphi(i*d)$$
$$S(n,m)=sum_{d|n} varphi(frac{n}{d}) S(d,lfloor frac{m}{d} floor)$$
记忆化一下即可。

而当$mu(n)=0$的时候就很方便了:
我们知道

$varphi(n)= n*Pi_{i=1}^{k} (1-(1/p_k)) $

其中p数组是n的全部质因数。我们已经可以算出每个质因数都只出现一次时的答案,那么对于剩下的质因数,直接乘进去就行。
也就是说,找到最大的x使得$|mu(frac{n}{x})|=1$,
$$S(n,m)=x*S(frac{n}{x},m)$$
后面这个S用最上面的方法求。
递归边界是$S(1,m)$,也就是 $mu $ 的前缀和,这个可以用杜教筛算。
最后,$ans=sum_{i=1}^{n}S(n,m)$

 1 #include<iostream>
 2 #include<algorithm>
 3 #include<cstring>
 4 #include<cstdio>
 5 #include<cmath>
 6 #include<map>
 7 #define LL long long
 8 using namespace std;
 9 const int mod=1e9+7;
10 const int mxn=5000010;
11 const int inv2=500000004;
12 int read(){
13     int x=0,f=1;char ch=getchar();
14     while(ch<'0' || ch>'9'){if(ch=='-')f=-1;ch=getchar();}
15     while(ch>='0' && ch<='9'){x=x*10+ch-'0';ch=getchar();}
16     return x*f;
17 }
18 int pri[mxn>>1],cnt=0;
19 int phi[mxn],mx[mxn];
20 int smm[mxn];
21 void init(){
22     phi[1]=1;mx[1]=1;
23     for(int i=2;i<mxn;i++){
24         if(!phi[i]){
25             pri[++cnt]=i;
26             phi[i]=i-1;
27             mx[i]=1;
28         }
29         for(int j=1;j<=cnt && pri[j]*i<mxn;j++){
30             int tmp=pri[j]*i;
31             if(i%pri[j]==0){
32                 mx[tmp]=mx[i]*pri[j];
33                 phi[tmp]=phi[i]*pri[j];
34                 break;
35             }
36             phi[tmp]=phi[i]*(pri[j]-1);
37             mx[tmp]=mx[i];
38         }
39     }
40     for(int i=1;i<mxn;i++)
41         smm[i]=((LL)smm[i-1]+phi[i])%mod;
42     return;
43 }
44 map<int,int>mphi;
45 int Calc_phi(int x){
46     if(x<mxn)return smm[x];
47     if(mphi.count(x))return mphi[x];
48     int res=(LL)x*(x+1)%mod*inv2%mod;
49     for(int i=2,pos;i<=x;i=pos+1){
50         int y=x/i; pos=x/y;
51         res=((LL)res-(pos-i+1)*Calc_phi(y))%mod;
52     }
53     res=(res+mod)%mod;
54     mphi[x]=res;
55     return res;
56 }
57 map<LL,int>mpS;
58 int solve(int n,int m){
59     if(!m)return 0;
60     if(n==1)return Calc_phi(m);
61     LL tar=(LL)n*mod+m;
62     if(mpS.count(tar))return mpS[tar];
63     int res=0;
64     for(int i=1;i*i<=n;i++){
65         if(n%i==0){
66             res=((LL)res+phi[n/i]*(LL)solve(i,m/i))%mod;
67             if(i*i!=n)res=((LL)res+phi[i]*(LL)solve(n/i,m/(n/i)))%mod;
68         }
69     }
70     mpS[tar]=res;
71     return res;
72 }
73 int n,m;
74 int main(){
75     int i,j;
76     init();
77     n=read();m=read();
78     int ans=0;
79     for(i=1;i<=n;i++)
80         ans=(ans+(LL)mx[i]*solve(i/mx[i],m))%mod;
81     ans=(ans+mod)%mod;
82     printf("%d
",ans);
83     return 0;
84 }
原文地址:https://www.cnblogs.com/SilverNebula/p/7020528.html