hdu 6069 Counting Divisors

Problem

acm.hdu.edu.cn/showproblem.php?pid=6069

Meaning

定义函数d(n) = n 的因子个数。
给定区间 [ l , r ] 和常数 k，求

$$(\ \sum^{r}_{i\ =\ l}d(i^{k})\ )\ mod\ 998244353$$

Analysis

朴素的想法，先把每一个 i 都质因数分解，$i = p_{1}^{e_{1}}…p_{n}^{e_{n}}$，于是 $i^{k}$ 对答案的贡献就是 $(e_{1}k+1)(e_{2}k+1)…(e_{n}k+1)$，最后答案就是：

$$(\ \sum_{i=l}^{r}\prod_{j=1}^{n}(e_{j}k+1)\ )\ mod\ 998244353$$
就先打一个 $10^{6}$ 的素数表（如果有大于 $10^{6}$ 的质因子，那么只能有一个，并且指数为 1），然后枚举 [ l , r ] 里的数，对每一个都做一次质数分解、统计答案，超时。
一种优化：对于某个 i，质因子分别为 $p_{1}、…、p_{n}$，那么在处理完 $i^{k}$ 的贡献后，可以一并处理掉 $(ie_{1})^{k}、…、(ie_{n})^{k}$ 的贡献，因为它们只间只有其中一个因式 $(e_{j}k+1)$ 变成 $((e_{j}+1)k+1)$ 这一个差别，这样就可以少做一些质因数分解（有点像素筛的过程），但还是超时。
可能主要是因为每次质因数分解都要扫一遍质数表，做了很多无用功。
反过来想，不从 i 分解质因子，而是用个数组装下区间[l , r]和它们的贡献（初始为 1），枚举 $10^{6}$ 内的质数，对每一个质数 p，在区间内找它的倍数来拆（分解），拆完之后更新这个数的贡献，最后再把贡献求和。

Code

#include <cstdio>
#include <cstring>
using namespace std;
const int N = 1000000, MOD = 998244353;

int prime[N+1], top; // prime table

void sieve()
{
    top = 0;
    memset(prime, -1, sizeof prime);
    for(int i = 2; i <= N; ++i)
    {
        if(prime[i])
            prime[top++] = i;
        for(int j = 0; j < top && prime[j] <= N / i; ++j)
        {
            prime[i * prime[j]] = 0;
            if(i % prime[j] == 0)
                break;
        }
    }
}

long long cntr[N+1]; // 区间内每个数的贡献
long long num[N+1]; // 区间内的数

int main()
{
    sieve();
    int T;
    scanf("%d", &T);
    while(T--)
    {
        long long l, r, k;
        scanf("%I64d%I64d%I64d", &l, &r, &k);
        for(int i = 0; i <= r - l; ++i)
        {
            num[i] = l + i;
            cntr[i] = 1;
        }
        for(int i = 0, p, beg; i < top; ++i)
        {
            p = prime[i];
            beg = (p - l % p) % p; // 区间里第一个是p倍数的数的下标
            // 枚举 p 的倍数
            for(int j = beg, cnt; j <= r - l; j += p)
            {
                // 质数分解
                for(cnt = 0; num[j] % p == 0; num[j] /= p)
                    ++cnt;
                // 更新贡献
                cntr[j] = cntr[j] * (cnt * k + 1) % MOD;
            }
        }
        // 补上大于 1e6 的质因数的贡献
        for(int i = 0; i <= r - l; ++i)
            if(num[i] > 1)
                cntr[i] = cntr[i] * (k + 1) % MOD;
        long long ans = 0;
        for(int i = 0; i <= r - l; ++i)
            ans = (ans + cntr[i]) % MOD;
        printf("%I64d\n", ans);
    }
    return 0;
}

TLE Code

#include <cstdio>
#include <cstring>
using namespace std;
const int N = 1000000, MOD = 998244353;

struct node
{
    int v; // 因子的值
    int num; // 因子的指数
} dv[N+10]; // 分解出来的因子

long long prime[N+10], top, d;
bool vis[N+10]; // 标记这个数有没有被处理过

void seive(int n) {
    top = 0;
    memset(prime, -1, sizeof prime);
    for(int i = 2; i <= N; ++i)
        if(prime[i])
        {
            prime[top++] = i;
            for(int j = i; j <= N; j+=i)
                prime[j] = 0;
        }
}

// 快速幂
long long fpow(long long a, long long b)
{
    long long ans = 1;
    for( ; b; b >>= 1, a = a * a % MOD)
        if(b & 1)
            ans = ans * a % MOD;
    return ans;
}

// 逆元
long long inv(long long a)
{
    return fpow(a, MOD - 2);
}

// 质因数分解
void dvide(long long x)
{
    d = 0;
    for(int i = 0; i < top && prime[i] * prime[i] <= x; ++i)
        if(x % prime[i] == 0)
        {
            dv[d].v = prime[i];
            dv[d].num = 0;
            while(x % prime[i] == 0)
            {
                ++dv[d].num;
                x /= prime[i];
            }
            ++d;
        }
    if(x != 1)
    {
        dv[d].v = x;
        dv[d++].num = 1;
    }
}

int main() {
    int t;
    scanf("%d", &t);
    seive(N);
    while(t--) {
        long long l, r;
        int k;
        scanf("%lld%lld%d", &l, &r, &k);
        memset(vis, false, sizeof vis);
        long long ans = 0;
        for(long long i = l, tmp; i <= r; ++i)
        {
            if(vis[i - l])
                continue;
            vis[i - l] = true;
            dvide(i); // 分解
            tmp = 1; // 贡献
            for(int j = 0; j < d; ++j)
                tmp = tmp * (dv[j].num * k + 1) % MOD;
            ans = (ans + tmp) % MOD;
            // 顺便处理 i * ej 的贡献
            for(int j = 0; j < d && i * dv[j].v <= r; ++j)
            {
                if(vis[i * dv[j].v - l])
                    continue;
                vis[i * dv[j].v - l] = true;
                ans = (ans + tmp + (tmp * inv(dv[j].num * k + 1) % MOD * k)) % MOD;
            }
        }
        printf("%lld\n", ans % MOD);
    }
    return 0;
}