bzoj3601 一个人的数论 (拉格朗日插值求系数)

题目链接:bzoj3601 一个人的数论

Description

有一天hjy96想到了一个数论问题:

对于一个非负整数d和一个正整数n,定义fa(n)为所有小于n且与n互质的正整数的d次方之和。如\(f_3(10) = 1^3+3^3+7^3+ 9^3\)

现给定d,n,求fa(n)的值。输出答案对\(10^9+7\)取模后的结果。

hjy96当然知道怎么做啦!但是他想考考你.....

Input

由于n可能很大,我们给出n的质因数分解式。

\[n=\prod\limits_{i=1}^wp_i^{\alpha_i} \]

第一行有用空格隔开的一个非负整数d,一个正整数w。

接下来w行,每行有两个用空格隔开的正整数\(p_i,\alpha_i\)。(保证\(p_i\)为素数且互不相同)

Output

一行,一个非负整数表示\(f_d(n)\)\(10^9+7\)取模后的结果。

Sample Input

3 2
2 1
5 1

Sample Output

1100

数据范围与约定

对于所有数据 \(1\le d\le 100,1 \le w \le 1000,2 \le p_i \le 10^9,1 \le \alpha_i \le 10^9.\)

Solution

由题可知

\[ans=\sum\limits_{i=1}^n[gcd(i,n)=1]i^d \]

推一下

\[ans=\sum\limits_{i=1}^n[gcd(i,n)=1]i^d\\ =\sum\limits_{i=1}^n\sum\limits_{t|x且t|n}\mu(t)i^d\\ =\sum\limits_{t|n}\mu(t)\sum\limits_{i=1}^{\frac{n}{t}}(i\times t)^d\\ =\sum\limits_{t|n}\mu(t)t^d\sum\limits_{i=1}^{\frac{n}{t}}i^d\\ \]

\(\sum\limits_{i=1}^{\frac{n}{t}}i^d可以表示为一个关于\frac{n}{t}的d+1次多项式\sum\limits_{i=0}^{d+1}f_i(\frac{n}{t})^i\)
\(f_i\)可以用拉格朗日插值或高斯消元求出(我用的是拉格朗日)​

\[\therefore ans=\sum\limits_{t|n}\mu(t)t^d\sum\limits_{i=0}^{d+1}f_i(\frac{n}{t})^i\\ =\sum\limits_{t|n}\mu(t)\sum\limits_{i=0}^{d+1}f_in^it^{d-i}\\ =\sum\limits_{i=0}^{d+1}\sum\limits_{t|n}\mu(t)f_in^it^{d-i}\\ =\sum\limits_{i=0}^{d+1}f_in^i\sum\limits_{t|n}\mu(t)t^{d-i}\\ \]

容易发现\(\sum\limits_{t|n}\mu(t)t^{d-i}\)是一个积性函数

\(F(n)=\sum\limits_{t|n}\mu(t)t^{d-i}\)

则有\(F(n)=\prod\limits_{i=1}^{w}F(p_i^{\alpha_i})\)

\(F(p_i^{\alpha_i})\)显然等于\(1-p_i^{d-i}\)

\[\therefore ans= \sum\limits_{i=0}^{d+1}f_in^i\prod\limits_{i=1}^{w}(1-p_i^{d-i})\\ \]

Code

#include<bits/stdc++.h>
const int N=1005;
const int mod=1e9+7;
int Pow(int x,int f=mod-2){
	int re=1;
	while(f){
		if(f&1)re=1ll*re*x%mod;
		f>>=1;
		x=1ll*x*x%mod;
	}
	return re;
}
namespace Lagrange{
	int inv[N],f_[N],f[N],y[N],tmp[N];
	void Solve(int n){
		for(int i=1;i<=n;i++){
			int div=1,lst=0;
			for(int j=1;j<=n;j++)
				if(i!=j)div=1ll*div*(mod+i-j)%mod;
			div=1ll*y[i]*Pow(div)%mod;
			for(int j=0;j<n;j++){
				tmp[j]=1ll*(lst+mod-f_[j])*inv[i]%mod;
				f[j]=(f[j]+(1ll*div*tmp[j]%mod))%mod;
				lst=tmp[j];
			}
		}
	}
	void Pre(int d){
		f_[0]=1;
		for(int i=1;i<=d+2;std::swap(f_,tmp),i++){
			tmp[0]=0,inv[i]=Pow(i);
			for(int j=1;j<=i;j++)tmp[j]=f_[j-1];
			for(int j=0;j<=i;j++)tmp[j]=(tmp[j]+(1ll*f_[j]*(mod-i)%mod))%mod;
		}
		for(int i=1;i<=d+2;i++)y[i]=(y[i-1]+Pow(i,d))%mod;
	}
}
int H(int t,int p){
	if(t<0)return (1+mod-Pow(Pow(p,-t)))%mod;
	return (1+mod-Pow(p,t))%mod;
}
int d,w,a[N],p[N],ans=0;
signed main(){
	using namespace Lagrange;
	scanf("%d%d",&d,&w);
	for(int i=1;i<=w;i++)scanf("%d%d",p+i,a+i);
	Pre(d);
	Solve(d+2);
	int pown=1;
	for(int i=1;i<=d+1;i++){
		int tmp=1;
		for(int j=1;j<=w;j++){
			pown=1ll*pown*Pow(p[j],a[j])%mod;
			tmp=1ll*tmp*H(d-i,p[j])%mod;
		}
		ans=(ans+(1ll*f[i]*pown%mod*tmp%mod))%mod;
	}
	printf("%d",ans);
	return 0;
}
posted @ 2019-05-09 15:56  The_KOG  阅读(403)  评论(0编辑  收藏  举报