home > さまざまな言語で数値計算 > 関数の近似 >

さまざまな言語で数値計算

Only Do What Only You Can Do

ネヴィル補間

与えられた4個の関数値 $ f(x_0), f(x_1), f(x_2), f(x_3) $ を通る3次式を求める場合, まず次のような表を考える.

このとき,

と定義し, 同様に

と定義すると, 与えられた $n$点を通る $n-1$次式は次のように表すことができる.

この式を使って, 与えられた点以外の点の値を求める.

例題として,

を近似する.

VBScript

Option Explicit

'データ点の数 - 1
Private Const N = 6

Dim x(): ReDim x(N)
Dim y(): ReDim y(N)

'1.5刻みで -4.5~4.5 まで, 7点だけ値をセット
Dim i
For i = 0 To N
    Dim d: d = i * 1.5 - 4.5
    x(i) = d
    y(i) = f(d)
Next

'0.5刻みで 与えられていない値を補間
For i = 0 To 18
    d = i * 0.5 - 4.5
    Dim d1: d1 = f(d)
    Dim d2: d2 = neville(d, x, y)

    '元の関数と比較
    WScript.StdOut.Write     Right(Space(5) & FormatNumber(d,       2, -1, 0, 0), 5) & vbTab
    WScript.StdOut.Write     Right(Space(8) & FormatNumber(d1,      5, -1, 0, 0), 8) & vbTab
    WScript.StdOut.Write     Right(Space(8) & FormatNumber(d2,      5, -1, 0, 0), 8) & vbTab
    WScript.StdOut.WriteLine Right(Space(8) & FormatNumber(d1 - d2, 5, -1, 0, 0), 8)
Next

'元の関数
Private Function f(ByVal x)
    f = x - (x ^ 3) / (3 * 2) + (x ^ 5) / (5 * 4 * 3 * 2)
End Function

'Neville (ネヴィル) 補間
Private Function neville(ByVal d, ByVal x(), ByVal y())
    Dim w(): ReDim w(N, N)
    Dim i
    For i = 0 To N
        w(0,i) = y(i)
    Next

    Dim j
    For j = 1 To N
        For i = 0 To N - j
            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))
        Next
    Next

    neville = w(N,0)
End Function
Z:\>cscript //nologo Z:\0702.vbs
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667    -0.00000
-3.50   -0.73099    -0.73099    -0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964    -0.00000
-2.00   -0.93333    -0.93333    -0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.47943     0.00000
 1.00    0.84167     0.84167    -0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984     0.00000

JScript

// データ点の数
var N = 7

var x = []
var y = []

// 1.5刻みで -4.5~4.5 まで, 7点だけ値をセット
for (var i = 0; i < N; i++)
{
    var d = i * 1.5 - 4.5
    x[i] = d
    y[i] = f(d)
}

// 0.5刻みで 与えられていない値を補間
for (var i = 0; i <= 18; i++)
{
    var d1  = i * 0.5 - 4.5
    var d2 = f(d1)
    var d3 = neville(d1, x, y)

    // 元の関数と比較
    WScript.StdOut.Write(("     "    + d1.toFixed(2)       ).slice(-5) + "\t")
    WScript.StdOut.Write(("        " + d2.toFixed(5)       ).slice(-8) + "\t")
    WScript.StdOut.Write(("        " + d3.toFixed(5)       ).slice(-8) + "\t")
    WScript.StdOut.Write(("        " + (d2 - d3).toFixed(5)).slice(-8) + "\n")
}

// 元の関数
function f(x)
{
    return x - Math.pow(x,3) / (3 * 2) + Math.pow(x,5) / (5 * 4 * 3 * 2);
}

// Neville (ネヴィル) 補間
function neville(d, x, y)
{
    w    = []
    w[0] = []
    for (var i = 0; i < N; i++)
        w[0][i] = y[i]

    for (var j = 1; j < N; j++)
    {
        w[j] = []
        for (var i = 0; i < N - j; i++)
            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])
    }

    return w[N-1][0]
}
Z:\>cscript //nologo Z:\0702.js
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667    -0.00000
-3.50   -0.73099    -0.73099    -0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964    -0.00000
-2.00   -0.93333    -0.93333    -0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.47943     0.00000
 1.00    0.84167     0.84167    -0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984     0.00000

