Analytical Computing 
by Laurent Bernardin


Listing One
> a[0] := 11/2; a[1] := 61/11;
> for i from 2 to 20 do
    a[i] := evalf( 111-1130/a[i-1]+3000/(a[i-1]*a[i-2]) )
  end do;
       a[2] := 5.59016393
       a[3] := 5.63343106
       a[4] := 5.67464821
       a[5] := 5.71332183
       a[6] := 5.74899458
       a[7] := 5.77961409
       a[8] := 5.77331611
       a[9] := 5.17967242
      a[10] := -6.8391055
      a[11] := 191.5387186
      a[12] := 102.8102519
      a[13] := 100.1612232
      a[14] := 100.0095189
      a[15] := 100.0005641
      a[16] := 100.0000335
      a[17] := 100.0000020
      a[18] := 100.0000001
      a[19] := 100.0000000
      a[20] := 100.0000000

Listing Two
> for i from 2 to 20 do
    a[i] := evalf( 111-1130/a[i-1]+3000/(a[i-1]*a[i-2]), 30 )
  end do;
       a[2] := 5.5901639344262295081967213114
       a[3] := 5.6334310850439882697947214080
       a[4] := 5.6746486205101509630400832965
       a[5] := 5.7133290523805155490321989960
       a[6] := 5.7491209197026380437051448349
       a[7] := 5.7818109204856155794683379136
       a[8] := 5.8113142382939957232038462063
       a[9] := 5.8376565489587119615631669961
      a[10] := 5.8609515225161319729296295759
      a[11] := 5.8813772158414186066348580484
      a[12] := 5.8991539057900653804656369740
      a[13] := 5.9145249506789843093265503806
      a[14] := 5.9277414077768100768697258775
      a[15] := 5.9390504854613682145398481297
      a[16] := 5.9486874924846266324223440920
      a[17] := 5.9568707319889872021880625846
      a[18] := 5.9637987220073017812321934030
      a[19] := 5.9696491639651381786739608575
      a[20] := 5.9745793622911322066401902108

Listing Three
> for i from 21 to 30 do
    a[i] := evalf( 111-1130/a[i-1]+3000/(a[i-1]*a[i-2]), 30 )
  end do;
      a[21] := 5.9787313055716646903432098926
      a[22] := 5.9823017419650597836243768506
      a[23] := 5.9866906078267980996399642244
      a[24] := 6.0136522709040277799769525928
      a[25] := 6.4232153924054698100400652873
      a[26] := 12.7415627706574382005210260605
      a[27] := 58.9699458767796995910794892852
      a[28] := 95.8304070556576619076782493096
      a[29] := 99.7392043815332393272183617503
      a[30] := 99.9843246386708478051625918056

Listing Four
> for i from 2 to 100 do
    a[i] := 111-1130/a[i-1]+3000/(a[i-1]*a[i-2]);
    printf("a[%d] := %a\n",i,evalf(a[i]))
  end do:
      a[2] := 5.590163934
      a[3] := 5.633431085
      a[4] := 5.674648621
      a[5] := 5.713329052
      a[6] := 5.749120920
      a[7] := 5.781810920
      a[8] := 5.811314238
      a[9] := 5.837656549
      a[10] := 5.860951523
      a[11] := 5.881377216
      a[12] := 5.899153906
      a[13] := 5.914524951
      a[14] := 5.927741408
      a[15] := 5.939050485
      a[16] := 5.948687492
      a[17] := 5.956870732
      a[18] := 5.963798721
      a[19] := 5.969649144
      a[20] := 5.974579029
            ...
      a[100] := 5.999999988

Listing Five
> g := codegen[optimize](f);
g := proc(x, y, z, t)
local t9, t18, t2, t1, t3, t6, t28;
    t1 := x^2; t2 := t1^2; t3 := y^2; t6 := z^2;
    t9 := t^2; t18 := t3^2; t28 := t6^2;
    t2 - 3*t1*t3 - 2*t1*t6 - t1*t9 + 4*x*t3*z + 4*y*z*x*t
       + t18 - 2*t3*y*t - 2*t3*t6 + t3*t9 - 2*t6*y*t + t28
end proc
> codegen[cost](g);
  7 storage + 7 assignments + 28 multiplications + 11 additions

Listing Six
> h := codegen[GRADIENT](g);
h := proc(x, y, z, t)
local df, t1, t18, t2, t28, t3, t6, t9;
    t1 := x^2; t2 := t1^2; t3 := y^2; t6 := z^2;
    t9 := t^2; t18 := t3^2; t28 := t6^2;
    df := array(1 .. 7);
    df[7] := 1; df[6] := 1; df[5] := -t1 + t3;
    df[4] := 2*df[7]*t6 - 2*t1 - 2*t3 - 2*y*t;
    df[3] := 2*df[6]*t3 - 3*t1 + 4*x*z - 2*y*t - 2*t6 + t9;
    df[2] := 1; df[1] := 2*df[2]*t1 - 3*t3 - 2*t6 - t9;
    (4*t3*z + 4*z*y*t + 2*df[1]*x,
     4*z*x*t - 2*t3*t - 2*t6*t + 2*df[3]*y,
     4*x*t3 + 4*x*y*t + 2*df[4]*z,
     4*y*x*z - 2*t3*y - 2*t6*y + 2*df[5]*t)
end proc

Listing Seven
Toeplits matrix
> codegen[C](h,ansi);
void h(double x,double y,double z,double t,double crea_par[4])
{
  double df[7];
  double t1;
  double t18;
  double t2;
  double t28;
  double t3;
  double t6;
  double t9;
  {
    t1 = x*x;
    t2 = t1*t1;
    t3 = y*y;
    t6 = z*z;
    t9 = t*t;
    t18 = t3*t3;
    t28 = t6*t6;
    df[6] = 1.0;
    df[5] = 1.0;
    df[4] = -t1+t3;
    df[3] = 2.0*df[6]*t6-2.0*t1-2.0*t3-2.0*y*t;
    df[2] = 2.0*df[5]*t3-3.0*t1+4.0*x*z-2.0*y*t-2.0*t6+t9;
    df[1] = 1.0;
    df[0] = 2.0*df[1]*t1-3.0*t3-2.0*t6-t9;
    crea_par[0] = 4.0*t3*z+4.0*z*y*t+2.0*df[0]*x;
    crea_par[1] = 4.0*z*x*t-2.0*t3*t-2.0*t6*t+2.0*df[2]*y;
    crea_par[2] = 4.0*x*t3+4.0*x*y*t+2.0*df[3]*z;
    crea_par[3] = 4.0*y*x*z-2.0*t3*y-2.0*y*t6+2.0*df[4]*t;
    return;
  }
}