#=== PLS functions (Recursive implamentation) ===#

P1(x,x1)=x1-x
P2(x,x1,x2)=P1(x,x1)*P1(x,x2)
P32(x,x1,x2)=P1(x,x2)*P1(x,x1)*P1(x,x1)
P3(x,x1,x2,x3)=P2(x,x1,x2)*P1(x,x3)
L0(x,x0,y0,a)=y0 + a*(x-x0)
L1(x,x0,y0,a0,a1)=L0(x,x0,y0,a0)*heav(x0-x)+heav(x-x0)*L0(x,x0,y0,a1)
L2(x,x0,y0,x1,y1,a0,a2)=L0(x,x0,y0,a0)*heav(x0-x)+heav(x-x0)*L1(x,x1,y1,(y1-y0)/(x1-x0),a2)
L3(x,x0,y0,x1,y1,x2,y2,a0,a3)=L0(x,x0,y0,a0)*heav(x0-x)+heav(x-x0)*L2(x,x1,y1,x2,y2,(y1-y0)/(x1-x0),a3)
S1(x,x0,y0,y1)=y0*heav(x0-x)+heav(x-x0)*y1
S2(x,x0,x1,y0,y1,y2)=S1(x,x0,y0,S1(x,x1,y1,y2) )
S3(x,x0,x1,x2,y0,y1,y2,y3)=S1(x,x0,y0,S2(x,x1,x2,y1,y2,y3))