program assignment2_2;

type points = array [0..3] of real;
{const xcal=[0,1,2,3];}
var xcal, ycal:points;
    x, y:real;
    error:boolean;

procedure interp(xcal, ycal:points; x:real; var y:real; var error:boolean);

var a, b, c, d:real;

begin

  if (x<xcal[0]) or (x>xcal[3]) then begin

    y:=0;
    error:=TRUE;

    writeln('Error: Value of x given outside of allowed range.');

  end else begin

    a:=ycal[0]/((xcal[0]-xcal[1])*(xcal[0]-xcal[2])*(xcal[0]-xcal[3]));
    b:=ycal[1]/((xcal[1]-xcal[0])*(xcal[1]-xcal[2])*(xcal[1]-xcal[3]));
    c:=ycal[2]/((xcal[2]-xcal[0])*(xcal[2]-xcal[1])*(xcal[2]-xcal[3]));
    d:=ycal[3]/((xcal[3]-xcal[0])*(xcal[3]-xcal[1])*(xcal[3]-xcal[2]));

    y:=a*(x-xcal[1])*(x-xcal[2])*(x-xcal[3])+b*(x-xcal[0])*(x-xcal[2])*(x-xcal[3])+c*(x-xcal[0])*(x-xcal[1])*(x-xcal[3])+d*(x-xcal[0])*(x-xcal[1])*(x-xcal[2]);
    error:=FALSE;

  end;

end;

begin

  xcal[0]:=0;
  xcal[1]:=1;
  xcal[2]:=2;
  xcal[3]:=3;

  ycal[0]:=0;
  ycal[1]:=3;
  ycal[2]:=14;
  ycal[3]:=39;

  writeln('What value of x would you like to interpolate for?');
  readln(x);

  interp(xcal, ycal, x, y, error);

  if error=FALSE then writeln('(',x:1:3,', ',y:1:3,')');

  readln;

end.