PowerShell

# データ点の数
set-variable -option constant -name N -value 7

# 元の関数
function f($x)
{
    $x - [Math]::Pow($x, 3) / (3 * 2) + [Math]::Pow($x, 5) / (5 * 4 * 3 * 2)
}

# Neville (ネヴィル) 補間
function neville($d, $x, $y)
{
    $w = New-Object "double[,]" $N,$N
    foreach ($i in 0..($N - 1))
    {
        $w[0,$i] = $y[$i]
    }
    foreach ($j in 1..($N - 1))
    {
        foreach ($i in 0..($N - $j - 1))
        {
            $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])
        }
    }
    $w[($N-1),0]
}

$x = New-Object double[] $N
$y = New-Object double[] $N

# 1.5刻みで -4.5~4.5 まで, 7点だけ値をセット
foreach ($i in 0..($N - 1))
{
    $d = $i * 1.5 - 4.5
    $x[$i] = $d
    $y[$i] = f($d)
}

# 0.5刻みで 与えられていない値を補間
foreach ($i in 0..18)
{
    $d  = $i * 0.5 - 4.5
    $d1 = f($d)
    $d2 = (neville $d $x $y)

    # 元の関数と比較
    Write-Host ([String]::Format("{0,5:F2}`t{1,8:F5}`t{2,8:F5}`t{3,8:F5}", $d, $d1, $d2, ($d1 - $d2)))
}
Z:\>powershell -file Z:\0702.ps1
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667     0.00000
-3.50   -0.73099    -0.73099     0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964     0.00000
-2.00   -0.93333    -0.93333     0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.47943     0.00000
 1.00    0.84167     0.84167     0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984     0.00000

Perl

# データ点の数 - 1
use constant N => 6;

my @x = ();
my @y = ();

# 1.5刻みで -4.5~4.5 まで, 7点だけ値をセット
for $i (0..N)
{
    my $d = $i * 1.5 - 4.5;
    $x[$i] = $d;
    $y[$i] = f($d);
}

# 0.5刻みで 与えられていない値を補間
for $i (0..18)
{
    my $d  = $i * 0.5 - 4.5;
    my $d1 = f($d);
    my $d2 = neville($d, \@x, \@y);

    # 元の関数と比較
    printf("%5.2f\t%8.5f\t%8.5f\t%8.5f\n", $d, $d1, $d2, $d1 - $d2);
}

# 元の関数
sub f
{
    my ($x) = @_;
    $x - ($x ** 3) / (3 * 2) + ($x ** 5) / (5 * 4 * 3 * 2);
}

# Neville (ネヴィル) 補間
sub neville
{
    my ($d, $x, $y) = @_;
    my @w = ();
    for $i (0..N)
    {
        $w[0][$i] = $$y[$i];
    }

    for $j (1..N)
    {
        for $i (0..(N - $j))
        {
            $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]);
        }
    }
    $w[N][0];
}
Z:\>perl Z:\0702.pl
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667    -0.00000
-3.50   -0.73099    -0.73099    -0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964    -0.00000
-2.00   -0.93333    -0.93333    -0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.47943     0.00000
 1.00    0.84167     0.84167    -0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984     0.00000

PHP

<?php
# データ点の数 - 1
define("N", 6);

$x = array();
$y = array();

# 1.5刻みで -4.5~4.5 まで, 7点だけ値をセット
foreach (range(0, N) as $i)
{
    $d = $i * 1.5 - 4.5;
    $x[$i] = $d;
    $y[$i] = f($d);
}

# 0.5刻みで 与えられていない値を補間
foreach (range(0, 18) as $i)
{
    $d  = $i * 0.5 - 4.5;
    $d1 = f($d);
    $d2 = neville($d, $x, $y);

    # 元の関数と比較
    printf("%5.2f\t%8.5f\t%8.5f\t%8.5f\n", $d, $d1, $d2, $d1 - $d2);
}

# 元の関数
function f($x)
{
    return $x - pow($x, 3) / (3 * 2) + pow($x, 5) / (5 * 4 * 3 * 2);
}

