fpc/packages/numlib/tests/intge2te.pas

program InitInfEx;

uses
  Typ,
  Spe,
  Int;

var
  num2: ArbInt;
  inte: ArbFloat;
const
  e = 2.71828182845905;

  function cx(x: ArbFloat): ArbFloat;
  begin
    cx := cos(x) / (sqr(x) + 1);
  end;

  function pcx(x: ArbFloat): ArbFloat;
  begin
    pcx := 2 * x * sin(x) / sqr(sqr(x) + 1);
  end;

  function ppcx(x: ArbFloat): ArbFloat;
  var
    s: ArbFloat;
  begin
    s    := sqr(x);
    ppcx := cos(x) * (2 - 6 * s) / ((s + 1) * sqr(s + 1));
  end;

  function cc2(x: ArbFloat): ArbFloat;
  var
    x2: ArbFloat;
  begin
    x2  := sqr(x);
    cc2 := cos(x) / (sqr(x2) + x2 + 1);
  end;

  function ucx(x: ArbFloat): ArbFloat;
  begin
    if x = 0 then
      ucx := 0
    else
      ucx := cx((1 - x) / x) / sqr(x);
  end;

  function ucc2(x: ArbFloat): ArbFloat;
  begin
    if x = 0 then
      ucc2 := 0
    else
      ucc2 := cc2((1 - x) / x) / sqr(x);
  end;

  function uz(x: ArbFloat): ArbFloat;
  begin
    uz := sin(x) * exp(-x);
  end;

  function ss1(x: ArbFloat): ArbFloat;    { f(«n(n-1))=0 (n=1,2,3,...) }
  var
    n, s, c: ArbFloat;                    { f(«ný)=2/(ný(n+1))  }
  begin                                   { overigens: f linear interpoleren}
    s := sqrt(2 * x);
    n := trunc(s);
    if n * (n + 1) / 2 <= x then
      n := n + 1;  { n z.d.d. «n(n-1) ó x ó «n(n-1) }
    c := 4 / (n * sqr(n) * (n + 1));
    if s < n then
      ss1 := c * (x - n * (n - 1) / 2)
    else
      ss1 := c * (n * (n + 1) / 2 - x);
  end;

  function ss2(x: ArbFloat): ArbFloat;   { als ss1 met f(«ný)=2/(n.2ü)  }
  var
    n, s, c: ArbFloat;
  begin
    s := sqrt(2 * x);
    n := trunc(s);
    if n * (n + 1) / 2 <= x then
      n := n + 1;
    c := spepow(2, 2 - n) / sqr(n);
    if s < n then
      ss2 := c * (x - n * (n - 1) / 2)
    else
      ss2 := c * (n * (n + 1) / 2 - x);
  end;

  function ss3(x: ArbFloat): ArbFloat;       { x  z.d.d. «.2ü ó x+1 ó 2ü  s=2ü}
  var
    s, c, f, x1: ArbFloat;
    n: ArbInt;        {n even: f(n)=-4/(n.2ü)}
  begin                                  {n oneven: f(n) = 4/(n.2ü)}
    n  := 0;
    s  := 1;
    x1 := x + 1;               { overigens: f lineair interpol.}
    repeat
      n := n + 1;
      s := s * 2
    until s > x1;
    c := 16 / (n * sqr(s));
    if x1 < 0.75 * s then
      f := c * (x1 - s / 2)
    else
      f := c * (s - x1);
    if odd(n) then
      ss3 := f
    else
      ss3 := -f;
  end;

  function ss4(x: ArbFloat): ArbFloat;        { 0 ó x ó 1}
  var
    y, h:  ArbFloat;                        { zij x = ä [1:ì] c(n)/3ü c(n)=0,1,2}
    ready: boolean;                         { zoek kleinste k met c(k)=1}
  begin                                     { dan geldt f(x)=u(y)/3k }
    y     := 3 * x;
    h     := 1 / 3;
    ready := False;         { met u(y)=|y-1«|, 1 ó y ó 2}
    repeat                                { en y =  ä [k:ì] (c(n)/3ü)*3k }
      if (y < 1) or (y > 2) then
      begin
        if y < 1 then
          y := 3 * y
        else
          y := 3 * (y - 2);
        h := h / 3;
        if h < macheps then
        begin
          ready := True;
          ss4   := 0;
        end;
      end
      else
      begin
        ready := True;
        if y < 1.5 then
          ss4 := h * (y - 1)
        else
          ss4 := h * (2 - y);
      end
    until ready;
  end;

  function ss5(x: ArbFloat): ArbFloat;   { uitbreiding ss4}
  var
    y, h:  ArbFloat;                   {functiewaarden op 'volgend' interval}
    ready: boolean;                    { [n, n+1] telkens halveren}
  begin
    y     := 3 * frac(x);
    h     := spepow(0.5, trunc(x)) / 3;
    ready := False;
    repeat
      if (y < 1) or (y > 2) then
      begin
        if y < 1 then
          y := 3 * y
        else
          y := 3 * (y - 2);
        h := h / 3;
        if h < macheps then
        begin
          ready := True;
          ss5   := 0;
        end;
      end
      else
      begin
        ready := True;
        if y < 1.5 then
          ss5 := h * (y - 1)
        else
          ss5 := h * (2 - y);
      end
    until ready;
  end;

  function ss6(x: ArbFloat): ArbFloat;   { 0 ó x ó 1, 'gladdere' variant van ss4}
  var
    y, h:  ArbFloat;
    ready: boolean;

    function f(y: ArbFloat): ArbFloat; { 1 ó y ó 2, 1 x cont. diff, symm. max in 1.5}
    begin
      if y > 1.5 then
        y := 3 - y;
      if y < 1.25 then
        f := sqr(y - 1)
      else
        f := 0.125 - sqr(1.5 - y);
    end;

  begin
    y     := 3 * x;
    h     := 1 / 3;
    ready := False;
    repeat
      if (y < 1) or (y > 2) then
      begin
        if y < 1 then
          y := 3 * y
        else
          y := 3 * (y - 2);
        h := h / 3;
        if h < macheps then
        begin
          ready := True;
          ss6   := 0;
        end;
      end
      else
      begin
        ready := True;
        ss6   := h * f(y);
      end
    until ready;
  end;

  function bb(x: ArbFloat): ArbFloat;
  begin
    bb := spepow(x, -x) * (ln(x) + 1);
  end;

