「VMware校园挑战赛」小V的和式

Description

给定 (n,m) ，求

[sumlimits_{x_1=1}^{n}sumlimits_{x_2=1}^{n}sumlimits_{y_1=1}^{m}sumlimits_{y_2=1}^{m}frac{(x_1y_2+x_2y_1) mod (max{x_1,x_2}max{ y_1,y_2})}{gcd(max{x_1,x_2},max{ y_1,y_2})} ]

数据范围 (1le n,mle 10^9)，答案对 (998244353) 取模。

Solution

约定

(S_k(n)=sumlimits_{i=1}^{n} i^k)

---

若 (n>m) ，则交换 (n,m) ，保证 (nle m)。

拆上下指标

先对上下指标做一下处理，得：

[[1le x_1le n][1le x_2le n][1le y_1le m][1le y_2le m] ]

[=([x_1<x_2]+[x_1=x_2]+[x_1>x_2])([y_1<y_2]+[y_1=y_2]+[y_1>y_2]) ]

[=2[x_1<x_2][y_1<y_2]+2[x_1>x_2][y_1<y_2]+2([x_1<x_2][y_1=y_2]+[x_1=x_2][y_1<y_2])+[x_1=x_2][y_1=y_2] ]

考虑拆开 ([x_1<x_2][y_1<y_2]) ：

[[1le x_1<x_2][1le y_1<y_2] ]

[=([0le x_1<x_2]-[0=x_1<x_2])([0le y_1<y_2]-[0=y_1<y_2]) ]

回带到原来式子，得到：

(2[0le x_1<x_2][0le y_1<y_2])

(-2([0=x_1<x_2][0le y_1<y_2]+[0le x_1<x_2][0=y_1<y_2]))

(+2[0=x_1<x_2][0=y_1<y_2])

(+2[x_1>x_2][y_1<y_2])

(+2([x_1<x_2][y_1=y_2]+[x_1=x_2][y_1<y_2]))

(+[x_1=x_2][y_1=y_2])

所以分别求出这 (6) 种情形即可，下面开始推导。

1. ([0le x_1<x_2][0le y_1<y_2])

即求

[f_1=sumlimits_{x_1=0}^{n}sumlimits_{x_2=x_1+1}^{n} sumlimits_{y_1=0}^{m} sumlimits_{y_2=y_1+1}^{m} frac{(x_1y_2+x_2y_1) mod (x_2y_2)}{gcd(x_2,y_2)} ]

[=sumlimits_{x_2=1}^{n} sumlimits_{y_2=1}^{m} frac{1}{gcd(x_2,y_2)} sumlimits_{x_1=0}^{x_2-1}sumlimits_{y_1=0}^{y_2-1} (x_1y_2+x_2y_1) mod (x_2y_2) ]

• 引理1

对于所有的 (a in [0,y), bin[0,x))，((ax+by) mod (xy)) 将集合 (S={ (a,b) | 0le a<y, 0le b<x}) 映射到 (T={ k imes gcd(a,b) | 0le k< frac{ab}{gcd(a,b)}}) ，并且对于 (T) 中的每一元素恰好有 (gcd(a,b)) 个原象。

感性理解一下即可，故：

[f_1=frac{1}{2}sumlimits_{x_2=1}^{n} sumlimits_{y_2=1}^{m} (frac{x_2^2y_2^2}{gcd(x_2,y_2)}-x_2y_2) ]

将前半部分拎出来：

[sumlimits_{x=1}^{n}sumlimits_{y=1}^{m} frac{x^2y^2}{gcd(x,y)} ]

[=sumlimits_{d=1}^{n} d^3 sumlimits_{i=1}^{n/d}sumlimits_{j=1}^{m/d}i^2j^2 [gcd(i,j)=1] ]

[=sumlimits_{d=1}^{n}d^3sumlimits_{k=1}^{n/d} mu(k)k^4 S_2(frac{n}{dk})S_2(frac{m}{dk}) ]