# Neville (ネヴィル) 補間
function neville($d, $x, $y)
{
    $w = array();
    foreach (range(0, N) as $i)
    {
        $w[0][$i] = $y[$i];
    }

    foreach (range(1, N) as $j)
    {
        foreach (range(0, N - $j) as $i)
        {
            $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]);
        }
    }
    return $w[N][0];
}
?>
Z:\>php Z:\0702.php
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667    -0.00000
-3.50   -0.73099    -0.73099    -0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964    -0.00000
-2.00   -0.93333    -0.93333    -0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.47943     0.00000
 1.00    0.84167     0.84167    -0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984     0.00000

Python

# coding: Shift_JIS

# データ点の数
N = 7

# 元の関数
def f(x):
    return x - (x ** 3) / (3 * 2) + (x ** 5) / (5 * 4 * 3 * 2)

# Neville (ネヴィル) 補間
def neville(d, x, y):
    w = [[0 for j in range(N)] for i in range(N)]
    for i in range(0, N):
        w[0][i] = y[i]

    for j in range(1, N):
        for i in range(0, N - j):
            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])

    return w[N-1][0]

x = [0 for i in range(N)]
y = [0 for i in range(N)]

# 1.5刻みで -4.5~4.5 まで, 7点だけ値をセット
for i in range(0, N):
    d = i * 1.5 - 4.5
    x[i] = d
    y[i] = f(d)

# 0.5刻みで 与えられていない値を補間
for i in range(0, 19):
    d  = i * 0.5 - 4.5
    d1 = f(d)
    d2 = neville(d, x, y)

    # 元の関数と比較
    print "%5.2f\t%8.5f\t%8.5f\t%8.5f" % (d, d1, d2, d1 - d2)
Z:\>python Z:\0702.py
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667    -0.00000
-3.50   -0.73099    -0.73099    -0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964    -0.00000
-2.00   -0.93333    -0.93333    -0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.47943     0.00000
 1.00    0.84167     0.84167    -0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984     0.00000

Ruby

# データ点の数 - 1
N = 6

# 元の関数
def f(x)
    x - (x ** 3) / (3 * 2) + (x ** 5) / (5 * 4 * 3 * 2)
end

# Neville (ネヴィル) 補間
def neville(d, x, y)
    w = Array.new(N+1) { Array.new(N) }
    (0..N).each do |i|
        w[0][i] = y[i]
    end

    (1..N).each do |j|
        (0..(N - j)).each do |i|
            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
    end

    w[N][0]
end

x = Array.new(N)
y = Array.new(N)

# 1.5刻みで -4.5~4.5 まで, 7点だけ値をセット
(0..N).each do |i|
    d = i * 1.5 - 4.5
    x[i] = d
    y[i] = f(d)
end

# 0.5刻みで 与えられていない値を補間
(0..18).each do |i|
    d  = i * 0.5 - 4.5
    d1 = f(d)
    d2 = neville(d, x, y)

    # 元の関数と比較
    printf("%5.2f\t%8.5f\t%8.5f\t%8.5f\n", d, d1, d2, d1 - d2)
end
Z:\>ruby Z:\0702.rb
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667    -0.00000
-3.50   -0.73099    -0.73099    -0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964    -0.00000
-2.00   -0.93333    -0.93333    -0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.47943     0.00000
 1.00    0.84167     0.84167    -0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984     0.00000

Groovy

Pascal

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 まで, 7点だけ値をセット
    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.
Z:\>fpc -v0 -l- Pas0702.pp

Z:\>Pas0702
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667     0.00000
-3.50   -0.73099    -0.73099     0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964     0.00000
-2.00   -0.93333    -0.93333     0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.47943     0.00000
 1.00    0.84167     0.84167     0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984     0.00000

Ada

with TEXT_IO, Ada.Long_Float_Text_IO;
use  TEXT_IO, Ada.Long_Float_Text_IO;

procedure Ada0702 is
    -- データ点の数 - 1
    N : Constant Integer := 6;

    type Long_Float_Array is array (0..N) of Long_Float;
    x : Long_Float_Array;
    y : Long_Float_Array;
    d, d1, d2 : Long_Float;

    -- 元の関数
    function f(x:Long_Float) return Long_Float is
    begin
        return x - Long_Float(x ** 3) / Long_Float(3 * 2) + Long_Float(x ** 5) / Long_Float(5 * 4 * 3 * 2);
    end f;

    -- Neville (ネヴィル) 補間
    function neville(d:Long_Float; x:Long_Float_Array; y:Long_Float_Array) return Long_Float is
        w : array (0..N, 0..N) of Long_Float;
    begin
        for i in y'Range loop
            w(0, i) := y(i);
        end loop;

        for j in 1 .. x'Last loop
            for i in x'First .. x'Last-j loop
                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 loop;
        end loop;
        return w(N,0);
    end;
