第六课 递推法

 递推就是逐步推导的过程。我们先看一个简单的问题。

例20：一个数列的第0项为0，第1项为1，以后每一项都是前两项的和，这个数列就是著名的裴波那契数列，求裴波那契数列的第N项。

分析：我们可以根据裴波那契数列的定义：从第2项开始，逐项推算，直到第N项。因此可以设计出如下算法：

F[0] := 1; F[1] := 2;

FOR I := 2 TO N DO

F[I] := F[I – 1] + F[I – 2];

从这个问题可以看出，在计算裴波那契数列的每一项目时，都可以由前两项推出。这样，相邻两项之间的变化有一定的规律性，我们可以将这种规律归纳成如下简捷的递推关系式：F_n=g(F_n-1)，这就在数的序列中，建立起后项和前项之间的关系。然后从初始条件（或是最终结果）入手，按递推关系式递推，直至求出最终结果（或初始值）。很多问题就是这样逐步求解的。

对一个试题，我们要是能找到后一项与前一项的关系并清楚其起始条件（或最终结果），问题就可以递推了，接下来便是让计算机一步步了。让高速的计算机从事这种重复运算，真正起到“物尽其用”的效果。

递推分倒推法和顺推法两种形式。算法流程如下： 

一、倒推法

所谓倒推法，就是在问题的解或目标是由初始值递推得到的问题中，已知解或目标，根据递推关系，采用倒推手段，一步步的倒推直至求得这个问题的初始陈述的方法。因为这类问题的运算过程是一一映射的，故可分析其递推公式。看看下面的例题。

 

例21：贮油点

一辆重型卡车欲穿过1000公里的沙漠，卡车耗汽油为1升/公里，卡车总载油能力为500公升。显然卡车装一次油是过不了沙漠的。因此司机必须设法在沿途建立若干个贮油点，使卡车能顺利穿过沙漠。试问司机如怎样建立这些贮油点？每一贮油点应存储多少汽油，才能使卡车以消耗最少汽油的代价通过沙漠？

编程计算及打印建立的贮油点序号，各贮油点距沙漠边沿出发的距离以及存油量。格式如下：

No.	Distance(k.m.)	Oil(litre)
1	× ×	× ×
2	× ×	× ×
…	… …	… …

 

分析：

设Way[I]——第I个贮油点到终点（I=0）的距离；

oil[I]——第I个贮油点的贮油量；

我们可以用倒推法来解决这个问题。从终点向始点倒推，逐一求出每个贮油点的位置及存油量。图19表示倒推时的返回点。

从贮油点I向贮油点I+1倒推的方法是：卡车在贮油点I和贮油点I+1间往返若干次。卡车每次返回I+1点时应该正好耗尽500公升汽油，而每次从I+1点出发时又必须装足500公升汽油。两点之间的距离必须满足在耗油最少的条件下，使I点贮足I*500公升汽油的要求（0≦I≦n-1）。具体来说，第一个贮油点I=1应距终点I=0处500km，且在该点贮藏500公升汽油，这样才能保证卡车能由I=1处到达终点I=0处，这就是说

Way[I]=500；oil[I]=500；

为了在I=1处贮藏500公升汽油，卡车至少从I=2处开两趟满载油的车至I=1处，所以I=2处至少贮有2*500公升汽油，即oil[2]=500*2=1000；另外，再加上从I=1返回至I=2处的一趟空载，合计往返3次。三次往返路程的耗油量按最省要求只能为500公升，即d₁₂=500/3km，Way[2]=Way[1]+d₁₂=Way[I]+500/3

此时的状况如图20所示。

为了在I=2处贮藏1000公升汽油，卡车至少从I=3处开三趟满载油的车至I=2处。所以I=3处至少贮有3*500公升汽油，即oil[3]=500*3=1500。加上I=2至I=3处的二趟返程空车，合计5次。路途耗油亦应500公升，即d₂₃=500/5，

Way[3]=Way[2]+d₂₃=Way[2]+500/5；

此时的状况如图21所示。

依次类推，为了在I=k处贮藏k*500公升汽油，卡车至少从I=k+1处开k趟满载车至I=k处，即oil[k+1]=(k+1)*500=oil[k]+500，加上从I=k返回I=k+1的k-1趟返程空间，合计2k-1次。这2k-1次总耗油量按最省要求为500公升，即d_k,k+1=500/(2k-1)，

