【luogu3768】简单的数学题 欧拉函数(欧拉反演)+杜教筛

$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}

\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)$$

$$\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$$

#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;
}


|转载请注明 [原文链接][作者] ，谢谢！

posted @ 2018-04-05 14:29  GXZlegend  阅读(559)  评论(0编辑  收藏  举报