begin
    -- 1.5刻みで -4.5~4.5 まで, 7点だけ値をセット
    for i in x'Range loop
        d    := Long_Float(i) * 1.5 - 4.5;
        x(i) := d;
        y(i) := f(d);
    end loop;

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

        -- 元の関数と比較
        Put(d,       Fore=>2, Aft=>2, Exp=>0);
        Put(Ascii.HT);
        Put(d1,      Fore=>3, Aft=>5, Exp=>0);
        Put(Ascii.HT);
        Put(d2,      Fore=>3, Aft=>5, Exp=>0);
        Put(Ascii.HT);
        Put(d1 - d2, Fore=>3, Aft=>5, Exp=>0);
        New_Line;
    end loop;
end Ada0702;
xxxxxx@yyyyyy /Z
$ gnatmake Ada0702.adb

xxxxxx@yyyyyy /Z
$ Ada0702
-4.50    -4.68984    -4.68984     0.00000
-4.00    -1.86667    -1.86667     0.00000
-3.50    -0.73099    -0.73099     0.00000
-3.00    -0.52500    -0.52500     0.00000
-2.50    -0.70964    -0.70964     0.00000
-2.00    -0.93333    -0.93333    -0.00000
-1.50    -1.00078    -1.00078     0.00000
-1.00    -0.84167    -0.84167     0.00000
-0.50    -0.47943    -0.47943     0.00000
 0.00     0.00000     0.00000     0.00000
 0.50     0.47943     0.47943     0.00000
 1.00     0.84167     0.84167     0.00000
 1.50     1.00078     1.00078     0.00000
 2.00     0.93333     0.93333     0.00000
 2.50     0.70964     0.70964    -0.00000
 3.00     0.52500     0.52500     0.00000
 3.50     0.73099     0.73099    -0.00000
 4.00     1.86667     1.86667     0.00000
 4.50     4.68984     4.68984     0.00000

VB.NET

Module VB0702
    'データ点の数 - 1
    Private Const N As Integer = 6

    Public Sub Main()
        Dim x(N) As Double
        Dim y(N) As Double

        '1.5刻みで -4.5~4.5 まで, 7点だけ値をセット
        For i As Integer = 0 To N
            Dim d As Double = i * 1.5 - 4.5
            x(i) = d
            y(i) = f(d)
        Next

        '0.5刻みで 与えられていない値を補間
        For i As Integer = 0 To 18
            Dim d  As Double = i * 0.5 - 4.5
            Dim d1 As Double = f(d)
            Dim d2 As Double = neville(d, x, y)

            '元の関数と比較
            Console.WriteLine(String.Format("{0,5:F2}{4}{1,8:F5}{4}{2,8:F5}{4}{3,8:F5}", d, d1, d2, d1 - d2, vbTab))
        Next
    End Sub

    '元の関数
    Private Function f(ByVal x As Double) As Double
        Return x - (x ^ 3) / (3 * 2) + (x ^ 5) / (5 * 4 * 3 * 2)
    End Function

    'Neville (ネヴィル) 補間
    Private Function neville(ByVal d As Double, ByVal x() As Double, ByVal y() As Double) As Double
        Dim w(N, N) As Double
        For i As Integer = 0 To N
            w(0,i) = y(i)
        Next

        For j As Integer = 1 To N
            For i As Integer = 0 To N - j
                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))
            Next
        Next

        Return w(N,0)
    End Function
End Module
Z:\>vbc -nologo VB0702.vb

Z:\>VB0702
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667     0.00000
-3.50   -0.73099    -0.73099     0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964     0.00000
-2.00   -0.93333    -0.93333     0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.47943     0.00000
 1.00    0.84167     0.84167     0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984     0.00000

C#

using System;

public class CS0702
{
    // データ点の数
    private const int N = 7;

