实数数列
原题：（权限限制没有链接（虽然我想应该是有的））
【问题描述】
一个实数数列共有n项，已知 a[i]=(a[i-1]-a[i+1])/2+d，(1<i<n) (n<＝60)
【输入文件】
输入第一行为n，m，d。第二行为，a[1],a[n]两个整数。
【输出文件】
输出a[m]，结果保留6位小数。
【样例输入】
3 2 2
1 －3
【样例输出】
4.000000
个人解法：
论矩阵乘法的不足，今天算是明白了。
曾经以为遇到数列就套矩阵乘法，不过还好做到了这个题。
显然一个通项公式是很容易搞出来的——
A[i+2]=A[i]-2*A[i+1]+2*d
如斐波那契数列一般，很容易发现有
接下来做一个矩阵快速幂。
之后我们就可以推出A[2]（因为A[2]可以列出一个只有A[2]和其他常数的方程）。
之后一个一个退到过去就可以了。
代码如下：
type
func=array[1..3,1..3] of extended;//矩阵
const
linit_func:func=((1,0,0),
(0,1,0),
(0,0,1));//初始矩阵（乘以任何一个矩阵A的答案仍然是矩阵A）
var
i,n,m:longint;
a1,an,ans,d:double;
A:array[1..62] of double;
handle_func:func;//操作矩阵——也就是上图中推导出来的那玩意
function multiply(a,b:func):func;//矩阵乘法
var
vi,vj,vk:longint;
begin
fillchar(multiply,sizeof(multiply),0);
for vi:=1 to 3 do
for vk:=1 to 3 do
for vj:=1 to 3 do
multiply[vi,vk]:=multiply[vi,vk]+a[vi,vj]*b[vj,vk];
end;
function power(a:func;b:longint):func;//矩阵快速幂
var
s,w:func;
begin
s:=linit_func;
w:=a;
while b<>0 do
begin
if b and 1<>0 then s:=multiply(s,w);
w:=multiply(w,w);
b:=b shr 1;
end;
exit(s);
end;
begin
read(n,m,d);//读入
read(a1,an);
if m=1 then//一大堆特判
begin
writeln(a1:0:6);
halt;
end;
if m=n then
begin
writeln(an:0:6);
halt;
end;
if n=3 then
begin
ans:=(a1-an)/2+d;
writeln(ans:0:6);
halt;
end;
handle_func[1,2]:=1;
handle_func[2,1]:=1;
handle_func[2,2]:=-2;
handle_func[3,2]:=2;
handle_func[3,3]:=1;
handle_func:=power(handle_func,n-2);//矩阵乘法快速幂
A[1]:=a1;
A[n]:=an;
A[2]:=(an-a1*Handle_func[1,2]-d*Handle_func[3,2])/Handle_func[2,2];//推导出A[2]
for i:=3 to m do
A[i]:=A[i-2]-2*A[i-1]+2*d;//递推
writeln(A[m]:0:6);//输出
end.
然而事实上似乎不太给力的样子。
考虑这个矩阵，最多做60次幂运算。之后考虑矩阵每一个格子的权值会相乘相加，事实上int64数储存不了这么大的数的。
显然——开extended或者double会直接变为科学计数法——精度就完全不能保证了。
怎么办？
考虑到我们 是直接以A[i]来推导的矩阵。是否可以换一个思考角度？
显然一个一个写过来，我们是可以得到每一项的系数——例如
设
A[i]=P[i]*A[2]+Q[i]*d+R[i]*A[1]
那么一个一个写过去就会发现
P[i]=P[i-2]-2*P[i-1]
Q[i]=Q[i-2]-2*Q[i-1]+2
R[i]=R[i-2]-2*R[i-1]
初始化就是
R[1]=1
P[2]=2
然后我们发现这个玩意——longint都不会超。
显然就比某矩阵乘法要强了。
（虽然这也是可以矩阵乘法的——但是这么简单的一个式子就不太愿意写矩阵乘法了）
还有一个注意点：
由于我们一个一个递推过去——用了除法，所以误差可能会有点大。
那么怎么解决呢？
我们知道——
A[i]=P[i]*A[2]+Q[i]*d+R[i]*A[1]
=P[i-1]*A[3]+Q[i-1]*d+R[i-1]*A[2]
=P[i-2]*A[4]+Q[i-2]*d+R[i-2]*A[3]
=……
=P[i-2+k]*A[k]+Q[i-2+k]*d+R[i-2+k]*A[k-1]
所以就有
A[k]=(A[n]-Q[n-k+2]*d+R[n-k+2]*A[k-1])/P[n-k+2]
由于P递减，所以最后误差会小很多
代码如下：
var
n,m,i:longint;
d:double;
A,P,Q,R:array[1..62] of double;
begin
read(n,m,d);
read(A[1],A[n]);//读入
R[1]:=1;
P[2]:=1;
for i:=3 to n do//递推
begin
P[i]:=P[i-2]-2*P[i-1];
Q[i]:=Q[i-2]-2*Q[i-1]+2;
R[i]:=R[i-2]-2*R[i-1];
end;
for i:=2 to m do
A[i]:=(A[n]-Q[n-i+2]*d-R[n-i+2]*A[i-1])/P[n-i+2];//计算
writeln(A[m]:0:6);//输出
end.
也算吃一堑长一智吧。