program matrices(input,output);
{ compute approx derivative N, and its inverse }
const
  hdim=125;
  edim=251;
  lmax=29;
  pdeg=100;
  vdim=130;

%include 'types.h';

var
  i,j: integer;
  cov,beta2,r,st: scalar;
  vt,vp,vc: vector;
  m,minv: matrix;

%include 'reps.i';
%include 'scalar.i';
%include 'vector.i';
%include 'linear.i';

procedure pconvinit(var c: scalar; var v: vector);
{ v:=exp(c.) truncated, used in pconv }
var
  i: integer;
begin { pconvinit }
  v[0]:=sone;
  for i:=1 to vdim do begin             { vf[i]=(c^i)/(i!)  }
    sprod(c,v[i-1],v[i]);
    isdiv(i,v[i]);
  end;
end { pconvinit };

begin { matrices }
  sinit;

  writeln('matrices: initializing convolution');
  sread('beta2.s',beta2);
  sread('cov.s',cov);
  pconvinit(cov,vc);
  writeln;

  writeln('matrices: reading and scaling approx fixed point');
  szero(r);
  vread('approx.v',r,vt);
  for i:=0 to pdeg do ismult(2,vt[i]);
  for i:=pdeg+1 to vdim do szero(vt[i]);
  vscale(beta2,vt,vp);
  writeln;

  write('matrices: computing columns of approx tangent map');
  for j:=0 to lmax do begin
    writeln;
    write(j:3);
    flush(output);
    pconv(lmax,vc,vp,vt);
    for i:=0 to lmax do minv[i,j]:=vt[i];
    for i:=vdim-1 downto 1 do sprod(beta2,vp[i-1],vp[i]);
    szero(vp[0]);
  end;
  writeln;

  writeln('matrices: computing and inverting the inverse of m');
  for i:=0 to lmax do begin
    for j:=0 to lmax do sneg(minv[i,j]);
    ssum(sone,minv[i,i],st);
    minv[i,i]:=st;
  end;
  mwrite('minv.m',minv);
  minverse(minv,m);
  mwrite('m.m',m);
  writeln;

  sdone('matrices');
end { matrices }.