Way[k+1]=Way[k]+d_k,k+1=Way[k]+500/(2k-1)；

此时的状况如图22所示。

最后，I=n至始点的距离为1000-Way[n],oil[n]=500*n。为了在I=n处取得n*500公升汽油，卡车至少从始点开n+1次满载车至I=n，加上从I=n返回始点的n趟返程空车，合计2n+1次，2n+1趟的总耗油量应正好为（1000-Way[n]）*(2n+1)，即始点藏油为oil[n]+(1000-Way[n])*(2n+1)。

程序设计如下：

program Oil_lib;

   var

     K: Integer;     {贮油点位置序号}

     D,                  {累计终点至当前贮油点的距离}

     D1: Real;         {I=n至终点的距离}

     Oil, Way: array [1 .. 10] of Real;

     i: Integer;       

   begin

     Writeln(‘No.’, ‘Distance’:30, ‘Oil’:80);

     K := 1;

     D := 500;       {从I=1处开始向终点倒推}

     Way[1] := 500;

     Oil[1] := 500;

     repeat

       K := K + 1;

       D := D + 500 / (2 * K – 1);

       Way[K] := D;

       Oil[K] := Oil[K – 1] + 500;

     until D >= 1000;

     Way[K] := 1000;                 {置始点到终点的距离值}

     D1 := 1000 – Way[K – 1];          {求I=n处至至点的距离}

     Oil[K] := D1 * (2 * k + 1) + Oil[K – 1];  {求始点贮油量}

     {由始点开始，逐一打印至当前贮油点的距离和贮油量}

     for i := 0 to K do  

       Writeln(i, 1000 – Way[K – i]:30, Oil[K – i]:80);

   end. 

二、顺推法

顺推法是从边界条件出发，通过递推关系式推出后项值，再由后项值按递推关系式推出再后项值……，依次类推，直至从问题初始陈述向前推进到这个问题的解为止。

看看下面的问题。 

例22昆虫繁殖

科学家在热带森林中发现了一种特殊的昆虫，这种昆虫的繁殖能力很强。每对成虫过x个月产y对卵，每对卵要过两个月长成成虫。假设每个成虫不死，第一个月只有一对成虫，且卵长成成虫后的第一个月不产卵(过X个月产卵)，问过Z个月以后，共有成虫多少对？x>=1,y>=1,z>=x

输入：x,y,z的数值

输出：成虫对数

事例：

    输入：x=1 y=2 z=8

    输出：37

分析：首先我们来看样例：每隔1个月产2对卵，求过8月（即第8+1=9月）的成虫个数

月份	1	2	3	4	5	6	7	8	9	…
新增卵	0	2	2	2	6	10	14	26	46	…
成虫	1	1	1	3	5	7	13	23	37	…

设数组A[i]表示第I月新增的成虫个数。

由于新成虫每过x个月产y对卵，则可对每个A[I]作如下操作:

A[i+k*x+2]:=A[i+k*x+2]+A[i]*y  （1<=k, I+k*x+2<=z+1）

因为A [i]的求得只与A[1]~A[i-1]有关，即可用递推求法。

则总共的成虫个数为：   

    程序如下：

program exam22;

var  x,y,z,i :integer;

ans :longint;

a    :array[1..60]of longint;

procedure add(i:integer);

  var j   :integer;

  begin

    j:=i+2+x; {新生成虫要过x 个月才开始产卵，即第I+2+x个月才出现第一群新成虫}

    repeat

      a[j]:=a[j]+a[i]*y;             {递推}

      j:=j+x

    until j>z+1

  end;

begin

  readln(x,y,z);

  a[1]:=1;                      {初始化}

  for i:=1 to z do add(i);          {对每个A[I]进行递推}   

  ans:=0;

  for i:=1 to z+1 do ans:=ans+a[i]; {累加总和}

  writeln(ans);

end.

 

例23：实数数列（NOI94第3题）

一个实数数列共有N项，已知

a_i=(a_i-1-a_i+1)/2+d，(1<I<N) (N<60)

键盘输入N，d，a₁，a_n，m，输出a_m。

输入数据均不需判错。

 

分析：

根据公式a_i=(a_i-1-a_i+1)/2+d 变形得，a_i+1=a_i-1-2a_i+2d，因此该数列的通项公式为：a_i=a_i-2-2a_i-1+2d，已知a₁，如果能求出a₂，这样就可以根据公式递推求出a_m