var
  integral, ae, err: ArbFloat;
  term:  ArbInt;
  intex: boolean;

  procedure Header;
  begin
    Write('int': num2, '': numdig - num2, ' ', 'err': 7, ' ': 4);
    if intex then
      Write('f': 6);
    writeln;
  end;

  procedure ShowResults;
  var
    f: ArbFloat;
  begin
    if intex then
      f := inte - integral;
    case term of
      1:
      begin
        Write(integral: numdig, ' ', err: 10, ' ');
        if intex then
          writeln(f: 10)
        else
          writeln;
      end;
      2:
      begin
        Write(integral: numdig, ' ', err: 10, ' ');
        if intex then
          writeln(f: 10)
        else
          writeln;
        Writeln('    process afgebroken, te hoge nauwkeurigheid?');
      end;
      3: Writeln('Verkeerde parameterwaarde (<=0) bij aanroep: ', ae: 8);
      4:
      begin
        Write(integral: numdig, ' ', err: 10, ' ');
        if intex then
          writeln(f: 10)
        else
          writeln;
        writeln('    process afgebroken, moeilijk, mogelijk divergent?');
      end;
    end;
  end;

begin
  num2 := numdig div 2;

  Writeln('     ì   ');
  Writeln('    ô  cos x           ã ');
  Writeln('    ³ ------- dx   =   -- ');
  Writeln('  0 õ  xý+ 1           2e ');
  writeln;
  ae := 1e-8;
  Writeln('   Gevraagde nauwkeurigheid', ae: 12);
  writeln;
  inte  := 0.5 * pi / e;
  intex := True;
  writeln(inte: numdig, ' is ''exacte'' oplossing');
  Int1fr(@cx, 0, infinity, ae, integral, err, term);
  Header;
  ShowResults;
  writeln;
  Writeln('berekend met Int1fr via transformatie x=(1-t)/t');
  writeln;
  Int1fr(@ucx, 0, 1, ae, integral, err, term);
  Header;
  ShowResults;
  Writeln('     ì   ');
  Writeln('    ô  2x sin x           ã ');
  Writeln('    ³ --------- dx   =    -- ');
  Writeln('  0 õ  (xý+ 1)ý           2e ');
  writeln;
  ae := 1e-8;
  Writeln('   Gevraagde nauwkeurigheid', ae: 12);
  writeln;
  Int1fr(@pcx, 0, infinity, ae, integral, err, term);
  Header;
  ShowResults;
  Writeln('     ì   ');
  Writeln('    ô (2-6xý)cos x            ã ');
  Writeln('    ³ ------------ dx   =    -- ');
  Writeln('  0 õ  (xý+ 1)3             2e ');
  writeln;
  ae := 1e-8;
  Writeln('   Gevraagde nauwkeurigheid', ae: 12);
  writeln;
  Int1fr(@ppcx, 0, infinity, ae, integral, err, term);
  Header;
  ShowResults;

  Writeln('     ì   ');
  Writeln('    ô     cos x ');
  Writeln('    ³ ------------ dx   =  (ã/û3) exp(-«û3) sin(ã/6+«) ');
  Writeln('  0 õ  (xý)ý+xý+ 1 ');
  writeln;
  writeln;
  ae := 1e-8;
  Writeln('   Gevraagde nauwkeurigheid', ae: 12);
  writeln;
  inte := (pi / sqrt(3)) * exp(-sqrt(0.75)) * sin(pi / 6 + 0.5);
  writeln(inte: numdig, ' is ''exacte'' oplossing');
  Int1fr(@cc2, 0, infinity, ae, integral, err, term);
  Header;
  ShowResults;
  writeln;
  Writeln('berekend met Int1fr via transformatie x=(1-t)/t');
  writeln;
  writeln;
  writeln(inte: numdig, ' is ''exacte'' oplossing');
  Int1fr(@ucc2, 0, 1, ae, integral, err, term);
  Header;
  ShowResults;

  Writeln('     ì   ');
  Writeln('    ô sin u                ');
  Writeln('    ³ ------ du   =  « ');
  Writeln('    õ exp(u)            ');
  writeln('   0 ');
  ae    := 1e-8;
  intex := True;
  inte  := 0.5;
  Writeln('   Gevraagde nauwkeurigheid', ae: 12);
  writeln;
  Int1fr(@uz, 0, infinity, ae, integral, err, term);
  Header;
  ShowResults;

  writeln(' functie ss1;  int = ä {1:ì}1/n(n+1)  =  1');
  ae    := 1e-8;
  intex := True;
  inte  := 1;
  Writeln('   Gevraagde nauwkeurigheid', ae: 12);
  int1fr(@ss1, 0, infinity, ae, integral, err, term);
  Header;
  Showresults;

  writeln(' functie ss2;  int = ä {1:ì} («)ü  =  1');
  ae    := 1e-8;
  intex := True;
  inte  := 1;
  Writeln('   Gevraagde nauwkeurigheid', ae: 12);
  int1fr(@ss2, 0, infinity, ae, integral, err, term);
  Header;
  Showresults;

  writeln(' functie ss3;  int = ä {1:ì} (-1)ü/n  =  ln(2)');
  ae    := 1e-8;
  intex := True;
  inte  := ln(2);
  Writeln('   Gevraagde nauwkeurigheid', ae: 12);
  int1fr(@ss3, 0, infinity, ae, integral, err, term);
  Header;
  Showresults;

  writeln(' functie ss4 (op [0,1]); int = 1/28 ');
  ae    := 1e-8;
  intex := True;
  inte  := 1 / 28;
  Writeln('   Gevraagde nauwkeurigheid', ae: 12);
  int1fr(@ss4, 0, 1, ae, integral, err, term);
  Header;
  Showresults;

  writeln(' functie ss5 (op [0,ì)); int = 1/14 ');
  ae    := 1e-8;
  intex := True;
  inte  := 1 / 14;
  Writeln('   Gevraagde nauwkeurigheid', ae: 12);
  int1fr(@ss5, 0, infinity, ae, integral, err, term);
  Header;
  Showresults;

  writeln(' functie ss6 (op [0,1]); int = 1/112 ');
  ae    := 1e-8;
  intex := True;
  inte  := 1 / 112;
  Writeln('   Gevraagde nauwkeurigheid', ae: 12);
  int1fr(@ss6, 0, 1, ae, integral, err, term);
  Header;
  Showresults;

  Writeln('     ì   ');
  Writeln('    ô  ln(x)+1         ');
  Writeln('    ³ ---------- dx   = 1  ');
  Writeln('    õ     xx            ');
  writeln('   1 ');
  ae    := 1e-8;
  intex := True;
  inte  := 1;
  Writeln('   Gevraagde nauwkeurigheid', ae: 12);
  writeln;
  Int1fr(@bb, 0, infinity, ae, integral, err, term);
  Header;
  ShowResults;
end.