P3768

发表于 更新于 分类于
简单的数学题

简单的数学题

找到了比较好的杜教筛练习题。

这个式子欧拉反演啥的应该相当简单。

$$
\begin{aligned}
& \sum_{i=1}^{n}\sum_{j=1}^{n}ij\gcd(ij)
\\
=&\sum_{i=1}^{n}i\sum_{j=1}^{n}j\sum_{k\mid i,k\mid j}\varphi(k)
\\
=& \sum_{k=1}\varphi(k)\sum_{i=1}^{\lfloor n/d\rfloor}ik\sum_{j=1}^{\lfloor n/d\rfloor}jk
\\
=& \sum_{k=1}k^2\varphi(k)\left(\sum_{i=1}^{\lfloor n/d\rfloor}i\right)^2
\end{aligned}
$$

设 $S(x)=\sum_{i=1}^{n}i$，那么原式就是：

$$
\sum_{k=1}k^2\varphi(k)S\left(\left\lfloor\dfrac{n}{d}\right\rfloor\right)
$$

后面的 $f$ 可以数论分块，前面的这个 $k^2\varphi(k)$ 需要求前缀和。

这里数据 $n\le 10^{10}$ 一看就不能暴力求，杜教筛看起来不错。

令 $f=\mathbf{id^2}\cdot \varphi$，$g=\mathbf{id^2}$，则 $f\ast g=\mathbf{id^3}$。

因为 $S_g(n)=n(n+1)(2n+1)/6$，$S_f(n)=S(n)^2$，都是 $O(1)$ 可求的，那么这个显然杜教筛可做。

时间复杂度 $O(n^{2/3}\log n)$。

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
#include<iostream>
#include<cstdio>
#include<map>
#define ll __int128
using namespace std;
namespace Ehnaev{
  inline ll read() {
    ll ret=0,f=1;char ch=getchar();
    while(ch<48||ch>57) {if(ch==45) f=-f;ch=getchar();}
    while(ch>=48&&ch<=57) {ret=(ret<<3)+(ret<<1)+ch-48;ch=getchar();}
    return ret*f;
  }
  inline void write(ll x) {
    static char buf[22];static ll len=-1;
    if(x>=0) {do{buf[++len]=x%10+48;x/=10;}while(x);}
    else {putchar(45);do{buf[++len]=-(x%10)+48;x/=10;}while(x);}
    while(len>=0) putchar(buf[len--]);
  }
}using Ehnaev::read;using Ehnaev::write;
inline void writeln(ll x) {write(x);putchar(10);}

const ll M=6e6;

ll mo,n,cnt,inv2,inv6;
ll prime[M+5],phi[M+5],sf1[M+5];
map<ll,ll> sf2;
bool f[M+5];

inline ll Pow(ll b,ll p) {
  ll r=1;while(p) {if(p&1) r=r*b%mo;b=b*b%mo;p>>=1;}return r;
}

inline ll F(ll x) {return (x*(x+1)%mo)*inv2%mo;}
inline ll G(ll x) {return (x*(x+1)%mo)*((2*x%mo+1)%mo)*inv6%mo;}

inline ll Sf(ll x) {
  if(x<=M) return sf1[x];
  if(sf2.find(x)!=sf2.end()) return sf2[x];sf2[x]=F(x)*F(x)%mo;
  for(ll i=2,j;i<=x;i=j+1) {
    j=x/(x/i);sf2[x]=(sf2[x]-((G(j)-G(i-1)+mo)%mo)*Sf(x/i)%mo+mo)%mo;
  }
  return sf2[x];
}

inline void Init() {
  f[1]=1;phi[1]=1;
  for(ll i=2;i<=M;i++) {
    if(!f[i]) {prime[++cnt]=i;phi[i]=i-1;}
    for(ll j=1;j<=cnt&&i*prime[j]<=M;j++) {
      f[i*prime[j]]=1;
      if(i%prime[j]==0) {
        phi[i*prime[j]]=phi[i]*prime[j];break;
      }
      phi[i*prime[j]]=phi[i]*(prime[j]-1);
    }
  }
  for(ll i=1;i<=M;i++) {sf1[i]=(sf1[i-1]+(phi[i]*i%mo)*i%mo)%mo;}
  inv2=Pow(2,mo-2);inv6=Pow(6,mo-2);
}

int main() {

  mo=read();n=read();
  Init();

  ll ans=0;
  for(ll i=1,j;i<=n;i=j+1) {
    j=n/(n/i);ans=(ans+(((Sf(j)-Sf(i-1)+mo)%mo)*F(n/i)%mo)*F(n/i)%mo)%mo;
  }

  write(ans);

  return 0;
}