51nod1222 最小公倍数计数

Posted 2020-10-26 Blue233333

求$\\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 }

View Code

方法零点五：还是求前缀和，换个方法。

$\\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 }

View Code

 

很好。

方法三：没关系我们东山再起。前面是把$\\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 }

View Code