SPOJ - DIVCNTK

题目链接:divcntk

Let σ0(n) σ 0 ( n ) be the number of positive divisors of n n .

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

Let

Sk(n)=i=1nσ0(ik). S k ( n ) = ∑ i = 1 n σ 0 ( i k ) .

Given n n and k, find Sk(n)mod264 S k ( n ) mod 2 64 .

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

The first line contains two integers n n and k ( 1n,k1010 1 ≤ n , k ≤ 10 10 ).

Output
For each test case, output a single line containing Sk(n)mod264 S k ( n ) mod 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: 1n10000 1 ≤ n ≤ 10000 , TL = 1s.
Input #2: 1T300, 1n107 1 ≤ T ≤ 300 ,   1 ≤ n ≤ 10 7 , TL = 5s.
Input #3: 1T75, 1n108 1 ≤ T ≤ 75 ,   1 ≤ n ≤ 10 8 , TL = 5s.
Input #4: 1T15, 1n109 1 ≤ T ≤ 15 ,   1 ≤ n ≤ 10 9 , TL = 5s.
Input #5: 1T5, 1n1010 1 ≤ T ≤ 5 ,   1 ≤ n ≤ 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)=σ0(nk) f ( n ) = σ 0 ( n k ) ,可得,用min_25筛搞下即可。

f(1)f(p)f(pc)=1=k+1=kc+1 { f ( 1 ) = 1 f ( p ) = k + 1 f ( p c ) = k c + 1

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

对于积性函数 f(x) f ( x ) ,若对于质数 p p ,f(p)可以表示为简单的多项式和
f(p)=ki=1cipi f ( p ) = ∑ i = 1 k c i p i ,则可以用min_25筛在 O(n34log) O ( n 3 4 l o g ) 的时间算出 f(x) f ( x ) 的前缀和

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

    • g(n,j)=xi=2ik[iPminp(i)>pj] g ( n , j ) = ∑ i = 2 x i k [ i ∈ P ⋃ m i n p ( i ) > p j ]
      p2j>n p j 2 > n ,则 g(n,j)=pkj[g(npj,j1)g(pj1,j1)] g ( n , j ) = p j k ∗ [ g ( n p j , j − 1 ) − g ( p j − 1 , j − 1 ) ]
    • pm p m <=n <= n 的最大质数,
      g(n,0) g ( n , 0 ) 初值为 ni=2ik ∑ i = 2 n i k
      所求为 g(n,m) g ( n , m ) .
    • 数论分块
      第一维有 O(n) O ( n ) 种取值,离散化后可以迭代处理
  • 积性函数 f(x) f ( x ) 求前缀和

    • 所求 S(n)=ni=1f(i) S ( n ) = ∑ i = 1 n f ( i )
      S(n)=f(1)+ni=2[iPminp(i)pm][iP]f(i) ⇒ S ( n ) = f ( 1 ) + ∑ i = 2 n [ i ∉ P ⋂ m i n p ( i ) ≤ p m ] ⋃ [ i ∈ P ] f ( i )

    • 定义 S(n,j)=ni=2f(x)[minp(i)pj] S ( n , j ) = ∑ i = 2 n f ( x ) [ m i n p ( i ) ≥ p j ]

    S(n)=g(n,j)j1i=1f(pi)+p2knk=jpe+1ke=1S(npek,k+1)f(pek)+f(pe+1k) ⇒ S ( n ) = g ( n , j ) − ∑ i = 1 j − 1 f ( p i ) + ∑ k = j p k 2 ≤ n ∑ e = 1 p k e + 1 S ( 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;
}
  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值