BZOJ 4916 神犇和蒟蒻(杜教筛)

Description

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

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

Input

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

Output

请你输出一个整数A=∑Ni=1μ(i2)A=∑i=1Nμ(i2);

请你输出一个整数B=∑Ni=1φ(i2)B=∑i=1Nφ(i2);

Sample Input

1

Sample Output

1

1

Solution

A=∑Ni=1μ(i2)A=∑i=1Nμ(i2)是拿来搞笑的。

B=∑Ni=1φ(i2)B=∑i=1Nφ(i2)化简一下，用杜教筛即可。

因为φ(n)=n⋅Πpi−1pi，φ(n2)=n2⋅Πpi−1piφ(n)=n·Πpi−1pi，φ(n2)=n2·Πpi−1pi

其中n=pk11pk22..pkmm,n2=p2k11p2k22..p2kmmn=p1k1p2k2..pmkm,n2=p12k1p22k2..pm2km

所以φ(n2)=n⋅φ(n)φ(n2)=n·φ(n)

题目就是求

∑i=1nφ(i2)=∑i=1ni⋅φ(i)∑i=1nφ(i2)=∑i=1ni·φ(i)

这个利用狄利克雷卷积知识就可以知道，(id⋅φ)∗id=id2(id·φ)∗id=id2，

即∑d|nd⋅φ(d)⋅nd=n∑d|nφ(d)=n2∑d|nd·φ(d)·nd=n∑d|nφ(d)=n2

我们设ϕ′(n)=∑ni=1i⋅φ(i)ϕ′(n)=∑i=1ni·φ(i)，利用平方前缀和公式，我们得到

n(n+1)(2n+1)6=∑i=1ni2=∑i=1n∑d|id⋅φ(d)⋅id=∑id=1id⋅∑d=1⌊nid⌋d⋅φ(d)=∑i=1ni⋅ϕ′(⌊ni⌋)n(n+1)(2n+1)6=∑i=1ni2=∑i=1n∑d|id·φ(d)·id=∑id=1id·∑d=1⌊nid⌋d·φ(d)=∑i=1ni·ϕ′(⌊ni⌋)

因此ϕ′(n)=n(n+1)(2n+1)6−∑ni=2i⋅ϕ′(⌊ni⌋)ϕ′(n)=n(n+1)(2n+1)6−∑i=2ni·ϕ′(⌊ni⌋)，预处理一下，用hash记忆化递归，时间就可以保证在O(n23)O(n23)了。

Code

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <cmath>
#include <algorithm>
#define N 1000010
#define M 2333333
#define Mod 1000000007
using namespace std;

typedef long long LL;
int n, cnt;
bool vis
;
int prime
;
LL phi
, hash_v[M], hash_n[M];

void Da(){
vis[1] = true;
phi[1] = 1ll;
for(int i = 2; i < N; i++){
if(!vis[i]){
prime[++cnt] = i;
phi[i] = (LL)(i-1);
}
for(int j = 1; j <= cnt && i * prime[j] < N; j++){
vis[i * prime[j]] = true;
if(i % prime[j] == 0){
phi[i * prime[j]] = phi[i] * (LL)prime[j];
break;
}
else  phi[i * prime[j]] = phi[i] * (LL)(prime[j]-1);
}
}
for(int i = 1; i < N; i++){
phi[i] = phi[i] % Mod * i % Mod;
phi[i] = (phi[i] + phi[i-1]) % Mod;
}
}

LL Find(LL x){
int p = x % M;
while(hash_n[p] && hash_n[p] != x)  p = (p+1) % M;
if(!hash_n[p])  return -1;
return hash_v[p];
}

void Push(LL v, LL x){
int p = x % M;
while(hash_n[p])  p = (p+1) % M;
hash_v[p] = v;
hash_n[p] = x;
}

LL Sum(LL x){
if(x < N)  return phi[x];
LL res = Find(x);
if(~ res)  return res;
else{
LL temp1 = x * (x + 1), temp2 = 2 * x + 1;
if(!(temp1 % 6))  temp1 /= 6;
else if(!(temp2 % 6))  temp2 /= 6;
else if(!(temp1 % 2) && !(temp2 % 3))  temp1 /= 2, temp2 /= 3;
else if(!(temp2 % 2) && !(temp1 % 3))  temp1 /= 3, temp2 /= 2;
res = ((temp1 % Mod) * (temp2 % Mod)) % Mod;
}
LL last;
for(LL i = 2; i <= x; i = last+1){
last = x/(x/i);
res = (res - (((last+i) * (last-i+1) / 2) % Mod * Sum(x/i)) % Mod + Mod) % Mod;
}
Push(res, x);
return res;
}

int main(){

freopen("bzoj4916.in", "r", stdin);
freopen("bzoj4916.out", "w", stdout);

Da();

scanf("%lld", &n);
printf("1\n%lld\n", Sum(n));

return 0;
}