    public static void Main()
    {
        double[] x = new double[N];
        double[] y = new double[N];

        // 1.5刻みで -4.5~4.5 まで, 7点だけ値をセット
        for (int i = 0; i < N; i++)
        {
            double d = i * 1.5 - 4.5;
            x[i] = d;
            y[i] = f(d);
        }

        // 0.5刻みで 与えられていない値を補間
        for (int i = 0; i <= 18; i++)
        {
            double d  = i * 0.5 - 4.5;
            double d1 = f(d);
            double d2 = neville(d, x, y);

            // 元の関数と比較
            Console.WriteLine(string.Format("{0,5:F2}\t{1,8:F5}\t{2,8:F5}\t{3,8:F5}", d, d1, d2, d1 - d2));
        }
    }

    // 元の関数
    private static double f(double x)
    {
        return x - Math.Pow(x,3) / (3 * 2) + Math.Pow(x,5) / (5 * 4 * 3 * 2);
    }

    // Neville (ネヴィル) 補間
    private static double neville(double d, double[] x, double[] y)
    {
        double[,] w = new double[N,N];
        for (int i = 0; i < N; i++)
            w[0,i] = y[i];

        for (int j = 1; j < N; j++)
        {
            for (int i = 0; i < N - j; i++)
                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]);
        }

        return w[N-1,0];
    }
}
Z:\>csc -nologo CS0702.cs

Z:\>CS0702
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667     0.00000
-3.50   -0.73099    -0.73099     0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964     0.00000
-2.00   -0.93333    -0.93333     0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.47943     0.00000
 1.00    0.84167     0.84167     0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984     0.00000

Java

public class Java0702 {

    // データ点の数
    private static final int N = 7;

    public static void main(String []args) {
        double[] x = new double[N];
        double[] y = new double[N];

        // 1.5刻みで -4.5~4.5 まで, 7点だけ値をセット
        for (int i = 0; i < N; i++) {
            double d = i * 1.5 - 4.5;
            x[i] = d;
            y[i] = f(d);
        }

        // 0.5刻みで 与えられていない値を補間
        for (int i = 0; i <= 18; i++) {
            double d  = i * 0.5 - 4.5;
            double d1 = f(d);
            double d2 = neville(d, x, y);

            // 元の関数と比較
            System.out.println(String.format("%5.2f\t%8.5f\t%8.5f\t%8.5f", d, d1, d2, d1 - d2));
        }
    }

    // 元の関数
    private static double f(double x) {
        return x - Math.pow(x,3) / (3 * 2) + Math.pow(x,5) / (5 * 4 * 3 * 2);
    }

    // Neville (ネヴィル) 補間
    private static double neville(double d, double[] x, double[] y) {
        double[][] w = new double[N][N];
        for (int i = 0; i < N; i++)
            w[0][i] = y[i];

        for (int j = 1; j < N; j++) {
            for (int i = 0; i < N - j; i++)
                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]);
        }

        return w[N-1][0];
    }
}
Z:\>javac Java0702.java

Z:\>java Java0702
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667    -0.00000
-3.50   -0.73099    -0.73099    -0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964    -0.00000
-2.00   -0.93333    -0.93333    -0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.47943     0.00000
 1.00    0.84167     0.84167    -0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984     0.00000

C++

#include <iostream>
#include <iomanip>
#include <math.h>
#include <stdio.h>
using namespace std;

// データ点の数
const int N = 7;

// 元の関数
double f(double x);

// Neville (ネヴィル) 補間
double neville(double d, double x[], double y[]);

int main()
{
    double x[N], y[N];

    // 1.5刻みで -4.5~4.5 まで, 7点だけ値をセット
    for (int i = 0; i < N; i++)
    {
        double d = i * 1.5 - 4.5;
        x[i] = d;
        y[i] = f(d);
    }

    // 0.5刻みで 与えられていない値を補間
    for (int i = 0; i <= 18; i++)
    {
        double d  = i * 0.5 - 4.5;
        double d1 = f(d);
        double d2 = neville(d, x, y);

        // 元の関数と比較
        cout << setw(5) << fixed << setprecision(2) << d       << '\t';
        cout << setw(8) << fixed << setprecision(5) << d1      << '\t';
        cout << setw(8) << fixed << setprecision(5) << d2      << '\t';
        cout << setw(8) << fixed << setprecision(5) << d1 - d2 << endl;
    }

   return 0;
}

