求\(\sum_{i=1}^{a}\sum_{j=1}^{b}\sum_{k=1}^{c}d(ijk)\),其中\(a,b,c\leq 10^5\),最多10组数据。
曾经考过这道题的弱化版:26.3.18 考试小记,但是这道题目中给出的标解是\(O(n^2log(n))\)的,在这道题中最多只能拿到30分,上午考试推了半小时推出了这个式子,然后写出来发现由于常数过大只能得到10分= =。之后我一直尝试从式子的角度去优化,但又推了两个多小时发现没有更好的办法。
标解给的很神奇,转化到了图论上,但是我并不会做,而且复杂度也很差。有的人去写了爆搜,好像效果还行,不过复杂度很玄学。
我们可以考虑从数学推导的角度接着搞。首先,无论用什么方法(从二维推到三维或者直接推),我们可以得到答案的式子:
$$\sum_{i=1}^{a}\sum_{j=1}^{b}\sum_{k=1}^{c} [(i,j)==1] [(j,k)=1] [(i,k)=1] \lfloor \frac{a}{i} \rfloor \lfloor \frac{b}{j} \rfloor\lfloor \frac{c}{k} \rfloor$$
我们容斥d=gcd(b,c),枚举a,这样还需要保证两组互质关系,写成式子的形式,就是我们最后要的(直接粘的上次写的):
$$\sum_{i=1}^{A} \lfloor \frac{A}{i} \rfloor \sum_{t=1} ^{min(B,C)} \mu(t) [gcd(t,i)=1] \; \sum_{j=1,(i,j)=1} ^{\lfloor \frac{B}{t} \rfloor } \lfloor \frac{B}{jt} \rfloor \; \sum_{k=1,(i,k)=1} ^{\lfloor \frac{C}{t} \rfloor } \lfloor \frac{C}{kt} \rfloor$$
我们考虑如何计算这个东西。最外层的\(i\)是必须直接枚举的,对于内层,显然如果我们能对\(t\)进行分块处理,那么复杂度才是可以接受的。但是这样就需要我们来搞一个前缀和,而前缀和又是与外层的\(i\)相关的,因此我们考虑快速计算这个前缀和。
首先我们定义函数 \( h(x,y)=\sum_{i=1}^{x} [gcd(i,y)=1] \lfloor \frac{x}{i} \rfloor \),我们的目的就是计算这个玩意。
接着,我们发现在分块过程中,最多只会用到sqrt(n)个h,那么我们可以考虑计算每个单点的h值,先把每个断点的函数值初始化成不管互质关系的前缀和。考虑对\(i\)分解质因数,也就是从小到大枚举i的每一个质因子,那么我们只要从大到小倒序枚举一个断点,令s[pos[x]]-=s[pos[x]/y]即可,这样,由于断点与商的值是等价的,我们将断点的值除掉y,也一定会出现在断点数组里面,就有了一个逐步晒掉所有不符合要求的质因子的递推过程。相似的,对于mu的处理,也是用这种思路,不过要用正序枚举,因为如果原来包含这个因子,现在还包含这个因子,对于mu而言会变成0,而我们的目的其实是去除影响,所以这里要特殊处理。
最后,我们就可以得到符合条件的mu的前缀和以及h值,就可以搞了,时间复杂度\(O(n^{1.5}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 |
#include <bits/stdc++.h> #pragma GCC optimize("-O3") using namespace std; #define ll long long #define mod 1000000007 #define mxn 100005 int prm[mxn],tot,mu[mxn],smu[mxn],sd[mxn],dv[mxn],mf[mxn],lf[mxn]; bool inp[mxn]; inline void gprm(){ register int i,j,k; mu[1]=1,mf[1]=lf[1]=1,dv[1]=1; for(i=2;i<=100000;++i){ if(!inp[i]) prm[++tot]=i,mu[i]=-1,mf[i]=i,lf[i]=1,dv[i]=2; for(j=1;j<=tot&&(k=i*prm[j])<=100000;++j){ inp[k]=1; if(i%prm[j]) mu[k]=-mu[i],mf[k]=prm[j],lf[k]=i,dv[k]=dv[i]<<1; else { mf[k]=mf[i],lf[k]=lf[i],dv[k]=(dv[i]/dv[i/lf[i]])*(dv[i/lf[i]]+1); break; } } } for(register int i=1;i<=100000;++i) smu[i]=smu[i-1]+mu[i],sd[i]=dv[i]+sd[i-1]; } int bo[mxn],cnt,a,b,c; int f[mxn],g[mxn],res[mxn],V1[mxn],V2[mxn]; bool vis[mxn]; inline int calc(int t){ register ll ret=0; register int v=1,tp=t,cur; while(tp>1) v=v*mf[tp],tp=lf[tp]; if(vis[v]) return res[v]; for(int i=1;i<=cnt;++i) f[bo[i]]=sd[bo[i]],g[bo[i]]=smu[bo[i]]; tp=t; while(tp>1) { cur=mf[tp],tp=lf[tp]; for(int i=1;i<=cnt;++i) g[bo[i]]+=g[bo[i]/cur]; for(int i=cnt;i>=1;--i) f[bo[i]]-=f[bo[i]/cur]; } for(register int i=1;i<=cnt&&bo[i]<=b;++i) ret+=1ll*(g[bo[i]]-g[bo[i-1]])*f[V1[i]]*f[V2[i]]; vis[v]=1; return res[v]=(ret+mod)%mod; } inline void work(){ scanf("%d%d%d",&a,&b,&c),memset(vis,0,sizeof(vis)); if(a>b) swap(a,b); if(a>c) swap(a,c); if(b>c) swap(b,c); cnt=0; for(int i=1,ls;i<=c;i=ls+1){ if(i<=b) ls=min(b/(b/i),c/(c/i)); else ls=c/(c/i); bo[++cnt]=ls,V1[cnt]=b/bo[cnt],V2[cnt]=c/bo[cnt]; } int ans=0; for(int i=1;i<=a;++i) ans+=1ll*(a/i)*calc(i)%mod,ans>=mod?ans-=mod:0; printf("%d\n",ans); } int main(){ //freopen("t2.in","r",stdin); gprm(); int cas; scanf("%d",&cas); while(cas--) work(); } |
叨叨几句... NOTHING