loj #6024. XLkxc

最新推荐文章于 2022-05-08 17:29:16 发布

原创最新推荐文章于 2022-05-08 17:29:16 发布 · 397 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#拉格朗日插值 #数论

拉格朗日插值专栏收录该内容

2 篇文章

订阅专栏

本文介绍了一种利用多项式技巧解决复杂求和问题的方法，通过将原问题转化为求多项式的值，实现了高效的计算。文章详细解释了如何通过递推公式构建多项式，并通过拉格朗日插值法求解最终结果。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

背景：

$\text{bzoj}$ 权限题？
没事 $\text{loj}$ 有原题。

题目传送门：

https://loj.ac/problem/6024

题意：

求： $\sum_{i=0}^{n}\sum_{j=1}^{a+i\cdot d}\sum_{l=1}^{j}l^k \mod p$ 。

思路：

做法很妙，实在想不到。
一般人会去化简式子，可是根本化不动，我太菜了。
巨佬们给的做法是：
设 $f_x=\sum_{i=1}^{x}i^k$ 。
我们发现 $f_x$ 是一个 $k + 1$ 次的多项式。
证明就是：
$\begin{aligned}f_x&=\sum_{i=1}^{x}i^k\\ &=\frac{1}{k+1}\sum_{i=1}^{k+1}C_{k+1}^{i}\cdot B_{k+1-i}\cdot (x+1)^i\\ &=\frac{1}{k+1}\sum_{i=1}^{k+1}C_{k+1}^{i}\cdot B_{k+1-i}\cdot \sum_{j=0}^{i}C(_{i}^{j})x^j1^{i-j}\\ &=\frac{1}{k+1}\sum_{i=1}^{k+1}C_{k+1}^{i}\cdot B_{k+1-i}\cdot \sum_{j=0}^{i}C(_{i}^{j})x^j\end{aligned}$
其中 $B$ 是伯努利数。
$j$ 的上限是 $i$ ， $i$ 的上限为 $k + 1$ ，故 $x^j$ 的最大次数为 $k + 1$ ，因此得证。

设 $g_x=\sum_{i=1}^{x}f_i$ 。
我们发现 $g$ 是一个 $k + 2$ 次的多项式。
证明就是：
$\begin{aligned}g_x&=\sum_{i=1}^{x}\sum_{j=1}^{j}j^k\\ &=\sum_{i=1}^{x}\frac{1}{k+1}\sum_{j=1}^{k+1}C_{k+1}^{j}\cdot B_{k+1-j}\cdot (i+1)^j\\ &=\frac{1}{k+1}\sum_{i=1}^{x}\sum_{j=1}^{k+1}C_{k+1}^{j}\cdot B_{k+1-j}\cdot \sum_{l=0}^{j}C_{j}^{l}\cdot i^l\cdot 1^{j-l}\\ &=\frac{1}{k+1}\sum_{i=1}^{x}\sum_{j=1}^{k+1}C_{k+1}^{j}\cdot B_{k+1-j}\cdot \sum_{l=0}^{j}C_{j}^{l}\cdot i^{l}\\ &=\frac{1}{k+1}\sum_{j=1}^{k+1}C_{k+1}^{j}\cdot B_{k+1-j}\cdot \sum_{l=0}^{j}C_{j}^{l}\cdot\sum_{i=1}^{x} i^{l}\end{aligned}$

上面我们已经证明了形如 $\sum_{i=1}^{x} i^{l}$ 是 $l + 1$ 次方的， $l$ 的上限是 $j$ ， $j$ 的上限是 $k + 1$ ，当 $l = j = k + 1$ 时，为 $k + 2$ 次方，得证。

设 $h_x=\sum_{i=1}^{x}g_{a+i\cdot d}$ 。
我们发现 $g$ 是一个 $k + 3$ 次的多项式。
证明就是：
在原来的基础上，再套一层即可。

剩下的按照次数拉格朗日插值即可。
时间复杂度好像有点大： $\Theta(n^3)$ 。
其实可以 $\Theta(n^2)$ （重心拉格朗日插值），当然我很懒，既然模板都没打，为什么要打呢。

代码：

#include<cstdio>
#include<cstring>
#include<algorithm>
#define LL long long
#define LD long double
#define mod 1234567891ll
using namespace std;
	int k;
	LL p,n,d;
	LL f[1010],g[1010],h[1010];
	struct node{LL x,y;} a[1010];
LL times(LL x,LL y)
{
	LD t=(LD)x*(LD)y-(LL)((LD)x*(LD)y/(LD)mod)*(LD)mod;
	LL tt=(LL)t;
	return tt>=0?tt%mod:(tt+mod)%mod;
}
LL dg(LL x,LL y)
{
	if(!y) return 1;
	LL op=dg(x,y>>1);
	if(y&1) return times(times(op,op),x); else return times(op,op);
}
LL inv(LL x)
{
	return dg(x,mod-2);
}
LL lagrange(LL x,int limit)
{
	LL sum=0;
	for(int i=0;i<=limit;i++)
	{
		LL op1=a[i].y,op2=1;
		for(int j=0;j<=limit;j++)
		{
			if(i==j) continue;
			op1=times(op1,(x-a[j].x+mod)%mod);
			op2=times(op2,(a[i].x-a[j].x+mod)%mod);
		}
		sum=(sum+times(op1,inv(op2)))%mod;
	}
	return sum;
}
int main()
{
	int T;
	scanf("%d",&T);
	while(T--)
	{
		scanf("%d %lld %lld %lld",&k,&p,&n,&d);
		for(int i=1;i<=k+3;i++)
			f[i]=(f[i-1]+dg(i,k))%mod;
		for(int i=1;i<=k+3;i++)
			g[i]=(g[i-1]+f[i])%mod;
		for(int i=0;i<=k+3;i++)
			a[i]=(node){(LL)i,g[i]};
		for(int i=0;i<=k+4;i++)
			h[i]=((!i?0:h[i-1])+lagrange((p+times(i,d))%mod,k+3))%mod;
		for(int i=0;i<=k+4;i++)
			a[i]=(node){(LL)i,h[i]};
		printf("%lld\n",lagrange(n,k+4));
	}
}