SPOJ：Counting Divisors (general)（min_25筛）

Let \sigma_0(n)σ​0​​(n) be the number of positive divisors of nn.

For example, \sigma_0(1) = 1σ​0​​(1)=1, \sigma_0(2) = 2σ​0​​(2)=2 and \sigma_0(6) = 4σ​0​​(6)=4.

LetS_k(n) = \sum _{i=1}^n \sigma_0(i^k).S​k​​(n)=​i=1​∑​n​​σ​0​​(i​k​​).

Given nn and kk, find S_k(n) \bmod 2^{64}S​k​​(n)mod2​64​​.

Input

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

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

Output

For each test case, output a single line containing S_k(n) \bmod 2^{64}S​k​​(n)mod2​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 100001≤n≤10000, TL = 1s.
• Input #2: 1 \le T \le 300,\ 1 \le n \le 10^71≤T≤300, 1≤n≤10​7​​, TL = 5s.
• Input #3: 1 \le T \le 75,\ 1 \le n \le 10^{8}1≤T≤75, 1≤n≤10​8​​, TL = 5s.
• Input #4: 1 \le T \le 15,\ 1 \le n \le 10^{9}1≤T≤15, 1≤n≤10​9​​, TL = 5s.
• Input #5: 1 \le T \le 5,\ 1 \le n \le 10^{10}1≤T≤5, 1≤n≤10​10​​ , TL = 5s.

My C++ solution runs in 5.6 sec. (total time)

题意：给出N，K，定义f(i)=i^k的因子数量，求1<=i<=N的f(i)之和。

思路：f是一个积性函数，且f(p^c) = kc+1，那么可以用min_25筛来解决，写了两个版本，递归版和非递归版。

关于min_25筛：

1、后者较前者速度慢一点，由于第i个素数的解通过第i+1个素数递推出来，因此可以用类似背包那样的滚动数组优化。

2、min_25筛用来求一类积性函数的和（不局限于这个用法），套路一般是先筛出1~N内的素数的F(i)值，F(i)可以是积性函数或者单纯是F(i)=i，还可以筛出1~N的质因子个数，这一步可以当模板用。

3、递归版f(x,i)表示2~X由大于等于i的质数组成的数（意味着包括合数）的函数和。

4、非递归版f(x,i)表示2~X由大于等于i的质数组成的数（意味着包括合数）的函数和 加上小于第i个质数的质数函数和。

递归版

# include <bits/stdc++.h>
using namespace std;
const int N = 1e5 + 1000;
typedef long long ll;
typedef unsigned long long ull;
ll n, P, id1[N],id2[N],w[2*N],m, prime[N];
ull h[N*2], k;
bool is[N];
void init(){
    for (ll i = 2; i <= 1e5; ++i) {
        if (!is[i]) prime[++prime[0]] = i;
        for (ll j = 1; j <= prime[0] && prime[j] * i <= 100000; ++j) {
            is[prime[j] * i] = 1;
            if (i % prime[j] == 0) break;
        }
    }
}
int getid(ll x) {
    if (x <= P) return id1[x]; else return id2[n / x];
}
void calc_gh() {
    m = 0;
    for (ll l=1; l<=n;){//根号分块，因为递推式都是下取整，真正用到的区间不超过2sqrt(n)个。
        ll v = n / l, r = n / v;
        if (v <= P) id1[v] = ++m; else id2[r] = ++m;//离散化区间
        w[m] = v;//区间长度
        h[m] = v - 1;//初始化假设全部都是素数
        l = r + 1;
    }
    for (ll i = 1; i <= prime[0]; ++i) {
        ll z = prime[i];
        for (ll j = 1; j <= m && z * z <= w[j]; ++j) {
            int op = getid(w[j]/z);
            h[j] = h[j] - (h[op] - (i - 1));
        }
    }
    //上面的h数组处理出区间内的素数个数。
    for (int i = 1; i <= m; ++i)
        h[i] = h[i] * (ull)(k+1);
}
ull ask(ll x,ll i) {
    if (x <= 1 || prime[i] > x) return 0;
    ull ans = h[getid(x)] - (k+1)*(i-1);//先算大于等于第i个质数的质数答案之和。
    for (ll j = i; j <= prime[0] && prime[j] * prime[j] <= x; ++j) {
        ll r = prime[j];
        for (ll e = 1; r * prime[j] <= x; e++, r *= prime[j])
            ans += ask(x / r, j + 1) * (k*e + 1) + (k*(e+1)+1);//后面加上仅由质数i组成的合数
    }
    return ans;
}
int main(){
    init();
    int T;
    for(scanf("%d",&T);T;--T){
        scanf("%lld%llu",&n,&k);
        P = sqrt(n);
        calc_gh();
        printf("%llu\n",ask(n,1) + 1);
    }
    return 0;
}

 

非递归版：

# include <bits/stdc++.h>
using namespace std;
const int N = 1e5 + 1000;
typedef long long ll;
typedef unsigned long long ull;
ll n, P, id1[N],id2[N],w[2*N],m, prime[N];
ull h[N*2], k;
bool is[N];
void init(){
    for (ll i = 2; i <= 1e5; ++i) {
        if (!is[i]) prime[++prime[0]] = i;
        for (ll j = 1; j <= prime[0] && prime[j] * i <= 100000; ++j) {
            is[prime[j] * i] = 1;
            if (i % prime[j] == 0) break;
        }
    }
}
int getid(ll x) {
    if (x <= P) return id1[x]; else return id2[n / x];
}
void calc_gh() {
    m = 0;
    for (ll l=1; l<=n;){
        ll v = n / l, r = n / v;
        if (v <= P) id1[v] = ++m; else id2[r] = ++m;
        w[m] = v;
        h[m] = v - 1;
        l = r + 1;
    }
    for (ll i = 1; i <= prime[0]; ++i) {
        ll z = prime[i];
        for (ll j = 1; j <= m && z * z <= w[j]; ++j) {
            int op = getid(w[j]/z);
            h[j] = h[j] - (h[op] - (i - 1));
        }
    }
    //上面的h数组处理出区间内的素数个数。
    for (int i = 1; i <= m; ++i)
        h[i] = h[i] * (ull)(k+1);
}
ull calc(ll n) {
    if (n <= 1) return 1;
    for (int j = prime[0]; j; j--) {
        if (prime[j] * prime[j] > n) continue;
        for (int i = 1; i <= m; ++i) {
            if (prime[j] * prime[j] > w[i]) break;
            ll l = prime[j];
            for (int e = 1; l * prime[j] <= w[i]; ++e,l*=prime[j]) {
                h[i] += (h[getid(w[i]/l)]-(k+1)*j) * (k*e + 1) + (k*(e+1)+1);
            }
        }
    }
    return h[1] + 1;
}
int main(){
    init();
    int T;
    for(scanf("%d",&T);T;--T){
        scanf("%lld%llu",&n,&k);
        P = sqrt(n);
        calc_gh();
        printf("%llu\n",calc(n));
    }
    return 0;
}