P1829 [国家集训队]Crash的数字表格 / JZPTAB

[国家集训队]Crash的数字表格 / JZPTAB

这个题可以低于线性，然后也可以杜教筛到 \(O(n^{2 / 3})\) 这个样子。

首先暴力推：

\[\begin{aligned} & \sum_{i=1}^{n}\sum_{j=1}^{m} \operatorname{lcm}(i,j) \\ =& \sum_{i=1}^{n}\sum_{j=1}^{m}\dfrac{ij}{\gcd(i,j)} \\ =& \sum_{d=1}\sum_{i=1}^{n}\sum_{j=1}^{m}\dfrac{ij[\gcd(i,j)=d]}{d} \\ =& \sum_{d=1}\sum_{i=1}^{\lfloor n/d\rfloor}\sum_{j=1}^{\lfloor m /d \rfloor}\dfrac{ijd^2}{d}[\gcd(i,j)=1] \\ =& \sum_{d=1}d\sum_{i=1}^{\lfloor n / d \rfloor}\sum_{j=1}^{\lfloor m / d \rfloor}ij\sum_{k\mid \gcd(i,j)}\mu(k) \\ =& \sum_{d=1}d\sum_{k=1}\mu(k)k^2\sum_{i=1}^{\lfloor n / kd\rfloor}i\sum_{j=1}^{\lfloor m / kd\rfloor}j \end{aligned} \]

设 \(h(n)=\sum_{i=1}^{n}i={n(n+1)} / {2}\)，于是这个就是：

\[\begin{aligned} =& \sum_{d=1}d\sum_{k=1}\mu(k)k^2 h\left(\left\lfloor\dfrac{n}{kd}\right\rfloor\right) h\left(\left\lfloor \dfrac{m}{kd}\right\rfloor\right) \end{aligned} \]

