{
    This file is part of the Numlib package.
    Copyright (c) 1986-2000 by
     Kees van Ginneken, Wil Kortsmit and Loek van Reij of the
     Computational centre of the Eindhoven University of Technology

    FPC port Code          by Marco van de Voort (marco@freepascal.org)
             documentation by Michael van Canneyt (Michael@freepascal.org)

    Interpolate and (curve) fitting.
    Slegpb in this unit patched parameters slightly. Units IPF and sle
    were not in the same revision in this numlib copy (which was a
    copy of the work directory of the author) .

    Contains two undocumented functions. If you recognize the algoritm,
    mail us.

    See the file COPYING.FPC, included in this distribution,
    for details about the copyright.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

 **********************************************************************}

{
 }
{$IFNDEF FPC_DOTTEDUNITS}
unit ipf;
{$ENDIF FPC_DOTTEDUNITS}
{$modeswitch exceptions}
{$I direct.inc}
interface

{$IFDEF FPC_DOTTEDUNITS}
uses NumLib.Typ, NumLib.Mdt, NumLib.Dsl, NumLib.Sle, NumLib.Spe;
{$ELSE FPC_DOTTEDUNITS}
uses typ, mdt, dsl, sle, spe;
{$ENDIF FPC_DOTTEDUNITS}

type
  THermiteSplineType = (
    hstMonotone // preserves monotonicity of the interpolated function by using
                // a Fritsch-Carlson algorithm
  );

{ Determine natural cubic spline "s" for data set (x,y), output to (a,d2a)
 term=1 success,
     =2 failure calculating "s"
     =3 wrong input (e.g. x,y is not sorted increasing on x)}
procedure ipffsn(n: ArbInt; var x, y, a, d2a: ArbFloat; var term: ArbInt);

{calculate d2s from x,y, which can be used to calculate s}
procedure ipfisn(n: ArbInt; var x, y, d2s: ArbFloat; var term: ArbInt);

{Calculate function value for dataset (x,y), with n.c. spline d2s for
x value t. Return (corrected) y value.
s calculated from x,y, with e.g. ipfisn}
function  ipfspn(n: ArbInt; var x, y, d2s: ArbFloat; t: ArbFloat;
                 var term: ArbInt): ArbFloat;

{Calculate minimum and maximum values for the n.c. spline d2s.
Does NOT take source points into account.}
procedure ipfsmm(n: ArbInt; var x, y, d2s, minv, maxv: ArbFloat; 
        var term: ArbInt);

{Calculates tangents for each data point (d1s), for a given array of input data
 points (x,y), by using a selected variant of a Hermite cubic spline interpolation.
 Inputs:
   hst - algorithm selection
   n - highest array index
   x[0..n] - array of X values (one value for each data point)
   y[0..n] - array of Y values (one value for each data point)
 Outputs:
   d1s[0..n] - array of tangent values (one value for each data point)
   term - status: 1 if function succeeded, 3 if less than two data points given
}
procedure ipfish(hst: THermiteSplineType; n: ArbInt; var x, y, d1s: ArbFloat; var term: ArbInt);

{Calculates interpolated function value for a given array of input data points
 (x,y) and tangents for each data point (d1s), for input value t, by using a
 Hermite cubic spline interpolation; d1s array can be obtained by calling the
 ipfish procedure.
 Inputs:
   n - highest array index
   x[0..n] - array of X values (one value for each data point)
   y[0..n] - array of Y values (one value for each data point)
   d1s[0..n] - array of tangent values (one value for each data point)
   t - input value X
 Outputs:
   term - status: 1 if function succeeded, 3 if less than two data points given
   result - interpolated function value Y
}
function ipfsph(n: ArbInt; var x, y, d1s: ArbFloat; t: ArbFloat; var term: ArbInt): ArbFloat;

{Calculate n-degree polynomal b for dataset (x,y) with m elements
 using the least squares method.}
procedure ipfpol(m, n: ArbInt; var x, y, b: ArbFloat; var term: ArbInt);


                {**** undocumented ****}

function spline(    n     : ArbInt;
                    x     : complex;
                var ac    : complex;
                var gammar: ArbFloat;
                    u1    : ArbFloat;
                    pf    : complex): ArbFloat;

                {**** undocumented ****}

procedure splineparameters
          (    n                 : ArbInt;
           var ac, alfadc        : complex;
           var lambda,
               gammar, u1,
               kwsom, energie    : ArbFloat;
           var pf                : complex);

implementation


procedure ipffsn(n: ArbInt; var x, y, a, d2a: ArbFloat; var term: ArbInt);

