题目大意：略 洛谷传送门
杜教筛入门题？
以下都是常规套路的变形，不再过多解释
$sumlimits_{i=1}^{N}sumlimits_{j=1}^{N}ijgcd(i,j)$
$sumlimits_{i=1}^{N}sumlimits_{j=1}^{N}ijsumlimits_{d|gcd(i,j)}varphi(d)$
$sumlimits_{d=1}^{N} varphi(d) sumlimits_{i=1}^{N}sumlimits_{j=1}^{N}ij$
$sumlimits_{d=1}^{N} varphi(d) d^{2} sumlimits_{i=1}^{left lfloor frac{N}{d} ight floor}sumlimits_{j=1}^{left lfloor frac{N}{d} ight floor}ij$
令$f(d)=varphi(d) d^{2}$，由于$varphi(d)$和$1(d)$都是积性函数，所以$f(d)$是积性函数
数据范围好像有点大？$n leq 10^{10}$？杜教筛！
考虑设计一个和$f$搭配且同为积性函数$g$，且$(f*g)$能很快计算出来，令$f(n)$前缀和是$S(n)$
给出杜教筛的推导式
$S(n)=sumlimits_{i=1}^{N}(f*g)(i)-sumlimits_{i=2}^{N}g(i)S(left lfloor frac{N}{i} ight floor)$
根据常见狄利克雷卷积形式$id=1*varphi$，容易推出比较优秀的$g(d)=d^{2}$，那么$(f*g)(n)=sumlimits_{d|n}f(d)*g(frac{n}{d})=varphi(d) d^{2}(frac{n}{d})^{2}=n^{3}$
$g$的前缀和就是平方和，$(f*g)$的前缀和就是立方和，可以$O(1)$求出来
一定要多取模，不然有些地方可能会爆$longlong$
1 #include <cstdio> 2 #include <cstring> 3 #include <algorithm> 4 #define N1 5000010 5 #define M1 3000010 6 #define ll long long 7 #define dd double 8 #define cll const long long 9 using namespace std; 10 11 const int maxn=5000000; 12 13 void exgcd(ll a,ll b,ll &x,ll &y) 14 { 15 if(!b){ x=1; y=0; return; } 16 exgcd(b,a%b,x,y); ll t=x; x=y; y=t-a/b*y; 17 } 18 struct Hsh{ 19 #define M 3000000 20 int head[M1],nxt[M1],val[M1],cte;ll to[M1]; 21 void ins(ll x,int w) 22 { 23 int j,y=x%M;ll v; 24 for(j=head[y];j;j=nxt[j]){ v=to[j]; if(v==x) return; } 25 cte++; to[cte]=x; val[cte]=w; nxt[cte]=head[y]; head[y]=cte; 26 } 27 int query(ll x) 28 { 29 int j,y=x%M;ll v; 30 for(j=head[y];j;j=nxt[j]){ v=to[j]; if(v==x) return val[j]; } 31 return -1; 32 } 33 #undef M 34 }h; 35 36 ll n,m,inv2,inv6; 37 int use[N1],pr[N1],cnt,phi[N1]; 38 ll f[N1],sf[N1]; 39 void get_mu(cll &jr) 40 { 41 int i,j; phi[1]=f[1]=sf[1]=use[1]=1; 42 for(i=2;i<=maxn;i++) 43 { 44 if(!use[i]){ pr[++cnt]=i; phi[i]=i-1; } 45 for(j=1;j<=cnt&&i*pr[j]<=maxn;j++) 46 { 47 use[i*pr[j]]=1; 48 if(i%pr[j]){ phi[i*pr[j]]=phi[i]*(pr[j]-1)%jr; } 49 else{ phi[i*pr[j]]=phi[i]*pr[j]%jr; break; } 50 } 51 f[i]=1ll*i*i%jr*phi[i]%jr; 52 sf[i]=(sf[i-1]+f[i])%jr; 53 } 54 } 55 int F(ll x,cll &jr) 56 { 57 if(x<=maxn) return sf[x]; 58 ll ans=h.query(x); if(ans!=-1) return ans; 59 ll i,la,pi,pla; ans=x%jr*(x+1)%jr*inv2%jr; ans=ans*ans%jr; 60 for(i=2;i<=x;i=la+1){ 61 la=x/(x/i); pi=i%jr; pla=la%jr; 62 ans-=(1ll*pla%jr*(pla+1)%jr*(2*pla+1)%jr-1ll*(pi-1)%jr*(pi)%jr*(2*pi-1)%jr)*inv6%jr*F(x/i,jr)%jr; 63 } 64 ans=(ans%jr+jr)%jr; 65 h.ins(x,ans); 66 return ans; 67 } 68 69 int main() 70 { 71 scanf("%lld%lld",&m,&n); 72 cll jr=m; get_mu(jr); ll i,la,x,y,ans=0; 73 exgcd(6,jr,inv6,y); exgcd(2,jr,inv2,y); inv6=(inv6%jr+jr)%jr; inv2=(inv2%jr+jr)%jr; 74 for(i=1;i<=n;i=la+1) 75 { 76 la=n/(n/i); 77 if((n/i)&1) x=(n/i+1)/2%jr*((n/i)%jr)%jr; 78 else x=(n/i)/2%jr*((n/i+1)%jr)%jr; 79 ans=(x*x%jr*(F(la,jr)-F(i-1,jr)+jr)+ans)%jr; 80 } 81 printf("%lld ",ans); 82 return 0; 83 }