[ACNOI2022]求和

题目

题目背景
蛋教筛终于取代了杜教筛，飘飘蛋名震四方！

“我有一颗蛋，足以慰风尘”。——《咏蛋》

题目描述
设 $x$ 的质因数分解为 ${\prod }_i{p}_i^{a_i}$，定义
$$f_k(x)=\prod _i(-1{)}^{a_i}[a_i⩽k]$$

不难发现 $f_1(x)=\mu (x)$ 。现在需求出
$$\sum _{i=1}^n\sum _{j=1}^n\sum _{d=1}^k{f}_d(\gcd (i,j))$$

答案对 $2^{30}$ 取模。

数据范围与提示
$n⩽10^{10}$ 而 $k⩽40$ 。

思路

不难发现答案是 $(2\phi -ϵ)\ast ({\sum }_i{f}_i)$ 的前缀和，其中 $\ast$ 代表狄利克雷卷积。

如果对每个 $f_i$ 单独求解，显然可以做到 $\mathcal{O}(\frac{n^{0.75}k}{\ln n})$ 。如果你跑得像兔子一样快，我不做评价。

观察 $f_k$ 的贝尔生成函数，不太能杜教筛。但是可以 $\text{powerful numbers}$ 。用什么拟合 $f_k$ 呢？要用 $\mu$，因为这样又有了一个美丽的巧合：$\frac{f_k}{\mu }$ 的每一项要么是 $0$ 要么是 $1$ 。因此，在搜索 $\text{powerful numbers}$ 的时候，可以同时用一个 $k\text{-bit}$ 的二进制数，记录哪些 $f_i(x)=1$ 。预处理这个简单的表，转移是 $\mathcal{O}(1)$ 的，且可以得到每个 $(2\phi -ϵ)\ast \mu$ 的前缀和的贡献系数。

显然 $\phi \ast \mu$ 可以杜教筛，时间复杂度 $\mathcal{O}(n^{2/3})$ 。

总结一下就是：一方面，$n=10^{10}$ 的时候 $n^{2/3}$ 算出来就是比 $\frac{n^{3/4}}{\ln n}$ 小嘛；另一方面，类似于 $\text{bitset}$ 的优化是没有减少计算量但是加快了运行速度的 😢

代码

#include <cstdio>
#include <algorithm> // Almighty XJX yyds!!
#include <cstring> // oracle: ZXY yydBUS!!!
#include <cctype> // Huge Egg Dog eat me!!!
#include <cstdint>
#include <cassert>
# define rep(i,a,b) for(int i=(a); i<=(b); ++i)
# define drep(i,a,b) for(int i=(a); i>=(b); --i)
# define rep0(i,a,b) for(int i=(a); i!=(b); ++i)
using llong = long long;
inline int readint(){
	int a = 0, c = getchar(), f = 1;
	for(; !isdigit(c); c=getchar()) if(c == '-') f = -f;
	for(; isdigit(c); c=getchar()) a = a*10+(c^48);
	return a*f;
}

const int MAXN = 5000000;
bool is_prime[MAXN];
unsigned primes[MAXN], primes_size, least[MAXN];
unsigned phi[MAXN], mu[MAXN], mp[MAXN];
void sieve(unsigned n = MAXN-1){
	memset(is_prime+2, true, n-1);
	phi[1] = mu[1] = mp[1] = 1;
	for(unsigned i=2,&len=primes_size=0; i<=n; ++i){
		if(is_prime[i]){
			primes[++ len] = least[i] = i;
			phi[i] = i-1, mu[i] = -1, mp[i] = i-2;
		}
		for(unsigned j=1; j<=len&&primes[j]*i<=n; ++j){
			const unsigned to = i*primes[j];
			is_prime[to] = false;
			if(!(i%primes[j])){
				least[to] = least[i]*primes[j];
				if(least[i] == i) // power of prime
					mp[to] = (phi[to]=to-i)-phi[i];
				else{
					phi[to] = phi[least[to]]*phi[i/least[i]];
					mp[to] = mp[least[to]]*mp[i/least[i]];
				}
				mu[to] = 0; break;
			}
			least[to] = primes[j], mu[to] = -mu[i];
			phi[to] = phi[i]*(primes[j]-1);
			mp[to] = mp[i]*(primes[j]-2); // J-xing
		}
	}
}

const int SQRTN = 100000;
llong w[SQRTN]; int haxi[SQRTN];
# define _id(x) haxi[n/(x)] // only for x > SQRTN
unsigned mu_sum[SQRTN], phi_sum[SQRTN], mp_sum[SQRTN];

const int LOGN = 33;
# define bitcnt(x) unsigned(__builtin_popcountll(x))
uint64_t choice[LOGN+1]; unsigned ans;
void dfs(const llong &n, int max_id, llong nowx, uint64_t nowv){
	if(!nowv) return; // no more chances
	{ // contribute to answer
		llong v = n/nowx; // what to multiply
		if(v < MAXN) ans += ((mp[v]<<1)-mu[v])*bitcnt(nowv);
		else v = _id(v), ans += (2*mp_sum[v]-mu_sum[v])*bitcnt(nowv);	
	}
	for(int i=1; i<=max_id; ++i){
		if(nowx*primes[i]*primes[i] > n) break;
		llong nxtx = nowx*primes[i];
		for(int j=2; true; ++j){
			if((nxtx *= primes[i]) > n) break;
			dfs(n, i-1, nxtx, nowv&choice[j]);
		}
	}
}

int main(){
	sieve(); // pre-compute
	rep0(i,2,MAXN){
		mu[i] += mu[i-1], phi[i] += phi[i-1];
		mp[i] += mp[i-1]; // all needed
	}
	llong n; scanf("%lld",&n);
	int k = readint(), tot = 0;
	for(llong l=1; true; l=n/w[tot]+1){
		if(n/l < MAXN) break; // done
		w[++tot] = n/l, _id(w[tot]) = tot;
	}
	for(int i=tot; i; --i){
		mu_sum[i] = 1; // mu * I = e
		phi_sum[i] = unsigned(w[i]*(w[i]+1)>>1);
		for(llong l=2,v,r; l!=w[i]+1; l=r){
			v = w[i]/l, r = w[i]/v+1;
			const unsigned len = unsigned(r-l);
			if(v < MAXN){ // transfer together
				mu_sum[i] -= len*mu[v];
				phi_sum[i] -= len*phi[v];
				mp_sum[i] -= len*mp[v];
			}
			else{
				v = _id(v); // uniform index
				mu_sum[i] -= len*mu_sum[v];
				phi_sum[i] -= len*phi_sum[v];
				mp_sum[i] -= len*mp_sum[v];
			}
		}
		mp_sum[i] += phi_sum[i];
	}
	rep(i,1,k){ // for each kind
		for(int j=2; j<=i&&j<=LOGN; j+=2) choice[j] |= 1ull<<i;
		if(!(i&1)) rep(j,i+1,LOGN) choice[j] |= 1ull<<i;
	}
	dfs(n, primes_size, 1, (2ull<<k)-2);
	ans &= (1u<<30)-1; // modulus
	printf("%u\n", ans);
	return 0;
}