【luogu3768】簡單的數學題 歐拉函數(歐拉反演)+杜教篩

題目描述

給出 $n$ 和 $p$ ，求 $(\sum\limits_{i=1}^n\sum\limits_{j=1}^nij\gcd(i,j))\mod p$ 。

$n\le 10^{10}$ 。

---

題解

歐拉函數(歐拉反演)+杜教篩

推式子：

$$\begin{align}&\sum\limits_{i=1}^n\sum\limits_{j=1}^nij\gcd(i,j)\\=&\sum\limits_{i=1}^n\sum\limits_{j=1}^nij\sum\limits_{d|\gcd(i,j)}\varphi(d)\\=&\sum\limits_{i=1}^n\sum\limits_{j=1}^nij\sum\limits_{d|i,d|j}\varphi(d)\\=&\sum\limits_{d=1}^n\varphi(d)\sum\limits_{i=1}^{\lfloor\frac nd\rfloor}id\sum\limits_{j=1}^{\lfloor\frac nd\rfloor}jd\\=&\sum\limits_{d=1}^nd^2\varphi(d)(\sum\limits_{i=1}^{\lfloor\frac nd\rfloor}i)^2\\=&\sum\limits_{d=1}^nd^2\varphi(d)(\frac{\lfloor\frac nd\rfloor(\lfloor\frac nd\rfloor+1)}2)^2\end{align}$$

對 $\lfloor\frac nd\rfloor$ 分塊處理，只須要求出 $f(n)=n^2\varphi(n)$ 的前綴和便可。

顯然這是一個積性函數，然而 $n$ 有 $10^{10}$ 之大，不能線性篩。

考慮杜教篩。設 $g(n)=n^2$ ，那麼有：

$$\begin{align}&(f·g)(n)\\=&\sum\limits_{d|n}f(d)g(\frac nd)\\=&\sum\limits_{d|n}d^2\varphi(d)·(\frac nd)^2\\=&n^2\sum\limits_{d|n}\varphi(d)\\=&n^3\end{align}$$

因此：

$$\begin{align}&\sum\limits_{i=1}^ni^3\\=&\sum\limits_{i=1}^n(f·g)(i)\\=&\sum\limits_{i=1}^n\sum\limits_{d|i}f(d)·g(\frac id)\\=&\sum\limits_{i=1}^n\sum\limits_{d|i}f(\frac id)g(d)\\=&\sum\limits_{d=1}^ng(d)\sum\limits_{i=1}^{\lfloor\frac nd\rfloor}f(i)\\=&\sum\limits_{d=1}^nd^2S(\lfloor\frac nd\rfloor)\end{align}$$

故有：

$$S(n)=\sum\limits_{i=1}^ni^3-\sum\limits_{i=2}^ni^2S(\lfloor\frac nd\rfloor)$$

線性篩預處理出 $n^{\frac 23}$ 之內的 $S(i)$ ，對超過 $n^{\frac 23}$ 的部分進行杜教篩便可。

可能須要用到的公式：

$$\sum\limits_{i=1}^ni^2=\frac{n(n+1)(2n+1)}6\\\sum\limits_{i=1}^ni^3=\frac{n^2(n+1)^2}4$$

時間複雜度 $O(n^{\frac 23})$ ，這裏偷懶使用map，複雜度多一個 $\log$ 。

#include <map>
#include <cstdio>
#define N 10000010
using namespace std;
typedef long long ll;
const int m = 10000000;
map<ll , ll> mp;
int phi[N] , prime[N] , tot , np[N];
ll sum[N] , p , inv4 , inv6;
void init()
{
    int i , j;
    phi[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] * phi[prime[j]];
        }
    }
    for(i = 1 ; i <= m ; i ++ ) sum[i] = (sum[i - 1] + 1ll * phi[i] * i % p * i) % p;
    inv4 = (p + 1) / 2;
    if(p % 3 == 1) inv6 = (2 * p + 1) / 3;
    else inv6 = (p + 1) / 3;
    inv6 = inv6 * inv4 % p , inv4 = inv4 * inv4 % p;
}
ll calc2(ll n) {n %= p; return n * (n + 1) % p * (2 * n + 1) % p * inv6 % p;}
ll calc3(ll n) {n %= p; return n * n % p * (n + 1) % p * (n + 1) % p * inv4 % p;}
ll solve(ll n)
{
    if(n <= m) return sum[n];
    if(mp.find(n) != mp.end()) return mp[n];
    ll i , last , ans = calc3(n);
    for(i = 2 ; i <= n ; i = last + 1) last = n / (n / i) , ans = (ans - (calc2(last) - calc2(i - 1) + p) * solve(n / i) % p + p) % p;
    return mp[n] = ans;
}
int main()
{
    ll n , i , last , ans = 0;
    scanf("%lld%lld" , &p , &n);
    init();
    for(i = 1 ; i <= n ; i = last + 1)
        last = n / (n / i) , ans = (ans + (solve(last) - solve(i - 1) + p) * calc3(n / i)) % p;
    printf("%lld\n" , ans);
    return 0;
}