var                       i, sr, n1s, ns1, ns2: ArbInt;
   s, lam, lam0, lam1, lambda, ey, ca, p, q, r: ArbFloat;
     px, py, pd, pa, pd2a,
  h, z, diagb, dinv, qty, qtdinvq, c, t, tl: ^arfloat1;
                                         ub: boolean;

  procedure solve; {n, py, qty, h, qtdinvq, dinv, lam, t, pa, pd2a, term}
  var i: ArbInt;
          p, q, r: ArbFloat;
             f, c: ^arfloat1;
               ca: ArbFloat = 0.0;
  begin
    getmem(f, 3*ns1); getmem(c, ns1);
    for i:=1 to n-1 do
      begin
        f^[3*i]:=qtdinvq^[3*i]+lam*t^[2*i];
        if i > 1
        then
          f^[3*i-1]:=qtdinvq^[3*i-1]+lam*t^[2*i-1];
        if i > 2
        then
          f^[3*i-2]:=qtdinvq^[3*i-2];
        if lam=0
        then
          c^[i]:=qty^[i]
        else
          c^[i]:=lam*qty^[i]
      end;
    slegpb(n-1, 2,{ 3,} f^[1], c^[1], pd2a^[1], ca, term);
    if term=2
    then
      begin
        freemem(f, 3*ns1); freemem(c, ns1);
        exit
      end;
    p:=1/h^[1];
    if lam=0
    then
      r:=1
    else
      r:=1/lam;
    q:=1/h^[2]; pa^[1]:=py^[1]-r*dinv^[1]*p*pd2a^[1];
    pa^[2]:=py^[2]-r*dinv^[2]*(pd2a^[2]*q-(p+q)*pd2a^[1]); p:=q;
    for i:=3 to n-1 do
      begin
        q:=1/h^[i];
        pa^[i]:=py^[i]-r*dinv^[i]*
                (p*pd2a^[i-2]-(p+q)*pd2a^[i-1]+q*pd2a^[i]);
        p:=q
      end;
    q:=1/h^[n];
    pa^[n]:=py^[n]-r*dinv^[n]*(p*pd2a^[n-2]-(p+q)*pd2a^[n-1]);
    pa^[n+1]:=py^[n+1]-r*dinv^[n+1]*q*pd2a^[n-1];
    if lam=0
    then
      for i:=1 to n-1 do
        pd2a^[i]:=0;
    freemem(f, 3*ns1); freemem(c, ns1);
  end; {solve}

    function e(var c, h: ArbFloat; n:ArbInt): ArbFloat;
    var i:ArbInt;
        s:ArbFloat;
        pc, ph: ^arfloat1;
    begin
      ph:=@h; pc:=@c;
      s:=ph^[1]*pc^[1]*pc^[1];
      for i:=1 to n-2 do
        s:=s+(pc^[i]*(pc^[i]+pc^[i+1])+pc^[i+1]*pc^[i+1])*ph^[i+1];
      e:=(s+pc^[n-1]*pc^[n-1]*ph^[n])/3
    end; {e}

    function cr(lambda: ArbFloat): ArbFloat;
    var s, crs: ArbFloat;
             i: ArbInt;
    begin
      cr:=0; lam:=lambda;
      solve; { n, py, qty, h, qtdinvq, dinv, lam, t, pa, pd2a, term }
      if term=2
      then
        exit;
      crs:=ey;
      if lam <> 0
      then
        begin
          crs:=crs+e(pd2a^[1], h^[1], n);
          s:=0;
          for i:=1 to n-1 do
            s:=s+pd2a^[i]*qty^[i];
          crs:=crs-2*s
        end;
      s:=0;
      for i:=1 to n+1 do
        s:=s+sqr(pa^[i]-py^[i])*diagb^[i];
      cr:=crs-s
    end; {cr}

    procedure roof1r(a, b, ae, re: ArbFloat; var x: ArbFloat);

    var fa, fb, c, fc, m, tol, w1, w2 : ArbFloat;
                                    k : ArbInt;
                                 stop : boolean;

    begin
      fa:=cr(a);
      if term=2
      then
        exit;
      fb:=cr(b);
      if term=2
      then
        exit;
      if abs(fb)>abs(fa)
      then
        begin
          c:=b; fc:=fb; x:=a; b:=a; fb:=fa; a:=c; fa:=fc
        end
      else
        begin
          c:=a; fc:=fa; x:=b
        end;
      k:=0;
      tol:=ae+re*spemax(abs(a), abs(b));
      w1:=abs(b-a); stop:=false;
      while (abs(b-a)>tol) and (fb<>0) and (not stop) do
        begin
          m:=(a+b)/2;
          if (k>=2) or (fb=fc)
          then
            x:=m
          else
            begin
              x:=(b*fc-c*fb)/(fc-fb);
              if abs(b-x)<tol
              then
                x:=b-tol*spesgn(b-a);
              if spesgn(x-m)=spesgn(x-b)
              then
                x:=m
            end;
          c:=b; fc:=fb; b:=x; fb:=cr(x);
          if term=2
          then
            exit;
          if spesgn(fa)*spesgn(fb)>0
          then
            begin
              a:=c; fa:=fc; k:=0
            end
          else
            k:=k+1;
          if abs(fb)>=abs(fa)
          then
            begin
              c:=b; fc:=fb; x:=a; b:=a; fb:=fa; a:=c; fa:=fc; k:=0
            end;
          tol:=ae+re*spemax(abs(a), abs(b));
          w2:=abs(b-a);
          if w2>=w1
          then
            stop:=true;
          w1:=w2
        end
    end; {roof1r}

