51nod1222 最小公倍数计数
求$\sum_{i=1}^{n}\sum_{j=1}^{n}[[i,j\leq n]$。n<1e11。
方法零：变不等为相等，直接反演，求其前缀和
$\sum_{i=1}^{n}\sum_{j=1}^{n}[[i,j]=n]$
闪一句：反演！！！
$=\sum_{d|n}\mu(d)\sum_{i=1}^n\sum_{j=1}^i[[i,j]|\frac{n}{d}]$
闪一句：$\sum_{i=1}^n\sum_{j=1}^i[[i,j]|\frac{n}{d}]=\sum_{i=1}^{n}\sum_{j=1}^{i}[i|\frac{n}{d}\wedge j|\frac{n}{d}]=\frac{(d(\frac{n}{d})+1)(d(\frac{n}{d}))}{2}$
$=\sum_{d|n}\mu(d)\frac{(d(\frac{n}{d})+1)(d(\frac{n}{d}))}{2}$
$=\frac{1}{2}(\sum_{d|n}\mu(d)d^2(\frac{n}{d})+\sum_{d|n}\mu(d)d(\frac{n}{d}))$
分别求括号里两坨东西的前缀和。
$\sum_{i=1}^{n}\sum_{d|i}\mu(d)d(\frac{i}{d})$
$=\sum_{d=1}^n\mu(d)\sum_{k=1}^{\left \lfloor \frac{n}{d} \right \rfloor}d(k)$
$=\sum_{d=1}^n\mu(d)\sum_{k=1}^{\left \lfloor \frac{n}{d} \right \rfloor}\left \lfloor \frac{n}{kd} \right \rfloor$
$\sum_{i=1}^{n}\sum_{d|i}\mu(d)d^2(\frac{i}{d})$
$=\sum_{d=1}^{n}\mu(d)\sum_{k=1}^{\left \lfloor \frac{n}{d} \right \rfloor}d^2(k)$
$\mu$的前缀和用杜教筛。整除函数的前缀和可用类似杜教筛的方法。约数个数的平方的前缀和可用洲阁筛，似乎最后可$\sqrt{n}$时间对每个$\frac{n}{i}$的取值的对应前缀和求出，我不太会就用了$n^{\frac{2}{3}}$的方法。
结果？妥妥的卡空间，空间开小卡时间。
1 //#include<iostream> 2 #include<cstring> 3 #include<cstdlib> 4 #include<cstdio> 5 //#include<map> 6 #include<math.h> 7 //#include<time.h> 8 //#include<complex> 9 #include<algorithm> 10 using namespace std; 11 12 #define LL long long 13 LL n,a,b; 14 15 LL md=3000000,mz=320000; 16 #define maxa 3000011 17 short miu[maxa],summiu[maxa],d[maxa],ci[maxa];unsigned int sumd2[maxa];int mi[maxa],sumd[maxa],prime[maxa/10],lp,lw; bool notprime[maxa]; 18 #define maxb 640011 19 LL g[maxb],ff[maxb],who[maxb],i0[maxb],sp[maxb]; 20 void pre() 21 { 22 miu[1]=summiu[1]=sumd[1]=sumd2[1]=d[1]=1; lp=0; 23 for (int i=2;i<=md;i++) 24 { 25 if (!notprime[i]) {prime[++lp]=i; miu[i]=-1; d[i]=2; mi[i]=i; ci[i]=1;} 26 summiu[i]=summiu[i-1]+miu[i]; 27 sumd[i]=sumd[i-1]+d[i]; 28 sumd2[i]=sumd2[i-1]+d[i]*d[i]; 29 for (int j=1,tmp;j<=lp && 1ll*i*prime[j]<=md;j++) 30 { 31 notprime[tmp=i*prime[j]]=1; 32 if (i%prime[j]) {miu[tmp]=-miu[i]; d[tmp]=d[i]*2; mi[tmp]=prime[j]; ci[tmp]=1;} 33 else {miu[tmp]=0; d[tmp]=d[i/mi[i]]*(ci[i]+2); mi[tmp]=mi[i]*prime[j]; ci[tmp]=ci[i]+1; break;} 34 } 35 } 36 int tmp=0; 37 for (int i=2;i<=mz;i++) {sp[i]=sp[i-1]+!notprime[i]; if (!notprime[i]) tmp++;} 38 lp=tmp-1; mz=prime[tmp]-1; 39 } 40 41 #define maxh 1000007 42 struct Hash 43 { 44 int first[maxh],le; struct Edge{LL to,v; int next;}edge[maxb]; 45 void clear(LL n) {le=2; for (LL i=1;i<=n;i=n/(n/i)+1) first[n/i%maxh]=0;} 46 void insert(LL x,LL y) {int z=x%maxh; Edge &e=edge[le]; e.to=x; e.v=y; e.next=first[z]; first[z]=le++;} 47 LL find(LL x) {int z=x%maxh; for (int i=first[z];i;i=edge[i].next) if (edge[i].to==x) return edge[i].v; return -1;} 48 }hm,hz; 49 50 LL dumiu(LL n) 51 { 52 if (n<=md) return summiu[n]; 53 LL tmp=hm.find(n); if (tmp!=-1) return tmp; 54 LL ans=1; 55 for (LL i=2,last;i<=n;i=last+1) 56 { 57 last=n/(n/i); 58 ans-=dumiu(n/i)*(last-i+1); 59 } 60 hm.insert(n,ans); return ans; 61 } 62 63 LL dud(LL n) 64 { 65 if (n<=md) return sumd[n]; 66 LL ans=0; 67 for (LL i=1,last;i<=n;i=last+1) 68 { 69 last=n/(n/i); 70 ans+=(last-i+1)*(n/i); 71 } 72 return ans; 73 } 74 75 void mg(LL n) 76 { 77 lw=0; hz.clear(n); for (LL i=1;i<=n;i=n/(n/i)+1) lw++,who[lw]=g[lw]=n/i,hz.insert(who[lw],lw); 78 memset(i0,0,sizeof(i0)); 79 for (int i=1;i<=lp;i++) 80 for (int j=1;j<=lw && (LL)prime[i]*prime[i]<=who[j];j++) 81 { 82 int k=hz.find(who[j]/prime[i]); 83 g[j]-=g[k]-(i-1-i0[k]); 84 i0[j]=i; 85 } 86 } 87 88 void mff(LL n) 89 { 90 for (int i=lp;i;i--) 91 for (int j=1;j<=lw && (LL)prime[i]*prime[i]<=who[j];j++) 92 { 93 if (prime[i+1]>who[j]) ff[j]=1; 94 else if ((LL)prime[i+1]*prime[i+1]>who[j]) ff[j]=(sp[min(mz,who[j])]-sp[prime[i+1]-1])*4+1; 95 LL tmp=prime[i],k2; 96 for (int c=1;tmp<=who[j];tmp*=prime[i],c++) 97 { 98 LL now=who[j]/tmp; 99 if (prime[i+1]>now) k2=1; 100 else if (prime[i+1]*(LL)prime[i+1]>now) k2=(sp[min(mz,now)]-sp[prime[i+1]-1])*4+1; 101 else k2=ff[hz.find(now)]; 102 ff[j]+=k2*(c+1)*(c+1); 103 } 104 } 105 } 106 107 LL list[maxb]; 108 void zh(LL n) 109 { 110 for (int i=1;i<=lw;i++) 111 { 112 if ((LL)4>who[i]) list[i]=1; 113 else list[i]=ff[i]; 114 if (who[i]>md) for (int j=1,last;j<=mz;j=last+1) 115 { 116 LL tmp=who[i]/j; last=min(mz,who[i]/tmp); 117 if (prime[lp+1]>tmp) continue; 118 int id=hz.find(tmp); 119 list[i]+=(sumd2[last]-sumd2[j-1])*(g[id]-1-(lp-i0[id]))*4; 120 } 121 else list[i]=sumd2[who[i]]; 122 } 123 } 124 125 LL calc(LL n) 126 { 127 mg(n); mff(n); zh(n); 128 LL ans=0; 129 hm.clear(n); 130 for (LL i=1,last;i<=n;i=last+1) 131 { 132 last=n/(n/i); 133 ans+=(dumiu(last)-dumiu(i-1))*dud(n/i); 134 } 135 for (LL i=1,last;i<=n;i=last+1) 136 { 137 last=n/(n/i); 138 ans+=(dumiu(last)-dumiu(i-1))*list[hz.find(n/i)]; 139 } 140 return ans>>1; 141 } 142 143 int main() 144 { 145 scanf("%lld%lld",&a,&b); 146 pre(); 147 LL ans=-calc(a-1); ans+=calc(b); 148 printf("%lld\n",ans); 149 return 0; 150 }
方法零点五：还是求前缀和，换个方法。
$\sum_{i=1}^{n}\sum_{j=1}^{i}[[i,j]=n]$
$=\sum_{d=1}^{n}\sum_{i=1}^{n}\sum_{j=1}^{i}[\frac{ij}{d}=n][(i,j)=d]$
$=\sum_{d=1}^n\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^i[ijd=n][(i,j)=1]$
闪一句：$d|n$是$[ijd=n]$的必要条件。
$=\sum_{d|n}\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^i[ij=\frac{n}{d}][(i,j)=1]$
在这里停住！
第一个$\sum$后面，与一个积性函数有关！！
把后面那坨东西记作$q(\frac{n}{d})$，设数$x$有$p(x)$个不同质因子，可得$q(x)=2^{p(x)-1}$。
而$(x,y)=1 \ \ -> \ \ p(xy)=p(x)+p(y)$，因此$2^{p(xy)}=2^{p(x)}*2^{p(y)}$，$q$的两倍是一个积性函数，可以洲阁筛筛其前缀和。
最后
$\sum_{i=1}^{n}\sum_{d|i}q(d)$
$=\sum_{d=1}^nq(d)\left \lfloor \frac{n}{d}\right \rfloor$
好的，洲阁筛在$\frac{n^{\frac{3}{4}}}{ln(n)}+n^{\frac{2}{3}}$时间内处理每个$\frac{n}{i}$的$q$的前缀和，最后记得$+1$，$\div 2$。
1 //#include<iostream> 2 #include<cstring> 3 #include<cstdlib> 4 #include<cstdio> 5 //#include<map> 6 #include<math.h> 7 #include<time.h> 8 //#include<complex> 9 #include<algorithm> 10 using namespace std; 11 12 #define LL long long 13 LL a,b,n,mz=320000,md=8000000; 14 #define maxm 8000011 15 LL sump[maxm]; short pp[maxm]; int prime[maxm/10],lp; bool notprime[maxm]; 16 #define maxn 640011 17 LL who[maxn],ff[maxn],g[maxn],sp[maxn];int i0[maxn],lw; 18 19 void pre() 20 { 21 pp[1]=0; lp=0; sump[1]=1; 22 for (int i=2;i<=md;i++) 23 { 24 if (!notprime[i]) {prime[++lp]=i; pp[i]=1;} 25 sump[i]=sump[i-1]+(1<<pp[i]); 26 for (int j=1,tmp;j<=lp && (LL)prime[j]*i<=md;j++) 27 { 28 notprime[tmp=i*prime[j]]=1; 29 if (i%prime[j]) {pp[tmp]=pp[i]+1;} 30 else {pp[tmp]=pp[i]; break;} 31 } 32 } 33 for (int i=2;i<=mz;i++) sp[i]=sp[i-1]+!notprime[i]; 34 lp=0; for (int i=2;i<=mz;i++) if (!notprime[i]) lp++; 35 mz=prime[lp]-1; lp--; 36 } 37 38 #define maxh 1000007 39 struct Hash 40 { 41 int first[maxh],le;struct Edge{LL to; int v,next;}edge[maxn]; 42 void clear(LL n) {le=2; for (LL i=1;i<=n;i=n/(n/i)+1) first[n/i%maxh]=0;} 43 void insert(LL x,int y) {int h=x%maxh; Edge &e=edge[le]; e.to=x; e.v=y; e.next=first[h]; first[h]=le++;} 44 int find(LL x) {int h=x%maxh; for (int i=first[h];i;i=edge[i].next) if (edge[i].to==x) return edge[i].v; return -1;} 45 }h; 46 47 void mg() 48 { 49 lw=0; h.clear(n); 50 for (LL i=1;i<=n;i=n/(n/i)+1) lw++,who[lw]=g[lw]=n/i,h.insert(who[lw],lw); 51 memset(i0,0,sizeof(i0)); 52 for (int i=1;i<=lp;i++) 53 for (int j=1;j<=lw && (LL)prime[i]*prime[i]<=who[j];j++) 54 { 55 int k=h.find(who[j]/prime[i]); 56 g[j]-=g[k]-(i-1-i0[k]); 57 i0[j]=i; 58 } 59 } 60 61 int cnt=0; 62 void mff() 63 { 64 for (int i=lp;i;i--) 65 for (int j=1;j<=lw && (LL)prime[i]*prime[i]<=who[j];j++) 66 { 67 if (prime[i+1]>who[j]) ff[j]=1; 68 else if ((LL)prime[i+1]*prime[i+1]>who[j]) ff[j]=(sp[min(mz,who[j])]-sp[prime[i+1]-1])*2+1; 69 LL tmp=prime[i],k2=0; 70 cnt++; 71 for (int c=1;tmp<=who[j];tmp*=prime[i],c++) 72 { 73 LL t=who[j]/tmp; 74 if (prime[i+1]>t) k2+=1; 75 else if ((LL)prime[i+1]*prime[i+1]>t) k2+=(sp[min(mz,t)]-sp[prime[i+1]-1])*2+1; 76 else k2+=ff[h.find(t)]; 77 } 78 ff[j]+=k2*2; 79 } 80 } 81 82 LL list[maxn]; 83 void zh() 84 { 85 for (int i=1;i<=lw;i++) 86 { 87 if (4>who[i]) list[i]=1; 88 else list[i]=ff[i]; 89 if (who[i]<=md) list[i]=sump[who[i]]; 90 else for (int j=1,last;j<=mz;j=last+1) 91 { 92 LL now=who[i]/j; last=min(mz,who[i]/now); 93 if (prime[lp+1]>now) continue; 94 int id=h.find(now); 95 list[i]+=(sump[last]-sump[j-1])*(g[id]-1-(lp-i0[id]))*2; 96 } 97 list[i]=(list[i]+1)>>1; 98 } 99 list[lw+1]=0; 100 } 101 102 LL calc(LL nn) 103 { 104 n=nn; mg(); mff(); zh(); 105 LL ans=0; 106 for (int i=1;i<=lw;i++) ans+=(n/who[i])*(list[i]-list[i+1]); 107 return ans; 108 } 109 110 int main() 111 { 112 scanf("%lld%lld",&a,&b); pre(); 113 printf("%lld\n",calc(b)-calc(a-1)); 114 return 0; 115 }
很好。
方法三：没关系我们东山再起。前面是把$\leq n$变$=n$然后算前缀和，这次直接算$\leq n$。为了方便，把$x \leq y$的限制删去，答案+b-a+1后/2。
$\sum_{i=1}^{n}\sum_{j=1}^{n}[[i,j]\leq n]$
$=\sum_{d=1}^{n}\sum_{i=1}^{n}\sum_{j=1}^{n}[\frac{ij}{d}\leq n][(i,j)=d]$
$=\sum_{d=1}^{n}\sum_{i=1}^{\left \lfloor \frac{n}{d} \right \rfloor}\sum_{j=1}^{\left \lfloor \frac{n}{d} \right \rfloor}[ijd \leq n][(i,j)=1]$
$=\sum_{d=1}^{n}\sum_{i=1}^{\left \lfloor \frac{n}{d} \right \rfloor}\sum_{j=1}^{\left \lfloor \frac{n}{d} \right \rfloor}[ijd \leq n]\sum_{t|i\wedge t|j}\mu(t)$
$=\sum_{t=1}^{n}\mu(t)\sum_{d=1}^{n}\sum_{i=1}^{\left \lfloor \frac{n}{dt^2} \right \rfloor}\sum_{j=1}^{\left \lfloor \frac{n}{dt^2} \right \rfloor}[ijd \leq \left \lfloor \frac{n}{t^2} \right \rfloor]$
$=\sum_{t=1}^{\sqrt{n}}\mu(t)\sum_{d=1}^{\left \lfloor \frac{n}{t^2} \right \rfloor}\sum_{i=1}^{\frac{n}{dt^2}}\sum_{j=1}^{\frac{n}{dt^2}}[ijd \leq \left \lfloor \frac{n}{t^2} \right \rfloor]$
嗯问题在右边那三个$\sum$。就是形如$abc \leq n$的东西。
令$a<b<c$，则$a$只要枚举到$n^{\frac{1}{3}}$，b枚举到$\sqrt{\frac{n}{a}}$，c直接得范围。
然后算三个不等、两个相等、三个相等的，乘对应组合数。
复杂度O（能过）
1 //#include<iostream> 2 #include<cstring> 3 #include<cstdlib> 4 #include<cstdio> 5 //#include<map> 6 #include<math.h> 7 #include<time.h> 8 //#include<complex> 9 #include<algorithm> 10 using namespace std; 11 12 #define LL long long 13 LL a,b,m=320000; 14 #define maxn 640011 15 short miu[maxn],summiu[maxn]; int prime[maxn],lp; bool notprime[maxn]; 16 void pre() 17 { 18 miu[1]=summiu[1]=1; 19 for (int i=2;i<=m;i++) 20 { 21 if (!notprime[i]) {prime[++lp]=i; miu[i]=-1;} 22 summiu[i]=summiu[i-1]+miu[i]; 23 for (int j=1,tmp;j<=lp && (LL)prime[j]*i<=m;j++) 24 { 25 notprime[tmp=i*prime[j]]=1; 26 if (i%prime[j]) {miu[tmp]=-miu[i];} 27 else {miu[tmp]=0; break;} 28 } 29 } 30 } 31 32 LL bbb(LL n) 33 { 34 LL A,B,C,D; A=B=C=D=0; 35 for (LL i=1;i*i*i<=n;i++) for (LL j=i+1,to=n/i;j*j<=to;j++) A+=to/j-j; A*=6; 36 for (LL i=1;i*i*i<=n;i++) B+=n/i/i-i; B*=3; 37 for (LL i=1;i*i*i<=n;i++) C+=(LL)(sqrt(n/i)+1e-6)-i; C*=3; 38 for (LL i=1;i*i*i<=n;i++) D++; 39 return A+B+C+D; 40 } 41 42 LL calc(LL n) 43 { 44 LL ans=0; 45 for (LL i=1,last;i*i<=n;i=last+1) 46 { 47 for (last=i;n/(i*i)==n/(last*last);last++); last--; 48 ans+=(summiu[last]-summiu[i-1])*bbb(n/(i*i)); 49 } 50 return ans; 51 } 52 53 int main() 54 { 55 scanf("%lld%lld",&a,&b); pre(); 56 printf("%lld\n",(calc(b)-calc(a-1)+b-a+1)>>1); 57 return 0; 58 }
