SPOJ - DIVCNTK

题目链接：divcntk

Let $\sigma_0(n)$ be the number of positive divisors of $n$.

For example, $\sigma_0(1) = 1$, $\sigma_0(2) = 2$ and $\sigma_0(6) = 4$.

Let

$$S_k(n) = \sum _{i=1}^n \sigma_0(i^k).$$

Given $n$ and $k$, find $S_k(n) \bmod 2^{64}$.

Input
There are multiple test cases. The first line of input contains an integer $T$ ($1 \le T \le 10000$), indicating the number of test cases. For each test case:

The first line contains two integers $n$ and $k$ ($1 \le n, k \le 10^{10}$).

Output
For each test case, output a single line containing $S_k(n) \bmod 2^{64}$.

Example
Input
5
1 3
2 3
3 3
10 3
100 3
Output
1
5
9
73
2302
Information
There are 5 Input files.

Input #1: $1 \le n \le 10000$, TL = 1s.
Input #2: $1 \le T \le 300,\ 1 \le n \le 10^7$, TL = 5s.
Input #3: $1 \le T \le 75,\ 1 \le n \le 10^{8}$, TL = 5s.
Input #4: $1 \le T \le 15,\ 1 \le n \le 10^{9}$, TL = 5s.
Input #5: $1 \le T \le 5,\ 1 \le n \le 10^{10}$ , TL = 5s.
My C++ solution runs in 5.6 sec. (total time)

Notes
This is general version of DIVCNT1, DIVCNT2 and DIVCNT3. You may want to solve these three problems first.

Sol:

设$f(n) =\sigma_0(n^k)$,可得,用min_25筛搞下即可。

$$\begin{cases} f(1)&=1\\ f(p)&=k+1\\ f(p^c)&=kc+1 \end{cases}$$

这里记录下min_25筛的惯用搞法

对于积性函数$f(x)$,若对于质数$p$,$f(p)$可以表示为简单的多项式和
$f(p) = \sum_{i=1}^{k} c_ip^i$,则可以用min_25筛在$O(\frac{n^\frac{3}{4}}{log})$ 的时间算出$f(x)$的前缀和

• 预处理:筛质数的$f(x)$ 之和

  • $g(n,j)=\sum_{i=2}^xi^k[i\in P \bigcup minp(i) > p_j]$
    若$p_j ^2 > n$,则$g(n,j) = p_j^k *[g(\frac{n}{p_j},j-1) -g(p_{j -1},j-1)]$。
  • 设$p_m$为 $<= \sqrt{n}$ 的最大质数,
    则$g(n,0)$初值为$\sum_{i=2}^{n} i^k$
    所求为 $g(n,m)$.
  • 数论分块
    第一维有$O(\sqrt{n})$种取值,离散化后可以迭代处理

• 积性函数 $f(x)$求前缀和

  • 所求$S(n) = \sum_{i=1}^{n} f(i)$
    $\Rightarrow S(n) = f(1) +\sum_{i=2}^{n}[i \notin P \bigcap minp(i) \le p_m] \bigcup[i\in P]f(i)$

  • 定义 $S(n,j) = \sum_{i=2}^{n}f(x)[minp(i) \ge p_j]$

  $\Rightarrow S(n) = g(n,j) - \sum_{i=1}^{j-1}f(p_i) + \sum_{k=j}^{p_k^2 \le n} \sum_{e = 1}^{p_k^{e+1}}S(\frac{n}{p_k^e},k+1)f(p_k^e) + f(p_k^{e+1})$

  • 然后就可以大力算前缀和了

code:

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;
typedef unsigned long long ull;
const int maxn = 1e6+10;
const int mod = 1e9+7;
const int inv2 = (mod>>1)+1;
ll Sq,n;
bool check[maxn];
ull pri[maxn];
int id[maxn][2];
int tot,cnt;
ull _k;
ull g[maxn],w[maxn];

void init(ll sq){
    tot = 0;
    check[1] = true;
    for(int i = 2;i<=sq;i++){
        if(!check[i]) {
            pri[++tot] = i;
        }
        for(int j = 1;j <= tot;j++){
            if(i * pri[j] > sq) break;
            check[i * pri[j]] = true;
            if(i % pri[j] == 0) break;
        }
    }
}

ull Sol(ll x,ll y){
    if(x<=1 || pri[y] > x) return 0;
    int k = (x <= Sq) ? id[x][0] : id[n/x][1];
    ull ret = (_k + 1) * (g[k] - y + 1);
    for(int i = y;i<=tot && pri[i] * pri[i] <= x;i++){
        ull t1 = pri[i],t2 = pri[i] * pri[i];
        for(int e = 1;t2 <= x;e++){
            ret += Sol(x/t1,i+1) * (_k * e + 1) + (_k * (e+1) + 1);
            t1 = t2;
            t2 *= pri[i];
        }
    }
    return ret;
}

int main(){
    int T;
    scanf("%d",&T);
    init(maxn-1);
    pri[++tot] = maxn;
    while(T--){
        scanf("%lld%llu",&n,&_k);
        Sq = sqrt(n)+0.5;
        ll r;
        cnt = 0;
        for(ll i = 1;i <= n;i = r + 1){
            r = n/(n/i);
            w[++cnt] = n/i;
            g[cnt] = w[cnt] - 1;
            if(w[cnt] <= Sq) id[w[cnt]][0] = cnt;
            else id[r][1] = cnt;
        }
        for(int j = 1;j<=tot;j++){
            for(int i = 1;i<=cnt && (ull)pri[j] * pri[j] <= w[i] ; ++ i){
                int k = (w[i]/pri[j]<=Sq) ? id[w[i]/pri[j]][0] : id[n/(w[i]/pri[j])][1];
                g[i] -= g[k] - (j - 1);
            }
        }
        ull ans = Sol(n,1) + 1;
        printf("%llu\n",ans);
    }
    return 0;
}