// local copy taken from module cvs2005 mtx07
// removing the spline functions not needed here
{$H+} // ansi strings
{$C+} // assertions on

unit mtx3;
// TODO -> there is a bug in dvd() with slices... ? DONE afaik
{ linklib j}   // link to objects in libj.a
{ L spline.o}  // link to local spline.o

interface

uses strutils;

type

// TODO assuming e==NIL is NOT true for local variables
// C style static mtx not possible use global to a unit instead
// global mtx are init with e==NIL ( but not local to a proc... )
// the fundamental reason is that one would need to alloc and free
// dyn. mem. on entry/exit of a proc... so rule is NO LOCAL mtx !
  mtx = record
    smx: boolean; // if it is a sub matrix dc<>col ou lgn if transpose
    nam: string;  //
    lgn,          // number of rows (ligne in french)
    col,          //        of columns
    max: longint; // e[1..max]: max+1 double allocated as e[0] is ignored
    dl: longint;  // index offset between lines
    dc: longint;  // index offset between columns
    t: boolean;   // true if transposed
    e: ^double;   // e[1] is elem[1,1]
    dyn: boolean; // if e is dynamic mem check of mtxfre
  end;

const
  ElmDelim=[' ',#9,','];
  RowDelim=[';',#10];
  EpsDouble = 1e-30;   // TODO see Forsythe & al.  to get a better value
  WordDelim=[#0..' ']; // assert ascii char set !
  BLANK=[#32,#9];      //

// interim mtx put on stdout all fields...
// for debug
procedure mpt(a:mtx);
procedure mpt(a:mtx; s:string);

function mem(var r:mtx; m:longint): boolean;
function mem(r:mtx): longint;
//
function mtystr(s:string):boolean;
// function read1st(var f:text; var s: string; var nam: string; var n:longint):boolean;
function mrd(var r: mtx; var f: text): boolean;
function mrd(var r:mtx):boolean;
function mrd(var r: mtx; var f: text; mnam:string): boolean;
function mrd(var r:mtx; mnam:string):boolean;
function mrd(var r:mtx; fnam:string; mnam:string):boolean;
//
function mwr(a:mtx; var f:text): boolean;
function mwr(a:mtx): boolean;
// write a in file s ( rewrite )
function mwr(a:mtx; s:string): boolean;

//
function elm(a:mtx; i, j: longint):double;
function elm(b:mtx; a:mtx):boolean;
function elm(a:mtx; i, j: longint; x:double):boolean;
function elm(b:mtx; x:double):boolean;
//
function trp(a:mtx):mtx;
//
function prd(var r: mtx; a, b: mtx):boolean; // TODO *** 1 ***
//
function smx(a:mtx; i1,j1:longint; lgn, col:longint):mtx;
//
function col(a:mtx; j1: longint): mtx;
function col(a:mtx; j1, j2: longint): mtx;
function col(a:mtx): longint;
//
function lgn(a:mtx; i1: longint): mtx;
function lgn(a:mtx; i1,i2: longint): mtx;
function lgn(a:mtx): longint;
function lgn(a:mtx; low, high:double):mtx;
// yet another variant of the lgn() selection tool
// remove rows of "a" if a[1] is not in one of the brackets
// defined by a min max pair in each line of r ( ROI ) mtx
// NOTE that size and content of a ARE MODIFIED
// TODO test on trp() or smx()... so far use only on PLAIN mtx
procedure lgn(var a:mtx; r:mtx);
//
function nam(p:mtx):string;
function nam(var p:mtx; name:string):boolean;
//
function cpy(var r:mtx; a:mtx):boolean; // r can be smx
//
function map(var x:double; lgn, col:longint):mtx;
function map(var a:mtx; lgn, col:longint):mtx;
function matinv(a:mtx; var  d: double): boolean;
//
function inv(var r:mtx; a:mtx; var d:double):boolean; // r must be a plain mtx
function inv(var r:mtx; a:mtx):boolean;               // r must be a plain mtx
function inv(r:mtx):boolean;                          // r must be a plain mtx
function inv(r:mtx; var d:double):boolean;            // r must be a plain mtx
//
function dim(var r:mtx; lgn, col:longint):boolean;
function dim(var r:mtx; lgn, col:longint; x:double):boolean;
function zer(var r:mtx; lgn, col: longint):boolean;
function zer(var r:mtx; a:mtx):boolean;
function one(var r:mtx; lgn, col: longint):boolean;
function one(var r:mtx; a:mtx):boolean;
function eye(var r:mtx; lgn, col: longint):boolean;
function eye(var r:mtx; a:mtx):boolean;
//
function ptr(a:mtx):pointer;
//
function mul(var r: mtx; a:mtx; b:mtx):boolean;
function mul(a: mtx; x: double):boolean;
//
function dvd(var r: mtx; a:mtx; b:mtx):boolean;
function dvd(a: mtx; x: double):boolean;

// TODO check trp() in arg. a and b ??? may be a bug in this case
// TODO same for all similar functions
function add(var r: mtx; a:mtx; b:mtx):boolean;
function add(a: mtx; x: double):boolean;
// elm/elm add by row (false) or col (true)
function add(var r: mtx; a:mtx; t:boolean):boolean;
//
function sub(var r: mtx; a:mtx; b:mtx):boolean; // TODO *** 1 ***
function sub(a: mtx; x: double):boolean;
//
procedure writx(var f: text; x: double);
procedure writx(x:double);
//
function chk(var r:mtx; a, b:mtx):boolean;
// TODO one() zer() eye() ***(r,a,true) summation/product of rows/columns
// TODO how about appending ? at least by lines cat(r,a,b,c...)
// TODO stat. may be mmm(r,a) calc. min max mean other moments in r ( /column )
// ? all summation like functions should operate on columns (use trp() to do by line)

{
// cubic spline interpolation must link with fortran code spline.f
function spl(var s:mtx; x, y: mtx):boolean;
function sev(var r:mtx; u:mtx; s:mtx):boolean;
function sev(u:double; s:mtx):double;
}

// scilab's linspace...  r = x0:dx:x1;
function lsp(var r:mtx; x0, dx, x1: double):boolean;

// defines n steps ( hence n+1 points ) between x0 and x1
// this should be the default variant better than the previous one
function lsp(var r:mtx; x0, x1: double; n:longint):boolean;

// make a diagonal mtx from a vector
function dia(var r:mtx; v: mtx):boolean;
// make a smx vector from the diagonal of a
function dia(a:mtx):mtx;

// basic linear fits...
// x[n,1] y[n,1] n>2
// xy[n,2] matrix n>2
procedure lin(x,y: mtx; var slope, y0, s2y: double);
procedure lin(xy: mtx; var slope, y0, s2y: double);

// largest/smallest element of a row or col vector
// assert lgn=1 or col=1
function max(a:mtx):double;
function min(a:mtx):double;

// TODO avg <- returns mean std of an array and/or mom(ents) kur(tosis) and the like

// convert a string to an mtx ( scilab like syntax )
// intended use to input mtx from a command line
// option. -m "1,2,3;4,5,6" with
// assert( mxs(a,optarg) );
function mxs(var r:mtx; s:string): boolean;

// Area below xy curve assert sorted by x
// 2d variant default kx=1 ky=2
function axy(a:mtx; kx, ky:longint): double;
function axy(a:mtx):double;

// mean and stdev of a vector
// TODO mtx variants col or row wise...
function avg(a:mtx; var dst:double):double;

// tool to write an identifier and a number ( longint or double )
// as a named (1,1) matrix...
procedure wnb(s:string; x:double; var f:text);
procedure wnb(s:string; n:longint; var f:text);
procedure wnb(s:string; x:double);
procedure wnb(s:string; n:longint);

implementation

// -------- wnb() avatars
var
  wnb0: mtx; // hidden mtx used by wnb avatars

procedure wnb(s:string; x:double; var f:text);
begin
  wnb0 := map(x,1,1);
  nam(wnb0,s);
  mwr(wnb0,f)
end;

procedure wnb(s:string; n:longint; var f:text);
var x: double;
begin
  x := double(n);
  wnb0 := map(x,1,1);
  nam(wnb0,s);
  mwr(wnb0,f);
end;

procedure wnb(s:string; x:double);
begin
  wnb(s,x,stdout)
end;

procedure wnb(s:string; n:longint);
begin
  wnb(s,n,stdout)
end;
// *** end of wnb() avatars ********

// interim stdout
procedure mpt(a:mtx);
begin
  writeln('# ====================================');
  writeln('# smx, trp, dyn: ',a.smx,' ',a.t,' ',a.dyn);
  writeln('# max, lgn, col, dl, dc: ',a.max,'   ',a.lgn,' ',a.col,'   ',a.dl,' ',a.dc);
  writeln('# ====================================');
  writeln;
  flush(output);
end;

// interim stdout avatar with 2d arg
procedure mpt(a:mtx; s:string);
begin
  writeln('# ========== mpt ====',s,'====');
  flush(output);
  mpt(a);
end;

// Prototypes of fortran routines spline_ and seval_
// source spline.f; compile with:  gfortran -c spline.f
// { linklib j}   // link to object in libj.a
// {$L spline.o}  // link to local spline.o

procedure spline_(var n:longint; x, y, b, c, d:pointer); cdecl; external;
function seval_(var n:longint; var u:double;  x, y, b, c, d:pointer):double; cdecl; external;

// Realloc mem. if needed and possible
function mem(var r:mtx; m:longint): boolean;
var p: pointer;
begin
  mem := false;
  //writeln(stderr,'**** in mem (',r.nam,') **** m=',m,'  r.e=', longint(r.e));
  assert(m>=0);
  if r.e = NIL then begin // WARNING FAILS WITH LOCAL MTX !!! VERY BAD
    //writeln(stderr,'***** mem ***** First alloc. *****');
    r.nam := '';
    r.lgn := 0;
    r.col := 0;
    r.t := false;
    r.dl := 0;
    r.dc := 0;
    r.max := m;
    r.dyn := true; // TODO check it should be reset in all smx functions
    r.smx := false;
    p := ReAllocMem(r.e, (m+1)*sizeof(double));
    assert( p <> NIL );
    r.e := p
  end else assert( not r.smx ); // I guess this make sense YES cf. zer one eye
  if ( r.max < m ) then begin
    //writeln(stderr,'***** mem ***** Realloc. *****');
    p := ReAllocMem(r.e, (m+1)*sizeof(double));
    assert( p <> NIL );
    r.max := m;
    r.e := p
  end;
  //writeln(stderr,'***** exit mem (',r.nam,') ***** r.e =', longint(r.e) );
  mem := true
end;

// return max. vector size
function mem(r:mtx): longint;
begin
  mem := r.max
end;

// replace buggy strutils.isEmptyStr()
function mtystr(s:string):boolean;
var i, n: longint;
begin
  mtystr := false;
  n := length(s);
  for i := 1 to n do begin
    if not ( s[i] in BLANK ) then exit;
  end;
  mtystr := true
end;


// new version to correct a pb with reading anonymous mtx after named ones...
// mnam must be in the arg. list... and a named mtx skipped if searching a noname
// prototype is different so I keep the same name.
// readln until it finds either a letter in col 1 ( first word is name of matrix )
// or a number (noname mtx i.e. name = '' )
// return true if found a match, nam is the name of the mtx or '' for a plain xy table
// n is the number of words in the first data line of the mtx/table s is this string
// which can be processed to retrieve the n numerical values
function read1st(mnam:string; var f:text; var s: string; var nam: string; var n:longint):boolean;
var s0: string; ok: boolean; p:longint; x:double; code:integer;
begin
  read1st := false;
  ok := false;
  repeat
    readln(f,s);
    if eof(f) then exit;
    if mtystr(s) then continue;
    if pos('#',s)>0 then continue;
    s0 := ExtractWordPos(1,s+' x',WordDelim,p); // hope this func. is ok ?
 // if mnam is neither '' nor '*' we can skip val() and set code=-1
    if (mnam='') or (mnam='*') then val(s0,x,code) else code:=-1;
    if  code=0 then begin
      // we assume that s is the first line of a plain noname xy table
      // and s0 the first number of the table
      nam := '';
      // *********** OK only if nam agree with mnam
      ok := (nam = mnam) or (mnam='*'); // ************
      // if name not ok then skip until empty string (or eof) TODO
      if not ok then while not mtystr(s) do readln(f,s)
    end else begin
      // s0 must be a name
      nam := s0;
      ok := (nam = mnam) or (mnam='*');
      if ok then readln(f,s) else while not mtystr(s) do readln(f,s)
      // s must be the first data line of mtx mnam
    end;
  until ok;
  // count the number of item in this first line
  // TODO achtung n returned is integer !!!!!!!!
  n := WordCount(s,WordDelim);
  // assert all words are valid numbers...
  read1st := n>0;
end;

function mrd(var r: mtx; var f: text; mnam:string): boolean;
var
  s, s0, nam: string;
  i, j, k, m, n: longint;
  x: double;
  code: integer;
begin
  mrd := false;
  nam := '';
  if not read1st(mnam,f,s,nam,n) then exit;
  assert( n>0 );
  // first default alloc is n ( 1 line then later 2n 4n 8n... as needed )
  assert( mem(r,n) ); // sanity checks made here including "not r.smx"
  assert( not (r.e=NIL) ); // hopefully
  r.nam := nam;
  k := 0;
  m := 0;
  j := n;
  while true do begin
    if pos('#',s)=0 then begin
        inc(m);
        //if  m > ( mx div n ) then exit;
        for i := 1 to n do begin
            s0 := ExtractWord(i,s,WordDelim);
            assert(s0<>'');
            //writeln(stderr,'---------------------> s0= "',s0,'"');
            val(s0,x,code);
            //writeln(stderr,'---------------------> x= ',x:15:5);
            //assert(code=0);
            if code<>0 then exit;
            inc(k);
            if k > r.max then begin
              //writeln(stderr,'--- mrd: k, r.max, r.e: ',k,', ',r.max,', ',longint(r.e));
              assert( mem(r, r.max*2) );
              //writeln(stderr,'After realloc...  r.e=',longint(r.e),'  r.max=',r.max);
            end;
            r.e[k] := x;
        end;
    end;
    {$I-}
    readln(f,s);
    {$I+}
    j := WordCount(s,WordDelim);
    if j=0 then break;
    if j<>n then exit;
  end;
  r.dl := 1;
  r.dc := n;
  r.col := n;
  r.lgn := m;
  mrd := true
end;

// read any matrix from a stream
function mrd(var r: mtx; var f: text): boolean;
begin
  mrd := mrd(r,f,'*')
end;

function mrd(var r:mtx):boolean;
begin
  mrd := mrd(r,input)
end;

function mrd(var r:mtx; mnam:string):boolean;
begin
  mrd := mrd(r,input,mnam)
end;

function mrd(var r:mtx; fnam:string; mnam:string):boolean;
var f: text;
begin
  mrd := false;
  assign(f,fnam);
  reset(f);
  if ioresult <> 0 then exit;
  mrd := mrd(r,f,mnam)
end;

// Set of functions to write a matrix to a text file
// Ugly kludge to make a nice generic double format
// TODO replace using sysutils:format function TODO TODO TODO
// I tried the ffgeneral format and found some very poor formatting
// afaik the the following kludge does a nicer job... may be slow?
procedure writx(var f: text; x: double);
const
  r=4.99E-7;
  wd=15;
  sf=5;
var
  y:double; p,pp,de,w,dc:longint;
  intg,n,ne:boolean;
begin
  if (abs(x)<=maxint) then intg:=abs(x-trunc(x))=0 else intg:=false;
  if intg then write(f,trunc(x):wd) else begin
    n := x<0;
    if n then y:=-x else y:=x;
    p:=0;
    while (y>=10) do begin y:=y/10; p:=succ(p) end;
    while (y<1)   do begin y:=y*10; p:=pred(p) end;
    if (p=0) then write(f,x:wd:sf-1) else begin
      ne:=p<0;
      if ne then pp:=-p else pp:=p;
      de:=2;
      while (pp>9) do begin pp:=pp div 10; de:=succ(de) end;
      if ne then de:=succ(de);
      y:=y+double(r);
      w:=wd-de;
      if (p>-(sf-2)) and (p<sf) then begin
        dc:=pred(sf-p);
        if (dc<1) then dc:=1;
        write(f,x:wd:dc)
      end else
        if n then write(f,-y:w:pred(sf),'E',p:1)
             else write(f,y:w:pred(sf),'E',p:1)
    end
  end
end;

procedure writx(x:double);
begin
  writx(stdout,x);
end;

// write mtx to a STREAM
function mwr(a:mtx; var f:text): boolean;
var
   i, j, ia, ia1: longint;
begin
  mwr := false;
  //writeln(f);
  if a.nam<>'' then writeln(f,a.nam);
  ia1 := 1;
  for i := 1 to a.lgn do begin
    ia := ia1;       // point to first elem of next line
    for j := 1 to a.col do begin
      writx(f,a.e[ia]);
      write(f,#9);
      ia := ia + a.dl; // point to next elem of current line
    end;
    ia1 := ia1 + a.dc; // update pointer to first elem of next line
    writeln(f);
  end;
  writeln(f);
  flush(f);
  mwr := true;
end;

// write to stdout
function mwr(a:mtx): boolean;
begin
  mwr := mwr(a,stdout);
end;

// (re)write to file name "s" wrt to stream to append
function mwr(a:mtx; s:string): boolean;
var
  f: text;
begin
  mwr := false;
  {$I-}
  assign(f,s);
  rewrite(f);
  if ioresult<>0 then exit;
  mwr(a,f);
  flush(f);
  close(f);
  mwr := true
  {$I+}
end;

// returns x := a[i,j]
function elm(a:mtx; i, j: longint):double;
var k: longint;
begin
  assert( (i>0) and (i<=a.lgn) );
  assert( (j>0) and (j<=a.col) );
  k := 1 + pred(i)*a.dc + pred(j)*a.dl;
  elm := a.e[k];
end;

// assigns a[i,j] := x
function elm(a:mtx; i, j: longint; x:double):boolean;
var k: longint;
begin
  elm := false;
  assert( (i>0) and (i<=a.lgn) );
  assert( (j>0) and (j<=a.col) );
  k := 1 + pred(i)*a.dc + pred(j)*a.dl;
  a.e[k] := x;
  elm := true
end;

// b := a  element/element copy of a into b.
// a and b  plain or smx with the SAME number of lgn and col.
// Therefore b MUST be dimensionned before:
//       dim(b,lgn(a),2);
//       elm(col(b,1),col(a,3));
//       elm(col(b,2),col(a,2));
function elm(b:mtx; a:mtx):boolean;
var i, j, lgn, col, ia, ib, ia1, ib1: longint;
begin
  elm := false;
  assert(a.lgn = b.lgn);
  assert(a.col = b.col);
  lgn := a.lgn;
  col := a.col;
  ia1 := 1;
  ib1 := 1;
  for i := 1 to lgn do begin
    ia := ia1;
    ib := ib1;
    for j := 1 to col do begin
       b.e[ib] := a.e[ia];
       ia := ia + a.dl;
       ib := ib + b.dl;
    end;
    ia1 := ia1 + a.dc;
    ib1 := ib1 + b.dc;
  end;
  elm := true
end;

// assign b[*,*] := x
function elm(b:mtx; x:double):boolean;
var i, j, lgn, col, ib, ib1: longint;
begin
  elm := false;
  lgn := b.lgn;
  col := b.col;
  ib1 := 1;
  for i := 1 to lgn do begin
    ib := ib1;
    for j := 1 to col do begin
       b.e[ib] := x;
       ib := ib + b.dl;
    end;
    ib1 := ib1 + b.dc;
  end;
  elm := true
end;

function zer(var r:mtx; lgn, col: longint):boolean;
begin
  zer := dim(r,lgn,col,0.0)
end;

function zer(var r:mtx; a:mtx):boolean;
begin
  zer := dim(r,lgn(a),col(a),0.0);
end;

function one(var r:mtx; lgn, col: longint):boolean;
begin
  one := dim(r,lgn,col,1.0)
end;

function one(var r:mtx; a:mtx):boolean;
begin
  one := dim(r,lgn(a),col(a),1.0);
end;

function eye(var r:mtx; lgn, col: longint):boolean;
var
  n, k: longint;
begin
  eye := dim(r,lgn,col,0.0);
  k := 1;
  n := col;
  if n > lgn then n := lgn;
  while n>0 do begin
    r.e[k] := 1.0;
    dec(n);
    k := k + col + 1
  end
end;

function eye(var r:mtx; a:mtx):boolean;
begin
  eye := eye(r,lgn(a),col(a));
end;

// transpose mtx access to a
function trp(a:mtx):mtx;
begin
  trp := a;
  trp.t := not a.t;
  trp.lgn := a.col;
  trp.col := a.lgn;
  if trp.t then begin
    trp.dl :=  a.dc;
    trp.dc :=  a.dl;
  end else begin
    trp.dl := a.dl;
    trp.dc := a.dc;
  end
end;

// If "r" is a submatrix ( field smx true ? TODO quid de trp() ) then
// returns true if product r = a*b is ok
// if "r" is a plain matrix
// redim and realloc dyn. mem if needed...
function chkprd(var r:mtx; a, b:mtx):boolean;
begin
  chkprd := false;
  if (a.col<>b.lgn) then exit;
  if r.e = NIL then r.smx := false;
  if (r.smx) then
    chkprd :=  (r.lgn=a.lgn) and (r.col=b.col)
  else
    chkprd := dim(r,a.lgn,b.col)
end;

// r = a*b matrix product dim of a and b are checked
// r is allocated and dim. as needed.
// r can be a slice
function prd(var r: mtx; a, b: mtx):boolean;

var
  m, n, p: longint;
  i, j, k: longint;
  ia, ib, ir: longint;
  ia1, ib1, ir1: longint;
  t0: double;

begin
  prd := false;
  assert( chkprd(r,a,b) ); // TODO overload assert if not -Sa option
  m := a.lgn;
  n := a.col;
  p := b.col;
  ia1 := 1;
  ib1 := 1;
  ir1 := 1;
  for i := 1 to m do begin
    ir := ir1;
    for j := 1 to p do begin
      t0 := 0;
      ia := ia1;
      ib := ib1;
      for k := 1 to n do begin
        t0 := t0 + a.e[ia]*b.e[ib];
        ia := ia + a.dl;
        ib := ib + b.dc;
      end;
      //writeln(stderr,i:5,j:5,ir:5);
      r.e[ir] := t0;
      ir := ir + r.dl;
      ib1 := ib1 + b.dl;
    end;
    ia1 := ia1 + a.dc;
    ib1 := 1;
    ir1 := ir1 + r.dc;
  end;
  prd := true
end;

// build a sub matrix (aka slice) to access "a" or die...
// main intended use through lgn() and col()
function smx(a:mtx; i1,j1:longint; lgn, col:longint):mtx;
var k: longint;
begin
  smx.smx := true;
  smx.max := 0; // ->freemem do nothing ? although .smx should be used
  smx.dyn := false;
  smx.lgn := lgn;
  smx.col := col;
  // calc. offset for this alias
  k := 1 + a.dc*(i1-1) + a.dl*(j1-1);
  //writeln(stderr,'k, e[k]: ',k:5,a.e[k]:10);
  smx.e  := @a.e[k] - sizeof(double); // 0 based vectors but only use e[1..max]
  smx.t  := a.t;
  smx.dl := a.dl;
  smx.dc := a.dc;
  smx.nam := ''
end;

// Returns mtx access to one or several adjacent columns of a
function col(a:mtx; j1: longint): mtx;
begin
  col := smx(a, 1, j1, a.lgn, 1)
end;
//
function col(a:mtx; j1, j2: longint): mtx;
begin
  col := smx(a, 1, j1, a.lgn, j2-j1+1)
end;
// Returns the number of columns of a
function col(a:mtx): longint;
begin
  col := a.col;
end;

// map one line of a
function lgn(a:mtx; i1: longint): mtx;
begin
  lgn := smx(a, i1, 1, 1, a.col)
end;
// map several contiguous lines in a sub matrix
function lgn(a:mtx; i1,i2: longint): mtx;
begin
  lgn := smx(a, i1, 1, i2-i1+1, a.col)
end;
// Returns the number of lines of a
function lgn(a:mtx): longint;
begin
  lgn := a.lgn;
end;
// select adjacent lines from "a" assuming a[:,1] sorted by increasing values
// with a[1,i]>low  and a[1,i]<high
// e.g. mwr( lgn(a,200.0,500.0); ACHTUNG use float constants ! "200.0" NOT "200"
// or you will get the wrong variant selecting a slice with line index i1 through i2
function lgn(a:mtx; low, high:double):mtx;
var
  ia, i, i1, i2: longint;
begin
  assert( low<=high );
  ia := 1;
  i1 := 1;
  i2 := 0;
  for i := 1 to a.lgn do begin
    if a.e[ia]<low then begin
        inc(i1);
    end;
    if a.e[ia]<high then begin
        inc(i2)
    end else break;
    ia := ia + a.dc
  end;
  lgn := lgn(a,i1,i2)
end;

// assign the name of a
function nam(var p:mtx; name:string):boolean;
begin
  p.nam := name;
  nam := true
end;
// returns the name of p
function nam(p:mtx):string;
begin
  nam := p.nam;
end;

// generic cpy function r := a
// r can be a plain mtx or or a smx but must be an identifier
// which could be assigned to... NOT a function call... unfortunately
// therefore I keep the functionally equivalent elm(b,a) to cpy
function cpy(var r:mtx; a:mtx):boolean;
var i, j, lgn, col, ia, ir, ia1, ir1: longint;
begin
  cpy := false;
  assert( chk(r,a,a) );
  lgn := a.lgn;
  col := a.col;
  ia1 := 1;
  ir1 := 1;
  for i := 1 to lgn do begin
    ia := ia1;
    ir := ir1;
    for j := 1 to col do begin
       r.e[ir] := a.e[ia];
       ia := ia + a.dl;
       ir := ir + r.dl;
    end;
    ia1 := ia1 + a.dc;
    ir1 := ir1 + r.dc;
  end;
  cpy := true
end;

// advanced feature to use static memory ( alias standard pascal array )
function map(var x:double; lgn, col:longint):mtx;
begin
  map.lgn := lgn;
  map.col := col;
  map.dl := 1;
  map.dc := col;
  map.dyn := false;
  map.max := lgn*col;
  map.smx := false; // TODO TODO either true or false depending on usage !
                    // TODO TODO fix and/or document usage
  map.e := @x - sizeof(double);
end;

function map(var a:mtx; lgn, col:longint):mtx;
begin
  map := map(a.e[1],lgn,col);
  map.smx := true; // TODO check logic of this and see above
  map.nam := a.nam
end;

// we dont care "t" flag because inv(A') = inv(A)' and det(a) = det(a')
// so insitu inversion is made: if t then we leave it like this
// so the inverse is also affected by the transpose flag
// legacy Gauss/Jordan max pivot in situ inversion...
function matinv(a:mtx; var  d: double): boolean;
const
  eps = EpsDouble; // TODO check for a good value ?
var
  n, i, j, k, l: longint;
  dd, amax: extended;
  save: double;
  ik, jk: ^longint;
begin
  matinv := false;
// only plain matrices are usable not a small part of a larger one
// so the col or lgn step must be a.col depending on the transpose
// flag... also inverse of transpose is transpose of inverse so
// we just keep the transpose flag as it is...
  if a.t then begin
    assert( ( a.dc = 1 ) and ( a.dl = a.col) )
  end else begin
    assert( ( a.dl = 1 ) and ( a.dc = a.col) )
  end;
  dd := 1.0;
  n := a.lgn;
  assert(a.col=n);
  GetMem(ik,(n+1)*sizeof(longint));
  GetMem(jk,(n+1)*sizeof(longint));
  for k := 1 to n do begin // k=1 to n
    amax := 0.0;
    for i := k to n do
      for j := k to n do begin
        save := a.e[(i-1)*n+j];
        if abs(save)>abs(amax) then begin
          amax := save; ik[k] := i; jk[k] := j
        end
      end;
    if abs(amax)<eps then begin
      d := 0;
      FreeMem(jk);
      FreeMem(ik);
      exit;
    end;
    i := ik[k];
    if i<>k then for j := 1 to n do begin
      save := a.e[(k-1)*n+j];
      a.e[(k-1)*n+j]:=a.e[(i-1)*n+j];
      a.e[(i-1)*n+j] := -save
    end;
    j := jk[k];
    if j<>k then for i := 1 to n do begin
      save := a.e[(i-1)*n+k];
      a.e[(i-1)*n+k]:=a.e[(i-1)*n+j];
      a.e[(i-1)*n+j]:=-save
    end;
    for i := 1 to n do
      if i<>k then
        a.e[(i-1)*n+k] := -a.e[(i-1)*n+k]/amax;
    for i := 1 to n do
      for j := 1 to n do
        if (i<>k)and(j<>k)then
          a.e[(i-1)*n+j] := a.e[(i-1)*n+j]+a.e[(i-1)*n+k]*a.e[(k-1)*n+j];
    for j := 1 to n do
      if j<>k then a.e[(k-1)*n+j]:=a.e[(k-1)*n+j]/amax;
    a.e[(k-1)*n+k] := 1.0/amax;
    dd := dd*amax
  end; // k=1 to n
  for l := 1 to n do begin // l=1 to n
    k := n-l+1;
    j := ik[k];
    if j>k then for i := 1 to n do begin
      save := a.e[(i-1)*n+k];
      a.e[(i-1)*n+k] := -a.e[(i-1)*n+j];
      a.e[(i-1)*n+j] := save
    end;
    i := jk[k];
    if i>k then for j := 1 to n do begin
      save := a.e[(k-1)*n+j];
      a.e[(k-1)*n+j] := -a.e[(i-1)*n+j];
      a.e[(i-1)*n+j] := save
    end
  end; // l=1 to n
  FreeMem(jk);
  FreeMem(ik);
  d := double(dd);
  matinv := true;
end;

// Gauss Jordan max. pivot  inversion after copy
// to r. Therefore a can be a sub matrix obtained with smx()
// and its sibblings col() and lgn()
function inv(var r:mtx; a:mtx; var d:double):boolean;
begin
  cpy(r,a);
  assert( not r.smx );
  inv := matinv(r,d);
end;

function inv(var r:mtx; a:mtx):boolean;
var dummy: double;
begin
  cpy(r,a);
  assert( not r.smx );
  inv := matinv(r,dummy);
end;

function inv(r:mtx):boolean;
var d: double;
begin
  assert( not r.smx );
  inv := matinv(r,d);
end;

function inv(r:mtx; var d:double):boolean;
begin
  assert( not r.smx );
  inv := matinv(r,d);
end;

// memory allocation (if needed) and dimensions of r
function dim(var r:mtx; lgn, col:longint):boolean;
begin
  dim := false;
  assert( mem(r,lgn*col) );
  assert( not r.smx );
  r.lgn := lgn;
  r.col := col;
  r.smx := false;
  r.dl := 1;
  r.dc := col;
  r.t := false;
  dim := true
end;

// dim and init. same purpose as scilab's ones and zeros...
function dim(var r:mtx; lgn, col:longint; x:double):boolean;
var k: longint;
begin
  dim := dim(r,lgn,col); // TODO check errors trapped in dim()
  for k := 1 to lgn*col do r.e[k] := x
end;

function ptr(a:mtx):pointer;
begin
  assert(a.e<>NIL);
  ptr := @a.e[1]
end;

// modified in the same way as cpy to allow smx for result
// but in this case check proper dim.
function mul(var r: mtx; a:mtx; b:mtx):boolean;
var
  i, j, ia, ia1, ib, ib1, ir, ir1: longint;
begin
  mul := false;
  assert( chk(r,a,b) ); // sanity check and mem. alloc. if needed
  ia1 := 1;
  ib1 := 1;
  ir1 := 1;
  for i := 1 to a.lgn do begin
    ia := ia1;
    ib := ib1;
    ir := ir1;
    for j := 1 to a.col do begin
       r.e[ir] := b.e[ib] * a.e[ia];
       ia := ia + a.dl;
       ib := ib + b.dl;
       ir := ir + r.dl
    end;
    ia1 := ia1 + a.dc;
    ib1 := ib1 + b.dc;
    ir1 := ir1 + r.dc
  end;
  mul := true
end;

function mul(a: mtx; x: double):boolean;
var
  i, j, ia, ia1: longint;
begin
  mul := false;
  ia1 := 1;
  for i := 1 to a.lgn do begin
    ia := ia1;
    for j := 1 to a.col do begin
       a.e[ia] := a.e[ia] * x;
       ia := ia + a.dl;
    end;
    ia1 := ia1 + a.dc;
  end;
  mul := true
end;

function dvd(var r: mtx; a:mtx; b:mtx):boolean;
var
  i, j, ia, ia1, ib, ib1, ir, ir1: longint;
begin
  dvd := false;
  assert( chk(r,a,b) ); // TODO check if chk() make sense here ? specially if a and b are slices...
  ia1 := 1;
  ib1 := 1;
  ir1 := 1;
  for i := 1 to a.lgn do begin
    ia := ia1;
    ib := ib1;
    ir := ir1;
    for j := 1 to a.col do begin
       r.e[ir] := a.e[ia] / b.e[ib]; // TODO exception ?
       ia := ia + a.dl;
       ib := ib + b.dl;
       ir := ir + r.dl
    end;
    ia1 := ia1 + a.dc;
    ib1 := ib1 + b.dc;
    ir1 := ir1 + r.dc
  end;
  dvd := true
end;

function dvd(a: mtx; x: double):boolean;
var
  i, j, ia, ia1: longint;
begin
  dvd := false;
  ia1 := 1;
  for i := 1 to a.lgn do begin
    ia := ia1;
    for j := 1 to a.col do begin
       a.e[ia] := a.e[ia] / x;
       ia := ia + a.dl;
    end;
    ia1 := ia1 + a.dc;
  end;
  dvd := true
end;

function add(var r: mtx; a:mtx; b:mtx):boolean;
var
  i, j, ia, ia1, ib, ib1, ir, ir1: longint;
begin
  add := false;
  assert( chk(r,a,b) ); // sanity check and mem. alloc. if needed
  ia1 := 1;
  ib1 := 1;
  ir1 := 1;
  for i := 1 to a.lgn do begin
    ia := ia1;
    ib := ib1;
    ir := ir1;
    for j := 1 to a.col do begin
       r.e[ir] := a.e[ia] + b.e[ib];
       ia := ia + a.dl;
       ib := ib + b.dl;
       ir := ir + r.dl
    end;
    ia1 := ia1 + a.dc;
    ib1 := ib1 + b.dc;
    ir1 := ir1 + r.dc
  end;
  add := true
end;

function add(a: mtx; x: double):boolean;
var
  i, j, ia, ia1: longint;
begin
  add := false;
  ia1 := 1;
  for i := 1 to a.lgn do begin
    ia := ia1;
    for j := 1 to a.col do begin
       a.e[ia] := a.e[ia] + x;
       ia := ia + a.dl;
    end;
    ia1 := ia1 + a.dc;
  end;
  add := true
end;

// check mem/dim for reduction by lgn or col ( mostly summation )
function chkred(var r:mtx; a:mtx; t:boolean):boolean;
var
  lgn, col: longint;
begin
  chkred := false;
  if t then begin
    lgn := 1;
    col := a.col;
  end else begin
    col := 1;
    lgn := a.lgn
  end;
  if r.smx then begin
    assert( (r.lgn=lgn) and (r.col=col) )
  end else begin
    assert( dim(r,lgn,col) )
  end;
  chkred := true
end;

function add(var r: mtx; a:mtx; t:boolean):boolean;
var
  i, j, ia, ia1, ir, ir1: longint;
  x: double;
  di, dj, dr, i0, j0: longint;
begin
  add := false;
  chkred(r,a,t); // sanity check and mem. alloc. if needed
//
  if t then begin
    i0 := col(a);
    j0 := lgn(a);
    di := a.dl;
    dj := a.dc;
    dr := r.dl;
  end else begin
    i0 := lgn(a);
    j0 := col(a);
    di := a.dc;
    dj := a.dl;
    dr := r.dc;
  end;
  ia1 := 1;
  ir1 := 1;
  for i := 1 to i0 do begin
    ia := ia1;
    ir := ir1;
    x := 0.0;
    for j := 1 to j0 do begin
      x := x + a.e[ia];
      ia := ia + dj
    end;
    r.e[ir] := x;
    ia1 := ia1 + di;
    ir1 := ir1 + dr
  end;
  add := true
end;

function sub(var r: mtx; a:mtx; b:mtx):boolean;
var
  i, j, ia, ia1, ib, ib1, ir, ir1: longint;
begin
  sub := false;
  assert( chk(r,a,b) ); // sanity check and mem. alloc. if needed
  ia1 := 1;
  ib1 := 1;
  ir1 := 1;
  for i := 1 to a.lgn do begin
    ia := ia1;
    ib := ib1;
    ir := ir1;
    for j := 1 to a.col do begin
       r.e[ir] := a.e[ia] - b.e[ib];
       ia := ia + a.dl;
       ib := ib + b.dl;
       ir := ir + r.dl
    end;
    ia1 := ia1 + a.dc;
    ib1 := ib1 + b.dc;
    ir1 := ir1 + r.dc
  end;
  sub := true
end;

function sub(a: mtx; x: double):boolean;
var
  i, j, ia, ia1: longint;
begin
  sub := false;
  ia1 := 1;
  for i := 1 to a.lgn do begin
    ia := ia1;
    for j := 1 to a.col do begin
       a.e[ia] := a.e[ia] - x;
       ia := ia + a.dl;
    end;
    ia1 := ia1 + a.dc;
  end;
  sub := true
end;

// If "r" is a submatrix ( field smx true ? TODO quid de trp() ) then
// returns true if it has the same dim as "a" else false
// if "r" is a plain matrix
// redim and realloc dyn. mem if needed...
// TODO use this function in all r = a oper b functions so that sub matrices
// can be used for r. In some cases r must be different from the operands in
// other cases it is not needed... this should be taken care of internally
// using temporary copies if needed...
function chk(var r:mtx; a, b:mtx):boolean;
begin
  chk := false;
  if (a.lgn<>b.lgn) or (a.col<>b.col) then exit;
  if r.e = NIL then r.smx := false;
  if (r.smx) then
    chk :=  (r.lgn=a.lgn) and (r.col=a.col)
  else
    chk := dim(r,a.lgn,a.col)
end;

{
// Try to make a nice wrapper for spline_ and seval_
function spl( var s:mtx; x, y: mtx):boolean;
var
  sx, sy, b, c, d: mtx;
  n: longint;
begin
  spl := false;
  assert( (x.lgn = y.lgn) and (x.col=y.col) );
  assert( (x.lgn=1) or (x.col=1) );
  assert( x.lgn<>x.col );
  if col(x)<lgn(x) then begin
    x := trp(x);
    y := trp(y);
  end;
  // sanity check of arguments... lgn or col vectors should be ok
  // so may be transpose if col vector..
  // sev() should also handle gracefully line/col vectors TODO
  // m := lgn(x);
  n := col(x);
  //assert( ( col(y)=n ) and (n>1)  );
  //assert( (lgn(y)=1) and (m=1) );
  // build s matrix with proper dim and copy x and y in first 2 lines
  dim(s, 5, n);
  elm(lgn(s, 1),x);
  elm(lgn(s, 2),y);
  // line 1 and 2 are copy of x and y and contiguous
  // so x and y arg. can be slices TODO check this
  // lines 3 4 5 will hold b, c and d coeff.
  sx := lgn(s,1);
  sy := lgn(s,2);
  b := lgn(s,3);
  c := lgn(s,4);
  d := lgn(s,5);
  spline_(n, ptr(sx), ptr(sy), ptr(b), ptr(c), ptr(d) );
  spl := true
end;

function sev(u:double; s:mtx):double;
var
  x, y, b, c, d: pointer;
  n: longint;
begin
  assert( (lgn(s)=5) and (col(s)>1 ) );
  x := ptr( lgn(s,1) );
  y := ptr( lgn(s,2) );
  b := ptr( lgn(s,3) );
  c := ptr( lgn(s,4) );
  d := ptr( lgn(s,5) );
  n := col(s);
  sev := seval_(n,u,x,y,b,c,d);
end;

//
function sev(var r:mtx; u:mtx; s:mtx):boolean;
var
  x, y, b, c, d: pointer;
  n: longint;
  k: longint;
  uk: double;
begin
  sev := false;
  assert( not u.smx ); // interim u is plain mtx
  assert( (lgn(s)=5) and (col(s)>1 ) );
  x := ptr( lgn(s,1) );
  y := ptr( lgn(s,2) );
  b := ptr( lgn(s,3) );
  c := ptr( lgn(s,4) );
  d := ptr( lgn(s,5) );
  n := col(s);
  assert(dim(r,lgn(u),col(u),0.0));
  for k := 1 to col(u)*lgn(u) do begin
    uk := u.e[k];
    r.e[k] := seval_(n,uk,x,y,b,c,d);
  end;
  sev := true
end;
}

// scilab linspace... TODO variants (mtx,double,double,longint) et (double,double)
// TODO TODO r should be zer() as it is a var arg...
function lsp(var r:mtx; x0, dx, x1: double):boolean;
var
  k, n: longint;
begin
  lsp := false;
  n := trunc( (x1-x0)/dx );
  zer(r,n+1,1); // enforces a col. vect.
  assert(n>0);
  assert( dim(r,1,n+1) );
  assert( (x0 + n*dx) <= x1 );
  for k := 0 to n do begin
    r.e[k+1] := x0+k*dx;
  end;
  r.lgn := 1;
  r.col := n+1;
  r.nam := '';
  lsp := true
end;

// scilab linspace... TODO variants (mtx,double,double,longint) et (double,double)
function lsp(var r:mtx; x0, x1: double; n:longint):boolean;
var
  k: longint;
begin
  lsp := false;
  zer(r,n+1,1);
  for k := 0 to n do begin
    r.e[k+1] := (x0*(n-k)+x1*k)/n;
  end;
  r.lgn := n+1;
  r.col := 1;
  r.nam := '';
  lsp := true
end;

// interim: r is plain... may be allow r smx
// hence check that size is ok...
// interim assert r plain created...
// check if ok with smx and transposed v
function dia(var r:mtx; v: mtx):boolean;
var
ir, iw, i, n: longint;
w: mtx;
begin
dia := false;
assert( not r.smx ); // interim
if lgn(v)=1 then
w := v
else if col(v)=1 then
w := trp(v)
else
assert(false);
// writeln('dia...');
n := col(w);
dim(r,n,n,0.0);
ir := 1;
iw := 1;
for i := 1 to n do begin
r.e[ir] := w.e[iw];
ir := ir + r.dl + r.dc;
iw := iw + w.dl;
end;
// mpt(r);
// writeln('...dia');
writeln;
dia := true
end;

// returns a line vector smx
function dia(a:mtx):mtx;
begin
assert( a.col=a.lgn );
assert( not a.smx );
dia := a;
dia.col := a.col;
dia.lgn := 1;
dia.dc := a.dc+a.dl;
dia.dl := a.dc+a.dl;
end;

// TODO corr. coeff.
procedure lin(x,y: mtx; var slope, y0, s2y: double);
var
  xi, yi, sumx, sumy, sumxy, sumx2 : double;
  i, m, n: longint;
  u, v: mtx;
begin
  //mpt(x,' x ');
  //mwr(x);
  //mpt(y,' y ');
  //mwr(y);
  m := lgn(x);
  n := col(x);
  //writeln('m:',m,'  n:',n);
  assert( (m>2) or (n>2) );
  assert( (lgn(y)=m) and (col(y)=n) );
  if m=1 then begin
    u := x;
    v := y;
  end else begin
    assert( n=1 );
    u := trp(x);
    v := trp(y);
    n := m;
    m := 1; // not needed
    //mpt(u,' u '); mwr(u);
    //mpt(v,' v '); mwr(v);
  end;
// here u and v must be row vectors
  sumx2 := 0;
  sumx := 0;
  sumxy := 0;
  sumy := 0;
  for i := 1 to n do begin
    xi := elm(u,1,i);
    yi := elm(v,1,i);
    //writeln('i=',i,' xi=',xi:0:3,' yi=',yi:0:3);
    sumx := sumx + xi;
    sumy := sumy + yi;
    sumx2 := sumx2 + sqr(xi);
    sumxy := sumxy + xi*yi;
  end;
  //writeln('sumx:',sumx);
  //writeln('sumy:',sumy);
  //writeln('sumx2:',sumx2);
  //writeln('sumxy:',sumxy);
  slope := (n*sumxy-sumx*sumy)/(n*sumx2-sqr(sumx));
  y0 := (sumx2*sumy-sumx*sumxy)/(n*sumx2-sqr(sumx));
  s2y := 0;
  for i := 1 to n do begin
    s2y := s2y + sqr(yi-y0-xi*slope);
  end;
  s2y := s2y/(n-2)
end;

procedure lin(xy: mtx; var slope, y0, s2y: double);
begin
  // we must know which dim is larger than 2 to use either rows or cols
  if ( col(xy)>2 ) then begin
    assert(lgn(xy)=2);
    lin(lgn(xy,1),lgn(xy,2),slope,y0,s2y);
  end else begin
    assert(col(xy)=2);
    lin(col(xy,1),col(xy,2),slope,y0,s2y);
  end
end;

function max(a:mtx):double;
var
  i: longint;
  r, s: double;
  b: mtx;
begin
  r := elm(a,1,1);
  if lgn(a)>col(a) then b := trp(a) else b := a;
  assert( lgn(b)=1 );
  r := elm(b,1,1);
  for i := 2 to col(b) do begin
    s := elm(b,1,i);
    if s>r then r:=s
  end;
  max := r
end;

function min(a:mtx):double;
var
  i: longint;
  r, s: double;
  b: mtx;
begin
  r := elm(a,1,1);
  if lgn(a)>col(a) then b := trp(a) else b := a;
  assert( lgn(b)=1 );
  r := elm(b,1,1);
  for i := 2 to col(b) do begin
    s := elm(b,1,i);
    if s<r then r:=s
  end;
  min := r
end;

function inrois(x:double; r:mtx):boolean;
var i: longint;
begin
  assert( col(r)=2 );
  inrois:= false;
  for i := 1 to lgn(r) do begin
    if (x>elm(r,i,1)) and (x<=elm(r,i,2)) then begin
      inrois := true;
      break
    end
  end
end;

// derived from roi3.pas unit...
procedure lgn(var a:mtx; r:mtx);
var i, j: longint;
begin
  assert( not a.smx ); // TODO interim before test...
  j := 0;
  for i := 1 to lgn(a) do begin
    if inrois(elm(a,i,1),r) then begin
      inc(j);
      if i<>j then elm(lgn(a,j),lgn(a,i));
    end
  end;
  dim(a,j,col(a))
end;

function mxs(var r:mtx; s:string): boolean;
var
  m: longint = 0;
  n: longint = 0;
  s0, s1: string;
  i, j: longint;
  x: double;
  code: integer;
begin
  mxs := false;
  m := WordCount(s,RowDelim);
  if m=0 then exit;
  for i := 1 to m do begin
    s0 := ExtractWord(i,s,RowDelim);
    if n=0 then begin
      // mem. alloc. and dim.
      n := WordCount(s0,ElmDelim);
      assert(n>0);
      assert( zer(r,m,n) );
    end else begin
      // sanity check
      assert(WordCount(s0,ElmDelim)=n);
    end;
    for j := 1 to n do begin
      s1 := ExtractWord(j,s0,ElmDelim);
      //write(' ',s1);
      val(s1,x,code);
      assert(code=0);
      elm(r,i,j,x);
    end;
    //writeln;
  end;
  //writeln;
  mxs := true
end;

// Area below xy curve assert sorted by x
function axy(a:mtx; kx, ky:longint): double;
var
  x1, y1, x2, y2: double;
  i: longint;
begin
  axy := 0.0;
  assert( lgn(a)>1 );
  assert( (kx>0) and (kx<=col(a)) );
  assert( (ky>0) and (ky<=col(a)) );
  x1 := elm(a,1,kx);
  y1 := elm(a,1,ky);
  for i := 2 to lgn(a) do begin
    x2 := elm(a,i,kx);
    y2 := elm(a,i,ky);
    axy := axy + (y1+y2)*(x2-x1);
    x1 := x2;
    y1 := y2
  end;
  axy := axy / 2.0
end;

// idem def. col 1 and 2
function axy(a:mtx):double;
begin
  axy := axy(a,1,2)
end;

// mean and stdev of a vector
// TODO mtx variants col or row wise...
function avg(a:mtx; var dst:double):double;
var
  a0: mtx;
  j, m: longint;
  sd2: double;
  sx1: double;
  x, x1, sqrx1: double;
  variance: double;
begin
  assert( (lgn(a)=1) or (col(a)=1) );
  if lgn(a)=1 then a0 := trp(a) else a0 := a;
  // now a0 is a column vector
  m := lgn(a0);
  assert( m>1 );
  x1 := elm(a0,1,1);
  sqrx1 := x1*x1;
  sd2 := 0.0;
  sx1 := 0.0;
  avg := x1;
  for j := 2 to m do begin
    x := elm(a0,j,1);
    sx1 := sx1 + x - x1;
    avg := x1 + sx1/j;
    sd2 := sd2 + x*x + sqrx1 - 2*x*x1;
  end;
  variance := sd2 - sx1*sx1/j;
  dst := sqrt(variance/(m-1))
end;

begin
  //writeln(stderr, '$Id: mtx3.pas,v 1.37 2008/04/11 11:33:27 jroux Exp $');
end.
