数论

【LOJ572】Misaka Network与求和

by AmanoKumiko

Description

一方通行成功接入了 Misaka Network。

现在他要使用超能力，自然计算式被送到了御坂网络进行处理。这次的计算式是这样子的：

\[\sum_{i=1}^{N}\sum_{i=1}^Nf(gcd(i,j))^k\ mod\ 2^{32} \]

其中\(f(x)\)表示\(x\)次大的质因数，重复的质因数计算多次，例如\(f(6)=2,f(4)=2\)。规定\(f(1)=0,f(p=1)\)，其中\(p\)为质数。

但是妹妹们都不会算这个式子……所以御坂 20001 号找到了你，希望你帮她算一下。

Input

一行两个正整数\(N\)和\(k\)

Output

一行一个整数表示答案

Sample Input

4 2

Sample Output

8

Data Constraint

对于所有数据\(1\le N,k \le 10^9\)

Solution

对\(Min25\)有了新的认识。。。

首先简单推推式子，变成了

\[\sum_{i=1}^N\lfloor\frac{N}{i}\rfloor^2\sum_{d|i}f(d)^kμ(\frac{i}{d}) \]

左边可以数论分块，右边可以考虑杜教筛

令\(S(n)=\sum_{i=1}^n\sum_{d|i}f(d)^kμ(\frac{i}{d})\)

要求\(S\)，由\(S*I=f*μ*I=f*(μ*I)=f\)，得

\[S(n)=\sum_{i=1}^nf(i)^k-\sum_{i=2}^nS(\lfloor\frac{n}{i}\rfloor) \]

现在关键是求\(f\)的前缀和

发现\(f\)在分解质因数的时候特别好求，而\(Min25\)的本质其实就是分解质因数

直接魔改一下

令\(F_m(n)\)表示上一个枚举的质数是\(p_m\)时的答案

\[F_m(n)=p_{m-1}^k(\sum_{i=1}^n[isprime(i)]-m+1)+\sum_{i>=m,p_i^2<=n}\sum_{c>=1,p_i^{c+1}<=n}F_{i+1}(\lfloor\frac{n}{p_i^c}\rfloor)+p_i^k \]

前一部分即在\(p_m\)之后再填一个质数

后一部分的\(p_i^{c+1}<=n\)是因为，当\(\frac{n}{p_i^c}<pi\)时，再递归一层就无法填比\(p_i\)大的质数

由于这样算实际上是在枚举合数，所以一开始还要加上质数的贡献

Code

#include<bits/stdc++.h>
#include<tr1/unordered_map>
using namespace std;
#define F(i,a,b) for(int i=a;i<=b;i++)
#define Fd(i,a,b) for(int i=a;i>=b;i--)
#define mo 4294967296
#define LL long long
#define N 3000010

tr1::unordered_map<LL,LL>val;
int tot,vis[N],sz,w1[N],w2[N],cnt;
LL n,k,h[N],fx[N],id[N],prime[N];

int get(int x){return x<=sz?w1[x]:w2[n/x];}

LL mi(LL x,LL y){
	if(y==1)return x;
	return y%2?x*mi(x*x%mo,y/2)%mo:mi(x*x%mo,y/2);
}

void pcnt(int x){
	vis[1]=1;
	F(i,2,x){
		if(!vis[i])prime[++tot]=i,fx[tot]=mi(i,k);
		for(int j=1;i*prime[j]<=x;j++){
			vis[i*prime[j]]=1;
			if(!(i%prime[j]))break;
		}
	}
}

LL Fk(int x,int y){
	if(x<=1||prime[y]>x)return 0;
	LL res=fx[y-1]*(h[get(x)]-y+1)%mo;
	for(int i=y;i<=tot&&prime[i]*prime[i]<=x;i++){
		LL P=prime[i];
		for(int j=1;P*prime[i]<=x;j++,P*=prime[i])res=(res+Fk(x/P,i+1)+fx[i])%mo;
	}
	return res;
}

void Fprime(){
	for(int l=1,r;l<=n;l=r+1){
		r=n/(n/l);
		id[++cnt]=n/l;
		h[cnt]=n/l-1;
		if(n/l<=sz)w1[n/l]=cnt;
		else w2[l]=cnt;
	}
	F(i,1,tot){
		for(int j=1;j<=cnt&&prime[i]*prime[i]<=id[j];j++){
			int pos=get(id[j]/prime[i]);
			h[j]=h[j]-h[pos]+i-1;
		}
	}
}

LL S(int x){
	if(val[x])return val[x];
	LL res=(Fk(x,1)+h[get(x)])%mo;
	for(LL l=2,r;l<=x;l=r+1){
		r=x/(x/l);
		res=(mo+res-S(x/l)*(r-l+1)%mo)%mo;
	}
	return (val[x]=res);
}

LL calc(){
	LL res=0;
	for(LL l=1,r;l<=n;l=r+1){
		r=n/(n/l);
		res=(res+(n/l)*(n/l)%mo*(S(r)-S(l-1)+mo)%mo)%mo;
	}
	return res;
}

int main(){
	scanf("%lld%lld",&n,&k);
	sz=sqrt(n);pcnt(sz);Fprime();
	printf("%lld",calc());
	return 0;
}