Program Pas0101(arg);
uses
    SysUtils, Math;
begin
    writeln(3 + 5);
    writeln(3 - 5);
    writeln(3 * 5);
    writeln(power(3, 5));
    writeln(5 / 3);
    writeln(5 div 3);
    writeln(5 mod 3);

    write(IntToStr(3 * 5) + #13#10);

    writeln(format('%3d',     [3 * 5]));
    writeln(format('%23.20f', [5 / 3]));
end.
Program Pas0102(arg);
uses
    SysUtils;
var
    i:integer = 3 * 5;
begin
    writeln(format('3 * 5 = %d', [i]));
end.
Program Pas0103(arg);
uses
    SysUtils;
var
    i:integer;
begin
    for i := 1 to 9 do
    begin
        write(format('%d, ', [i]));
    end;
    writeln();
end.
Program Pas0104(arg);
uses
    SysUtils;
var
    i:integer;
begin
    for i := 1 to 9 do
    begin
        if (i mod 3 = 0) then
            write(format('%d, ', [i]));
    end;
    writeln();
end.
Program Pas0105(arg);
var
    i:integer;
    sum:integer = 0;
begin
    for i := 1 to 99 do
    begin
        if (i mod 3 = 0) then
            sum := sum + i;
    end;
    writeln(sum);
end.
Program Pas0301(arg);
{$MODE delphi}

// 初項:a, 公差:a で, 上限:lim の数列の総和を返す関数
function sn(a:Integer; lim:Integer):Integer;
var
    n, l:Integer;
begin
    n := lim div a;               // 項数:n  =  上限:lim / 公差:a
    l := n * a;                   // 末項:l  =  項数:n   * 公差:a
    result := (a + l) * n div 2;  // 総和:sn = (初項:a   + 末項:l) * 項数:n / 2
end;

begin
    // 3 の倍数の合計
    writeln(sn(3, 999));
end.
Program Pas0302(arg);
var
    n:Integer;
begin
    // 10000 までの 自然数の和
    // 項目数 n = 10000
    n := 10000;
    writeln( n * (n + 1) div 2 );
end.
Program Pas0303(arg);
var
    n:Integer;
begin
    // 10000 までの 偶数の和
    // 項目数 n = 5000
    n := 10000 div 2;
    writeln( n * (n + 1) );
end.
Program Pas0304(arg);
uses
    SysUtils, Math;
var
    n:Integer;
begin
    // 10000 までの 奇数の和
    // 項目数 n = 5000
    n := 10000 div 2;
    writeln( format('%g', [power(n, 2)]) ); // power は Extended 型
end.
Program Pas0305(arg);
var
    n:Integer;
begin
    // 1000 までの 自然数の２乗の和
    n := 1000;
    writeln( n * (n + 1) * (2 * n + 1) div 6 );
end.
Program Pas0306(arg);
uses
    SysUtils, Math;
var
    n:Integer;
begin
    // 100 までの 自然数の３乗の和
    n := 100;
    writeln( format('%g', [(power(n, 2) * power((n + 1), 2) / 4)]) );
end.
Program Pas0303(arg);
uses
    SysUtils, Math;
var
    n:Integer;
    a:Integer;
    r:Integer;
begin
    // 初項 2, 公比 3, 項数 10 の等比数列の和
    n := 10;
    a := 2;
    r := 3;
    writeln( format('%g', [(a * (power(r, n) - 1)) / (r - 1)]) );
end.
Program Pas0401(arg);
const
    a:Integer = 5;  // 初項 5
    d:Integer = 3;  // 公差 3
    n:Integer = 10; // 項数 10
var
    p:Int64   = 1;  // 積
    m:Integer;
    i:Integer;
begin
    for i := 1 to n do
    begin
        m := a + (d * (i - 1));
        p := p * m;
    end;
    writeln(p);
end.
Program Pas0402(arg);
{$MODE delphi}

function product(m:Integer; d:Integer; n:Integer):Int64;
begin
    if n = 0 then
        result := 1
    else
        result := m * product(m + d, d, n - 1);
end;

begin
    // 初項 5, 公差 3, 項数 10 の数列の積を表示する
    writeln(product(5, 3, 10));
end.
Program Pas0403(arg);
    // 階乗を求める関数
    function Fact(n: Integer): Longint;
    begin
        if n <= 1 then
            Fact := 1
        else
            Fact := n * Fact(n - 1);
    end;

begin
    // 10の階乗
    writeln(Fact(10));
    writeln(10 * 9 * 8 * 7 * 6 * 5 * 4 * 3 * 2 * 1);
end.
Program Pas0404(arg);

    // 下降階乗冪
    function FallingFact(x: Integer; n: Integer): Longint;
    begin
        if n <= 1 then
            FallingFact := x
        else
            FallingFact := x * FallingFact(x - 1, n - 1);
    end;

begin
    // 10 から 6 までの 総乗
    writeln(FallingFact(10, 5));
    writeln(10 * 9 * 8 * 7 * 6);
end.
Program Pas0405(arg);

    // 上昇階乗冪
    function RisingFact(x: Integer; n: Integer): Longint;
    begin
        if n <= 1 then
            RisingFact := x
        else
            RisingFact := x * RisingFact(x + 1, n - 1);
    end;

begin
    // 10 から 14 までの 総乗
    writeln(RisingFact(10, 5));
    writeln(10 * 11 * 12 * 13 * 14);
end.
Program Pas0406(arg);
uses
    SysUtils;

    // 階乗
    function Fact(n: Integer): Longint;
    begin
        if n <= 1 then
            Fact := 1
        else
            Fact := n * Fact(n - 1);
    end;

    // 下降階乗冪
    function FallingFact(x: Integer; n: Integer): Longint;
    begin
        if n <= 1 then
            FallingFact := x
        else
            FallingFact := x * FallingFact(x - 1, n - 1);
    end;

var
    n: Integer;
    r: Integer;
begin
    // 順列 (異なる 10 個のものから 5 個取ってできる順列の総数)
    n := 10;
    r := 5;
    writeln(format('%g', [Fact(n) / Fact(n - r)]));
    writeln(FallingFact(n, r));
end.
Program Pas0407(arg);
uses
    SysUtils, Math;
var
    n: Integer;
    r: Integer;
begin
    // 重複順列 (異なる 10 個のものから重複を許して 5 個取ってできる順列の総数)
    n := 10;
    r := 5;
    writeln(format('%g', [power(n, r)]));
end.
Program Pas040101(arg);

    // 組合せ
    function Comb(n: Integer; r: Integer): Longint;
    begin
        if (r = 0) or (r = n) then
            Comb := 1
        else if r = 1 then
            Comb := n
        else
            Comb := Comb(n - 1, r - 1) + Comb(n - 1, r);
    end;

var
    n: Integer;
    r: Integer;
begin
    // 組合せ (異なる 10 個のものから 5 個取ってできる組合せの総数)
    n := 10;
    r := 5;
    writeln(Comb(n, r));
end.
Program Pas0409(arg);

    // 組合せ
    function Comb(n: Integer; r: Integer): Longint;
    begin
        if (r = 0) or (r = n) then
            Comb := 1
        else if r = 1 then
            Comb := n
        else
            Comb := Comb(n - 1, r - 1) + Comb(n - 1, r);
    end;

var
    n: Integer;
    r: Integer;
begin
    // 重複組合せ (異なる 10 個のものから重複を許して 5 個とる組合せの総数)
    n := 10;
    r := 5;
    writeln(Comb(n + r - 1, r));
end.
Program Pas0501(arg);
{$MODE delphi}

uses
    SysUtils, Math;

// 自作の正弦関数
function mySin(x:Double; n:Integer; nega:Boolean; numerator:Double; denominator:Double; y:Double):Double;
var
    m: Integer;
    a: Double;
begin
    m           := 2 * n;
    denominator := denominator * (m + 1) * m;
    numerator   := numerator   * x * x;
    a           := numerator / denominator;
    // 十分な精度になったら処理を抜ける
    if (a <= 0.00000000001) then
        result := y
    else
    begin
        if not nega then a := -a;
        result := y + mySin(x, n + 1, not nega, numerator, denominator, a);
    end;
end;

var
    i:      Integer;
    degree: Integer;
    radian: Double;
    d1:     Double;
    d2:     Double;
begin
    for i := 0 to 24 do
    begin
        degree := i * 15;
        if (degree mod 30 = 0) or (degree mod 45 = 0) then
        begin
            radian := DegToRad(degree);
            // 自作の正弦関数
            d1 := mySin(radian, 1, false, radian, 1.0, radian);
            // 標準の正弦関数
            d2 := Sin(radian);
            // 標準関数との差異
            writeln(format('%3d : %13.10f - %13.10f = %13.10f', [degree, d1, d2, d1 - d2]));
        end;
    end;
end.
Program Pas0502(arg);
{$MODE delphi}

uses
    SysUtils, Math;

// 自作の余弦関数
function myCos(x:Double; n:Integer; nega:Boolean; numerator:Double; denominator:Double; y:Double):Double;
var
    m: Integer;
    a: Double;
begin
    m           := 2 * n;
    denominator := denominator * m * (m - 1);
    numerator   := numerator   * x * x;
    a           := numerator / denominator;
    // 十分な精度になったら処理を抜ける
    if (a <= 0.00000000001) then
        result := y
    else
    begin
        if not nega then a := -a;
        result := y + myCos(x, n + 1, not nega, numerator, denominator, a);
    end;
end;

var
    i:      Integer;
    degree: Integer;
    radian: Double;
    d1:     Double;
    d2:     Double;
begin
    for i := 0 to 24 do
    begin
        degree := i * 15;
        if (degree mod 30 = 0) or (degree mod 45 = 0) then
        begin
            radian := DegToRad(degree);
            // 自作の余弦関数
            d1 := myCos(radian, 1, false, 1.0, 1.0, 1.0);
            // 標準の余弦関数
            d2 := Cos(radian);
            // 標準関数との差異
            writeln(format('%3d : %13.10f - %13.10f = %13.10f', [degree, d1, d2, d1 - d2]));
        end;
    end;
end.
Program Pas0503(arg);
{$MODE delphi}

uses
    SysUtils, Math;

// 自作の正接関数
function myTan(x:Double; x2:Double; n:Integer; t:Double):Double;
begin
    t := x2 / (n - t);
    n := n - 2;
    if (n <= 1) then
        result := x / (1 - t)
    else
        result := myTan(x, x2, n, t);
end;

var
    i:      Integer;
    degree: Integer;
    radian: Double;
    x2:     Double;
    d1:     Double;
    d2:     Double;
begin
    for i := 0 to 12 do
    begin
        if (i * 15 mod 180 <> 0) then
        begin
            degree := i * 15 - 90;
            radian := DegToRad(degree);
            // 自作の正接関数
            x2 := radian * radian;
            d1 := myTan(radian, x2, 15, 0.0); // 15:必要な精度が得られる十分大きな奇数
            // 標準の正接関数
            d2 := Tan(radian);
            // 標準関数との差異
            writeln(format('%3d : %13.10f - %13.10f = %13.10f', [degree, d1, d2, d1 - d2]));
        end;
    end;
end.
Program Pas0504(arg);
{$MODE delphi}

uses
    SysUtils, Math;

// 自作の指数関数
function myExp(x:Double; n:Integer; numerator:Double; denominator:Double; y:Double):Double;
var
    a: Double;
begin
    denominator := denominator * n;
    numerator   := numerator   * x;
    a           := numerator / denominator;
    // 十分な精度になったら処理を抜ける
    if (Abs(a) <= 0.00000000001) then
        result := y
    else
        result := y + myExp(x, n + 1, numerator, denominator, a);
end;

var
    i:  Integer;
    x:  Double;
    d1: Double;
    d2: Double;
begin
    for i := 0 to 20 do
    begin
        x  := (i - 10) / 4.0;
        // 標準の指数関数
        d1 := Exp(x);
        // 自作の指数関数
        d2 := myExp(x, 1, 1.0, 1.0, 1.0);
        // 標準関数との差異
        writeln(format('%5.2f : %13.10f - %13.10f = %13.10f', [x, d1, d2, d1 - d2]));
    end;
end.
Program Pas0505(arg);
{$MODE delphi}

uses
    SysUtils, Math;

// 自作の指数関数
function myExp(x:Double; x2:Double; n:Integer; t:Double):Double;
begin
    t := x2 / (n + t);
    n := n - 4;

    if (n < 6) then
        result := 1 + ((2 * x) / (2 - x + t))
    else
        result := myExp(x, x2, n, t);
end;

var
    i:  Integer;
    x:  Double;
    x2: Double;
    d1: Double;
    d2: Double;
begin
    for i := 0 to 20 do
    begin
        x  := (i - 10) / 4.0;
        // 標準の指数関数
        d1 := Exp(x);
        // 自作の指数関数
        x2 := x * x;
        d2 := myExp(x, x2, 30, 0.0); // 30:必要な精度が得られるよう, 6から始めて4ずつ増加させる
        // 標準関数との差異
        writeln(format('%5.2f : %13.10f - %13.10f = %13.10f', [x, d1, d2, d1 - d2]));
    end;
end.
Program Pas0506(arg);
{$MODE delphi}

uses
    SysUtils, Math;

// 自作の対数関数
function myLog(x2:Double; numerator:Double; denominator:Double; y:Double):Double;
var
    a :Double;
begin
    denominator := denominator + 2;
    numerator   := numerator   * x2 * x2;
    a           := numerator / denominator;
    // 十分な精度になったら処理を抜ける
    if (Abs(a) <= 0.00000000001) then
        result := y
    else
        result := y + myLog(x2, numerator, denominator, a);
end;

var
    i:  Integer;
    x:  Double;
    x2: Double;
    d1: Double;
    d2: Double;
begin
    for i := 1 to 20 do
    begin
        x  := i / 5.0;
        // 標準の対数関数
        d1 := Ln(x);
        // 自作の対数関数
        x2 := (x - 1) / (x + 1);
        d2 := 2 * myLog(x2, x2, 1.0, x2);
        // 標準関数との差異
        writeln(format('%5.2f : %13.10f - %13.10f = %13.10f', [x, d1, d2, d1 - d2]));
    end;
end.
Program Pas0507(arg);
{$MODE delphi}

uses
    SysUtils, Math;

// 自作の対数関数
function myLog(x:Double; n:Integer; t:Double):Double;
var
    n2 :Integer;
    x2 :Double;
begin
    n2 := n;
    x2 := x;
    if (n > 3) then
    begin
        if (n Mod 2 = 0) then
            n2 := 2;
        x2 := x * (n Div 2);
    end;
    t := x2 / (n2 + t);

    if (n <= 2) then
        result := x / (1 + t)
    else
        result := myLog(x, n - 1, t);
end;

var
    i:  Integer;
    x:  Double;
    d1: Double;
    d2: Double;
begin
    for i := 1 to 20 do
    begin
        x  := i / 5.0;
        // 標準の対数関数
        d1 := Ln(x);
        // 自作の対数関数
        d2 := myLog(x - 1, 27, 0.0); // 27:必要な精度が得られる十分大きな奇数
        // 標準関数との差異
        writeln(format('%5.2f : %13.10f - %13.10f = %13.10f', [x, d1, d2, d1 - d2]));
    end;
end.
Program Pas0508(arg);
{$MODE delphi}

uses
    SysUtils, Math;

// 自作の双曲線正弦関数
function mySinh(x:Double; n:Integer; numerator:Double; denominator:Double; y:Double):Double;
var
    m: Integer;
    a: Double;
begin
    m           := 2 * n;
    denominator := denominator * (m + 1) * m;
    numerator   := numerator   * x * x;
    a           := numerator / denominator;
    // 十分な精度になったら処理を抜ける
    if (Abs(a) <= 0.00000000001) then
        result := y
    else
        result := y + mySinh(x, n + 1, numerator, denominator, a);
end;

var
    i:  Integer;
    x:  Integer;
    d1: Double;
    d2: Double;
begin
    for i := 0 to 20 do
    begin
        x  := i - 10;
        // 自作の双曲線正弦関数
        d1 := mySinh(x, 1, x, 1.0, x);
        // 標準の双曲線正弦関数
        d2 := Sinh(x);
        // 標準関数との差異
        writeln(format('%3d : %17.10f - %17.10f = %13.10f', [x, d1, d2, d1 - d2]));
    end;
end.
Program Pas0509(arg);
{$MODE delphi}

uses
    SysUtils, Math;

// 自作の双曲線余弦関数
function myCosh(x:Double; n:Integer; numerator:Double; denominator:Double; y:Double):Double;
var
    m: Integer;
    a: Double;
begin
    m           := 2 * n;
    denominator := denominator * m * (m - 1);
    numerator   := numerator   * x * x;
    a           := numerator / denominator;
    // 十分な精度になったら処理を抜ける
    if (Abs(a) <= 0.00000000001) then
        result := y
    else
        result := y + myCosh(x, n + 1, numerator, denominator, a);
end;

var
    i:  Integer;
    x:  Integer;
    d1: Double;
    d2: Double;
begin
    for i := 0 to 20 do
    begin
        x  := i - 10;
        // 自作の双曲線余弦関数
        d1 := myCosh(x, 1, 1.0, 1.0, 1.0);
        // 標準の双曲線余弦関数
        d2 := Cosh(x);
        // 標準関数との差異
        writeln(format('%3d : %17.10f - %17.10f = %13.10f', [x, d1, d2, d1 - d2]));
    end;
end.
Program Pas0601(arg);
{$MODE delphi}

uses
    SysUtils, Math;

function f(x:Double):Double;
begin
    result := 4 / (1 + x * x);
end;

const
    a:Double = 0;
    b:Double = 1;
var
    n, i, j:Integer;
    h, s, x:Double;
begin
    // 台形則で積分
    n := 2;
    for j := 1 to 10 do
    begin
        h := (b - a) / n;
        s := 0;
        x := a;
        for i := 1 to n - 1 do
        begin
            x := x + h;
            s := s + f(x);
        end;
        s := h * ((f(a) + f(b)) / 2 + s);
        n := n * 2;

        // 結果を π と比較
        writeln(format('%2d : %13.10f, %13.10f', [j, s, s - PI]));
    end;
end.
Program Pas0602(arg);
{$MODE delphi}

uses
    SysUtils, Math;

function f(x:Double):Double;
begin
    result := 4 / (1 + x * x);
end;

const
    a:Double = 0;
    b:Double = 1;
var
    n, i, j:Integer;
    h, s, x:Double;
begin
    // 中点則で積分
    n := 2;
    for j := 1 to 10 do
    begin
        h := (b - a) / n;
        s := 0;
        x := a + (h / 2);
        for i := 1 to n do
        begin
            s := s + f(x);
            x := x + h;
        end;
        s := h * s;
        n := n * 2;

        // 結果を π と比較
        writeln(format('%2d : %13.10f, %13.10f', [j, s, s - PI]));
    end;
end.
Program Pas0603(arg);
{$MODE delphi}

uses
    SysUtils, Math;

function f(x:Double):Double;
begin
    result := 4 / (1 + x * x);
end;

const
    a:Double = 0;
    b:Double = 1;
var
    n, i, j:Integer;
    h, s, x:Double;
    s2, s4 :Double;
begin
    // Simpson則で積分
    n := 2;
    for j := 1 to 5 do
    begin
        h  := (b - a) / n;
        s2 := 0;
        s4 := 0;
        x  := a + h;
        for i := 1 to (n div 2) do
        begin
            s4 := s4 + f(x);
            x  := x  + h;
            s2 := s2 + f(x);
            x  := x  + h;
        end;
        s2 := (s2 - f(b)) * 2 + f(a) + f(b);
        s4 := s4 * 4;
        s  := (s2 + s4) * h / 3;
        n  := n * 2;

        // 結果を π と比較
        writeln(format('%2d : %13.10f, %13.10f', [j, s, s - PI]));
    end;
end.
Program Pas0604(arg);
{$MODE delphi}

uses
    SysUtils, Math;

function f(x:Double):Double;
begin
    result := 4 / (1 + x * x);
end;

const
    a:Double = 0;
    b:Double = 1;
var
    n, i, j:Integer;
    h, s, x:Double;
    t:array [1..6, 1..6] of Double;
begin
    // 台形則で積分
    n := 2;
    for i := 1 to 6 do
    begin
        h := (b - a) / n;
        s := 0;
        x := a;
        for j := 1 to n - 1 do
        begin
            x := x + h;
            s := s + f(x);
        end;
        // 結果を保存
        t[i,1] := h * ((f(a) + f(b)) / 2 + s);
        n      := n * 2;
    end;

    // Richardsonの補外法
    n := 4;
    for j := 2 to 6 do
    begin
        for i := j to 6 do
        begin
            t[i,j] := t[i, j - 1] + (t[i, j - 1] - t[i - 1, j - 1]) / (n - 1);
            if i = j then
            begin
                // 結果を π と比較
                writeln(format('%2d : %13.10f, %13.10f', [j, t[i,j], t[i,j] - PI]));
            end;
        end;
        n := n * 4;
    end;
end.
Program Pas0701(arg);
{$MODE delphi}

uses
    SysUtils, Math;

// 元の関数
function f(x:Double):Double;
begin
    result := x - power(x,3) / (3 * 2) + power(x,5) / (5 * 4 * 3 * 2);
end;

// Lagrange (ラグランジュ) 補間
function lagrange(d:Double; x:array of Double; y:array of Double):Double;
var
    sum, prod :Double;
    i, j      :Integer;
begin
    sum := 0;
    for i := Low(x) to High(x) do
    begin
        prod := y[i];
        for j := Low(x) to High(x) do
        begin
            if j <> i then
                prod := prod * (d - x[j]) / (x[i] - x[j]);
        end;
        sum := sum + prod;
    end;
    result := sum;
end;

const
    // データ点の数 - 1
    N = 6;
var
    i :Integer;
    x :array [0..N] of Double;
    y :array [0..N] of Double;
    d, d1, d2 :Double;
begin
    // 1.5刻みで -4.5〜4.5 まで, ７点だけ値をセット
    for i := Low(x) to High(x) do
    begin
        d    := (i - 1) * 1.5 - 4.5;
        x[i] := d;
        y[i] := f(d);
    end;

    // 0.5刻みで 与えられていない値を補間
    for i := 0 to 18 do
    begin
        d  := i * 0.5 - 4.5;
        d1 := f(d);
        d2 := lagrange(d, x, y);

        // 元の関数と比較
        writeln(format('%5.2f'#9'%8.5f'#9'%8.5f'#9'%8.5f', [d, d1, d2, d1 - d2]));
    end;
end.
Program Pas0702(arg);
{$MODE delphi}

uses
    SysUtils, Math;

const
    // データ点の数 - 1
    N = 6;

// 元の関数
function f(x:Double):Double;
begin
    result := x - power(x,3) / (3 * 2) + power(x,5) / (5 * 4 * 3 * 2);
end;

// Neville (ネヴィル) 補間
function neville(d:Double; x:array of Double; y:array of Double):Double;
var
    w :array [0..N, 0..N] of Double;
    i, j      :Integer;
begin
    for i := Low(y) to High(y) do
        w[0,i] := y[i];

    for j := 1 to High(x) do
    begin
        for i := Low(x) to (High(x) - j) do
            w[j,i] := w[j-1,i+1] + (w[j-1,i+1] - w[j-1,i]) * (d - x[i+j]) / (x[i+j] - x[i]);
    end;
    result := w[N,0];
end;

var
    i :Integer;
    x :array [0..N] of Double;
    y :array [0..N] of Double;
    d, d1, d2 :Double;
begin
    // 1.5刻みで -4.5〜4.5 まで, ７点だけ値をセット
    for i := Low(x) to High(x) do
    begin
        d    := i * 1.5 - 4.5;
        x[i] := d;
        y[i] := f(d);
    end;

    // 0.5刻みで 与えられていない値を補間
    for i := 0 to 18 do
    begin
        d  := i * 0.5 - 4.5;
        d1 := f(d);
        d2 := neville(d, x, y);

        // 元の関数と比較
        writeln(format('%5.2f'#9'%8.5f'#9'%8.5f'#9'%8.5f', [d, d1, d2, d1 - d2]));
    end;
end.
Program Pas0703(arg);
{$MODE delphi}

uses
    SysUtils, Math;

const
    // データ点の数 - 1
    N = 6;

// 元の関数
function f(x:Double):Double;
begin
    result := x - power(x,3) / (3 * 2) + power(x,5) / (5 * 4 * 3 * 2);
end;

// Newton (ニュートン) 補間
function newton(d:Double; x:array of Double; a:array of Double):Double;
var
    sum, prod :Double;
    i, j      :Integer;
begin
    sum := a[0];
    for i := 1 to High(x) do
    begin
        prod := a[i];
        for j := Low(x) to i-1 do
            prod := prod * (d - x[j]);
        sum := sum + prod;
    end;
    result := sum;
end;

var
    i, j :Integer;
    x :array [0..N] of Double;
    y :array [0..N] of Double;
    a :array [0..N] of Double;
    d :array [0..N, 0..N] of Double;
    d1, d2, d3 :Double;
begin
    // 1.5刻みで -4.5〜4.5 まで, ７点だけ値をセット
    for i := Low(x) to High(x) do
    begin
        d1   := i * 1.5 - 4.5;
        x[i] := d1;
        y[i] := f(d1);
    end;

    // 差分商の表を作る
    for j := Low(y) to High(y) do
        d[0,j] := y[j];

    for i := 1 to High(x) do
    begin
        for j := Low(x) to High(x)-i do
            d[i,j] := (d[i-1,j+1] - d[i-1,j]) / (x[j+i] - x[j]);
    end;

    // ｎ階差分商
    for j := Low(a) to High(a) do
        a[j] := d[j,0];

    // 0.5刻みで 与えられていない値を補間
    for i := 0 to 18 do
    begin
        d1 := i * 0.5 - 4.5;
        d2 := f(d1);
        d3 := newton(d1, x, a);

        // 元の関数と比較
        writeln(format('%5.2f'#9'%8.5f'#9'%8.5f'#9'%8.5f', [d1, d2, d3, d2 - d3]));
    end;
end.
Program Pas0704(arg);
{$MODE delphi}

uses
    SysUtils, Math;

const
    // データ点の数 - 1
    N   =  6;
    Nx2 = 13;

// 元の関数
function f(x:Double):Double;
begin
    result := x - power(x,3) / (3 * 2) + power(x,5) / (5 * 4 * 3 * 2);
end;

// 導関数
function fd(x:Double):Double;
begin
    result := 1 - power(x,2) / 2 + power(x,4) / (4 * 3 * 2);
end;

// Hermite (エルミート) 補間
function hermite(d:Double; z:array of Double; a:array of Double):Double;
var
    sum, prod :Double;
    i, j      :Integer;
begin
    sum := a[0];
    for i := 1 to High(a) do
    begin
        prod := a[i];
        for j := Low(z) to i-1 do
            prod := prod * (d - z[j]);
        sum := sum + prod;
    end;
    result := sum;
end;

var
    i, j :Integer;
    x  :array [0..N] of Double;
    y  :array [0..N] of Double;
    yd :array [0..N] of Double;
    z  :array [0..Nx2] of Double;
    a  :array [0..Nx2] of Double;
    d  :array [0..Nx2, 0..Nx2] of Double;
    d1, d2, d3 :Double;
begin
    // 1.5刻みで -4.5〜4.5 まで, ７点だけ値をセット
    for i := Low(x) to High(x) do
    begin
        d1    := i * 1.5 - 4.5;
        x[i]  := d1;
        y[i]  := f(d1);
        yd[i] := fd(d1);
    end;

    // 差分商の表を作る
    for i := Low(z) to High(z) do
    begin
        j      := i div 2;
        z[i]   := x[j];
        d[0,i] := y[j];
    end;

    for i := 1 to High(z) do
    begin
        for j := Low(z) to High(z)-i do
        begin
            if (i = 1) and (j mod 2 = 0) then
                d[i,j] := yd[j div 2]
            else
                d[i,j] := (d[i-1,j+1] - d[i-1,j]) / (z[j+i] - z[j]);
        end;
    end;

    // ｎ階差分商
    for j := Low(a) to High(a) do
        a[j] := d[j,0];

    // 0.5刻みで 与えられていない値を補間
    for i := 0 to 18 do
    begin
        d1 := i * 0.5 - 4.5;
        d2 := f(d1);
        d3 := hermite(d1, z, a);

        // 元の関数と比較
        writeln(format('%5.2f'#9'%8.5f'#9'%8.5f'#9'%8.5f', [d1, d2, d3, d2 - d3]));
    end;
end.
Program Pas0705(arg);
{$MODE delphi}

uses
    SysUtils, Math;

const
    // データ点の数 - 1
    N = 6;

// 元の関数
function f(x:Double):Double;
begin
    result := x - power(x,3) / (3 * 2) + power(x,5) / (5 * 4 * 3 * 2);
end;

// Spline (スプライン) 補間
function spline(d:Double; x:array of Double; y:array of Double; z:array of Double):Double;
var
    d1, d2, d3 :Double;
    i, k :Integer;
begin
    // 補間関数値がどの区間にあるか
    k := -1;
    for i := 1 to High(x) do
    begin
        if d <= x[i] then
        begin
            k := i - 1;
            break;
        end;
    end;
    if k < 0 then
        k := i - 1;

    d1     := x[k+1] - d;
    d2     := d      - x[k];
    d3     := x[k+1] - x[k];
    result := (z[k] * power(d1,3) + z[k+1] * power(d2,3)) / (6.0 * d3)
            + (y[k]   / d3 - z[k]   * d3 / 6.0) * d1
            + (y[k+1] / d3 - z[k+1] * d3 / 6.0) * d2;
end;

var
    i, j :Integer;
    x  :array [0..N] of Double;
    y  :array [0..N] of Double;
    a  :array [0..N] of Double;
    b  :array [0..N] of Double;
    c  :array [0..N] of Double;
    d  :array [0..N] of Double;
    g  :array [0..N] of Double;
    s  :array [0..N] of Double;
    z  :array [0..N] of Double;
    d1, d2, d3 :Double;
begin
    // 1.5刻みで -4.5〜4.5 まで, ７点だけ値をセット
    for i := Low(x) to High(x) do
    begin
        d1   := i * 1.5 - 4.5;
        x[i] := d1;
        y[i] := f(d1);
    end;

    // ３項方程式の係数の表を作る
    for i := 1 to High(x)-1 do
    begin
        a[i] :=         x[i]   - x[i-1];
        b[i] := 2.0 *  (x[i+1] - x[i-1]);
        c[i] :=         x[i+1] - x[i];
        d[i] := 6.0 * ((y[i+1] - y[i]) / (x[i+1] - x[i]) - (y[i] - y[i-1]) / (x[i] - x[i-1]));
    end;
    // ３項方程式を解く (ト−マス法)
    g[1] := b[1];
    s[1] := d[1];
    for i := 2 to High(a)-1 do
    begin
        g[i] := b[i] - a[i] * c[i-1] / g[i-1];
        s[i] := d[i] - a[i] * s[i-1] / g[i-1];
    end;
    z[0]   := 0;
    z[N]   := 0;
    z[N-1] := s[N-1] / g[N-1];
    for i := N-2 downto 1 do
    begin
        z[i] := (s[i] - c[i] * z[i+1]) / g[i];
    end;

    // 0.5刻みで 与えられていない値を補間
    for i := 0 to 18 do
    begin
        d1 := i * 0.5 - 4.5;
        d2 := f(d1);
        d3 := spline(d1, x, y, z);

        // 元の関数と比較
        writeln(format('%5.2f'#9'%8.5f'#9'%8.5f'#9'%8.5f', [d1, d2, d3, d2 - d3]));
    end;
end.
Program Pas0801(arg);
{$MODE delphi}

uses
    SysUtils, Math;

const
    // 重力加速度
    g = -9.8;
    // 空気抵抗係数
    k = -0.01;
    // 時間間隔(秒)
    h = 0.01;
    // 角度
    degree = 45;

// 空気抵抗による水平方向の減速分
function fx(vx:Double; vy:Double):Double;
begin
    result := k * Sqrt(vx * vx + vy * vy) * vx;
end;

// 重力と空気抵抗による鉛直方向の減速分
function fy(vx:Double; vy:Double):Double;
begin
    result := g + (k * Sqrt(vx * vx + vy * vy) * vy);
end;

var
    // 角度
    radian:Double;
    // 初速
    v:Double;
    // 水平方向の速度
    vx:array [0..1] of Double;
    // 鉛直方向の速度
    vy:array [0..1] of Double;
    // 経過秒数
    t:Double = 0.0;
    // 位置
    x:Double = 0.0;
    y:Double = 0.0;

    i:Integer;
begin
    // 角度
    radian := degree * PI / 180.0;
    // 初速 250 km/h -> 秒速に変換
    v := 250 * 1000 div 3600;
    // 水平方向の速度
    vx[0] := v * Cos(radian);
    // 鉛直方向の速度
    vy[0] := v * Sin(radian);

    // Euler法
    i := 1;
    while y >= 0.0 do
    begin
        // 経過秒数
        t := i * h;

        // 位置
        x := x + h * vx[0];
        y := y + h * vy[0];

        writeln(format('%4.2f'#9'%8.5f'#9'%9.5f'#9'%9.5f'#9'%8.5f', [t, vx[0], vy[0], x, y]));

        // 速度
        vx[1] := vx[0] + h * fx(vx[0], vy[0]);
        vy[1] := vy[0] + h * fy(vx[0], vy[0]);
        vx[0] := vx[1];
        vy[0] := vy[1];

        inc(i);
    end;
end.
Program Pas0802(arg);
{$MODE delphi}

uses
    SysUtils, Math;

const
    // 重力加速度
    g = -9.8;
    // 空気抵抗係数
    k = -0.01;
    // 時間間隔(秒)
    h = 0.01;
    // 角度
    degree = 45;

// 空気抵抗による水平方向の減速分
function fx(vx:Double; vy:Double):Double;
begin
    result := k * Sqrt(vx * vx + vy * vy) * vx;
end;

// 重力と空気抵抗による鉛直方向の減速分
function fy(vx:Double; vy:Double):Double;
begin
    result := g + (k * Sqrt(vx * vx + vy * vy) * vy);
end;

var
    // 角度
    radian:Double;
    // 初速
    v:Double;
    // 水平方向の速度
    vx:array [0..2] of Double;
    // 鉛直方向の速度
    vy:array [0..2] of Double;
    // 経過秒数
    t:Double = 0.0;
    // 位置
    x:array [0..2] of Double;
    y:array [0..2] of Double;

    i:Integer;
begin
    // 角度
    radian := degree * PI / 180.0;
    // 初速 250 km/h -> 秒速に変換
    v := 250 * 1000 div 3600;
    // 水平方向の速度
    vx[0] := v * Cos(radian);
    // 鉛直方向の速度
    vy[0] := v * Sin(radian);
    // 位置
    x[0] := 0.0;
    y[0] := 0.0;

    // Heun法
    i := 1;
    while y[0] >= 0.0 do
    begin
        // 経過秒数
        t := i * h;

        // 位置・速度
        x[1]  :=  x[0] + h *    vx[0];
        y[1]  :=  y[0] + h *    vy[0];
        vx[1] := vx[0] + h * fx(vx[0], vy[0]);
        vy[1] := vy[0] + h * fy(vx[0], vy[0]);

        x[2]  :=  x[0] + h * (  vx[0]          +    vx[1]        ) / 2;
        y[2]  :=  y[0] + h * (  vy[0]          +    vy[1]        ) / 2;
        vx[2] := vx[0] + h * (fx(vx[0], vy[0]) + fx(vx[1], vy[1])) / 2;
        vy[2] := vy[0] + h * (fy(vx[0], vy[0]) + fy(vx[1], vy[1])) / 2;

        x[0]  :=  x[2];
        y[0]  :=  y[2];
        vx[0] := vx[2];
        vy[0] := vy[2];

        writeln(format('%4.2f'#9'%8.5f'#9'%9.5f'#9'%9.5f'#9'%8.5f', [t, vx[0], vy[0], x[0], y[0]]));

        inc(i);
    end;
end.
Program Pas0803(arg);
{$MODE delphi}

uses
    SysUtils, Math;

const
    // 重力加速度
    g = -9.8;
    // 空気抵抗係数
    k = -0.01;
    // 時間間隔(秒)
    h = 0.01;
    // 角度
    degree = 45;

// 空気抵抗による水平方向の減速分
function fx(vx:Double; vy:Double):Double;
begin
    result := k * Sqrt(vx * vx + vy * vy) * vx;
end;

// 重力と空気抵抗による鉛直方向の減速分
function fy(vx:Double; vy:Double):Double;
begin
    result := g + (k * Sqrt(vx * vx + vy * vy) * vy);
end;

var
    // 角度
    radian:Double;
    // 初速
    v:Double;
    // 水平方向の速度
    vx:array [0..1] of Double;
    wx:Double;
    // 鉛直方向の速度
    vy:array [0..1] of Double;
    wy:Double;
    // 経過秒数
    t:Double = 0.0;
    // 位置
    x:array [0..1] of Double;
    y:array [0..1] of Double;

    i:Integer;
begin
    // 角度
    radian := degree * PI / 180.0;
    // 初速 250 km/h -> 秒速に変換
    v := 250 * 1000 div 3600;
    // 水平方向の速度
    vx[0] := v * Cos(radian);
    // 鉛直方向の速度
    vy[0] := v * Sin(radian);
    // 位置
    x[0] := 0.0;
    y[0] := 0.0;

    // 中点法
    i := 1;
    while y[0] >= 0.0 do
    begin
        // 経過秒数
        t := i * h;

        // 位置・速度
        vx[1] := h * fx(vx[0], vy[0]);
        vy[1] := h * fy(vx[0], vy[0]);

        wx    := vx[0] + vx[1] / 2;
        wy    := vy[0] + vy[1] / 2;
        vx[0] := vx[0] + h * fx(wx, wy);
        vy[0] := vy[0] + h * fy(wx, wy);
        x[0]  := x[0]  + h *    wx;
        y[0]  := y[0]  + h *    wy;

        writeln(format('%4.2f'#9'%8.5f'#9'%9.5f'#9'%9.5f'#9'%9.5f', [t, vx[0], vy[0], x[0], y[0]]));

        inc(i);
    end;
end.
Program Pas0804(arg);
{$MODE delphi}

uses
    SysUtils, Math;

const
    // 重力加速度
    g = -9.8;
    // 空気抵抗係数
    k = -0.01;
    // 時間間隔(秒)
    h = 0.01;
    // 角度
    degree = 45;

// 空気抵抗による水平方向の減速分
function fx(vx:Double; vy:Double):Double;
begin
    result := k * Sqrt(vx * vx + vy * vy) * vx;
end;

// 重力と空気抵抗による鉛直方向の減速分
function fy(vx:Double; vy:Double):Double;
begin
    result := g + (k * Sqrt(vx * vx + vy * vy) * vy);
end;

var
    // 角度
    radian:Double;
    // 初速
    v:Double;
    // 水平方向の速度
    vx:array [0..4] of Double;
    wx:Double;
    // 鉛直方向の速度
    vy:array [0..4] of Double;
    wy:Double;
    // 経過秒数
    t:Double = 0.0;
    // 位置
    x:array [0..4] of Double;
    y:array [0..4] of Double;

    i:Integer;
begin
    // 角度
    radian := degree * PI / 180.0;
    // 初速 250 km/h -> 秒速に変換
    v := 250 * 1000 div 3600;
    // 水平方向の速度
    vx[0] := v * Cos(radian);
    // 鉛直方向の速度
    vy[0] := v * Sin(radian);
    // 位置
    x[0] := 0.0;
    y[0] := 0.0;

    // Runge-Kutta法
    i := 1;
    while y[0] >= 0.0 do
    begin
        // 経過秒数
        t := i * h;

        // 位置・速度
        x[1]  := h *    vx[0];
        y[1]  := h *    vy[0];
        vx[1] := h * fx(vx[0], vy[0]);
        vy[1] := h * fy(vx[0], vy[0]);

        wx    := vx[0] + vx[1] / 2;
        wy    := vy[0] + vy[1] / 2;
        x[2]  := h *    wx;
        y[2]  := h *    wy;
        vx[2] := h * fx(wx, wy);
        vy[2] := h * fy(wx, wy);

        wx    := vx[0] + vx[2] / 2;
        wy    := vy[0] + vy[2] / 2;
        x[3]  := h *    wx;
        y[3]  := h *    wy;
        vx[3] := h * fx(wx, wy);
        vy[3] := h * fy(wx, wy);

        wx    := vx[0] + vx[3];
        wy    := vy[0] + vy[3];
        x[4]  := h *    wx;
        y[4]  := h *    wy;
        vx[4] := h * fx(wx, wy);
        vy[4] := h * fy(wx, wy);

        x[0]  := x[0]  + ( x[1] +  x[2] * 2 +  x[3] * 2 +  x[4]) / 6;
        y[0]  := y[0]  + ( y[1] +  y[2] * 2 +  y[3] * 2 +  y[4]) / 6;
        vx[0] := vx[0] + (vx[1] + vx[2] * 2 + vx[3] * 2 + vx[4]) / 6;
        vy[0] := vy[0] + (vy[1] + vy[2] * 2 + vy[3] * 2 + vy[4]) / 6;

        writeln(format('%4.2f'#9'%8.5f'#9'%9.5f'#9'%9.5f'#9'%9.5f', [t, vx[0], vy[0], x[0], y[0]]));

        inc(i);
    end;
end.
Program Pas0805(arg);
{$MODE delphi}

uses
    SysUtils, Math;

const
    // 重力加速度
    g = -9.8;
    // 空気抵抗係数
    k = -0.01;
    // 時間間隔(秒)
    h = 0.01;
    // 角度
    degree = 45;

// 空気抵抗による水平方向の減速分
function fx(vx:Double; vy:Double):Double;
begin
    result := k * Sqrt(vx * vx + vy * vy) * vx;
end;

// 重力と空気抵抗による鉛直方向の減速分
function fy(vx:Double; vy:Double):Double;
begin
    result := g + (k * Sqrt(vx * vx + vy * vy) * vy);
end;

var
    // 角度
    radian:Double;
    // 初速
    v:Double;
    // 水平方向の速度
    vx:array [0..4] of Double;
    wx:Double;
    // 鉛直方向の速度
    vy:array [0..4] of Double;
    wy:Double;
    // 経過秒数
    t:Double = 0.0;
    // 位置
    x:array [0..4] of Double;
    y:array [0..4] of Double;

    i:Integer;
begin
    // 角度
    radian := degree * PI / 180.0;
    // 初速 250 km/h -> 秒速に変換
    v := 250 * 1000 div 3600;
    // 水平方向の速度
    vx[0] := v * Cos(radian);
    // 鉛直方向の速度
    vy[0] := v * Sin(radian);
    // 位置
    x[0] := 0.0;
    y[0] := 0.0;

    // Runge-Kutta-Gill法
    i := 1;
    while y[0] >= 0.0 do
    begin
        // 経過秒数
        t := i * h;

        // 位置・速度
        x[1]  := h *    vx[0];
        y[1]  := h *    vy[0];
        vx[1] := h * fx(vx[0], vy[0]);
        vy[1] := h * fy(vx[0], vy[0]);

        wx    := vx[0] + vx[1] / 2;
        wy    := vy[0] + vy[1] / 2;
        x[2]  := h *    wx;
        y[2]  := h *    wy;
        vx[2] := h * fx(wx, wy);
        vy[2] := h * fy(wx, wy);

        wx    := vx[0] + vx[1] * ((Sqrt(2.0) - 1) / 2) + vx[2] * (1 - 1 / Sqrt(2.0));
        wy    := vy[0] + vy[1] * ((Sqrt(2.0) - 1) / 2) + vy[2] * (1 - 1 / Sqrt(2.0));
        x[3]  := h *    wx;
        y[3]  := h *    wy;
        vx[3] := h * fx(wx, wy);
        vy[3] := h * fy(wx, wy);

        wx    := vx[0] - vx[2] / Sqrt(2.0) + vx[3] * (1 + 1 / Sqrt(2.0));
        wy    := vy[0] - vy[2] / Sqrt(2.0) + vy[3] * (1 + 1 / Sqrt(2.0));
        x[4]  := h *    wx;
        y[4]  := h *    wy;
        vx[4] := h * fx(wx, wy);
        vy[4] := h * fy(wx, wy);

        x[0]  := x[0]  + ( x[1] +  x[2] * (2 - Sqrt(2.0)) +  x[3] * (2 + Sqrt(2.0)) +  x[4]) / 6;
        y[0]  := y[0]  + ( y[1] +  y[2] * (2 - Sqrt(2.0)) +  y[3] * (2 + Sqrt(2.0)) +  y[4]) / 6;
        vx[0] := vx[0] + (vx[1] + vx[2] * (2 - Sqrt(2.0)) + vx[3] * (2 + Sqrt(2.0)) + vx[4]) / 6;
        vy[0] := vy[0] + (vy[1] + vy[2] * (2 - Sqrt(2.0)) + vy[3] * (2 + Sqrt(2.0)) + vy[4]) / 6;

        writeln(format('%4.2f'#9'%8.5f'#9'%9.5f'#9'%9.5f'#9'%9.5f', [t, vx[0], vy[0], x[0], y[0]]));

        inc(i);
    end;
end.
Program Pas0901(arg);
{$MODE delphi}

uses
    SysUtils, Math;

function f(x:Double):Double;
begin
    result := x * x - 2.0;
end;

function bisection(a:Double; b:Double):Double;
var
    c:  Double;
    fc: Double;
begin
    while true do
    begin
        { 区間 (a, b) の中点 c = (a + b) / 2 }
        c := (a + b) / 2;
        writeln(format('%12.10f'#9'%13.10f', [c, c - sqrt(2)]));

        fc := f(c);
        If Abs(fc) < 0.0000000001 then break;

        if fc < 0 then
        begin
            { f(c) < 0 であれば, 解は区間 (c, b) の中に存在 }
            a := c;
        end
        else
        begin
            { f(c) > 0 であれば, 解は区間 (a, c) の中に存在 }
            b := c;
        end;
    end;

    result := c;
end;

var
    a:  Double = 1.0;
    b:  Double = 2.0;
begin
    writeln(format('%12.10f', [bisection(a, b)]));
end.
Program Pas0902(arg);
{$MODE delphi}

uses
    SysUtils, Math;

function f(x:Double):Double;
begin
    result := x * x - 2.0;
end;

function falseposition(a:Double; b:Double):Double;
var
    c:  Double;
    fc: Double;
begin
    while true do
    begin
        { 点 (a,f(a)) と 点 (b,f(b)) を結ぶ直線と x軸の交点 }
        c := (a * f(b) - b * f(a)) / (f(b) - f(a));
        writeln(format('%12.10f'#9'%13.10f', [c, c - sqrt(2)]));

        fc := f(c);
        If Abs(fc) < 0.0000000001 then break;

        if fc < 0 then
        begin
            { f(c) < 0 であれば, 解は区間 (c, b) の中に存在 }
            a := c;
        end
        else
        begin
            { f(c) > 0 であれば, 解は区間 (a, c) の中に存在 }
            b := c;
        end;
    end;

    result := c;
end;

var
    a:  Double = 1.0;
    b:  Double = 2.0;
begin
    writeln(format('%12.10f', [falseposition(a, b)]));
end.
Program Pas0903(arg);
{$MODE delphi}

uses
    SysUtils, Math;

function g(x:Double):Double;
begin
    result := (x / 2) + (1 / x);
end;

function iterative(x0:Double):Double;
var
    x1: Double;
begin
    while true do
    begin
        x1 := g(x0);
        writeln(format('%12.10f'#9'%13.10f', [x1, x1 - sqrt(2)]));

        If Abs(x1 - x0) < 0.0000000001 then break;
        x0 := x1;
    end;

    result := x1;
end;

var
    x:  Double = 1.0;
begin
    writeln(format('%12.10f', [iterative(x)]));
end.
Program Pas0904(arg);
{$MODE delphi}

uses
    SysUtils, Math;

function f0(x:Double):Double;
begin
    result := x * x - 2;
end;

function f1(x:Double):Double;
begin
    result := 2 * x;
end;

function newton(x0:Double):Double;
var
    x1: Double;
begin
    while true do
    begin
        x1 := x0 - (f0(x0) / f1(x0));
        writeln(format('%12.10f'#9'%13.10f', [x1, x1 - sqrt(2)]));

        If Abs(x1 - x0) < 0.0000000001 then break;
        x0 := x1;
    end;

    result := x1;
end;

var
    x:  Double = 2.0;
begin
    writeln(format('%12.10f', [newton(x)]));
end.
Program Pas0905(arg);
{$MODE delphi}

uses
    SysUtils, Math;

function f0(x:Double):Double;
begin
    result := x * x - 2;
end;

function f1(x:Double):Double;
begin
    result := 2 * x;
end;

function f2(x:Double):Double;
begin
    result := 2;
end;

function bailey(x0:Double):Double;
var
    x1: Double;
begin
    while true do
    begin
        x1 := x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2 * f1(x0)))));
        writeln(format('%12.10f'#9'%13.10f', [x1, x1 - sqrt(2)]));

        If Abs(x1 - x0) < 0.0000000001 then break;
        x0 := x1;
    end;

    result := x1;
end;

var
    x:  Double = 2.0;
begin
    writeln(format('%12.10f', [bailey(x)]));
end.
Program Pas0906(arg);
{$MODE delphi}

uses
    SysUtils, Math;

function f(x:Double):Double;
begin
    result := x * x - 2;
end;

function secant(x0:Double; x1:Double):Double;
var
    x2: Double;
begin
    while true do
    begin
        x2 := x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0));
        writeln(format('%12.10f'#9'%13.10f', [x1, x1 - sqrt(2)]));

        If Abs(x1 - x0) < 0.0000000001 then break;
        x0 := x1;
        x1 := x2;
    end;

    result := x2;
end;

var
    x0:  Double = 1.0;
    x1:  Double = 2.0;
begin
    writeln(format('%12.10f', [secant(x0, x1)]));
end.