procedure NoodGreep;
var I, j: ArbInt;
begin
  i:=1;
  while i <= n do
    begin
      if (pd^[i] <= 0) or (px^[i+1] <= px^[i])
      then
        begin
          term:=3;
          exit
        end;
      i:=i+1
    end;
  if pd^[n+1] <= 0
  then
    begin
      term:=3;
      exit
    end;
  for i:=1 to n+1 do
    dinv^[i]:=1/pd^[i];
  for i:=1 to n do
    h^[i]:=px^[i+1]-px^[i];
  t^[2]:=(h^[1]+h^[2])/3;
  for i:=2 to n-1 do
    begin
      t^[2*i]:=(h^[i]+h^[i+1])/3; t^[2*i-1]:=h^[i]/6
    end;
  move(t^[1], tl^[1], ns2);
  mdtgpb(n-1, 1, 2, tl^[1], ca, term);
  if term=2
  then
    exit;
  z^[1]:=1/(h^[1]*tl^[2]);
  for j:=2 to n-1 do
    z^[j]:=-(tl^[2*j-1]*z^[j-1])/tl^[2*j];
  s:=0;
  for j:=1 to n-1 do
    s:=s+sqr(z^[j]);
  diagb^[1]:=s;
  z^[1]:=(-1/h^[1]-1/h^[2])/tl^[2];
  if n>2
  then
    z^[2]:=(1/h^[2]-tl^[3]*z^[1])/tl^[4];
  for j:=3 to n-1 do
    z^[j]:=-tl^[2*j-1]*z^[j-1]/tl^[2*j];
  s:=0;
  for j:=1 to n-1 do
    s:=s+sqr(z^[j]);
  diagb^[2]:=s;
  for i:=2 to n-2 do
    begin
      z^[i-1]:=1/(h^[i]*tl^[2*(i-1)]);
      z^[i]:=(-1/h^[i]-1/h^[i+1]-tl^[2*i-1]*z^[i-1])/tl^[2*i];
      z^[i+1]:=(1/h^[i+1]-tl^[2*i+1]*z^[i])/tl^[2*(i+1)];
      for j:=i+2 to n-1 do
        z^[j]:=-tl^[2*j-1]*z^[j-1]/tl^[2*j];
      s:=0;
      for j:=i-1 to n-1 do
        s:=s+sqr(z^[j]);
      diagb^[i+1]:=s
    end;
  z^[n-2]:=1/(h^[n-1]*tl^[2*(n-2)]);
  z^[n-1]:=(-1/h^[n-1]-1/h^[n]-tl^[2*n-3]*z^[n-2])/tl^[2*(n-1)];
  s:=0;
  for j:=n-2 to n-1 do
    s:=s+sqr(z^[j]);
  diagb^[n]:=s;
  diagb^[n+1]:=1/sqr(h^[n]*tl^[2*(n-1)]);
  p:=1/h^[1];
  for i:=2 to n do
    begin
      q:=1/h^[i]; qty^[i-1]:=py^[i+1]*q-py^[i]*(p+q)+py^[i-1]*p;
      p:=q
    end;
  p:=1/h^[1]; q:=1/h^[2]; r:=1/h^[3];
  qtdinvq^[3]:=dinv^[1]*p*p+dinv^[2]*(p+q)*(p+q)+dinv^[3]*q*q;
  if n>2
  then
    begin
      qtdinvq^[6]:=dinv^[2]*q*q+dinv^[3]*(q+r)*(q+r)+dinv^[4]*r*r;
      qtdinvq^[5]:=-(dinv^[2]*(p+q)+dinv^[3]*(q+r))*q;
      p:=q; q:=r;
      for i:=3 to n-1 do
        begin
          r:=1/h^[i+1];
          qtdinvq^[3*i]:=dinv^[i]*q*q+dinv^[i+1]*(q+r)*(q+r)+dinv^[i+2]*r*r;
          qtdinvq^[3*i-1]:=-(dinv^[i]*(p+q)+dinv^[i+1]*(q+r))*q;
          qtdinvq^[3*i-2]:=dinv^[i]*p*q;
          p:=q; q:=r
        end
    end;
  dslgpb(n-1, 1, 2, tl^[1], qty^[1], c^[1], term);
  if term=2
  then
    exit;
  ey:=e(c^[1], h^[1], n);
  lam0:=0;
  s:=cr(lam0);
  if term=2
  then
    exit;
  if s >= 0
  then
    begin
      lambda:=0; term:=4
    end
  else
    begin
      lam1:=1e-8; ub:=false;
      while (not ub) and (lam1<=1.1e8) do
        begin
          s:=cr(lam1);
          if term=2
          then
            exit;
          if s  >= 0
          then
            ub:=true
          else
            begin
              lam0:=lam1; lam1:=10*lam1
            end
        end;
      if not ub
      then
        begin
          term:=4; lambda:=lam0
        end
      else
        roof1r(lam0, lam1, 0, 1e-6, lambda);
      if term=2
      then
        exit
    end;

end;

