hdu6607(min25筛+杜教筛+拉格朗日插值)

|

题目链接

http://acm.hdu.edu.cn/showproblem.php?pid=6607

题意

给定 $n$ 和 $k$

$n\le10^{10},k\le100$

题解

感觉就是板子大合集。。

但是公式不难推。。

显然对 $\lfloor\frac{n}{d} \rfloor$ 分块即可，这样我们需要求出 $\sum i^2\varphi(i)$ 和 $\sum_{d\in prime}d^{k+1}$

前面的直接杜教筛

后面的上 min25 筛即可，只需要用到前半部分的预处理。。


然后问题就是如何快速求 $k$ 次幂和了。。

这个需要拉格朗日插值。。

如果知道了多项式的点值 $(x_i,y_i)$ ，可以通过以下方法构造

对 $k$ 幂和，我们知道其为 $k+1$ 次多项式，所以我们需要 $k+2$ 个点来求解 ，最简单的就是 $1,2..,k+2$ ，然后对每个要求的 $x$ ，对分子分母用前缀和预处理，然后直接求解即可。。复杂度为 $O(k\sqrt n)$




代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
/**
 *　　　　　　  　　┏┓　　 　┏┓
 * 　　　　　  　　┏┛┗━━━━━━━┛┗━━━┓
 * 　　　　　  　　┃　　　　　　　┃ 　
 * 　　　　　  　　┃　　　━　　 　┃
 * 　　　　　  　　┃　＞　　　＜　┃
 * 　　　　　  　　┃　　　　　　　┃
 * 　　　　　  　　┃...　⌒　... 　┃
 * 　　　　　　  　┃              ┃
 * 　　　　　　  　┗━┓          ┏━┛
 * 　　　　　　　　　┃          ┃　Code is far away from bug with the animal protecting　　　　　　　　　　
 * 　　　　　　　　　┃          ┃   神兽保佑,代码无bug
 * 　　　　　　　　　┃          ┃　　　　　　　　　　　
 * 　　　　　　　　　┃          ┃  　　　　　　
 * 　　　　　　　　　┃          ┃
 * 　　　　　　　　　┃          ┃　　　　　　　　　　　
 * 　　　　　　　　　┃          ┗━━━┓
 * 　　　　　　　　　┃              ┣┓
 * 　　　　　　　　　┃              ┏┛
 * 　　　　　　　　　┗┓┓┏━━━━━━━━┳┓┏┛
 * 　　　　　　　　　　┃┫┫       ┃┫┫
 * 　　　　　　　　　　┗┻┛       ┗┻┛
 */ 
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<iostream>
#include<queue>
#include<map>
#include<stack>
#include<set>
#include<bitset>
#include<cmath>
#include<unordered_map>
#define inc(i,l,r) for(int i=l;i<=r;i++)
#define dec(i,l,r) for(int i=l;i>=r;i--)
#define link(x) for(edge *j=h[x];j;j=j->next)
#define mem(a) memset(a,0,sizeof(a))
#define ll long long
#define eps 1e-8
#define succ(x) (1<<x)
#define lowbit(x) (x&(-x))
#define sqr(x) ((x)*(x))
#define mid (x+y>>1)
#define NM 200005
#define nm 10000005
#define pi 3.1415926535897931
using namespace std;
const int inf=1e9+7;
ll read(){
    ll x=0,f=1;char ch=getchar();
    while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
    while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
    return f*x;
}
 
 


inline void reduce(ll&x){x+=x>>63&inf;}
ll n,ans,f[NM],w[NM],pre[NM],prime[nm],phi[nm];
ll inv[105],invp[105];
int _k,m,tot,cnt;
bool v[nm];
ll qpow(ll x,ll t){
    ll s=1;
    for(;t;t>>=1,x=x*x%inf)if(t&1)s=s*x%inf;
    return s;
}
inline int id(ll x){return x<=m/2?m-x+1:n/x;}

void init(){
    cnt=1e7;phi[1]=1;
    inc(i,2,cnt){
	if(!v[i])prime[++tot]=i,phi[i]=i-1;
	inc(j,1,tot){
	    if(i*prime[j]>cnt)break;
	    v[i*prime[j]]++;
	    if(i%prime[j]==0){
		phi[i*prime[j]]=phi[i]*prime[j];
		break;
	    }
	    phi[i*prime[j]]=phi[i]*phi[prime[j]];
	}
    }
    inc(i,1,cnt)phi[i]=(phi[i-1]+phi[i]*i%inf*i%inf)%inf;
    inv[1]=invp[0]=invp[1]=1;
    inc(i,2,103)inv[i]=inv[inf%i]*(inf-inf/i)%inf,invp[i]=invp[i-1]*inv[i]%inf;
}
void init(ll n){
    m=sqrt(n);
    for(tot=1;prime[tot]<=m;tot++)
	reduce(pre[tot]=pre[tot-1]+qpow(prime[tot],_k)-inf);
    tot--;
    inc(i,1,m)w[i]=n/i;
    while(w[m]>1)w[m+1]=w[m]-1,m++;
}

unordered_map<ll,int>mp;

ll p2(ll n){return n*(n+1)%inf*(2*n+1)%inf*(inf+1)/6%inf;}
ll PH(ll n){
    if(n<=cnt)return phi[n];
    if(mp.count(n))return mp[n];
    ll ans=(n+1)%inf*(n%inf)%inf*(inf+1)/2%inf;
    ans=ans*ans%inf;
    for(ll x=2,y;x<=n;x=y+1){
	y=n/(n/x);
	reduce(ans-=PH(n/x)*(p2(y%inf)-p2((x-1)%inf)+inf)%inf);
    }
    return mp[n]=ans;
}


ll sum[105];
ll fun(ll x){
    ll pre[105],suc[105];
    int n=_k+2;
    mem(pre);mem(suc);
    pre[0]=1;
    inc(i,1,n)pre[i]=pre[i-1]*(x-i)%inf;
    suc[n+1]=1;
    dec(i,n,1)suc[i]=suc[i+1]*(x-i)%inf;
    ll ans=0;
    inc(i,1,n)
	if((n-i)&1)ans+=inf-pre[i-1]*suc[i+1]%inf*invp[i-1]%inf*invp[n-i]%inf*sum[i]%inf,ans%=inf;
	else ans+=inf+pre[i-1]*suc[i+1]%inf*invp[i-1]%inf*invp[n-i]%inf*sum[i]%inf,ans%=inf;
    return ans;
}


int main(){
    init();
    int _=read();while(_--){
	mp.clear();
	n=read();_k=read()+1;
	init(n);ans=0;
	inc(i,1,_k+2)sum[i]=(sum[i-1]+qpow(i,_k))%inf;
	inc(i,1,m)f[i]=fun(w[i]%inf)-1;
	//inc(i,1,m)printf("%lld ",w[i]);putchar('\n');
	//inc(i,1,m)printf("%lld ",f[i]);putchar('\n');
	inc(j,1,tot){
	    ll t=qpow(prime[j],_k);
	    inc(i,1,m)if(prime[j]*prime[j]<=w[i]){
		reduce(f[i]-=t*(f[id(w[i]/prime[j])]-pre[j-1]+inf)%inf);
	    }else break;
	}
	//inc(i,1,m)printf("%lld ",f[i]);putchar('\n');
	f[m+1]=0;
	for(ll x=1,y;x<=n;x=y+1){
	    y=n/(n/x);
	    ans+=PH(n/x)*(f[id(y)]-f[id(x-1)]+inf)%inf;
	    ans%=inf;
	}
	printf("%lld\n",ans);
    }
    return 0;
}