矩阵乘法

前一阵子FYH和XYY连讲了两周的矩阵乘法，不过单学一个矩乘顶多是在学模拟对吧，所以我们要知道这个算法怎么应用，这里的话主要讲怎样优化递推。

首先，我们要知道怎样操作矩乘，公式如下：

C[i,j]=∑a[i,k]*b[k,j]

前提是a的列数和b的行数一致。

而且我们还要知道，一切矩阵的0次方（即单位矩阵）是这个样子的

众所周知，一般来讲不少递推的递推式是这样的

F[i]=c[1]*f[i-1]+c[2]*f[i-2]+…+c[k]*f[i-k]+q

然后我们可以构造两个矩阵来优化这一递推式

A矩阵：

[c₁,c₂,c₃,…,c_k-1,c_k,q

 1 ,0 ,0,0,…,0,0,0

 0 ,1 ,0,0,…,0,0,0

 ………………

 0 ,0 ,0,0,…,1,0,0

0 ,0 ,0 ,0,…,0,0,1]

B矩阵

[f_i-1

 f_i-2

 f_i-3

 ……

 f₀

1]

然后A*B=

[f_i

 f_i-1

 f_i-2

 ……

 f₁

1]

当然如果想要推导第p项的话，我们可以得到这么一条式子

A^p*B=

[f_p

 f_p-1

 f_p-2

 ^……

 f_p-k

 1]

好了，然后我们来围观一下POJ的一道题3744

题目点击此处

我们可以很轻易地知道踩在某一格子上的几率是这样的

F[i]=p*f[i-1]+(1-p)*f[i-2]

然后再进行分段计算即可，分段的依据是地雷所在位置。

对于任何地雷所在位置i，跨过它的几率是f[i-1]*(1-p)

把这些几率相乘就得到结论

注意：

1.     f[1]=1,f[0]=0（任何一段都是一样的）

2.     把地雷按其所在位置排序后要先进行特判，如果起点就是地雷，挂；如果连续两格是地雷，挂

3.     注意每一段的起点的取值

好吧，上代码……

type matrix=array[1..2,1..2]of extended;
     line=array[1..2]of extended;
var bas,a:matrix;
    res:line;
    ans:extended;
    n,i,j:longint;
    p,q:extended;
    mine:array[1..100]of longint;
procedure sort(l,r: longint);//把地雷位置排序（不要问我为什么开快排，我当时有点无聊）
var
i,j,x,y: longint;
begin
  i:=l;
  j:=r;
  x:=mine[(l+r) div 2];
  repeat
    while mine[i]<x do
    inc(i);
    while x<mine[j] do
    dec(j);
    if not(i>j) then
    begin
      y:=mine[i];
      mine[i]:=mine[j];
      mine[j]:=y;
      inc(i);
      j:=j-1;
    end;
  until i>j;
  if l<j then
  sort(l,j);
  if i<r then
  sort(i,r);
end;
function mul(a,b:matrix):matrix;//矩阵乘法
var t:matrix;i,j,k:longint;
begin
  for i:=1 to 2 do
    for j:=1 to 2 do
    begin
      t[i,j]:=0;
      for k:=1 to 2 do
      t[i,j]:=t[i,j]+a[i,k]*b[k,j];
    end;
  exit(t);
end;
function solve(a:line;b:matrix):line;//都是矩阵乘法，不过这里操作的是只有一列的矩阵
var t:line;i,j:longint;
begin
  fillchar(t,sizeof(t),0);
  for i:=1 to 2 do
    for j:=1 to 2 do
    t[j]:=t[j]+a[i]*b[i,j];
  exit(t);
end;
function calc(n:longint):extended;//分段计算过程
var res:matrix;a:line;
begin
  bas[1,1]:=p;
  bas[1,2]:=q;
  bas[2,1]:=1;
  bas[2,2]:=0;
  res[1,1]:=1;
  res[1,2]:=0;
  res[2,1]:=0;
  res[2,2]:=1;//任何矩阵的零次方都是单位矩阵
  a[1]:=0;
  a[2]:=1;//f[1]=1,f[0]=0
  dec(n);
  while n>0 do
  begin
    if n mod 2=1 then res:=mul(bas,res);
    n:=n shr 1;
    bas:=mul(bas,bas);
  end;
  a:=solve(a,res);//a=a*bas^n
  exit(a[1]*q);
end;
procedure main;
var last:longint;
begin
  for i:=1 to n do
  read(mine[i]);
  readln;
  sort(1,n);
  if mine[1]=1 then
  begin
    writeln('0.0000000');
    exit;
  end;
  for i:=2 to n do
  if mine[i]-mine[i-1]=1 then
  begin
    writeln('0.0000000');
    exit;
  end;
  ans:=1;
  last:=1;
  q:=1-p;
  for i:=1 to n do
  begin
    ans:=ans*calc(mine[i]-last+1);
    last:=mine[i]+1;
  end;
  writeln(ans:0:7);
end;
begin
  while not eof do
  begin
    readln(n,p); 
    main;
  end;
end.