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

1002题

---

杜教筛（min_25筛，洲阁筛应该都可以），数论分块

没本事，做数学一般先打表，k这个幂比较麻烦，所以先考虑k=1的情况，打出n从1到5的表如下
$$[1]$$
$$[1,1,1,3]$$
$$[1,1,1,1,3,1,1,1,7]$$
$$[1,1,1,1,1,3,1,3,1,1,7,1,1,3,1,15]$$
$$[1,1,1,1,1,1,3,1,3,1,1,1,7,1,1,1,3,1,15,1,1,1,1,1,31]$$
看不出什么东西啊，只能看出1很多，不妨将这些数再打个表，得到
$$1:1$$
$$1:3,3:1$$
$$1:7,3:1,7:1$$
$$1:11,3:3,7:1,15:1$$
$$1:19,3:3,7:1,15:1,31:1$$

明确一下我们的目标，我们需要求所有的值，即上面序列中个数和值相乘的和，显然，这个值一看就是$$2^i-1$$这种样子的，考虑后面的系数，不妨再打个表，得到

[0, 0, 1, 1, 3, 3, 7, 7, 11, 11, 19, 19, 23, 23, 35, 35, 43, 43, 55, 55, 63, 63, 83, 83, 91, 91, 115, 115, 127, 127]

这种东西

那么这个一看就很有规律吧，考虑如何找到这其中的规律，我当然不会告诉你我是上oeis上面搜的，给出一种比较规范的做法。

不妨将题目中给的式子做一下划分，显然我们已经知道值是2的幂-1，
那么原式变为
$$\sum _{d=1}^n\sum _{i=1}^n\sum _{j=1}^n[gcd(i,j)==d](2^d-1{)}^k$$
是不是这边要求的系数即是
$$\sum _{i=1}^n\sum _{j=1}^n[gcd(i,j)==d]$$
的个数

如果看过oi-wiki上莫比乌斯反演的内容，会发现，这个式子就是习题一的样子，化简式子得到

$$\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{d}\rfloor }I[gcd(i,j)==1]$$

此处与wiki上推导不同，上式中I表示示性函数，条件成立为1，否则为0.

考虑i和j的对称性，可以推导得到

$$\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }(I(gcd(i,i)==1)+2\sum _{j=1}^{i-1}I(gcd(i,j)==1))$$

这个式子就是将i和j分开，因为想想因子的gcd，i与j互换不影响结果，如果j等于i那么无需乘2.

显然

$$gcd(i,i)=1$$

得到

$$\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }(-I(gcd(i,i)==1)+2\sum _{j=1}^{i}I(gcd(i,j)==1))$$

就是将前面的一个1放到后面去，方便下面操作

应该可以看出，2乘后面的式子就是欧拉函数的定义

那么就可以得到

$$\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }(-I(i==1)+2\varphi (i))$$

显然其中那个i不受外面的求和符号影响，故可以得到

$$(-1+2\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\varphi (i))$$

这个式子的意思就是欧拉函数的前缀和，去oeis上面也是这个式子。

对于欧拉函数的前缀和，一般是采用线性筛来做的，但是此处数据量达到1e9，显然不可用，考虑数论分块，其实上面的推导是没啥用的，因为杜教筛求积性函数前缀和的过程中，我们推导这个式子，可以往莫比乌斯函数的方向推导，得到

$$\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{d}\rfloor }(gcd(i,j)==1)=\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{d|i,d|j}\mu (d)=\sum _{d=1}^{\lfloor \frac{n}{d}\rfloor }\mu (d)\lfloor \frac{n}{d}{\rfloor }^2$$

这个式子的右半部分只需要转变求和顺序即可，$$1\sim \lfloor \frac{n}{d}\rfloor$$中1的倍数就是这个值。

由此我们只需要求到莫比乌斯函数的前缀和。
$$ϵ=[n=1]=\mu \ast 1=\sum _{d|n}\mu (d)$$
$$S_1(n)=\sum _{i=1}^nϵ(i)-\sum _{i=2}^n{S}_1(\lfloor \frac{n}{i}\rfloor )=1-\sum _{i=2}^n{S}_1(\lfloor \frac{n}{i}\rfloor )$$

这个式子是可以使用杜教筛的，由此可以得解。
还用直接用杜教筛来直接求解欧拉函数的前缀和的，可以去oi-wiki上看。

min_25筛和州阁筛也是可以做的，但我没有写出来。

上面那个式子只是为了更方便理解系数的意义，如果看到了那个式子，肯定会想到如何求欧拉函数前缀和，就有方向了，当然，一般能推出来的也不需要知道这个式子，这只是方便我这种打表用户。

由此，解决了题目的一大难题，下面是我认为的第二大难题，对$$2^i-1$$的求和。