[=sumlimits_{T=1}^{n} S_2(frac{n}{T})S_2(frac{m}{T})G_1(T) ]

其中 (G_1(n)=n^3 sumlimits_{d|n} mu(d)d)

因此

[f_1=frac{1}{2}sumlimits_{T=1}^{n} S_2(frac{n}{T})S_2(frac{m}{T})G_1(T)-frac{n(n+1)m(m+1)}{8} ]

2. ([0= x_1<x_2][0le y_1<y_2])

即求

[f_2=sumlimits_{x_2=1}^{n} sumlimits_{y_1=0}^{m}sumlimits_{y_2=y_1+1}^{m} frac{(x_2y_1) mod (x_2y_2)}{gcd(x_2,y_2)} ]

[=sumlimits_{x_2=1}^{n} sumlimits_{y_1=0}^{m}sumlimits_{y_2=y_1+1}^{m} frac{x_2y_1}{gcd(x_2,y_2)} ]

其实还有对称的 ([0le x_1<x_2][0=y_1<y_2]) ，推导类似。

3. ([0=x_1<x_2][0=y_1<y_2])

[f_3=sumlimits_{x_2=1}^{n}sumlimits_{y_2=1}^{m}frac{0}{gcd(n,m)}=0 ]

4. ([x_1>x_2][y_1<y_2])

[f_4=sumlimits_{x_2=1}^{n}sumlimits_{x_1=x_2+1}^{n}sumlimits_{y_1=1}^{m}sumlimits_{y_2=y_1+1}^{m} frac{(x_1y_2+x_2y_1) mod (x_1y_2)}{gcd(x_1,y_2)} ]

[=sumlimits_{x_1=1}^{n}sumlimits_{y_2=1}^{m}frac{1}{gcd(x_1,y_2)} sumlimits_{x_2=1}^{x_1-1}sumlimits_{y_1=1}^{y_2-1}x_2y_1 ]

[=sumlimits_{x_1=1}^{n}sumlimits_{y_2=1}^{m}frac{S_1(x_1-1)S_1(y_2-1)}{gcd(x_1,y_2)} ]

[=frac{1}{4}sumlimits_{d=1}^{n}frac{1}{d} sumlimits_{x_1=1}^{n}sumlimits_{y_2=1}^{m}x_1(x_1-1)y_2(y_2-1)[(x_1,y_2)=d] ]

[=frac{1}{4}sumlimits_{d=1}^{n} frac{1}{d}sumlimits_{x_1=1}^{n/d}sumlimits_{y_2=1}^{m/d}(x_1^2y_2^2d^3-x_1y_2^2d^2-x_1^2y_2d^2+x_1y_2d)[gcd(x_1,y_2=1)] ]

[=frac{1}{4}sumlimits_{d=1}^{n}sumlimits_{t=1}^{n/d}mu(t)sumlimits_{x_1=1}^{n/dt}sumlimits_{y_2=1}^{m/dt}(x_1^2y_2^2d^3t^4-x_1y_2^2d^2t^3-x_1^2y_2d^2t^3+x_1y_1dt^2) ]

[=frac{1}{4}sumlimits_{T=1}^{n}(S_2(frac{n}{T})S_2(frac{m}{T})G_1(T)-(S_1(frac{n}{T})S_2(frac{m}{T})+S_2(frac{n}{T})S_1(frac{m}{T}))G_2(T)+S_1(frac{n}{T})S_1(frac{m}{T})G_3(T)) ]

其中 (G_2(n)=n^2sumlimits_{d|n}mu(d)d)

(G_3(n)=nsumlimits_{d|n}mu(d)d)

5. ([x_1=x_2][y_1<y_2])

[f_5=sumlimits_{x_2=1}^{n}sumlimits_{y_1=1}^{m}sumlimits_{y_2=y_1+1}^{m} frac{x_2(y_1+y_2) mod (x_2y_2)}{gcd(x2,y2)} ]