// 元の関数
double f(double x)
{
    return x - pow(x,3) / (3 * 2) + pow(x,5) / (5 * 4 * 3 * 2);
}

// Neville (ネヴィル) 補間
double neville(double d, double x[], double y[])
{
    double w[N][N];
    for (int i = 0; i < N; i++)
        w[0][i] = y[i];

    for (int j = 1; j < N; j++)
    {
        for (int i = 0; i < N - j; i++)
              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]);
    }

    return w[N-1][0];
}
Z:\>bcc32 -q CP0702.cpp
cp0702.cpp:

Z:\>CP0702
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667     0.00000
-3.50   -0.73099    -0.73099     0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964     0.00000
-2.00   -0.93333    -0.93333    -0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.47943     0.00000
 1.00    0.84167     0.84167     0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964    -0.00000
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     0.73099    -0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984     0.00000

Objective-C

#import <Foundation/Foundation.h>
#import <math.h>

// データ点の数
const int N = 7;

// 元の関数
double f(double x);

// Neville (ネヴィル) 補間
double neville(double d, double x[], double y[]);

int main()
{
    double x[N], y[N];

    // 1.5刻みで -4.5~4.5 まで, 7点だけ値をセット
    int i;
    for (i = 0; i < N; i++)
    {
        double d = i * 1.5 - 4.5;
        x[i] = d;
        y[i] = f(d);
    }

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

        // 元の関数と比較
        printf("%5.2f\t%8.5f\t%8.5f\t%8.5f\n", d, d1, d2, d1 - d2);
    }

   return 0;
}

// 元の関数
double f(double x)
{
    return x - pow(x,3) / (3 * 2) + pow(x,5) / (5 * 4 * 3 * 2);
}

// Neville (ネヴィル) 補間
double neville(double d, double x[], double y[])
{
    int i, j;
    double w[N][N];
    for (i = 0; i < N; i++)
        w[0][i] = y[i];

    for (j = 1; j < N; j++)
    {
        for (i = 0; i < N - j; i++)
              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]);
    }

    return w[N-1][0];
}
xxxxxx@yyyyyy /Z
$ gcc -o OC0702 OC0702.m -lobjc -lgnustep-base -I $INCLUDE -L $LIB $CFLAGS

xxxxxx@yyyyyy /Z
$ OC0702
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667     0.00000
-3.50   -0.73099    -0.73099     0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964     0.00000
-2.00   -0.93333    -0.93333    -0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.47943     0.00000
 1.00    0.84167     0.84167     0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964    -0.00000
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     0.73099    -0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984     0.00000

D

import std.stdio;
import std.math;

// データ点の数
const int N = 7;

void main(string[] args)
{
    double x[N];
    double y[N];

    // 1.5刻みで -4.5~4.5 まで, 7点だけ値をセット
    for (int i = 0; i < N; i++)
    {
        double d = i * 1.5 - 4.5;
        x[i] = d;
        y[i] = f(d);
    }

    // 0.5刻みで 与えられていない値を補間
    for (int i = 0; i <= 18; i++)
    {
        double d  = i * 0.5 - 4.5;
        double d1 = f(d);
        double d2 = neville(d, x, y);

        // 元の関数と比較
        writefln("%5.2f\t%8.5f\t%8.5f\t%8.5f", d, d1, d2, d1 - d2);
    }
}

// 元の関数
double f(double x)
{
    return x - pow(x,3) / (3 * 2) + pow(x,5) / (5 * 4 * 3 * 2);
}

// Neville (ネヴィル) 補間
double neville(double d, double x[], double y[])
{
    double w[N][N];
    for (int i = 0; i < N; i++)
        w[0][i] = y[i];

    for (int j = 1; j < N; j++)
    {
        for (int i = 0; i < N - j; i++)
              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]);
    }

    return w[N-1][0];
}
Z:\>dmd D0702.d

Z:\>D0702
-4.50   -4.68984    -4.68984    -0.00000
-4.00   -1.86667    -1.86667     0.00000
-3.50   -0.73099    -0.73099     0.00000
-3.00   -0.52500    -0.52500    -0.00000
-2.50   -0.70964    -0.70964     0.00000
-2.00   -0.93333    -0.93333     0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000     0.00000    -0.00000
 0.50    0.47943     0.47943     0.00000
 1.00    0.84167     0.84167     0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964    -0.00000
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     0.73099    -0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984     0.00000

