拉格朗日插值法

题面：

统计K维空间中与原点切比雪夫距离不超过N的点与原点的切比雪夫距离和

solution：

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstring>
#include<vector>
#include<queue>
#include<set>
#include<map>
#include<stack>
#include<bitset>
#include<ext/pb_ds/priority_queue.hpp>
using namespace std;

const int maxn = 1E6 + 10;
typedef long long LL;
const LL mo = 1000000007;

int T,N,K,tot,pri[maxn*2];
bool not_pri[maxn*2];
LL Ans,F[maxn],fac[maxn],_fac[maxn],pre[maxn],suf[maxn],mi[maxn*2];

int Mul(const LL &x,const LL &y) {return x * y % mo;}
int Add(const int &x,const int &y) {return (x + y) % mo;}
int Dec(const int &x,const int &y) {return (x - y + mo) % mo;}

int ksm(int x,int y)
{
	int ret = 1;
	for (; y; y >>= 1)
	{
		if (y & 1) ret = Mul(ret,x);
		x = Mul(x,x);
	}
	return ret;
}

int main()
{
	#ifdef DMC
		freopen("DMC.txt","r",stdin);
	#else
		freopen("distance.in","r",stdin);
		freopen("distance.out","w",stdout);
	#endif
	
	cin >> T; fac[0] = _fac[0] = 1;
	for (int i = 1; i < maxn; i++)
		fac[i] = Mul(fac[i-1],i),_fac[i] = Mul(_fac[i-1],Dec(0,i));
	fac[maxn-1] = ksm(fac[maxn-1],mo - 2);
	for (int i = maxn - 2; i >= 0; i--)
		fac[i] = Mul(fac[i+1],i + 1);
	_fac[maxn-1] = ksm(_fac[maxn-1],mo - 2);
	for (int i = maxn - 2; i >= 0; i--)
		_fac[i] = Mul(_fac[i+1],Dec(0,i + 1));
	
	while (T--)
	{
		cin >> K >> N; int G = 2*(K+2) + 1; mi[1] = 1;
		for (int i = 2; i <= G; i++)
		{
			if (!not_pri[i])
			{
				pri[++tot] = i;
				mi[i] = ksm(i,K);
			}
			for (int j = 1; j <= tot; j++)
			{
				LL Nex = 1LL*pri[j]*i;
				if (Nex > G) break;
				not_pri[Nex] = 1;
				mi[Nex] = Mul(mi[i],mi[pri[j]]);
				if (i % pri[j] == 0) break;
			}
		}
		for (int i = 2; i <= G; i++) not_pri[i] = 0; tot = 0;
		for (int i = 1; i <= K + 2; i++)
			F[i] = Add(F[i-1],Mul(i,Dec(mi[2*i + 1],mi[2*i - 1])));
		Ans = 0; pre[0] = suf[K+3] = 1;
		for (int i = 1; i <= K + 2; i++)
			pre[i] = Mul(pre[i-1],Dec(N,i));
		for (int i = K + 2; i >= 0; i--)
			suf[i] = Mul(suf[i+1],Dec(N,i));
		for (int i = 1; i <= K + 2; i++)
		{
			LL A = Mul(pre[i - 1],suf[i + 1]);
			LL B = Mul(fac[i - 1],_fac[K + 2 - i]);
			Ans = Add(Ans,Mul(Mul(A,B),F[i]));
		}
		cout << Ans << endl;
	}
	return 0;
}