hdu5275 (插值法)
2015-06-30 00:26
441 查看
首先要明确,所谓插值,就是多项式的点值表达式转系数表达式。
对于n个x不相同的点(xi,yi),可以唯一确定一个阶为n的多项式y=c0+c1∗x+c2∗x2+c3∗x3+...+cn−1∗xn−1 y=c_0+c_1*x+c_2*x^2+c_3*x^3+...+c_{n-1}*x^{n-1}
一般的,有两种常用的插值方法:
1.拉格朗日插值法
直接构造多项式y=∑k=0n−1yk∏j≠k(xj−xk)∏j≠k(x−xj)y=\sum^{n-1}_{k=0}\frac{y_k}{\prod_{j\ne k}(x_j-x_k)}\prod_{j\ne k}(x-x_j)
简单粗暴,证明只需要把x0,x1,..xn−1x_0,x_1,..x_{n-1}带入,看对应y值是否与给定的相同就可以了
2.牛顿插值
虽然拉格朗日插值法比较简单也便于记忆,但对于一些要对子串进行插值的题目,牛顿插值显然更加方便
令y=c0+c1∗(x−x0)+c2∗(x−x1)∗(x−x0)+...cn−1∗(x−x0)∗...∗(x−xn−2)y=c_0+c_1*(x-x_0)+c_2*(x-x_1)*(x-x_0)+...c_{n-1}*(x-x_0)*...*(x-x_{n-2})
注意对于x0,x1..xn−1x_0,x_1..x_{n-1}进行插值,上式并不包含xn−1x_{n-1},而计算c0,c1...cn−1c_0,c_1...c_{n-1}的方法就是分别把(x0,y0),(x1,y1)...(xn−1,yn−1)(x_0,y_0),(x_1,y_1)...(x_{n-1},y_{n-1})代入,看上去很麻烦,但是引入差商表之后,cic_i的计算显得异常容易。
定义xi的0阶差商为yi,xi的k阶差商f(i,k)=f(i+1,k−1)−f(i,k−1)xi+k−xix_i的0阶差商为y_i,x_i的k阶差商f(i,k)=\frac{f(i+1,k-1)-f(i,k-1)}{x_{i+k}-x_i}
根据定义进行计算可以发现c0,c1...cic_0,c_1...c_i其实就是f(0,0),f(0,1)...f(0,i)f(0,0),f(0,1)...f(0,i)
更进一步的,假如你想对(xl,yl),(xl+1,yl+1)...(xr,yr)(x_l,y_l),(x_{l+1},y_{l+1})...(x_{r},y_{r})
对应插值所得的系数就是f(l,0),f(l,1)...f(l,r−l)f(l,0),f(l,1)...f(l,r-l)
对于上面插值的理解,可以参考我hdu5275的代码
首先是很适合给子串插值的牛顿插值
然后是拉格朗日插值法,对于本题需要一些观察和技巧
对于n个x不相同的点(xi,yi),可以唯一确定一个阶为n的多项式y=c0+c1∗x+c2∗x2+c3∗x3+...+cn−1∗xn−1 y=c_0+c_1*x+c_2*x^2+c_3*x^3+...+c_{n-1}*x^{n-1}
一般的,有两种常用的插值方法:
1.拉格朗日插值法
直接构造多项式y=∑k=0n−1yk∏j≠k(xj−xk)∏j≠k(x−xj)y=\sum^{n-1}_{k=0}\frac{y_k}{\prod_{j\ne k}(x_j-x_k)}\prod_{j\ne k}(x-x_j)
简单粗暴,证明只需要把x0,x1,..xn−1x_0,x_1,..x_{n-1}带入,看对应y值是否与给定的相同就可以了
2.牛顿插值
虽然拉格朗日插值法比较简单也便于记忆,但对于一些要对子串进行插值的题目,牛顿插值显然更加方便
令y=c0+c1∗(x−x0)+c2∗(x−x1)∗(x−x0)+...cn−1∗(x−x0)∗...∗(x−xn−2)y=c_0+c_1*(x-x_0)+c_2*(x-x_1)*(x-x_0)+...c_{n-1}*(x-x_0)*...*(x-x_{n-2})
注意对于x0,x1..xn−1x_0,x_1..x_{n-1}进行插值,上式并不包含xn−1x_{n-1},而计算c0,c1...cn−1c_0,c_1...c_{n-1}的方法就是分别把(x0,y0),(x1,y1)...(xn−1,yn−1)(x_0,y_0),(x_1,y_1)...(x_{n-1},y_{n-1})代入,看上去很麻烦,但是引入差商表之后,cic_i的计算显得异常容易。
定义xi的0阶差商为yi,xi的k阶差商f(i,k)=f(i+1,k−1)−f(i,k−1)xi+k−xix_i的0阶差商为y_i,x_i的k阶差商f(i,k)=\frac{f(i+1,k-1)-f(i,k-1)}{x_{i+k}-x_i}
根据定义进行计算可以发现c0,c1...cic_0,c_1...c_i其实就是f(0,0),f(0,1)...f(0,i)f(0,0),f(0,1)...f(0,i)
更进一步的,假如你想对(xl,yl),(xl+1,yl+1)...(xr,yr)(x_l,y_l),(x_{l+1},y_{l+1})...(x_{r},y_{r})
对应插值所得的系数就是f(l,0),f(l,1)...f(l,r−l)f(l,0),f(l,1)...f(l,r-l)
对于上面插值的理解,可以参考我hdu5275的代码
首先是很适合给子串插值的牛顿插值
#include<bits/stdc++.h> using namespace std; typedef long long Int; const int M=1e9+7; int rev[250002]; int dp[3002][3002]; int x[3002],y[3002]; inline int getrev(int x){return x<0?-rev[-x]:rev[x];} int main() { rev[1]=1; for(int i=2;i<=250000;i++)rev[i]=(M-M/i)*(Int)rev[M%i]%M; int n; while(scanf("%d",&n)!=EOF) { for(int i=1;i<=n;i++)scanf("%d%d",x+i,y+i); for(int i=1;i<=n;i++)dp[0][i]=y[i]; for(int i=1;i<=n;i++) for(int j=1;j+i<=n;j++) dp[i][j]=((dp[i-1][j+1]-dp[i-1][j])*(Int)getrev(x[i+j]-x[j])%M+M)%M; int m;scanf("%d",&m); while(m--) { int l,r,q;scanf("%d%d%d",&l,&r,&q); int ans=0,cur=1; for(int i=0;i<=r-l;i++) { ans+=dp[i][l]*(Int)cur%M; if(ans>=M)ans-=M; cur=(cur*(Int)(q-x[l+i])%M+M)%M; }printf("%d\n",ans); } } }
然后是拉格朗日插值法,对于本题需要一些观察和技巧
#include<bits/stdc++.h> using namespace std; const int Maxn=3002,M=1e9+7; typedef long long Int; int x[Maxn],y[Maxn]; int n,m; int dp[Maxn][Maxn]; int rev[250002],frev[250002]; int powmod(int x,int y) { int t=x,ret=1; while(y) { if(y&1)ret=ret*(Int)t%M; y>>=1; t=t*(Int)t%M; } return ret; } inline int getrev(int x) { if(x<0)return frev[-x]; return rev[x]; } int main() { for(int i=1;i<=250000;i++)rev[i]=powmod(i,M-2),frev[i]=powmod(M-i,M-2); while(scanf("%d",&n)!=EOF) { for(int i=1;i<=n;i++)scanf("%d%d",x+i,y+i); scanf("%d",&m); for(int i=1;i<=n;i++) { dp[i][i]=1; for(int j=i-1;j>=1;j--)dp[i][j]=dp[i][j+1]*(Int)getrev(x[i]-x[j])%M; for(int j=i+1;j<=n;j++)dp[i][j]=dp[i][j-1]*(Int)getrev(x[i]-x[j])%M; } while(m--) { int l,r,q; scanf("%d%d%d",&l,&r,&q); int cs=-1,ans=0; int tot=1; for(int i=l;i<=r;i++) { if(x[i]!=q)tot=((q-x[i])*(Int)tot%M+M)%M; else{cs=i;} } if(cs!=-1) { ans=y[cs]*(Int)dp[cs][l]%M*dp[cs][r]%M*tot%M; } else { for(int k=l;k<=r;k++) { ans+=y[k]*(Int)dp[k][l]%M*dp[k][r]%M*tot%M*getrev(q-x[k])%M; if(ans>=M)ans-=M; } }printf("%d\n",ans); } } }
相关文章推荐
- SQLServer 订阅过期解决方法
- 【NGUI】NGUI血条制作,当人物移出屏幕后不显示血条,优化代码
- c3p0的DriverManagerDataSource源代码看单例
- UVa-512(经典栈)
- ImageLoader 线程池 单例 使用模板
- LeetCode Simplify Path
- 2015年6月29日 课设日志
- win7-时间更新
- 2015062910 - 一周搞定Effective java
- android searchView加ListView共同使用
- mysql定时器
- 定义一个功能,记录程序运行的次数。满足5次后,给出提示,试用次数已到,请注册
- 输入框只能输入数字 以及其他限制的总结
- AFNetworking获取服务端plist文件,XML文件
- POJ 3592 强连通缩点+spfa最长路
- 公安项目项目感想
- 注解开发spring-aop 入门
- 2015062909 - 黑鹰坠落
- .Net 自定义事件
- .Net 自定义事件