@51nod - 1847@ 奇怪的数学题


@description@

给出 N, K, 请计算下面这个式子:

\[\sum_{i=1}^{N}\sum_{j=1}^{N}sgcd(i, j)^k \]

其中, sgcd(i, j)表示(i, j)的所有公约数中第二大的,特殊地,如果gcd(i, j) = 1, 那么sgcd(i, j) = 0。

考虑到答案太大,请输出答案对2^32取模的结果.

原题传送门。

@solution@

记 mnp[i] 表示 i 的最小质因子。则:

\[\begin{aligned} \sum_{i=1}^{N}\sum_{j=1}^{N}sgcd(i, j)^k =& \sum_{d=1}^{N}\sum_{i=1}^{N}\sum_{j=1}^{N}[\gcd(i, j) = d]\times (\frac{d}{mnp[d]})^{k}\\ =& \sum_{d=1}^{N}(\frac{d}{mnp[d]})^{k}\sum_{i=1}^{\lfloor \frac{N}{d}\rfloor}\sum_{j=1}^{\lfloor \frac{N}{d}\rfloor}[\gcd(i, j) = 1]\\ =& \sum_{d=1}^{N}(\frac{d}{mnp[d]})^{k}(\sum_{i=1}^{\lfloor \frac{N}{d}\rfloor}\phi(i) - 1) \end{aligned}\]

后面那个 \(\phi\) 前缀和,杜教筛基础,不讲了。前面整除分块,接下来讲讲具体怎么求。

\(S(N) = \sum_{d=1}^{N}(\frac{d}{mnp[d]})^{k}\),考虑怎么求 \(S\)。分成两部分:d 为质数,贡献为 \(\pi(N)\)(N 以内质数个数);d 为合数。
注意到这种涉及到最小质因子的可以类比 min-25 筛的思想搞(其实就是埃氏筛法)。

我们做 min-25 筛前半部分时,有一个辅助数组 \(f(N, p)\) 表示 N 以内的 质数 或 最小质因子 > 第 p 个质数的数 的 K 次幂之和。
利用 \(f(N, p)\) 就可以得到 N 以内最小质因子 >= 第 p 个质数的数的 K 次幂之和 \(g(N, p)\),这一点很容易做到。

然后有:

\[S(N) = \pi(N) + \sum_{i=1}^{pcnt}g(\lfloor \frac{N}{prime_i}\rfloor, i) \]

自然数幂和用斯特林数处理好些,因为模数不是质数,不一定有逆元。

@accepted code@

#include <cstdio>
#include <algorithm>
using namespace std;

const unsigned long long MOD = (1LL << 32);

typedef unsigned int ui;

const int MAXN = 1000000;

ui S[55][55], inv[55];
void init() {
	for(int i=0;i<=50;i++) {
		S[i][0] = (i == 0);
		for(int j=1;j<=i;j++)
			S[i][j] = S[i-1][j-1] + S[i-1][j] * ui(j);
	}
	inv[1] = 1;
	for(int i=3;i<=51;i+=2) {
		for(int j=1;j<=i;j++)
			if( (j*MOD + 1) % i == 0 ) {
				inv[i] = (j*MOD + 1) / i % MOD;
				break;
			}
	}
}
ui sum(ui n, int k) {
	int cnt = 0; ui ans = 0, tmp = 1; n++;
	for(int i=0;i<=k;i++) {
		ui del = (n - i), pw = 0;
		if( del == 0 ) break;
		while( del % 2 == 0 ) del /= 2, cnt++;
		tmp *= del;
		
		del = i + 1;
		while( del % 2 == 0 ) del /= 2, pw++;
		
		if( (cnt - pw) < 60 ) {
			ui p = ui(1LL << (cnt - pw));
			ans += S[k][i]*tmp*inv[del]*p;
		}
	}
	
	return ans;
}

int prm[MAXN + 5], pcnt;
bool nprm[MAXN + 5]; ui phi[MAXN + 5];
void sieve() {
	phi[1] = 1;
	for(int i=2;i<=MAXN;i++) {
		if( !nprm[i] ) prm[++pcnt] = i, phi[i] = i - 1;
		for(int j=1;i*prm[j]<=MAXN;j++) {
			nprm[i*prm[j]] = true;
			if( i % prm[j] == 0 ) {
				phi[i*prm[j]] = phi[i]*prm[j];
				break;
			}
			else phi[i*prm[j]] = phi[i]*phi[prm[j]];
		}
	}
	for(int i=1;i<=MAXN;i++)
		phi[i] += phi[i - 1];
}

int n, k;
ui powk(ui b) {
	ui ret = 1;
	for(int i=k;i;i>>=1,b*=b)
		if( i & 1 ) ret *= b;
	return ret;
}

int a[MAXN + 5], id1[MAXN + 5], id2[MAXN + 5], cnt;
int id(int x) {return x <= MAXN ? id1[x] : id2[n / x];}

bool tag[MAXN + 5]; ui sp[MAXN + 5];
ui phisum(int n) {
	if( n <= MAXN ) return sp[id(n)] = phi[n];
	if( tag[id(n)] ) return sp[id(n)];
	ui ans = sum(n, 1);
	for(int i=2;i<=n;i++) {
		int p = n / i, j = n / p;
		ans -= phisum(p) * (sum(j, 0) - sum(i - 1, 0));
		i = j;
	}
	tag[id(n)] = true; return sp[id(n)] = ans;
}
ui f[MAXN + 5], g[MAXN + 5], h[MAXN + 5];
void solve() {
	for(int i=1;i<=cnt;i++)
		f[i] = a[i] - 1, g[i] = sum(a[i], k) - 1, h[i] = 0;
	ui tmp = 0;
	for(int i=1;i<=pcnt;i++) {
		long long p = 1LL*prm[i]*prm[i]; ui pw = powk(prm[i]);
		if( p > n ) break;
		for(int j=1;j<=cnt;j++) {
			if( p > a[j] ) break;
			h[j] += (g[id(a[j]/prm[i])] - tmp);
			f[j] -= f[id(a[j]/prm[i])] - (i - 1);
			g[j] -= (g[id(a[j]/prm[i])] - tmp)*pw;
		}
		tmp += pw;
	}
	for(int i=1;i<=cnt;i++) h[i] += f[i];
}
void get_id() {
	for(int i=1;i<=n;i=n/(n/i)+1) {
		int p = n / i;
		if( p <= MAXN ) id1[p] = (++cnt);
		else id2[n / p] = (++cnt);
		a[cnt] = p;
	}
}
int main() {
	init();
	sieve(), scanf("%d%d", &n, &k), get_id(), solve();
	
	ui ans = 0;
	for(int i=2;i<=n;i++) {
		int p = n / i, j = n / p;
		ans += (2*phisum(p) - 1)*(h[id(j)] - h[id(i - 1)]);
		i = j;
	}
	printf("%lld\n", (long long)ans);
}

@details@

妄想用线性求逆元然后发现这种方法并不适用于模数为合数的情况。当时我人就傻在那里。

然后不想打 exgcd,于是突然脑洞大开写了暴力解不定方程求逆元。

posted @ 2020-03-19 22:05  Tiw_Air_OAO  阅读(218)  评论(0编辑  收藏  举报