[=sumlimits_{x_2=1}^{n}sumlimits_{y_1=1}^{m}sumlimits_{y_2=y_1+1}^{m} frac{x_2y_1}{gcd(x_2,y_2)} ]

可以发现 (f_5) 和 (f_2) 的式子完全一致，且符号相反，因此可以抵消。

显然对于另一半对称的式子也是如此。

6. ([x_1=x_2][y_1=y_2])

[f_6=sumlimits_{x_1=1}^{n}sumlimits_{y_1=1}^{m} frac{0}{gcd(x_2,y_2)}=0 ]

合并

答案为

[2f_1-2f_2+2f_3+2f_4+2f_5+f_6 ]

[=2f_1+2f_4 ]

因此我们只需要快速求出 (G_1(n),G_2(n),G_3(n)) 的前缀和即可。

杜教筛

考虑杜教筛，以 (G_1(n)) 为例。

(f(n)=n^3sumlimits_{d|n}mu(d)d)

记 (S(n)=sumlimits_{i=1}^{n}f(i))

则 (S(n)=sumlimits_{i=1}^{n}i^3sumlimits_{d|i}mu(d)d)

(=sumlimits_{d=1}^{n}mu(d)d^4S_3(lfloor frac{n}{d} floor))

所以我们需要杜教筛求 (f'(i)=mu(i)i^4) 的前缀和。

考虑 (g=ID^4) ，则 (h(n)=f'*g=[n==1])

故 (S(n)=1-sumlimits_{i=2}^{n}i^4 imes S(lfloor frac{n}{i} floor))

实现

时间复杂度：(O(n^{frac{5}{6}}))

空间复杂度：(O(n^{frac{2}{3}}))

#pragma GCC optimize(2)
#pragma GCC optimize(3)
#pragma GCC optimize("Ofast")
#pragma GCC optimize("inline")
#pragma GCC optimize("-fgcse")
#pragma GCC optimize("-fgcse-lm")
#pragma GCC optimize("-fipa-sra")
#pragma GCC optimize("-ftree-pre")
#pragma GCC optimize("-ftree-vrp")
#pragma GCC optimize("-fpeephole2")
#pragma GCC optimize("-ffast-math")
#pragma GCC optimize("-fsched-spec")
#pragma GCC optimize("unroll-loops")
#pragma GCC optimize("-falign-jumps")
#pragma GCC optimize("-falign-loops")
#pragma GCC optimize("-falign-labels")
#pragma GCC optimize("-fdevirtualize")
#pragma GCC optimize("-fcaller-saves")
#pragma GCC optimize("-fcrossjumping")
#pragma GCC optimize("-fthread-jumps")
#pragma GCC optimize("-funroll-loops")
#pragma GCC optimize("-fwhole-program")
#pragma GCC optimize("-freorder-blocks")
#pragma GCC optimize("-fschedule-insns")
#pragma GCC optimize("inline-functions")
#pragma GCC optimize("-ftree-tail-merge")
#pragma GCC optimize("-fschedule-insns2")
#pragma GCC optimize("-fstrict-aliasing")
#pragma GCC optimize("-fstrict-overflow")
#pragma GCC optimize("-falign-functions")
#pragma GCC optimize("-fcse-skip-blocks")
#pragma GCC optimize("-fcse-follow-jumps")
#pragma GCC optimize("-fsched-interblock")
#pragma GCC optimize("-fpartial-inlining")
#pragma GCC optimize("no-stack-protector")
#pragma GCC optimize("-freorder-functions")
#pragma GCC optimize("-findirect-inlining")
#pragma GCC optimize("-fhoist-adjacent-loads")
#pragma GCC optimize("-frerun-cse-after-loop")
#pragma GCC optimize("inline-small-functions")
#pragma GCC optimize("-finline-small-functions")
#pragma GCC optimize("-ftree-switch-conversion")
#pragma GCC optimize("-foptimize-sibling-calls")
#pragma GCC optimize("-fexpensive-optimizations")
#pragma GCC optimize("-funsafe-loop-optimizations")
#pragma GCC optimize("inline-functions-called-once")
#pragma GCC optimize("-fdelete-null-pointer-checks")
#include<bits/stdc++.h>
using namespace std;
#define rint register int
#define rep(i,l,r) for(rint i=l;i<=r;i++)
#define per(i,l,r) for(rint i=l;i>=r;i--)
#define ll long long
#define ull unsigned long long
#define pii pair<int,int>
#define pll pair<ll,ll>
#define pb push_back
#define fir first
#define sec second
#define mset(s,t) memset(s,t,sizeof(s))
template<typename T1,typename T2>void ckmin(T1 &a,T2 b){if(a>b)a=b;}
template<typename T1,typename T2>void ckmax(T1 &a,T2 b){if(a<b)a=b;}
template<typename T>T gcd(T a,T b){return b?gcd(b,a%b):a;}
int read(){
  int x=0,f=0;
  char ch=getchar();
  while(!isdigit(ch))f|=ch=='-',ch=getchar();
  while(isdigit(ch))x=10*x+ch-'0',ch=getchar();
  return f?-x:x;
}
const int N=10000005;
const int mod=998244353;
ll qpow(ll a,ll b=mod-2){
  ll res=1;
  a%=mod;
  while(b>0){
    if(b&1)res=res*a%mod;
    a=a*a%mod;
    b>>=1;
  }
  return res;
}
const ll inv2=qpow(2);
const ll inv4=qpow(4);
const ll inv6=qpow(6);
const ll inv8=qpow(8);
const ll inv30=qpow(30);
int mu[N],g1[N],g2[N],g3[N];
bool vis[N];
int pr[N>>3],len,L=1e7;
void init(int n){
  mu[1]=1;
  for(register int i=2;i<=n;i++){
    if(!vis[i])pr[len++]=i,mu[i]=-1;
    for(register int j=0;j<len&&pr[j]*i<=n;j++){
      vis[pr[j]*i]=1;
      if(i%pr[j]==0)break;
      mu[pr[j]*i]=-mu[i];
    }
  }
  rep(i,1,n){
    g1[i]=(g1[i-1]+1ll*i*i%mod*i%mod*i%mod*mu[i])%mod;
    g2[i]=(g2[i-1]+1ll*i*i%mod*i%mod*mu[i])%mod;
    g3[i]=(g3[i-1]+1ll*i*i%mod*mu[i])%mod;
  }
}
ll S1(ll x){
  x%=mod;
  return x*(x+1)%mod*inv2%mod;
}
ll S2(ll x){
  x%=mod;
  return x*(x+1)%mod*(2*x+1)%mod*inv6%mod;
}
ll S3(ll x){
  x%=mod;
  return S1(x)*S1(x)%mod;
}
ll S4(ll x){
  x%=mod;
  return x*(x+1)%mod*(2*x+1)%mod*(3*x*x%mod+3*x-1)%mod*inv30%mod;
}
int n,m;

unordered_map<int,ll>Map1,Map2,Map3;
ll GG1(int n){
  if(n<=L)return g1[n];
  if(Map1[n])return Map1[n];
  ll ans=1;
  for(int i=2,j;i<=n;i=j+1){
    j=n/(n/i);
    ans=(ans-GG1(n/i)*(S4(j)-S4(i-1)+mod))%mod;
  }
  return Map1[n]=(ans%mod+mod)%mod;
}
ll GG2(ll n){
  if(n<=L)return g2[n];
  if(Map2[n])return Map2[n];
  ll ans=1;
  for(int i=2,j;i<=n;i=j+1){
    j=n/(n/i);
    ans=(ans-GG2(n/i)*(S3(j)-S3(i-1)+mod))%mod;
  }
  return Map2[n]=(ans%mod+mod)%mod;
}
ll GG3(ll n){
  if(n<=L)return g3[n];
  if(Map3[n])return Map3[n];
  ll ans=1;
  for(int i=2,j;i<=n;i=j+1){
    j=n/(n/i);
    ans=(ans-GG3(n/i)*(S2(j)-S2(i-1)+mod))%mod;
  }
  return Map3[n]=(ans%mod+mod)%mod;
}
unordered_map<int,ll>map4,map5,map6;
ll G1(int n){
  if(map4[n])return map4[n];
  ll ans=0;
  for(int i=1,j;i<=n;i=j+1){
    j=n/(n/i);
    ans=(ans+(GG1(j)-GG1(i-1)+mod)*S3(n/i))%mod;
  }
  return map4[n]=ans;
}
ll G2(int n){
  if(map5[n])return map5[n];
  ll ans=0;
  for(int i=1,j;i<=n;i=j+1){
    j=n/(n/i);
    ans=(ans+(GG2(j)-GG2(i-1)+mod)*S2(n/i))%mod;
  }
  return map5[n]=ans;
}
ll G3(int n){
  if(map6[n])return map6[n];
  ll ans=0;
  for(int i=1,j;i<=n;i=j+1){
    j=n/(n/i);
    ans=(ans+(GG3(j)-GG3(i-1)+mod)*S1(n/i))%mod;
  }
  return map6[n]=ans;
}
ll f1(){
  ll ans=0;
  for(int i=1,j;i<=n;i=j+1){
    j=min(n/(n/i),m/(m/i));
    ans=(ans+S2(n/i)*S2(m/i)%mod*(G1(j)-G1(i-1)+mod))%mod;
  }
  ans=ans*inv2%mod;
  ans=(ans+mod-S1(n)*S1(m)%mod*inv2%mod)%mod;
  return (ans+mod)%mod;
}
ll f4(){
  ll ans=0;
  for(int i=1,j;i<=n;i=j+1){
    j=min(n/(n/i),m/(m/i));
    ans=(ans+S2(n/i)*S2(m/i)%mod*(G1(j)-G1(i-1)+mod))%mod;
    ans=(ans-(S1(n/i)*S2(m/i)+S2(n/i)*S1(m/i))%mod*(G2(j)-G2(i-1)+mod))%mod;
    ans=(ans+S1(n/i)*S1(m/i)%mod*(G3(j)-G3(i-1)+mod))%mod;
  }
  ans=ans*inv4%mod;
  return (ans+mod)%mod;
}
int main(){
  init(L);
  //n=1e9,m=1e9;
  scanf("%d%d",&n,&m);
  if(n>m)swap(n,m);
  printf("%lld
",2ll*(f1()+f4())%mod);
  return 0;
}

---

关于对 ([0= x_1<x_2][0le y_1<y_2]) 的彩(推)蛋(导)

qwq 其实是因为我一开始没意料到能相互抵消，所以推了不少，这里就保留一下好了。

[=sumlimits_{d=1}^{n} sumlimits_{x_2=1}^{n/d}sumlimits_{y_2=1}^{m/d}x_2S_1(y_2d-1)[gcd(x_2,y_2)=1] ]

[=sumlimits_{d=1}^{n}sumlimits_{k=1}^{n/d}mu(k)kS_1(frac{n}{dk})sumlimits_{y_2=1}^{m/dk}S_1(y_2dk-1) ]

右边 (y_2) 的式子，令 (t=m/dk) ，经化简得

[frac{1}{2}(d^2k^2S_2(t)-dkS_1(t)) ]

因此：

[f_2=frac{1}{2}sumlimits_{d=1}^{n}sumlimits_{k=1}^{n/d}mu(k)kS_1(frac{n}{dk})(d^2k^2S_2(frac{m}{dk})-dkS_1(frac{m}{dk})) ]

然后令 (T=dk) 化简一下就行了。。。