【题解】51nod1575 LCM and GCD (min25筛)

【题解】51nod1575 LCM and GCD (min25筛)

http://www.51nod.com/Challenge/Problem.html#problemId=1575

设\(f(n)=\sum_\limits{i,j} [(i,n),(n,j)]\)，打出以下的表，发现是个积性函数

1 :: 1
2 :: 7
3 :: 19
4 :: 42
5 :: 61
6 :: 133
7 :: 127
8 :: 228
9 :: 273
10 :: 427
11 :: 331
12 :: 798
13 :: 469
14 :: 889
15 :: 1159
16 :: 1160
17 :: 817
18 :: 1911
19 :: 1027
20 :: 2562

什么，你要问证明?

\[f(n)={\sum_{d_1|n}{\sum_{d_2|n}{\text{lcm}(d_1,d_2)\varphi(\frac{n}{d_1})\varphi(\frac{n}{d_2})}}}\\={\sum_{d_1|n}{\sum_{d_2|n}{\sum_{d|d_1,d|d_2}{\sum_{d'|\frac{d_1}{d},d'|\frac{d_2}{d}}{\mu(d')d'^2\frac{d_1}{dd'}\frac{d_2}{dd'}d}}\varphi(\frac{n}{\frac{d_1}{dd'}dd'})\varphi(\frac{n}{\frac{d_2}{dd'}dd'})}}}\\={\sum_{d|n}{\sum_{d'|\frac{n}{d}}{d\mu(d')d'^2\sum_{\frac{d_1}{dd'}|\frac{n}{dd'}}{\frac{d_1}{dd'}\varphi(\frac{\frac{n}{dd'}}{\frac{d_1}{dd'}})}\sum_{\frac{d_2}{dd'}|\frac{n}{dd'}}{\frac{d_2}{dd'}\varphi(\frac{\frac{n}{dd'}}{\frac{d_2}{dd'}})}}}}\\={\sum_{d|n}{\sum_{d'|\frac{n}{d}}{d\mu(d')d'^2(\sum_{\frac{d''}{dd'}|\frac{n}{dd'}}{\frac{d''}{dd'}\varphi(\frac{\frac{n}{dd'}}{\frac{d''}{dd'}})})^2}}}\\={\sum_{xyz=i}{x\mu(y)y^2(\sum_{d|z}{d\varphi(\frac{z}{d})})^2}}\\={(id\ast id_2\mu\ast(id\ast \varphi)^2)(n)} \]

抄自这

考虑下min25筛，考虑求质数次幂时的值

\[f(p^k)=\sum_i \sum_j [(p^k,i),(p^k,j)] \\ =\sum_{p|i} \sum_{p|j} [(p^k,i),(p^k,j)] +2\sum_{p\not \mid i}\sum_{p|j} [(p^k,i),(p^k,j)]+\sum_{p\not \mid i}\sum_{p\not \mid j}[(p^k,i),(p^k,j)] \\ =2 \sum_i^k ip^i -\sum_i^kp^i+2\varphi(p^k)\sum_i^k p^i +\varphi(p^k)^2 \\ =(2k+1)p^{2k}-(2k+1)p^{2k-1}+p^{k-1} \]

而

\[f(p)=3p^2-3p+1 \]

Min25 核tao心lu式子

\[g(n,j)=g(n,j-1)-p^k\big(g(\lfloor {n\over p}\rfloor,j-1)-{\rm sump} (j-1)\big) \]

//@winlere
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>

using namespace std;  typedef long long ll;
typedef unsigned long long ull;
typedef unsigned int uint;
inline int qr(){
	int ret=0,f=0,c=getchar();
	while(!isdigit(c)) f|=c==45,c=getchar();
	while( isdigit(c)) ret=ret*10+c-48,c=getchar();
	return f?-ret:ret;
}
const int maxn=1e5+5;
uint k2[maxn],k1[maxn],k0[maxn],tot[maxn],sump2[maxn],sump1[maxn],sump0[maxn],sump[maxn];
int id1[maxn],id2[maxn],val[maxn],pr[maxn],usd[maxn],cnt,n,N,id;
uint ksm(uint base,uint p){
	uint ret=1;
	while(p){
		if(p&1) ret*=base;
		base*=base;
		p>>=1;
	}
	return ret;
}
int&getId(int val){return val>N?id2[n/val]:id1[val];}
uint sum1(uint x){return 1ull*x*(x+1)>>1;}
uint sum2(ull x){
 	int k=x%6;
 	if(k==0) return x*(x+1)/6*(2*x+1);
 	if(k==2) return x*(x+1)/6*(2*x+1);
 	if(k==3) return x*(x+1)/6*(2*x+1);
 	if(k==5) return (x+1)*(2*x+1)/6*x;
 	if(k==1) return (x+1)*(2*x+1)/6*x;
	return (2*x+1)*x/6*(x+1);
}

uint fPrime(uint v,int k){
	return (2u*k+1)*(ksm(v,2*k)-ksm(v,2*k-1))+ksm(v,k-1);
}

void seive(int n){
	for(int t=2;t<=n;++t){
		if(!usd[t]) pr[++cnt]=t,sump[cnt]=sump[cnt-1]+3*t*t-3*t+1,sump2[cnt]=sump2[cnt-1]+1*t*t,sump1[cnt]=sump1[cnt-1]+1*t,sump0[cnt]=sump0[cnt-1]+1;
		for(int i=1,v;i<=cnt;++i){
			v=pr[i]*t;
			if(v>n) break;
			usd[v]=1;
			if(t%pr[i]==0) break;
		}
	}
}

uint S(int n,int j){
	if(n==0||pr[j]>n) return 0;
	int k=getId(n);
	uint ret=tot[k]-sump[j-1];
	for(int t=j;t<=cnt&&pr[t]*pr[t]<=n;++t)
		for(ll p=pr[t],e=1;p*pr[t]<=n;++e,p*=pr[t])
			ret+=fPrime(pr[t],e+1)+fPrime(pr[t],e)*S(n/p,t+1);
	return ret;
}

void init(){
  	id=0;
  	for(int l=1,r=1;r<=n;l=++r){
  		r=n/(n/l);
  		getId(n/l)=++id;
  		k2[id]=sum2(n/l)-1,k1[id]=sum1(n/l)-1,k0[id]=n/l-1;
  		val[id]=n/l;
  	}
  	for(int t=1;t<=cnt&&pr[t]*pr[t]<=val[1];++t)
  		for(int i=1;i<=id&&pr[t]*pr[t]<=val[i];++i){
  			int k=getId(val[i]/pr[t]);
  			k2[i]-=pr[t]*pr[t]*(k2[k]-sump2[t-1]);
  			k1[i]-=      pr[t]*(k1[k]-sump1[t-1]);
  			k0[i]-=            (k0[k]-sump0[t-1]);
		}
  	for(int t=1;t<=id;++t) tot[t]=3u*k2[t]-3u*k1[t]+k0[t];
	printf("%u\n",S(n,1)+1u);
}

int main(){
	int T=qr();
	seive(maxn-2);
	while(T--) n=qr(),N=sqrt(n),init();
	return 0;
}