begin
  term:=1;
  if n < 2
  then
    begin
      term:=3; exit
    end;
  sr:=sizeof(ArbFloat);
  n1s:=(n+1)*sr;
  ns2:=2*(n-1)*sr;
  ns1:=(n-1)*sr;
  getmem(dinv, n1s);
  getmem(h, n*sr);
  getmem(t, ns2);
  getmem(tl, ns2);
  getmem(z, ns1);
  getmem(diagb, n1s);
  getmem(qtdinvq, 3*ns1);
  getmem(c, ns1);
  getmem(qty, ns1);

   getmem(pd, n1s);
 { pd:=@d; }
  px:=@x;
  py:=@y;
  pa:=@a;
  pd2a:=@d2a;
  { de gewichten van de punten worden op 1 gezet}
  for i:=1 to n+1 do
    pd^[i]:=1;

  NoodGreep;

  freemem(dinv, n1s);
  freemem(h, n*sr);
  freemem(t, ns2);
  freemem(tl, ns2);
  freemem(z, ns1);
  freemem(diagb, n1s);
  freemem(qtdinvq, 3*ns1);
  freemem(c, ns1);
  freemem(qty, ns1);

  freemem(pd, n1s);
end; {ipffsn}

procedure ortpol(m, n: ArbInt; var x, alfa, beta: ArbFloat);
// this function used to use mark/release.
var
                             i, j, ms : ArbInt;
    xppn1, ppn1, ppn, p, alfaj, betaj : ArbFloat;
               px, pal, pbe, pn, pn1 : ^arfloat1;
begin
  px:=@x; pal:=@alfa; pbe:=@beta; ms:=m*sizeof(ArbFloat);
  getmem(pn, ms); getmem(pn1, ms);
  xppn1:=0; ppn1:=m;
  for i:=1 to m do
    begin
      pn^[i]:=0; pn1^[i]:=1; xppn1:=xppn1+px^[i]
    end;
  pal^[1]:=xppn1/ppn1; pbe^[1]:=0;
  for j:=2 to n do
    begin
      alfaj:=pal^[j-1]; betaj:=pbe^[j-1];
      ppn:=ppn1; ppn1:=0; xppn1:=0;
      for i:=1 to m do
        begin
          p:=(px^[i]-alfaj)*pn1^[i]-betaj*pn^[i];
          pn^[i]:=pn1^[i]; pn1^[i]:=p; p:=p*p;
          ppn1:=ppn1+p; xppn1:=xppn1+px^[i]*p
        end; {i}
      pal^[j]:=xppn1/ppn1; pbe^[j]:=ppn1/ppn
    end; {j}
    freemem(pn); freemem(pn1);
end; {ortpol}

procedure ortcoe(m, n: ArbInt; var x, y, alfa, beta, a: ArbFloat);
// this function used to use mark/release.
var                        i, j, mr : ArbInt;
         fpn, ppn, p, alphaj, betaj : ArbFloat;
    px, py, pal, pbe, pa, pn, pn1 : ^arfloat1;

begin
  mr:=m*sizeof(ArbFloat);
  px:=@x; py:=@y; pal:=@alfa; pbe:=@beta; pa:=@a;
  getmem(pn, mr); getmem(pn1, mr);
  fpn:=0;
  for i:=1 to m do
    begin
      pn^[i]:=0; pn1^[i]:=1; fpn:=fpn+py^[i]
    end; {i}
  pa^[1]:=fpn/m;
  for j:=1 to n do
    begin
      fpn:=0; ppn:=0; alphaj:=pal^[j]; betaj:=pbe^[j];
      for i:=1 to m do
        begin
          p:=(px^[i]-alphaj)*pn1^[i]-betaj*pn^[i];
          pn^[i]:=pn1^[i]; pn1^[i]:=p;
          fpn:=fpn+py^[i]*p; ppn:=ppn+p*p
        end; {i}
      pa^[j+1]:=fpn/ppn
    end; {j}
    freemem(pn); freemem(pn1);  
end; {ortcoe}

procedure polcoe(n:ArbInt; var alfa, beta, a, b: ArbFloat);

var            k, j : ArbInt;
           pal, pbe : ^arfloat1;
            pa, pb  : ^arfloat0;

begin
  pal:=@alfa; pbe:=@beta; pa:=@a; pb:=@b;
  move(pa^[0], pb^[0], (n+1)*sizeof(ArbFloat));
  for j:=0 to n-1 do
    for k:=n-j-1 downto 0 do
      begin
        pb^[k+j]:=pb^[k+j]-pal^[k+1]*pb^[k+j+1];
        if k+j<>n-1
        then
          pb^[k+j]:=pb^[k+j]-pbe^[k+2]*pb^[k+j+2]
      end
end; {polcoe}

procedure ipfpol(m, n: ArbInt; var x, y, b: ArbFloat; var term: ArbInt);

var                      i, ns: ArbInt;
                          fsum: ArbFloat;
                py, alfa, beta: ^arfloat1;
                         pb, a: ^arfloat0;
