2023“钉耙编程”中国大学生算法设计超级联赛（5）1002 GCD Magic

终于有一道赛时写出来的数学题了，赶紧写写想法

题意

求
$$\sum _{i=1}^n\sum _{j=1}^n[gcd(2^i-1,2^j-1){]}^K$$

前置芝士

莫比乌斯反演、二项式展开、杜教筛、数论分块

切入

首先需要想到一个转化$gcd(2^i-1,2^j-1)=2^{gcd(i,j)}-1$

通过$2^i-1$的二进制是$(\mathrm{111...11}{)}_2$感性理解

这里我的想法是一个二进制位全是1的数只能被另一个二进制位全是1的数整除，而且这另一个二进制位全是1的数的长度也一定是原来那个数的长度的约数。

所以最后会转化成二进制全是$1$的一个二进制位最长的数也就是长度为$gcd(i,j)$

严谨的数学推理这里用题解的数学推理：
$$(2^i-1,2^j-1)[j>i]=(2^i-1,2^j-2^i)=(2^i-1,2^i(2^{j-i}-1))=(2^i-1,2^{j-i}-1)$$
这个步骤就是辗转相除法的步骤

题解

$$\begin{gathered} \sum _{i=1}^n\sum _{j=1}^n[gcd(2^i-1,2^j-1){]}^K \\ \sum _{i=1}^n\sum _{j=1}^n(2^{gcd(i,j)}-1{)}^K \\ =\sum _{d=1}^n\sum _{i=1}^n\sum _{j=1}^n(2^d-1{)}^K[gcd(i,j)=d] \\ 这里是最基础的一步，几乎所有莫比乌斯反演题都是从这里开始 \\ =\sum _{d=1}^n\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{d}\rfloor }(2^d-1{)}^K[gcd(i,j)=1] \\ =\sum _{d=1}^n(2^d-1{)}^K\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{d}\rfloor }[gcd(i,j)=1] \\ 使f(n)=\sum _{i=1}^n\sum _{j=1}^n[gcd(i,j)=1] \\ 有f(n)=\sum _{d=1}^n\mu (d)\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{d}\rfloor }1 \\ 即f(n)=\sum _{d=1}^n\mu (d)\lfloor \frac{n}{d}{\rfloor }^2 \\ 此时原式可以写作 \\ \sum _{d=1}^n(2^d-1{)}^{K}f(\lfloor \frac{n}{d}\rfloor ) \\ 展开 \\ \sum _{d=1}^n(2^d-1{)}^K\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\mu (i)\lfloor \frac{n}{d\ast i}{\rfloor }^2 \\ 我们知道后面这个f(n)可以使用杜教筛与数论分块的方式处理出来 \\ 而前面这个式子无法进行前缀和优化,所以需要一个转化 \\ 将(2^d-1{)}^K拆开可以得到 \\ \sum _{p=0}^K{C}_K^p{2}^{p\ast d}(-1{)}^{p-d} \\ 转移到原式中并进行sigma交换 \\ \sum _{p=0}^K{C}_K^p(-1{)}^{K-p}[\sum _{d=1}^n{2}^{d\ast p}\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\mu (i)\lfloor \frac{n}{d\ast i}{\rfloor }^2] \\ 这样就可以用等差数列求和去计算2的幂次求和 \end{gathered}$$

关于数论分块的一些解析

我们知道形如${\sum }_{d=1}^{n}f(d){\sum }_{i=1}^{\lfloor \frac{n}{d}\rfloor }g(i)$可以进行数论分块的基础就是可以快速求出${\sum }_{i=l}^{r}f(i)$

这样在两个指针跳动的时候才可以求出整块的答案

代码

#include <bits/stdc++.h>

const int N = 1e7 + 6, mod = 998244353;
int n, k, C[20][20], INV[20];

inline int add(int a, int b) { return (a += b) >= mod ? a - mod : a; }
inline int sub(int a, int b) { return (a -= b) < 0 ? a + mod : a; }
inline int mul(int a, int b) { return 1LL * a * b % mod; }
inline void Add(int &a, int b) { a = add(a, b); }
inline void Sub(int &a, int b) { a = sub(a, b); }
inline void Mul(int &a, int b) { a = mul(a, b); }

inline int quickpow(int a, int b) {
    int re = 1;
    while(b) {
        if(b & 1) Mul(re, a);
        b >>= 1;
        Mul(a, a);
    }
    return re;
}

int prime[N], mu[N], m;
bool vis[N + 1];
inline void Prime(int x) {
	mu[1] = 1;
	for(int i = 2; i <= x; ++i) {
		if(!vis[i])	prime[++m] = i, mu[i] = mod - 1;
		for(int j = 1; j <= m && 1LL * prime[j] * i <= x; ++j) {
			vis[prime[j] * i] = 1;
			if(!(i % prime[j]))	break;
			mu[i * prime[j]] = sub(mod, mu[i]);
		}
	}
    for(int i = 1; i <= x; ++i) mu[i] = add(mu[i], mu[i - 1]);
}

//从1到n的和
std::unordered_map < int , int > ansmu, anscal1;
//杜教筛莫比乌斯函数（模板）其实也是数论分块来的
inline int ans_mu(int x) {
	if(x <= 1000000) return mu[x];
	if(ansmu[x]) return ansmu[x];
	int re = 0;
	for(int l = 2, r = 0; l <= x; l = r + 1)
		r = x / (x / l), Add(re, mul((r - l + 1), ans_mu(x / l)));
	return ansmu[x] = sub(1, re);
}

inline int calc(int n) {//计算f(n)
    if(anscal1[n]) return anscal1[n];
    int re = 0;
    for(int l = 1, r = 0; l <= n; l = r + 1) {
        r = n / (n / l);
        int cnt = n / l;
        Add(re, mul(mul(cnt, cnt), sub(ans_mu(r), ans_mu(l - 1))));
    }
    return anscal1[n] = re;
}

inline int calc2(int n, int p) {
    int re = 0;
    for(int l = 1, r = 0; l <= n; l = r + 1) {
        r = n / (n / l);
        if(p == 0) Add(re, mul(r - l + 1, calc(n / l)));//特别的
        else Add(re, mul(mul(sub(quickpow(2, 1LL * (r + 1) * p % (mod - 1)), quickpow(2, 1LL * l * p % (mod - 1))), INV[p]), calc(n / l)));
    }
    return re;
}

void solve() {
    std::cin >> n >> k;
    int ans = 0;
    for(int p = 0; p <= k; ++p) {
        int pos = C[k][p];
        if((k - p) & 1) pos = sub(mod, pos);
        Add(ans, mul(calc2(n, p), pos));
    }
    std::cout << ans << '\n';
}

signed main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);
    Prime(1000000);
    for(int i = 0; i <= 10; ++i) {
        C[i][0] = 1;
        for(int j = 1; j <= i; ++j) {
            C[i][j] = add(C[i - 1][j], C[i - 1][j - 1]);
        }
    }
    for(int i = 15; i >= 1; --i) INV[i] = quickpow(sub(quickpow(2, i), 1), mod - 2);
    int T = 1;
    std::cin >> T;
    while(T--) solve();
    return 0;
}

tips

这道题有点卡常所以需要在计算$f(n)$的时候用一个桶记录下来答案