这个东西其实是很简单(贱)的，考虑二项展开

$$(Ax+B{)}^k=\sum _{i=0}^k{A}^i{B}^{k-i}{x}^i$$

你明白我把$$2^i-1$$理解成一个新数然后后推的痛苦吗！

这个式子一看就知道$$A=2^0,B=-1,x=2^i$$

那么就不用多说了，直接遍历k就行，注意当k=0时，这个等比数列求和是0，要把常数项拿出来处理，然后就没有坑点了。

时间复杂度:$$O(k{n}^{\frac{2}{3}}t)$$(其中t远小于$$2\sqrt{n}$$)
空间复杂度:$$O(\sqrt{n})$$

我代码是超时的，但是我感觉这个方法是对的，可以指导我一下（求求了）

#include <cmath>
#include <algorithm>
#include <cstdio>
#include <cstring>
#include <map>
#include <iostream>
#include <bits/stdc++.h>
using namespace std;
const int maxn = 2000010;
const long long mod = 998244353;
long long T, n, k,res,temp,temp1,temp2,temp3,pri[maxn], cur, mu[maxn], sum_mu[maxn];
bool vis[maxn];


map<long long, long long> mp_mu;

long long S_mu(long long x) {  // 求mu的前缀和
    if (x < maxn) return sum_mu[x];
    if (mp_mu[x]) return mp_mu[x];  // 如果map中已有该大小的mu值，则可直接返回
    long long ret = (long long)1;
    for (long long i = 2, j; i <= x; i = j + 1) {
        j = x / (x / i);
        ret -= S_mu(x / i) * (j - i + 1);
    }
    return mp_mu[x] = ret;  // 路径压缩，方便下次计算
}

long long S_phi(long long x) {  // 求phi的前缀和
    long long ret = (long long)0;
    long long j;
    for (long long i = 1; i <= x; i = j + 1) {
        j = x / (x / i);
        ret += (S_mu(j) - S_mu(i - 1)) * (x / i) * (x / i);
    }
    return (ret - 1) / 2 + 1;
}
long long gg(long long x){
    return 2 * S_phi(x) - 1;

}
long long  int qmi1(long long int  a, long long int  b, long long int  p)
{
    long long res1 = 1;
    while(b)
    {
        if(b & 1) res1 = (long long ) res1 * a % p;
        b >>= 1;
        a = (long long ) a * a % p;
    }
    return res1;
}
int c[20][20];
int main() {
    for(int i=0;i<11;i++)
    {
        for(int j = 0;j<=i;j++)
        {
            if(!j)
                c[i][j]=1;
            else
            {
                c[i][j]=(c[i-1][j]+c[i-1][j-1])%mod;
            }
        }
    }
    mu[1] = 1;
    for (int i = 2; i < maxn; i++) {  // 线性筛预处理mu数组
        if (!vis[i]) {
            pri[++cur] = i;
            mu[i] = -1;
        }
        for (int j = 1; j <= cur && i * pri[j] < maxn; j++) {
            vis[i * pri[j]] = true;
            if (i % pri[j])
                mu[i * pri[j]] = -mu[i];
            else {
                mu[i * pri[j]] = 0;
                break;
            }
        }
    }
    for (int i = 1; i < maxn; i++){
        sum_mu[i] = sum_mu[i - 1] + mu[i];  // 求mu数组前缀和
    }
    scanf("%lld",&T);
    while(T--){
        res = 0;
        scanf("%lld %lld",&n,&k);
        temp = (long long)sqrt(n)+1;
        res = 0;
        for(int i = 1 ;i < temp;i++){
            temp1 = n / i;
            (res += gg(temp1)*(qmi1((long long)2,(long long)(i),mod)-(long long)1))%=mod;
        }
        for(int i = 1 ;i < temp;i++){
            temp1 = n/i;
            temp2 = n/(i+1);
            temp1 = temp1 - temp2;
            temp3 = 0;
            for(int l = k;l>=0;l--){
                long long tt = qmi1(2,temp2 + 1,mod);
                (temp3 += (qmi1(tt,(long long)l,mod) * c[k][l] * ((k-l)&1 ? -1 : 1) * max(qmi1(qmi1(2,l,mod),temp1,mod)-1,(long long )1) * max((qmi1((qmi1(2,l,mod)-1),mod-2,mod)),(long long )1)+mod))%=mod;
                if (l == 0){
                    temp3 += qmi1(tt,(long long)l,mod) * c[k][l] * ((k-l)&1 ? -1 : 1)*max((long long )0,n-2);
                }
            }
            (res += temp3 * gg(i))%=mod;
        }
        printf("%lld\n",res);
    }
    return 0;
}