begin
  if (n<0) or (m<1)
  then
    begin
      term:=3; exit
    end;
  term:=1;
  if n = 0
  then
    begin
      py:=@y; pb:=@b;
      fsum:=0;
      for i:=1 to m do
        fsum:=fsum+py^[i];
      pb^[0]:=fsum/m
    end
  else
    begin
      if n>m-1
      then
        begin
          pb:=@b;
          fillchar(pb^[m], (n-m+1)*sizeof(ArbFloat), 0);
          n:=m-1
        end;
      ns:=n*sizeof(ArbFloat);
      getmem(alfa, ns); getmem(beta, ns);
      getmem(a, (n+1)*sizeof(ArbFloat));
      ortpol(m, n, x, alfa^[1], beta^[1]);
      ortcoe(m, n, x, y, alfa^[1], beta^[1], a^[0]);
      polcoe(n, alfa^[1], beta^[1], a^[0], b);
      freemem(alfa, ns); freemem(beta, ns);
      freemem(a, (n+1)*sizeof(ArbFloat));
    end
end; {ipfpol}

procedure ipfisn(n: ArbInt; var x, y, d2s: ArbFloat; var term: ArbInt);

var
                s, i, L : ArbInt;
                   p, q : ArbFloat;
        px, py, h, b, t : ^arfloat0;
                   pd2s : ^arfloat1;
                     ca : ArbFloat = 0.0;
begin
  term:=1;
  if n < 1
  then
    begin
      term:=3; exit
    end; {n<1}
  if n = 1 then
    exit;

  px:=@x; py:=@y; pd2s:=@d2s;
  s:=sizeof(ArbFloat);
  getmem(h, n*s);
  getmem(b, (n-1)*s);
  getmem(t, 2*(n-1)*s);
  for i:=0 to n-1 do
    h^[i]:=px^[i+1]-px^[i];
  q:=1/6; p:=2*q;
  t^[1]:=p*(h^[0]+h^[1]);
  for i:=2 to n-1 do
    begin
      t^[2*i-1]:=p*(h^[i-1]+h^[i]); t^[2*i-2]:=q*h^[i-1]
    end; {i}
  p:=1/h^[0];
  for i:=2 to n do
    begin
      q:=1/h^[i-1]; b^[i-2]:=py^[i]*q-py^[i-1]*(p+q)+py^[i-2]*p; p:=q
    end;
  if n > 2 then L := 1 else L := 0;
  slegpb(n-1, L, {2,} t^[1], b^[0], pd2s^[1], ca, term);
  freemem(h, n*s);
  freemem(b, (n-1)*s);
  freemem(t, 2*(n-1)*s);
end; {ipfisn}

function ipfspn(n: ArbInt; var x, y, d2s: ArbFloat; t:ArbFloat;
                var term: ArbInt): ArbFloat;

var
   px, py       : ^arfloat0;
   pd2s         : ^arfloat1;
   i, j, m      : ArbInt;
   d, s3, h, dy : ArbFloat;
begin
  term:=1;
  if n<1
  then
    begin
      term:=3; exit
    end; {n<1}
  px:=@x; py:=@y; pd2s:=@d2s;
  if n = 1
  then
    begin
      h:=px^[1]-px^[0];
      dy:=(py^[1]-py^[0])/h;
      ipfspn:=py^[0]+(t-px^[0])*dy
    end { n = 1 }
  else
  if t <= px^[0]
  then
    begin
      h:=px^[1]-px^[0];
      dy:=(py^[1]-py^[0])/h-h*pd2s^[1]/6;
      ipfspn:=py^[0]+(t-px^[0])*dy
    end { t <= x[0] }
  else
  if t >= px^[n]
  then
    begin
      h:=px^[n]-px^[n-1];
      dy:=(py^[n]-py^[n-1])/h+h*pd2s^[n-1]/6;
      ipfspn:=py^[n]+(t-px^[n])*dy
    end { t >= x[n] }
  else
    begin
      i:=0; j:=n;
      while j <> i+1 do
        begin
          m:=(i+j) div 2;
          if t>=px^[m]
          then
            i:=m
          else
            j:=m
        end; {j}
      h:=px^[i+1]-px^[i];
      d:=t-px^[i];
      if i=0
      then
        begin
          s3:=pd2s^[1]/h;
          dy:=(py^[1]-py^[0])/h-h*pd2s^[1]/6;
          ipfspn:=py^[0]+d*(dy+d*d*s3/6)
        end
      else
      if i=n-1
      then
        begin
          s3:=-pd2s^[n-1]/h;
          dy:=(py^[n]-py^[n-1])/h-h*pd2s^[n-1]/3;
          ipfspn:=py^[n-1]+d*(dy+d*(pd2s^[n-1]/2+d*s3/6))
        end
      else
        begin
          s3:=(pd2s^[i+1]-pd2s^[i])/h;
          dy:=(py^[i+1]-py^[i])/h-h*(2*pd2s^[i]+pd2s^[i+1])/6;
          ipfspn:=py^[i]+d*(dy+d*(pd2s^[i]/2+d*s3/6))
        end
    end { x[0] < t < x[n] }
end; {ipfspn}

procedure ipfsmm(
  n: ArbInt; var x, y, d2s, minv, maxv: ArbFloat; var term: ArbInt);

