BZOJ4916: 神犇和蒟蒻(杜教筛)

Description

很久很久以前，有一只神犇叫yzy;

很久很久之后，有一只蒟蒻叫lty;

Input

请你读入一个整数N;1<=N<=1E9,A、B模1E9+7;

Output

请你输出一个整数$A=sum_{i=1}^N{mu (i^2)}$

请你输出一个整数$B=sum_{i=1}^N{varphi (i^2)}$

Sample Input

Sample Output

1
1

解题思路：

第一问考察$mu$的定义所以答案为1。

第二问才是重头戏，主要应用性质$phi(i^2)=iphi(i)$

已知$sum_{i|N}phi(i)=N$

很显然要用杜教筛，这里的g最好用id()

根据杜教筛算法：设f(x)=$xphi(x)$

$S(N)=(f*g)(N)-sum_{i=2}^{N}g(i)*S(left lfloor frac{N}{i} ight floor)$

前面那个就等于$sum_{i=1}^{N}i^2=frac{n(n+1)(2n+1)}{6}$

后面就是$iS(leftlfloorfrac N i ight floor)$

然后就可以愉快的结束了。

代码：

 1 #include<map>
 2 #include<cstdio>
 3 #include<cstring>
 4 #include<algorithm>
 5 typedef long long lnt;
 6 const lnt mod=(lnt)(1e9+7);
 7 const int N=4000000;
 8 struct JDr_map{
 9     std::map<int,lnt>A;
10     lnt a[N];
11     void Insert(int pos,lnt v)
12     {
13         if(pos<N)
14             a[pos]=v;
15         else
16             A[pos]=v;
17         return ;
18     }
19     bool hav(int pos)
20     {
21         if(pos<N)
22             return true;
23         return A.find(pos)!=A.end();
24     }
25     lnt val(int pos)
26     {
27         if(pos<N)
28             return a[pos];
29         return A[pos];
30     }
31 }s;
32 int prime[N];
33 int phi[N];
34 bool vis[N];
35 int cnt;
36 void gtp(void)
37 {
38     phi[1]=1;
39     for(int i=2;i<N;i++)
40     {
41         if(!vis[i])
42         {
43             prime[++cnt]=i;
44             phi[i]=i-1;
45         }
46         for(int j=1;j<=cnt&&prime[j]*i<N;j++)
47         {
48             vis[i*prime[j]]=true;
49             if(i%prime[j]==0)
50             {
51                 phi[i*prime[j]]=phi[i]*prime[j];
52                 break;
53             }
54             phi[i*prime[j]]=phi[i]*(prime[j]-1);
55         }
56     }
57     s.Insert(0,0);
58     for(int i=1;i<N;i++)
59         s.Insert(i,(s.val(i-1)+1ll*i*phi[i])%mod);
60     return ;
61 }
62 lnt SLP(int pos)
63 {
64     if(s.hav(pos))
65         return s.val(pos);
66     lnt tmp=0;
67     for(int i=2,j;i<=pos;i=j+1)
68     {
69         j=pos/(pos/i);
70         (tmp+=(1ll*(j-i+1)*(j+i)/2)%mod*SLP(pos/i))%=mod;
71     }
72     tmp=((1ll*(pos)*(2ll*pos+1)%mod*(pos+1))%mod*166666668ll%mod-tmp)%mod;
73     s.Insert(pos,tmp);
74     return tmp;
75 }
76 int main()
77 {
78     int n;
79     gtp();
80     scanf("%d",&n);
81     printf("%d
%lld
",1,(SLP(n)%mod+mod)%mod);
82     return 0;
83 }