题解 P6060 【[加油武汉]传染病研究】
sochiji
2020-02-10 14:52:02
### 蒟蒻来抛砖引玉写一发50分的题解。
## 定义
函数 $D(x),x\in\mathbb{Z_+}$表示$x$的约数个数。
>例：$D(6)=4,\quad D(7)=2$
## 单个函数值的一般计算方法
设$x=p_1^{\alpha _1}\cdot p_2^{\alpha _2}\cdot ... \cdot p_k^{\alpha _k}$，其中$p_i$为素数且互不相同，$\alpha _i$为正整数。
对于$x$的任意正因数$d=p_1^{\beta _1}\cdot p_2^{\beta _2}\cdot ... \cdot p_k^{\beta _k}$和任意$1\leq i\leq k$，有$0\leq\beta _i\leq\alpha_i$;
根据乘法原理，$\displaystyle D(x)=\prod^k_{i=1}(\alpha _i+1)$
## 性质
1. 若$x_1,x_2$互素，有$D(x_1\cdot x_2)=D(x_1)\cdot D(x_2)$；
>证明：
>设$\displaystyle x_1=\prod^k_{i=1}p_i^{\alpha_i},\quad x_2=\prod^{k+l}_{i=k+1}p_i^{\alpha_i}$；
>由于对于任意$1\leq i\leq k,\quad k+1\leq j\leq k+l$，有$p_i\neq p_j$，
>所以$\displaystyle D(x_1\cdot x_2)=\prod^{k+l}_{i=1}(\alpha _i+1)=\prod^{k}_{i=1}(\alpha _i+1)\cdot \prod^{k+l}_{i=k+1}(\alpha _i+1)=D(x_1)\cdot D(x_2)$；
>**性质1**表明$D(x)$是**积性函数**。
2. 若$x=p^\alpha$且$p$为素数，则：
- $D(x)=\alpha+1$；
- $D(x^n)=\alpha\cdot n+1=n\cdot D(x)-n+1$；
>证明：
>第一条由上可知显然。
>第二条：$D(x^n)=D(p^{n\alpha})=\alpha\cdot n+1$
>$D(x^n)=n(D(x)-1)+1=n\cdot D(x)-n+1$
3. 若$x=p^\alpha \cdot c$，$p$为素数且$p$与$c$互素，则$D(x^n)=(n\alpha +1)\cdot D(c^n)$；
>证明：
>$D(x^n)=D(p^{n\alpha}\cdot c^n)$；
>因为$p$与$c$互素，所以$p^{n\alpha }$与$c^n$互素，进而有$D(p^{n\alpha}\cdot c^n)=D(p^{n\alpha})\cdot D(c^n)=(n\alpha +1)\cdot D(c^n)$；
---
## 使用线性欧拉筛法生成$k=1$时的函数值表
>阅读之前请确保已经对[与素数及欧拉筛法相关的知识](https://zhuanlan.zhihu.com/p/74477008)有简单的了解；
欧拉筛法使得每个合数仅能通过“它的**最小真因子**乘上它的**最大真因子**”这种方法被构造出，而素数不会被筛掉。我们构造出一个新的合数$x$时，有下面几种情况：
1. $x$本身是素数：
>这种情况下$D(x)=2$
2. $x$的最大真因子$p_2$是素数，且与最小真因子$p_1$**不相等**：
>这种情况可以表示为$x=p_1\cdot p_2$，那么$D(x)=D(p_1\cdot p_2)=D(p_1)\cdot D(p_2)$
3. $x$的最大真因子$p$是素数，且与最小真因子**相等**：
>这种情况可以表示为$x=p\cdot p$，那么$D(x)=D(p^2)=3$
4. $x$的最大真因子$c$是合数，且与最小真因子$p$**互素**：
>这种情况可以表示为$x=p\cdot c$，那么$D(x)=D(p\cdot c)=D(p)\cdot D(c)=2\cdot D(c)$
5. $x$的最大真因子$m$是合数，且$m$是$x$的最小真因子$p$的整数次方的**倍数**，即$m=p^\alpha \cdot c$：
>这种情况可以表示为$x=p\cdot m=p\cdot p^{\alpha}\cdot c=p^{\alpha +1} \cdot c$，其中$p$与$c$互素即$\alpha =\mathrm{ord}_{p}m$，那么$D(x)=D(p^{\alpha +1}\cdot c)=D(p^{\alpha +1})\cdot D(c)=(\alpha +2)\cdot D(c)$
为了使欧拉筛法可以递推$D(x)$，我们不仅得记录每个数的$D(x)$，还有必要记录它的：
1. 最小素因子$p$，代码中表示为`minfac`；
2. 最大阶$\alpha=\mathrm{ord}_px$，代码中表示为`order`；
3. 与$p$互素的`c`；
4. 函数值$D(x)$，代码中表示为`d`
```cpp
struct rec
{
int minfac;
int order;
int c;
int d;
} dlist[10000001];
```
在欧拉筛法的代码中加入计算`dlist[]`的语句即可生成$k=1$时的$D(x^k)$函数值表。
由于此题的测试点的数据较多，可以通过记录$k=1$时函数值表的前$n$项和$\displaystyle\sum_{i=1}^{n}D(i)$来优化。
当$x$是$p$的整数次幂时有$c=1$，因此要注意初值条件：
```
dlist[1].c = 1;
dlist[1].d = 1;
```
线性筛代码：
```cpp
//primeMap[]是素性映射表，primeList[]是素数表，sumList[]是k==1时函数值前缀和表
int getPrimeList(int n, bool primeMap[], int primeList[], int sumList[])
{
//初始化
int count = 0;
for (int i = 2; i <= n; i++)
primeMap[i] = true;
primeMap[0] = false;
primeMap[1] = false;
dlist[1].c = 1;
dlist[1].d = 1;
int half = n / 2;
for (int i = 2; i <= half; i++)
if (primeMap[i])
{
dlist[i] = {i, 1, 1, 2}; //第1种情况
primeList[count] = i;
count++;
for (int j = 0; j < count && i >= primeList[j] && i * primeList[j] <= n; j++)
{
primeMap[i * primeList[j]] = false;
if (i == primeList[j])
{
dlist[i * primeList[j]] = {primeList[j], 2, 1, 3}; //第3种情况
break;
}
else
dlist[i * primeList[j]] = {primeList[j], 1, i, dlist[i].d * dlist[primeList[j]].d}; //第2种情况
}
}
else
for (int j = 0; j < count && i * primeList[j] <= n; j++)
{
primeMap[i * primeList[j]] = false;
if (dlist[i].minfac == primeList[j]) //由于存储了i的最小素因子 故写法与普通线性筛不同
{
dlist[i * primeList[j]] = {primeList[j], dlist[i].order + 1, dlist[i].c, (dlist[i].order + 2) * dlist[dlist[i].c].d}; //第5种情况
break;
}
else
dlist[i * primeList[j]] = {primeList[j], 1, i, 2 * dlist[i].d}; //第4种情况
}
for (int i = half + 1; i <= n; i++)
if (primeMap[i])
{
dlist[i] = {i, 1, 1, 2}; //第1种情况
primeList[count] = i;
count++;
}
for (int i = 2; i <= n; i++)
sumList[i] = sumList[i - 1] + dlist[i].d; //生成前缀和
return count;
}
```
>注意：在非初始化处对`struct`进行整体赋值的语法在`C++11`标准之前是非法的。
据题意，我们应当在`main()`中调用上面的函数，并给参数`n`赋值$10^7$。之后我们就能从`sumList[]`立即输出任何$k=1$时的$\displaystyle\sum_{i=1}^{n}D(i)$。
---
## 根据$k=1$时的$D(x)$函数值表生成$k$取任意值的$D(x^k)$函数值表
根据**性质3**，我们可以在$O(n)$的时间内递推出$\displaystyle\sum_{i=1}^{n}D(i^k)$全部项并求和。代码如下：
```cpp
const int MODNUM = 998244353;
long long dpow[10000001];//存放n^k的函数值，以n为索引
int getSum(int n, int k)
{
long long sum = 0;
sum = 1;
for (int i = 2; i <= n; i++)
{
dpow[i] = (k * dlist[i].order + 1) * dpow[dlist[i].c;
sum += dpow[i];
sum %= MODNUM;
}
return (int)sum;
}
```
在`main()`中如下调用上面的函数并注意初值条件：
```cpp
dpow[1] = 1;
sumlist[1] = 1;
int t;
scanf("%d", &t);
for (int i = 1; i <= t; i++)
{
int n, k;
scanf("%d%d", &n, &k);
printf("%d\n", (k == 1) ? sumlist[n] : getSum(n, k));
}
```
不限制$k=1$时，总的时间复杂度为$O(T\cdot n)$，上面的方法最多需要进行$10^{11}$次加法运算，会导致5个测试点超时。
要AC此题，应该得以$O(T)$的时间复杂度解决整个问题，但本蒟蒻实在找不出$\displaystyle\sum_{i=1}^{n}D(i)$与$\displaystyle\sum_{i=1}^{n}D(i^k)$的关系（应该是一个高次多项式函数），所以这篇题解只值50分。
## 感谢阅读。