unit Unit_func;

interface
uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
  StdCtrls;
{
////////////////////////////////////////////////////////////////////////////////////
// y=1/(1+x^2)
//區間為：[-5,5]
const PN=100;
const N=11; //節點
const X:array[1..N] of double=(-5,-4,-3,-2,-1,0,1,2,3,4,5);
const Y:array[1..N] of double=();
////////////////////////////////////////////////////////////////////////////////////
}
function fun1(a, b:real; n:integer):real;
procedure POLINT(XA, YA:array of real; N:integer; X:real; var Y, DY:real);
implementation

function fun1(a, b:real; n:integer):real;
var
  XA, YA:array[0..10] of real;
  X, Y, FF, DY:real;
  F:TextFile;
  I, NFUNC:Integer;
const
  N2 = 10;
  x_left = -5; //左右邊界
  y_right = 5;
  // PI = 3.1415926;
  s1 = '%12.6f';
  s2 = '.0000E+00';
  //---------y=1/(1+x^2)------------------
  function cal(x:real):real;
  begin
    result := 1 / (1 + x * x);
  end;
  //---------------------------
begin
  //輸出計算結果到文件
  AssignFile(F, 'func1.txt');
  Rewrite(F);
  Writeln(F, 'Generation of interpolation tables');
  Writeln(F, '... y=1/(1+x^2)    -5<x<5');
  Writeln(F, 'How many entries go in these tables?(Note: N<10)');
  Writeln(F, N2);
  for I := 1 to N2 do
  begin
    XA[I] := I * (y_right - x_left) / N2; //N等份后 ，每等份的值
    YA[I] := 1 / (1 + XA[I] * XA[I]); //Sin(XA[I]);
  end;
  Writeln(F, '       x           f(x)    interpolated   error');
  for I := 0 to 10 do
  begin
    //[-5,5]
     // X := (-0.05 + I / 10) * PI;
    X := -5 + i;
    FF := 1 / (1 + X * X);
    POLINT(XA, YA, N2, X, Y, DY);
    Writeln(F, Format(s1, [X]), Format(s1, [FF]), Format(s1, [Y]),
      '    ', FormatFloat(s2, DY));
    //end;

    //屏幕顯示計算結果
  end;
  CloseFile(F);
end;

procedure POLINT(XA, YA:array of real; N:integer; X:real; var Y, DY:real);
var
  C, D:array[0..10] of real;
  DIF, DIFT, HO, HP, W, DEN:real;
  NS, I, M:integer;
begin
  NS := 1;
  DIF := Abs(X - XA[1]);
  for I := 1 to N do
  begin
    DIFT := Abs(X - XA[I]);
    if DIFT < DIF then
    begin
      NS := I;
      DIF := DIFT;
    end;
    C[I] := YA[I];
    D[I] := YA[I];
  end;
  Y := YA[NS];
  NS := NS - 1;
  for M := 1 to N - 1 do
  begin
    for I := 1 to N - M do
    begin
      HO := XA[I] - X;
      HP := XA[I + M] - X;
      W := C[I + 1] - D[I];
      DEN := HO - HP;
      if DEN = 0 then
      begin
        ShowMessage('PAUSE');
        Exit;
      end;
      DEN := W / DEN;
      D[I] := HP * DEN;
      C[I] := HO * DEN;
    end;
    if 2 * NS < N - M then
      DY := C[NS + 1]
    else
    begin
      DY := D[NS];
      NS := NS - 1;
    end;
    Y := Y + DY;
  end;
end;
end.