好了，现在我们设 \(h'(d)=\sum_{k=1}\mu(k)k^2h(\lfloor n / kd\rfloor)h(\lfloor m / kd\rfloor)\)。

\[\begin{aligned}=& \sum_{d=1}dh'(d)\end{aligned} \]

这个显然可以数论分块 \(O(\sqrt{n})\) 做，然后 \(h'(d)\) 又是 \(O(\sqrt{n})\) 数论分块做。加上线性筛之类的东西，这个复杂度就是 \(O(n)\) 的了。

这个题已经能做出来了。然而线性还不够 nb，我们还能有低于线性的解法。

首先设 \(T=kd\)，则 \(k= T / d\)，于是：

\[\begin{aligned} =& \sum_{T=1}\sum_{d\mid T}d\mu \left(\dfrac{T}{d}\right) \left(\dfrac{T}{d}\right)^2 h\left(\left\lfloor\dfrac{n}{T}\right\rfloor\right) h\left(\left\lfloor\dfrac{m}{T}\right\rfloor\right) \end{aligned} \]

我们来研究这个比较怪异的卷积，看一看能不能杜教筛一类，先设：

\[f=(\mathbf{id}^2\cdot \mu)*\mathbf{id} \]

很容易能得到它的贝尔级数：

\[f _ p ( x ) = \dfrac { 1 - p ^ 2 x } { 1 - p x } \]

然后我们肯定取 \(g=\mathbf{id}^2\)，然后卷上去：

\[(\mathbf{id}^2\cdot \mu)* \mathbf{id} * \mathbf{id^2}=\mathbf{id} \]

那么 \(g\) 和 \(f*g\) 的前缀和都非常好求。\(g\) 的前缀和随便给一下：

\[S_g(n)=\sum_{i=1}^ng(i)=\dfrac{n(n+1)(2n+1)}{6} \]

然后我们杜教筛直接套就完了。

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

杜教筛不熟的我可以写个式子：

\[g(1)S_f(n)=\sum_{i=1}^n(f \ast g)(i)-\sum_{y=2}^{n}g(y)S_f\left (\left\lfloor\dfrac{n}{y}\right\rfloor\right) \]

然后这里的记忆化讲实话确实难搞，所以我直接躺平，用 Hash 随便取个质数模一下就完了。然后这里我取的模数是 \(786433\)。当然用 map 也是可以的。

因为 \(2\) 和 \(6\) 都是与 \(20101009\) 互质的，所以直接求个逆元。

然而因为我预处理直接躺平，所以最后跑得不如一开始提出的线性算法。

（现在把博客内容改了，但是我还是不想重写杜教筛的代码，你们看得懂就好了。）

代码（杜教筛 \(O(n^{2 / 3 })\)）：

#include<iostream>
#include<cstdio>
#include<cstring>
#define ll long long
using namespace std;

const ll N=1e7,M=1e6,mo=20101009,inv2=10050505,inv6=16750841,hmo=786433;

ll n,m,ans,cnt;

ll prime[M+5],mu[M+5],Sf1[M+5],Sf2[M+5];
bool f[M+5],vis[M+5];

inline void Init() {
  mu[1]=1;f[1]=1;
  for(ll i=2;i<=M;i++) {
    if(!f[i]) {prime[++cnt]=i;mu[i]=-1;}
    for(ll j=1;j<=cnt&&i*prime[j]<=M;j++) {
      f[i*prime[j]]=1;
      if(i%prime[j]==0) {mu[i*prime[j]]=0;break;}
      mu[i*prime[j]]=-mu[i];
    }
  }
  for(ll i=1;i<=M;i++) {
    for(ll j=i;j<=M;j+=i) {
      Sf1[j]=(Sf1[j]+(mu[i]*i*i+mo)%mo*(j/i)%mo)%mo;
    }
  }
  for(ll i=1;i<=M;i++) Sf1[i]=Sf1[i-1]+Sf1[i];
}

inline ll H(ll x) {return x%hmo;}

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

inline ll Sf(ll x) {
  if(x<=M) return Sf1[x];ll tmp=H(x);
  if(vis[tmp]) return Sf2[tmp];
  vis[tmp]=1;Sf2[tmp]=S1(x);
  for(ll i=2,j;i<=x;i=j+1) {
    j=x/(x/i);
    Sf2[tmp]=(Sf2[tmp]-(S2(j)-S2(i-1)+mo)%mo*Sf(x/i)%mo+mo)%mo;
  }
  return Sf2[tmp];
}

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--]);
}

int main() {

  n=read();m=read();if(n<m) swap(n,m);
  Init();

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

  write(ans);

  return 0;
}

代码（\(O(n)\)）：

#include<iostream>
#include<cstdio>
#include<cstring>
#define ll long long
using namespace std;

const ll N=1e7,mo=20101009;

ll n,m,ans,cnt;

ll prime[N+5],mu[N+5],s[N+5];
bool f[N+5];

inline void Init() {
  mu[1]=1;f[1]=1;
  for(ll i=2;i<=n;i++) {
    if(!f[i]) {prime[++cnt]=i;mu[i]=-1;}
    for(ll j=1;j<=cnt&&i*prime[j]<=n;j++) {
      f[i*prime[j]]=1;
      if(i%prime[j]==0) {mu[i*prime[j]]=0;break;}
      mu[i*prime[j]]=-mu[i];
    }
  }
  for(ll i=1;i<=n;i++) {
    s[i]=(s[i-1]+mu[i]*i%mo*i%mo)%mo;
  }
}

inline ll S(ll x) {return (x*(x+1)/2)%mo;}

inline ll F(ll x) {
  ll res=0;
  for(ll i=1,j;i<=m/x;i=j+1) {
    j=min((n/x)/(n/x/i),(m/x)/(m/x/i));
    res=(res+(s[j]-s[i-1]+mo)*S(n/x/i)%mo*S(m/x/i)%mo)%mo;
  }
  return res;
}

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--]);
}

inline void writeln(ll x) {write(x);putchar(10);}

int main() {

  n=read();m=read();
  if(n<m) swap(n,m);Init();

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

  write(ans);

  return 0;
}