Go

package main

import "fmt"
import "math"

// データ点の数
const N = 7

func main() {
    var x [N]float64
    var y [N]float64

    // 1.5刻みで -4.5~4.5 まで, 7点だけ値をセット
    for i := 0; i < N; i++ {
        var d float64 = float64(i) * 1.5 - 4.5
        x[i] = d
        y[i] = f(d)
    }

    // 0.5刻みで 与えられていない値を補間
    for i := 0; i <= 18; i++ {
        var d  float64 = float64(i) * 0.5 - 4.5
        var d1 float64 = f(d)
        var d2 float64 = neville(d, x[:], y[:])

        // 元の関数と比較
        fmt.Printf("%5.2f\t%8.5f\t%8.5f\t%8.5f\n", d, d1, d2, d1 - d2)
    }
}

// 元の関数
func f(x float64) float64 {
    return x - math.Pow(x,3) / (3 * 2) + math.Pow(x,5) / (5 * 4 * 3 * 2)
}

// Neville (ネヴィル) 補間
func neville(d float64, x []float64, y []float64) float64 {
    var w[N][N] float64
    for i := 0; i < N; i++ {
        w[0][i] = y[i]
    }

    for j := 1; j < N; j++ {
        for i := 0; i < N - j; i++ {
              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])
        }
    }

    return w[N-1][0]
}
Z:\>8g GO0702.go

Z:\>8l -o GO0702.exe GO0702.8

Z:\>GO0702
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667    -0.00000
-3.50   -0.73099    -0.73099    -0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964    -0.00000
-2.00   -0.93333    -0.93333    -0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.47943     0.00000
 1.00    0.84167     0.84167    -0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984     0.00000

Scala

object Scala0702 {
    // データ点の数 - 1
    val N = 6

    def main(args: Array[String]) {
        // 1.5刻みで -4.5~4.5 まで, 7点だけ値をセット
        val x = (0 to N).map(_ * 1.5 - 4.5)
        val y = x.map(f)

        // 0.5刻みで 与えられていない値を補間
        val d1 = (0 to 18).map(_ * 0.5 - 4.5)
        val d2 = d1.map(f)
        val d3 = d1.map(neville(_, x, y))

        (d1 zip d2 zip d3).foreach {
            case ((d1, d2), d3) =>
                println("%5.2f\t%8.5f\t%8.5f\t%8.5f".format(d1, d2, d3, d2 - d3))
        }
    }

    // 元の関数
    def f(x:Double) = {
        x - Math.pow(x,3) / (3 * 2) + Math.pow(x,5) / (5 * 4 * 3 * 2)
    }

    // Neville (ネヴィル) 補間
    def neville(d:Double, x:IndexedSeq[Double], y:IndexedSeq[Double]) = {
        val w = Array.ofDim[Double](N + 1, N + 1)
        for (i <- 0 to N)
            w(0)(i) = y(i)

        for (j <- 1 to N) {
            for (i <- 0 to N - j)
                  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))
        }

        w(N)(0)
    }
}
Z:\>scala Scala0702.scala
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667    -0.00000
-3.50   -0.73099    -0.73099    -0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964    -0.00000
-2.00   -0.93333    -0.93333    -0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.47943     0.00000
 1.00    0.84167     0.84167    -0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984     0.00000

F#

module Fs0702

open System

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

// 元の関数
let f (x:double):double =
    x - Math.Pow(x,3.0) / (float (3 * 2)) + Math.Pow(x,5.0) / (float (5 * 4 * 3 * 2))

// Neville (ネヴィル) 補間
let neville(d:double) (x:double list) (y:double list) =
    let w = Array2D.zeroCreate<double> (N+1) (N+1)
    for i in [0..N] do
        w.[0,i] <- y.[i]
    for j in [1..N] do
        for i in [0..N-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])
    w.[N,0]

// 1.5刻みで -4.5~4.5 まで, 7点だけ値をセット
let x = [0..N] |> List.map(fun i -> (float i) * 1.5 - 4.5)
let y = x |> List.map(f)

// 0.5刻みで 与えられていない値を補間
let d1 = [0..18] |> List.map(fun i -> (float i) * 0.5 - 4.5)
let d2 = d1 |> List.map(f)
let d3 = d1 |> List.map(fun d -> (neville d x y))

