【bzoj4916】神犇和蒟蒻 杜教筛

题目描述

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

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

输入

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

输出

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

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

样例输入

1

样例输出

1
1

---

题解

杜教筛

第一问的答案毫无疑问肯定是1（wow~ ⊙o⊙），毕竟除了i=1以外质因子的幂次一定大于等于2。

对于第二问欧拉函数Φ(i^2)=iΦ(i)，因为将i分解质因数，质因子的幂次一定大于等于2且为偶数。我们计算欧拉函数时，如果出现p^a，那么是算进答案(p-1)*p^(a-1)。如果把p^(a/2)的部分取出，剩下的是Φ(i)的部分，而取出的是i的部分，因此答案是∑iΦ(i)。

现在只要求∑iΦ(i)即可。

我们把它与id求卷积，发现形式就是n^2。

所以有：

其中重点是因子与乘积的转化，在公式第2、3行将枚举乘积i转化为枚举2行的i/d和d，即3行的i和j。

先预处理出n^2/3的部分，超过的部分对于每个n/i进行递归，记忆化搜索，使用map储存。

总时间复杂度O(n^2/3logn)。

这里有一个小细节：除的数固定并且很小时不需要求逆元，直接将模数扩大相应倍数，正常取模，最后再模下去。

#include <cstdio>
#include <map>
#define N 1000010
#define mod 6000000042ll
using namespace std;
typedef long long ll;
map<ll , ll> f;
map<ll , ll>::iterator it;
ll m = 1000000 , phi[N] , prime[N] , tot , sum[N];
bool np[N];
ll s1(ll l , ll r)
{
    return (l + r) * (r - l + 1) % mod / 2;
}
ll s2(ll x)
{
    return x * (x + 1) % mod * (2 * x + 1) % mod / 6;
}
ll query(ll n)
{
    if(n <= m) return sum[n];
    it = f.find(n);
    if(it != f.end()) return it->second;
    ll ans = s2(n) , i , last;
    for(i = 2 ; i <= n ; i = last + 1) last = n / (n / i) , ans = (ans - s1(i , last) * query(n / i) % mod + mod) % mod;
    f[n] = ans;
    return ans;
}
int main()
{
    ll i , j , n;
    phi[1] = sum[1] = 1;
    for(i = 2 ; i <= m ; i ++ )
    {
        if(!np[i]) phi[i] = i - 1 , prime[++tot] = i;
        for(j = 1 ; j <= tot && i * prime[j] <= m ; j ++ )
        {
            np[i * prime[j]] = 1;
            if(i % prime[j] == 0)
            {
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }
            else phi[i * prime[j]] = phi[i] * (prime[j] - 1);
        }
        sum[i] = (sum[i - 1] + i * phi[i]) % mod;
    }
    scanf("%lld" , &n);
    printf("1
%lld
" , query(n) % 1000000007);
    return 0;
}