∵  a_i=a_i-2-2a_i-1+2d                 ……①

=a_i-2-2（a_i-3-2a_i-2+2d）+2d

=-2a_i-3+5（a_i-4-2a_i-3+2d）-2d

=5a_i-4-12a_i-3+8d

……

一直迭代下去，直到最后，可以建立a_i和a₁与a₂的关系式。

设a_i=P_ia₂+Q_id+R_ia₁，我们来寻求P_i,Q_i,R_i的变化规律。

∵  a_i=a_i-2-2a_i-1+2d

∴  a_i=P_i-2a₂+Q_i-2d+R_i-2a₁-2（P_i-1a₂+Q_i-1d+R_i-1a₁）+2d

      =(P_i-2-2P_i-1)a₂+(Q_i-2-2Q_i-1+2)d+(R_i-2-2R_i-1)a₁

∴  P_i=P_i-2-2P_i-1                                   ……②

Q_i=Q_i-2-2Q_i-1+2                                   ……③

R_i=R_i-2-2R_i-1                                ……④

显然，P₁=0 Q₁=0 R₁=1  （i=1）

P₂=1 Q₂=0 R₂=0   （i=2）

将初值P₁、Q₁、R₁和P₂、Q₂、R₂代入②③④可以求出P_n、Q_n、R_n

∵ a_n=P_na₂+Q_nd+R_na₁

∴ a₂=(a_n-Q_nd+R_na₁)/P_n

然后根据公式①递推求出a_m，问题解决。

但仔细分析，上述算法有一个明显的缺陷：在求由于在求a₂要运用除法，因此会存在实数误差，这个误差在以后递推求a_m的过程又不断的扩大。在实际中，当m超过30时，求出的a_m就明显偏离正确值。显然，这种算法虽简单但不可靠。

为了减少误差，我们可设计如下算法：

∵ a_i=P_ia₂+Q_id+R_ia₁

         =P_i-1a₃+Q_i-1d+R_i-1a₂

         =P_i-2a₄+Q_i-2d+R_i-2a₃

       ……

     =P_i-2+ka_k+Q_i-2+kd+R_i-2+ka_k-1

∴ a_n=P_n-k+2a_k+Q_n-k+2d+R_n-k+2a_k-1

a_k=(a_n-Q_n-k+2d+R_n-k+2a_k-1)/P_n-k+2         ……⑤

根据公式⑤，可以顺推a₂、a₃、…、a_M。虽然仍然存在实数误差，但由于P_n-k+2递减，因此最后得出的a_m要比直接利用公式①精确得多。

 

程序如下：

program NOI94_3;

 const

   MaxN = 60;

 var

   N, M, i: Integer;

   D: Real;

   A: array [1 .. MaxN] of Real;       

   F: array [1 .. MaxN, 1 .. 3] of Real;   

{F[i,1]：对应P_i；F[i,2]：对应Q_i；F[i,3]：对应R_i}

 

 procedure Init;

   begin

     Write(‘N, M, D =’);

     Readln(N, M, D);               {输入项数、输出项序号和常数}

     Write(‘A1, A’, N, ‘ =’);

     Readln(A[1], A[N]);            {输入a₁和a_n}

   end;

 

 procedure Solve;

 {根据公式P_ißP_i-2-2*P_i-1,Q_ißQ_i-2-2*Q_i-1,R_ißR_i-2-2*R_i-1求P_i、Q_i、R_i }

  begin

     F[1, 1] := 0; F[1, 2] := 0; F[1, 3]:= 1;

     F[2, 1] := 1; F[2, 2] := 0; F[2, 3] := 0;

     for i := 3 to N do

begin

         F[i, 1] := F[i – 2, 1] – 2 * F[i – 1, 1];

         F[i, 2] := F[i – 2, 2] – 2 * F[i – 1, 2] + 2;

         F[i, 3] := F[i – 2, 3] – 2 * F[i – 1, 3];

      end;

  end;

 

procedure Main;

  begin

    Solve;     

    {递推A₂…A_m}

    for i := 2 to M do

      A[i]:=(A[N]–F[N–i+2,2]*D–F[N–i+2,3]*A[i–1])/F[N–i+2,1];

    Writeln(‘a’, m, ‘ =’, A[M]:20 :10);    

  end;

 

 begin

   Init;           

   Main;           

 end.