(List.zip (List.zip d1 d2) d3)
|> List.iter (fun ((d1, d2), d3) -> printfn "%5.2f\t%8.5f\t%8.5f\t%8.5f" d1 d2 d3 (d2 - d3))

exit 0
Z:\>fsi  --nologo --quiet Fs0702.fs
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667     0.00000
-3.50   -0.73099    -0.73099     0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964     0.00000
-2.00   -0.93333    -0.93333     0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.47943     0.00000
 1.00    0.84167     0.84167     0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984     0.00000

Clojure

(def N 7)

; 元の関数
(defn f[x]
    (+ (- x (/ (Math/pow x 3.0) (* 3 2))) (/ (Math/pow x 5.0) (* 5 (* 4 (* 3 2))))))

; Neville (ネヴィル) 補間
(defn neville [d x y]
    (def w (cons y nil))
    (doseq [j (range 1 N)]
        (def t_list (list))
        (doseq [i (range 0 (- N j))]
            (def t (nth w 0))
            (def t_list
                (cons
                    (+ (nth t (+ i 1))
                       (/ (* (- (nth t (+ i 1)) (nth t i))
                             (- d (nth x (+ i j))))
                          (- (nth x (+ i j)) (nth x i))))
                    t_list
                )
            )
        )
        (def w (cons (reverse t_list) w))
    )
    (nth (nth w 0) 0)
)

; 1.5刻みで -4.5~4.5 まで, 7点だけ値をセット
(def x (map #(- (* % 1.5) 4.5) (range 0 N)))
(def y (map #(f %) x))

; 0.5刻みで 与えられていない値を補間
(def d1 (map #(- (* % 0.5) 4.5) (range 0 19)))
(def d2 (map #(f %) d1))
(def d3 (map #(neville % x y) d1))

(doseq [d (map list d1 d2 d3)]
    (def d1 (nth d 0))
    (def d2 (nth d 1))
    (def d3 (nth d 2))
    (println (format "%5.2f\t%8.5f\t%8.5f\t%8.5f" d1 d2 d3 (- d2 d3))))
Z:\>java -cp C:\ProgramFiles\clojure-1.5.1\clojure-1.5.1.jar clojure.main Clj0702.clj
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667    -0.00000
-3.50   -0.73099    -0.73099    -0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964    -0.00000
-2.00   -0.93333    -0.93333    -0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.47943     0.00000
 1.00    0.84167     0.84167    -0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984     0.00000

Haskell

import Text.Printf
import Control.Monad

-- データ点の数 - 1
n = 6 :: Int

-- 元の関数
f::Double->Double
f x = x - (x ^ 3) / (fromIntegral (3 * 2)) + (x ^ 5) / (fromIntegral (5 * 4 * 3 * 2))

-- Neville (ネヴィル) 補間
neville::Double->[Double]->[Double]->Int->Double
neville d x w j =
    let
        t = map(\i -> do
                (w!!(i+1) + (w!!(i+1) - w!!i) * (d - x!!(i+j)) / (x!!(i+j) - x!!i))
            ) [0..(n-j)::Int]
    in
        if j == n then t!!0
                  else (neville d x t (j + 1))

main = do
    -- 1.5刻みで -4.5~4.5 まで, 7点だけ値をセット
    let x = map(\i -> (fromIntegral i) * 1.5 - 4.5) [0..n]
    let y = map(\i -> f(i)) x

    -- 0.5刻みで 与えられていない値を補間
    let d1 = map(\i -> (fromIntegral i) * 0.5 - 4.5) [0..18]
    let d2 = map(\i -> (f i)) d1
    let d3 = map(\i -> (neville i x y 1)) d1

    forM_ (zip (zip d1 d2) d3) $ \((d1, d2), d3) -> do
        printf "%5.2f\t%8.5f\t%8.5f\t%8.5f\n" d1 d2 d3 (d2 - d3)
Z:\>runghc Hs0702.hs
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667    -0.00000
-3.50   -0.73099    -0.73099    -0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964    -0.00000
-2.00   -0.93333    -0.93333    -0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.47943     0.00000
 1.00    0.84167     0.84167    -0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984     0.00000
inserted by FC2 system