var
  i: ArbInt;
  h: ArbFloat;
  px, py: ^arfloat0;
  pd2s: ^arfloat1;

  procedure UpdateMinMax(v: ArbFloat);
  begin
    if (0 >= v) or (v >= h) then exit;
    v := ipfspn(n, x, y, d2s, px^[i]+v, term);
    if v < minv then
      minv := v;
    if v > maxv then
      maxv := v;
  end;

  procedure MinMaxOnSegment;
  var
    a, b, c: ArbFloat;
    d: ArbFloat;
  begin
    h:=px^[i+1]-px^[i];
    if i=0
    then
      begin
        a:=pd2s^[1]/h/2;
        b:=0;
        c:=(py^[1]-py^[0])/h-h*pd2s^[1]/6;
      end
    else
    if i=n-1
    then
      begin
        a:=-pd2s^[n-1]/h/2;
        b:=pd2s^[n-1];
        c:=(py^[n]-py^[n-1])/h-h*pd2s^[n-1]/3;
      end
    else
      begin
        a:=(pd2s^[i+1]-pd2s^[i])/h/2;
        b:=pd2s^[i];
        c:=(py^[i+1]-py^[i])/h-h*(2*pd2s^[i]+pd2s^[i+1])/6;
      end;
    if a=0 then exit;
    d := b*b-4*a*c;
    if d<0 then exit;
    d:=Sqrt(d);
    UpdateMinMax((-b+d)/(2*a));
    UpdateMinMax((-b-d)/(2*a));
  end;

begin
  term:=1;
  if n<1 then begin
    term:=3;
    exit;
  end;
  if n = 1 then
    exit;
  px:=@x; py:=@y; pd2s:=@d2s;
  for i:=0 to n-1 do
    MinMaxOnSegment;
end;

procedure ipfish(hst: THermiteSplineType; n: ArbInt; var x, y, d1s: ArbFloat; var term: ArbInt);
var
  px, py, pd1s : ^arfloat0;
  i : ArbInt;
  dks : array of ArbFloat;
