BeginPackage["alginterpolation`"]

ainterpolation::usage =
   "ainterpolation[F,N] gives a N-valued algebraic function y(x)\n
which is a set of roots of the polynomial equation\n
1 + a1 x + a2 y = 0   (N=1)\n
1 + a1 x + a2 y + b1 x^2 + b2 x y + b3 y^2 = 0   (N=2)\n
1 + a1 x + a2 y + b1 x^2 + b2 x y + b3 y^2\n
         + c1 x^3 + c2 x^2 y + ... = 0   (N=3) ...\n
where the coefficients a1, a2, b1, b2, b3, c1, c2, ... are determined\n
from the set of M = (N+1)(N+2)/2-1 linear equations at\n
{x=x1=F(1,1),y=y1=F(1,2)}, {x=x2=F(2,1),y=y2=F(2,2)}, ..., {x=x_M=F(M,1),y=y_M=F(M,2)}.\n\n
N should be a positive integer.\n\n
F should be a double-indexed array {{x1,y1},{x2,y2},...,{x_M,y_M},...(the rest is not used)...}\n
   of length at least (M,2).\n\n
Returnes the pure function.\n\n
This is the basic function that is called from the interpolation subroutine alginterpolation"

alginterpolation::usage =
   "alginterpolation[F,N,xval] fits data points {{X1,Y1},{X2,Y2},...}\n
by a branch of the algebraic function y(x)\n
which is a root of the polynomial equation\n
1 + a1 x + a2 y = 0   (N=1)\n
1 + a1 x + a2 y + b1 x^2 + b2 x y + b3 y^2 = 0   (N=2)\n
1 + a1 x + a2 y + b1 x^2 + b2 x y + b3 y^2\n
         + c1 x^3 + c2 x^2 y + ... = 0   (N=3) ...\n
where the coefficients a1, a2, b1, b2, b3, c1, c2, ... are determined\n
from the set of M = (N+1)(N+2)/2-1 linear equations at\n
{x=x1=F(1,1),y=y1=F(1,2)}, {x=x2=F(2,1),y=y2=F(2,2)}, ..., {x=x_M=F(M,1),y=y_M=F(M,2)}.\n
Here, points x1,...,x_M are M the closest to x0 points from the set {X1,X2,X3,...}\n\n
N should be a non-negative integer (if N=0, returnes the constant y1).\n
F should be a double-indexed array {{x1,y1},{x2,y2},...,{x_M,y_M},...}\n
   of length at least (M,2).\n
xval should be a number.\n
Returnes y(xval).\n
An option intermediatePoints (9 by default) is the number of intermediate points\n
   where y(x) is calculated to trace the correct continuous branch"

ipoints=intermediatePoints

Begin["`private`"]

ainterpolation[xy_List,nn_Integer]:=(
mm=(nn+1)(nn+2)/2-1; (* number of variables *)
eq=1+Sum[Sum[ny=n-nx;c[nx,ny] x^nx y^ny,{nx,0,n}],{n,nn}]==0;
varlist=Table[Table[ny=n-nx;c[nx,ny],{nx,0,n}],{n,nn}]//Flatten;
eqs=Table[eq/.{x->xy[[m,1]],y->xy[[m,2]]},{m,mm}];
scf=Solve[eqs,varlist];
Table[y/.(Solve[(eq/.x->#)/.scf,y][[n]]),{n,nn}]&
)

Options[alginterpolation]={ipoints->9}

alginterpolation[xy_List,nn_Integer,xval_,opts___]:=(
xysort=Sort[xy,Abs[xval-#1[[1]]]<Abs[xval-#2[[1]]]&];
{xnearest,ynearest}=xysort//First;
If[nn==0||xval==xnearest,ynearest,
   f=ainterpolation[xysort,nn];
   If[nn==1,f[xval][[1]], (* If N=1, then tracing of the branch is trivial *)
      pts=ipoints/.Flatten[{opts,Options[alginterpolation]}];
      mp=pts+1;
      y1=ynearest;
      Do[x1=xnearest+np*(xval-xnearest)/mp;
         y1=Sort[f[x1],Abs[y1-#1]<Abs[y1-#2]&]//First
        ,{np,mp}];
      y1
     ]
  ]
)

End[]

EndPackage[]