拉格朗日插值法

拉格朗日

公式
f ( x ) = ∑ i = 1 k + 1 y i ∏ j = 1 , i ≠ j k + 1 x − x j x i − x j f(x)=\sum_{i=1}^{k+1}y_i\prod_{j=1,i\neq j}^{k+1} \frac{x-x_j}{x_i-x_j} f(x)=i=1k+1yij=1,i=jk+1xixjxxj

插值基函数

公式
l i ( x ) = ∏ j = 1 , i ≠ j k + 1 x − x j x i − x j l_i(x)=\prod_{j=1,i\neq j}^{k+1} \frac{x-x_j}{x_i-x_j} li(x)=j=1,i=jk+1xixjxxj

1、无论x取和值,它们的和都为1
2、插值基函数是线性无关的

模板

对k次多项式进行插值,需要k+1个点,dy中有k+1项,为d[1]……d[k+1],x为当前所求函数值,返回值即为f(x)

ll Lagrange(int k,ll dy[],ll x)
{
	if(x<=k)
		return dy[x];
	pref[0]=1;
	for(int j=1;j<=k+1;++j)
		pref[j]=1ll*pref[j-1]*(x-j)%mod;
	
	suff[k+2]=1;
	for(int j=k+1;j>=1;--j)
		suff[j]=1ll*suff[j+1]*(x-j)%mod;
	
	ll ans=0;
	for(int i=1;i<=k+1;++i)
	{
		ll tmp=1ll*dy[i]*pref[i-1]%mod*suff[i+1]%mod*finv[i-1]%mod*finv[k+1-i]%mod;	
		if(k-i+1&1)
			ans=(ans-tmp+mod)%mod;
		else
			ans=(ans+tmp)%mod;	
	}
	return ans;
} 

多项式次数的判断
首先插值成一个数组,然后做多次差分,如果第k次差分为0,那么该式子就是k-1次多项式

示例:对 f ( x ) = ∑ i = 1 x i k f(x)=\sum_{i=1}^x i^k f(x)=i=1xik,作差分

#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int mod=1e9+7;

ll qpow(ll base,ll n,ll mod)
{
	ll ret=1;
	while(n)
	{
		if(n&1)
			ret=ret*base%mod;
		base=base*base%mod;
		n>>=1;
	}
	return ret;
}

ll f[20],d[20];

int main()
{
	int k;
	while(cin>>k)
	{
		f[0]=0;
		for(int i=1;i<=10;++i)
			f[i]=f[i-1]+qpow(i,k,mod);

		for(int i=0;i<=10;++i)
		{
			for(int j=i+1;j<=10;++j)
			{
				d[j]=f[j]-f[j-1];
				cout<<d[j]<<" ";
			}
			for(int j=i+1;j<=10;++j)
				f[j]=d[j];				
			cout<<endl;
		}	
	}  
	return 0;
}
/*结果
4
1 16 81 256 625 1296 2401 4096 6561 10000
15 65 175 369 671 1105 1695 2465 3439
50 110 194 302 434 590 770 974
60 84 108 132 156 180 204
24 24 24 24 24 24
0 0 0 0 0
0 0 0 0
0 0 0
0 0
0
*/

F. The Sum of the k-th Powers

题意:给定n、k,求解 ∑ i = 1 n i k \sum_{i=1}^n i^k i=1nik, ( 1 ≤ n ≤ 1 0 9 , 0 ≤ k ≤ 1 0 6 ) (1\le n\le 10^9,0\le k\le 10^6) (1n109,0k106)
思路:k+1次多项式,需要插入k+2个点

#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int mod=1e9+7;

ll qpow(ll base,ll n,ll mod)
{
	ll ret=1;
	while(n)
	{
		if(n&1)
			ret=ret*base%mod;
		base=base*base%mod;
		n>>=1;
	}
	return ret;
}

const int K=1e6+10;
int pref[K+10],suff[K+10];
int fac[K+10],finv[K+10];

void init()
{
	fac[0]=fac[1]=1;
	for(int i=2;i<=K;++i)
		fac[i]=1ll*fac[i-1]*i%mod;
	finv[K]=qpow(fac[K],mod-2,mod);
	for(int i=K-1;i>=0;--i)
		finv[i]=1ll*finv[i+1]*(i+1)%mod;
}

//函数说明:对k次多项式进行插值,dy中有k+1项,为d[1]……d[k+1],x为当前所求函数值,返回值即为f(x)
ll Lagrange(int k,ll dy[],ll x)
{
	if(x<=k)
		return dy[x];
	pref[0]=1;
	for(int j=1;j<=k+1;++j)
		pref[j]=1ll*pref[j-1]*(x-j)%mod;
	
	suff[k+2]=1;
	for(int j=k+1;j>=1;--j)
		suff[j]=1ll*suff[j+1]*(x-j)%mod;
	
	ll ans=0;
	for(int i=1;i<=k+1;++i)
	{
		ll tmp=1ll*dy[i]*pref[i-1]%mod*suff[i+1]%mod*finv[i-1]%mod*finv[k+1-i]%mod;	
		if(k-i+1&1)
			ans=(ans-tmp+mod)%mod;
		else
			ans=(ans+tmp)%mod;	
	}
	return ans;
} 

int n,k;
ll dy[K+10];

int main()
{
	init();
	scanf("%d%d",&n,&k);
	for(int i=1;i<=k+2;++i)
		dy[i]=qpow(i,k,mod);
	for(int i=1;i<=k+2;++i)
		dy[i]=(dy[i]+dy[i-1])%mod;
	
	printf("%lld\n",Lagrange(k+1,dy,n));
	return 0;
}

「BZOJ3453」「Tyvj1858」 XLkxc

题意
给定 k,a,n,d,p, f ( i ) = 1 k + 2 k + 3 k + ⋯ + i k f(i)=1^k+2^k+3^k+\dots+i^k f(i)=1k+2k+3k++ik g ( x ) = f ( 1 ) + f ( 2 ) + f ( 3 ) + . . . . + f ( x ) g(x)=f(1)+f(2)+f(3)+....+f(x) g(x)=f(1)+f(2)+f(3)+....+f(x)
( g ( a ) + g ( a + d ) + g ( a + 2 d ) + . . . . . . + g ( a + n d ) ) m o d   p (g(a)+g(a+d)+g(a+2d)+......+g(a+nd)) mod\ p (g(a)+g(a+d)+g(a+2d)+......+g(a+nd))mod p
1 ≤ k ≤ 123 , 0 ≤ a , n , d ≤ 123456789 , p = 1234567891 1≤k≤123,0≤a,n,d≤123456789 ,p=1234567891 1k123,0a,n,d123456789,p=1234567891
思路:
g ( x ) = ∑ j = 1 x ∑ l = 1 j l k g(x)=\sum_{j=1}^x\sum_{l=1}^j l^k g(x)=j=1xl=1jlk,g(x)为k+2次多项式

f ( x ) = ∑ i = 0 x g ( a + i d ) f(x)=\sum_{i=0}^xg(a+id) f(x)=i=0xg(a+id),f(x)为k+4次多项式

先暴力求出g(x)的前k+3项,再通过得到的函数,暴力求出f(x)的前k+5项

参考博客:
树形dp
https://blog.csdn.net/weixin_43973966/article/details/85338458

几道例题:
BZOJ4559
HDU4059
HDU6239

  • 0
    点赞
  • 2
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值