begin
  term:=1;
  if n < 1 then
  begin
    term:=3;
    exit;
  end;
  px:=@x;
  py:=@y;
  pd1s:=@d1s;

  {Monotone cubic Hermite interpolation}
  {See: https://en.wikipedia.org/wiki/Monotone_cubic_interpolation
   and: https://en.wikipedia.org/wiki/Cubic_Hermite_spline}

  {For each two adjacent data points, calculate tangent of the segment between them}
  SetLength(dks,n);
  for i:=0 to n-1 do
    dks[i]:=(py^[i+1]-py^[i])/(px^[i+1]-px^[i]);

  {As proposed by Fritsch and Carlson: For each data point - except the first and
   the last one - assign point's tangent (stored in a "d1s" array) as an average
   of tangents of the two adjacent segments (this is called 3PD, three-point
   difference) - but only if both tangents are either positive (segments are
   raising) or negative (segments are falling); in all other cases there is a local
   extremum at the data point, or a non-monotonic range begins/continues/ends there,
   so spline at this point must be flat to preserve monotonicity - so assign point's
   tangent as zero}
  for i:=0 to n-2 do
  if ((dks[i] > 0) and (dks[i+1] > 0)) or ((dks[i] < 0) and (dks[i+1] < 0)) then
    pd1s^[i+1]:=0.5*(dks[i]+dks[i+1])
  else
    pd1s^[i+1]:=0;

  {For the first and the last data point, assign point's tangent as a tangent of
   the adjacent segment (this is called one-sided difference)}
  pd1s^[0]:=dks[0];
  pd1s^[n]:=dks[n-1];

  {As proposed by Fritsch and Carlson: Reduce point's tangent if needed, to prevent
   overshoot}
  for i:=0 to n-1 do
  if dks[i] <> 0 then
  try
    if pd1s^[i]/dks[i] > 3 then
      pd1s^[i]:=3*dks[i];
    if pd1s^[i+1]/dks[i] > 3 then
      pd1s^[i+1]:=3*dks[i];
  except
    {There may be an exception for dks[i] values that are very close to zero}
    pd1s^[i]:=0;
    pd1s^[i+1]:=0;
  end;

  {Addition to the original algorithm: For the first and the last data point,
   modify point's tangent in such a way that the cubic Hermite interpolation
   polynomial has its inflection point exactly at the data point - so there
   will be a smooth transition to the extrapolated part of the graph}
  pd1s^[0]:=1.5*dks[0]-0.5*pd1s^[1];
  pd1s^[n]:=1.5*dks[n-1]-0.5*pd1s^[n-1];
end; {ipfish}

function ipfsph(n: ArbInt; var x, y, d1s: ArbFloat; t: ArbFloat; var term: ArbInt): ArbFloat;
var
   px, py, pd1s : ^arfloat0;
   i, j, m : ArbInt;
   h : ArbFloat;
begin
  term:=1;
  if n < 1 then
  begin
    term:=3;
    exit;
  end;
  px:=@x;
  py:=@y;
  pd1s:=@d1s;
  if t <= px^[0] then
    ipfsph:=py^[0]+(t-px^[0])*pd1s^[0]
  else
  if t >= px^[n] then
    ipfsph:=py^[n]+(t-px^[n])*pd1s^[n]
  else
  begin
    i:=0;
    j:=n;
    while j <> i+1 do
    begin
      m:=(i+j) div 2;
      if t>=px^[m] then
        i:=m
      else
        j:=m;
    end; {j}
    h:=px^[i+1]-px^[i];
    t:=(t-px^[i])/h;
    ipfsph:= py^[i]*(1+2*t)*Sqr(1-t) + h*pd1s^[i]*t*Sqr(1-t) + py^[i+1]*Sqr(t)*(3-2*t) + h*pd1s^[i+1]*Sqr(t)*(t-1);
  end;
end; {ipfsph}

function p(x, a, z:complex): ArbFloat;
begin
      x.sub(a);
      p:=x.Inp(z)
end;

function e(x, y: complex): ArbFloat;
const c1: ArbFloat = 0.01989436788646;
var s: ArbFloat;
begin x.sub(y);
      s := x.norm;
      if s=0 then e:=0 else e:=c1*s*ln(s)
end;

function spline(    n     : ArbInt;
                    x     : complex;
                var ac    : complex;
                var gammar: ArbFloat;
                    u1    : ArbFloat;
                    pf    : complex): ArbFloat;
var i     : ArbInt;
    s     : ArbFloat;
    a     : arcomp0 absolute ac;
    gamma : arfloat0 absolute gammar;
begin
    s := u1 + p(x, a[n-2], pf);
    for i:=0 to n do s := s + gamma[i]*e(x,a[i]);
    spline := s
end;

procedure splineparameters
          (    n                 : ArbInt;
           var ac, alfadc        : complex;
           var lambda,
               gammar, u1,
               kwsom, energie    : ArbFloat;
           var pf                : complex);

   procedure SwapC(var v, w: complex);
   var x: complex;
   begin
       x := v; v := w; w := x
   end;

   procedure pxpy(a, b, c: complex; var p:complex);
   var det: ArbFloat;
   begin
        b.sub(a); c.sub(a); det := b.xreal*c.imag-b.imag*c.xreal;
        b.sub(c); p.Init(b.imag/det, -b.xreal/det)
   end;

   procedure pfxpfy(a, b, c: complex; f: vector; var pf: complex);
   begin
      b.sub(a); c.sub(a);
      f.j := f.j-f.i; f.k := f.k-f.i;
      pf.init(f.j*c.imag - f.k*b.imag, -f.j*c.xreal + f.k*b.xreal);
      pf.scale(1/(b.xreal*c.imag - b.imag*c.xreal))
   end;

   function InpV(n: ArbInt; var v1, v2: ArbFloat): ArbFloat;
   var i: ArbInt;
       a1: arfloat0 absolute v1;
       a2: arfloat0 absolute v2;
       s : ArbFloat;
   begin
       s := 0;
       for i:=0 to n-1 do s := s + a1[i]*a2[i];
       InpV := s
   end;

   PROCEDURE SPDSOL(    N  : INTEGER;
                    VAR AP : pointer;
                    VAR B  : ArbFloat);
   VAR I, J, K : INTEGER;
       H       : ArbFloat;
       a       : ^ar2dr absolute ap;
       bx      : arfloat0 absolute b;
   BEGIN
      for k:=0 to n do
      BEGIN
          h := sqrt(a^[k]^[k]-InpV(k, a^[k]^[0], a^[k]^[0]));
          a^[k]^[k] := h;
          FOR I:=K+1 TO N do a^[i]^[k] := (a^[i]^[k] - InpV(k, a^[k]^[0], a^[i]^[0]))/h;
          BX[K] := (bx[k] - InpV(k, a^[k]^[0], bx[0]))/h
      END;
      FOR I:=N DOWNTO 0 do
      BEGIN
          H := BX[I];
          FOR J:=I+1 TO N DO H := H - A^[J]^[I]*BX[J];
          BX[I] := H/A^[I]^[I]
      END
   END;

var i, j, i1 : ArbInt;
    x, h,
    absdet,
    absdetmax,
    s, s1, ca: ArbFloat;
    alfa, dv, hulp,
    u, v, w  : vector;
    e22      : array[0..2] of vector;
    e21, b   : ^arvect0;
    k, c     : ^ar2dr;
    gamma    : arfloat0 absolute gammar;
    an2, an1, an, z,
    vz, wz   : complex;
    a        : arcomp0 absolute ac;
    alfad    : arcomp0 absolute alfadc;

begin

  i1:=0;
  x:=a[0].xreal;
  for i:=1 to n do
  begin
       h:=a[i].xreal;
       if h<x then begin i1:=i; x:=h end
  end;
  SwapC(a[n-2], a[i1]);
  SwapC(alfad[n-2], alfad[i1]);

  x:=a[0].xreal;
  i1 := 0;
  for i:=1 to n do
  begin
       h:=a[i].xreal;
       if h>x then begin i1:=i; x:=h end
  end;
  SwapC(a[n-1], a[i1]);
  SwapC(alfad[n-1], alfad[i1]);

  vz := a[n-2]; vz.sub(a[n-1]);

  absdetmax := -1;
  for i:=0 to n do
  begin
    wz := a[i]; wz.sub(a[n-2]);
    absdet := abs(wz.imag*vz.xreal-wz.xreal*vz.imag);
    if absdet>absdetmax then begin i1:=i; absdetmax:=absdet end
  end;
  SwapC(a[n], a[i1]);
  SwapC(alfad[n], alfad[i1]);

  an2 := a[n-2]; an1 := a[n-1]; an := a[n];
  alfa.i := alfad[n-2].xreal; dv.i := alfad[n-2].imag;
  alfa.j := alfad[n-1].xreal; dv.j := alfad[n-1].imag;
  alfa.k := alfad[n  ].xreal; dv.k := alfad[n  ].imag;

  n := n - 3;

  GetMem(k, (n+1)*SizeOf(pointer));
  for j:=0 to n do GetMem(k^[j], (j+1)*SizeOf(ArbFloat));

  GetMem(e21, (n+1)*SizeOf(vector));
  GetMem(b, (n+1)*SizeOf(vector));

  pxpy(an2,an1,an,z); for i:=0 to n do b^[i].i:=1+p(a[i],an2,z);
  pxpy(an1,an,an2,z); for i:=0 to n do b^[i].j:=1+p(a[i],an1,z);
  pxpy(an,an2,an1,z); for i:=0 to n do b^[i].k:=1+p(a[i],an,z);

  e22[0].init(0,e(an1,an2),e(an,an2));
  e22[1].init(e(an1,an2),0,e(an,an1));
  e22[2].init(e(an,an2),e(an,an1),0);

  for j:=0 to n do e21^[j].init(e(an2,a[j]),e(an1,a[j]),e(an,a[j]));

  GetMem(c, (n+1)*SizeOf(pointer));
  for j:=0 to n do GetMem(c^[j], (j+1)*SizeOf(ArbFloat));

  for i:=0 to n do
  for j:=0 to i do
  begin
    if j=i then s:=0 else s:=e(a[i],a[j]);
    hulp.init(b^[j].Inprod(e22[0]), b^[j].Inprod(e22[1]), b^[j].Inprod(e22[2]));
    hulp.sub(e21^[j]);
    k^[i]^[j] := s+b^[i].InProd(hulp)-b^[j].Inprod(e21^[i]);
    if j=i then s:=1/alfad[i].imag else s:=0;
    hulp.init(b^[j].i/dv.i, b^[j].j/dv.j, b^[j].k/dv.k);
    c^[i]^[j] := k^[i]^[j] + (s + b^[i].Inprod(hulp))/lambda
  end;

  for i:=0 to n do gamma[i]:=alfad[i].xreal - b^[i].Inprod(alfa);

  SpdSol(n, pointer(c), gamma[0]);

  for j:=n downto 0 do FreeMem(c^[j], (j+1)*SizeOf(ArbFloat));
  FreeMem(c, (n+1)*SizeOf(pointer));

  s:=0; for j:=0 to n do s:=s+b^[j].i*gamma[j]; w.i:=s; gamma[n+1] := -s;
  s:=0; for j:=0 to n do s:=s+b^[j].j*gamma[j]; w.j:=s; gamma[n+2] := -s;
  s:=0; for j:=0 to n do s:=s+b^[j].k*gamma[j]; w.k:=s; gamma[n+3] := -s;
  FreeMem(b, (n+1)*SizeOf(vector));

  u.init(w.i/dv.i, w.j/dv.j, w.k/dv.k);
  u.scale(1/lambda);
  u.add(alfa);

  s:=0; for j:=0 to n do s:=s+e21^[j].i*gamma[j]; v.i := e22[0].inprod(w)-s;
  s:=0; for j:=0 to n do s:=s+e21^[j].j*gamma[j]; v.j := e22[1].inprod(w)-s;
  s:=0; for j:=0 to n do s:=s+e21^[j].k*gamma[j]; v.k := e22[2].inprod(w)-s;
  FreeMem(e21, (n+1)*SizeOf(vector));

  u.add(v);

  pfxpfy(an2, an1, an, u, pf); u1:=u.i;

  kwsom := 0; for j:=0 to n do kwsom:=kwsom+sqr(gamma[j])/alfad[j].imag;
  kwsom := kwsom+sqr(w.i)/dv.i+sqr(w.j)/dv.j+sqr(w.k)/dv.k;
  kwsom := kwsom/sqr(lambda);

  s:=0;
  for i:=0 to n do
  begin s1:=0;
        for j:=0 to i do s1:=s1+k^[i]^[j]*gamma[j];
        for j:=i+1 to n do s1:=s1+k^[j]^[i]*gamma[j];
        s := gamma[i]*s1+s
  end;
  for j:=n downto 0 do FreeMem(k^[j], (j+1)*SizeOf(ArbFloat));
  FreeMem(k, (n+1)*SizeOf(pointer));
  energie := s

end {splineparameters};

end.