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

ニュートン補間

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

このとき,

と定義する.
これを $ x_0 $ と $ x_1 $ の第１差分商という.
同様に

と定義し, これを $ x_0, x_1, x_2 $ の第２差分商という.
第$n$差分商を

で表すと, 与えられた $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 まで, ７点だけ値をセット
Dim i
For i = 0 To N
    Dim d1: d1 = i * 1.5 - 4.5
    x(i) = d1
    y(i) = f(d1)
Next

'差分商の表を作る
Dim d(): ReDim d(N, N)
Dim j
For j = 0 To N
    d(0,j) = y(j)
Next

For i = 1 To N
    For j = 0 To (N - i)
        d(i,j) = (d(i-1,j+1) - d(i-1,j)) / (x(j+i) - x(j))
    Next
Next

'ｎ階差分商
Dim a(): ReDim a(N)
For j = 0 To N
    a(j) = d(j,0)
Next

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

    '元の関数と比較
    WScript.StdOut.Write     Right(Space(5) & FormatNumber(d1,      2, -1, 0, 0), 5) & vbTab
    WScript.StdOut.Write     Right(Space(8) & FormatNumber(d2,      5, -1, 0, 0), 8) & vbTab
    WScript.StdOut.Write     Right(Space(8) & FormatNumber(d3,      5, -1, 0, 0), 8) & vbTab
    WScript.StdOut.WriteLine Right(Space(8) & FormatNumber(d2 - d3, 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

'Newton (ニュートン) 補間
Private Function newton(ByVal d, ByVal x(), ByVal a())
    Dim sum: sum = a(0)
    Dim i, j
    For i = 1 To N
        Dim prod: prod = a(i)
        For j = 0 To (i - 1)
            If j <> i Then
                prod = prod * (d - x(j))
            End If
        Next
        sum = sum + prod
    Next
    newton = sum
End Function

Z:\>cscript //nologo Z:\0703.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 まで, ７点だけ値をセット
for (var i = 0; i < N; i++)
{
    var d = i * 1.5 - 4.5
    x[i] = d
    y[i] = f(d)
}

// 差分商の表を作る
var d = []
d[0]  = []
for (var j = 0; j < N; j++)
    d[0][j] = y[j]

for (var i = 1; i < N; i++)
{
    d[i]  = []
    for (var j = 0; j < N - i; j++)
        d[i][j] = (d[i-1][j+1] - d[i-1][j]) / (x[j+i] - x[j])
}
// ｎ階差分商
var a = []
for (var j = 0; j < N; j++)
    a[j] = d[j][0]

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

    // 元の関数と比較
    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)
}

// Newton (ニュートン) 補間
function newton(d, x, a)
{
    var sum = a[0]
    for (var i = 1; i < N; i++)
    {
        var prod = a[i]
        for (var j = 0; j < i; j++)
            prod *= (d - x[j])
        sum += prod
    }

    return sum
}

Z:\>cscript //nologo Z:\0703.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)
}

# Newton (ニュートン) 補間
function newton($d, $x, $a)
{
    $sum = $a[0]
    foreach ($i in 1..($N - 1))
    {
        $prod = $a[$i]
        foreach ($j in 0..($i - 1))
        {
            $prod *= ($d - $x[$j])
        }
        $sum += $prod
    }
    $sum
}

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

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

# 差分商の表を作る
$d = New-Object "double[,]" $N,$N
foreach ($j in 0..($N - 1))
{
    $d[0,$j] = $y[$j]
}
foreach ($i in 1..($N - 1))
{
    foreach ($j in 0..($N - $i - 1))
    {
        $d[$i,$j] = ($d[($i-1),($j+1)] - $d[($i-1),$j]) / ($x[($j+$i)] - $x[$j])
    }
}
# ｎ階差分商
$a = New-Object double[] $N
foreach ($j in 0..($N - 1))
{
    $a[$j] = $d[$j,0]
}

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

    # 元の関数と比較
    Write-Host ([String]::Format("{0,5:F2}`t{1,8:F5}`t{2,8:F5}`t{3,8:F5}", $d1, $d2, $d3, ($d2 - $d3)))
}

Z:\>powershell -file Z:\0703.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 まで, ７点だけ値をセット
for $i (0..N)
{
    my $d1 = $i * 1.5 - 4.5;
    $x[$i] = $d1;
    $y[$i] = f($d1);
}

# 差分商の表を作る
my @d = ();
for $j (0..N)
{
    $d[0][$j] = $y[$j];
}
for $i (1..N)
{
    for $j (0..(N - $i))
    {
        $d[$i][$j] = ($d[$i-1][$j+1] - $d[$i-1][$j]) / ($x[$j+$i] - $x[$j]);
    }
}
# ｎ階差分商
my @a = ();
for $j (0..N)
{
    $a[$j] = $d[$j][0];
}

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

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

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

# Newton (ニュートン) 補間
sub newton
{
    my ($d, $x, $a) = @_;
    my $sum = $$a[0];
    for $i (1..N)
    {
        my $prod = $$a[$i];
        for $j (0..($i - 1))
        {
            $prod *= ($d - $$x[$j]);
        }
        $sum += $prod;
    }
    $sum;
}

Z:\>perl Z:\0703.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 まで, ７点だけ値をセット
foreach (range(0, N) as $i)
{
    $d1 = $i * 1.5 - 4.5;
    $x[$i] = $d1;
    $y[$i] = f($d1);
}

# 差分商の表を作る
$d = array();
foreach (range(0, N) as $j)
{
    $d[0][$j] = $y[$j];
}
foreach (range(1, N) as $i)
{
    foreach (range(0, N - $i) as $j)
    {
        $d[$i][$j] = ($d[$i-1][$j+1] - $d[$i-1][$j]) / ($x[$j+$i] - $x[$j]);
    }
}
# ｎ階差分商
$a = array();
foreach (range(0, N) as $j)
{
    $a[$j] = $d[$j][0];
}

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

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

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

# Newton (ニュートン) 補間
function newton($d, $x, $a)
{
    $sum = $a[0];
    foreach (range(1, N) as $i)
    {
        $prod = $a[$i];
        foreach (range(0, $i - 1) as $j)
        {
            $prod *= ($d - $x[$j]);
        }
        $sum += $prod;
    }
    return $sum;
}
?>

Z:\>php Z:\0703.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)

# Newton (ニュートン) 補間
def newton(d, x, a):
    sum = a[0]
    for i in range(1, N):
        prod = a[i]
        for j in range(0, i):
            prod *= (d - x[j])
        sum += prod
    return sum

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

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

# 差分商の表を作る
d = [[0 for j in range(N)] for i in range(N)]
for j in range(0, N):
    d[0][j] = y[j]

for i in range(1, N):
    for j in range(0, N - i):
        d[i][j] = (d[i-1][j+1] - d[i-1][j]) / (x[j+i] - x[j])

# ｎ階差分商
a = [0 for i in range(N)]
for j in range(0, N):
    a[j] = d[j][0]

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

    # 元の関数と比較
    print "%5.2f\t%8.5f\t%8.5f\t%8.5f" % (d1, d2, d3, d2 - d3)

Z:\>python Z:\0703.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

# Newton (ニュートン) 補間
def newton(d, x, a)
    sum = a[0]
    (1..N).each do |i|
        prod = a[i]
        (0..(i - 1)).each do |j|
            if j != i
                prod *= (d - x[j])
            end
        end
        sum += prod
    end
    sum
end

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

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

# 差分商の表を作る
d = Array.new(N+1) { Array.new(N) }
(0..N).each do |j|
    d[0][j] = y[j]
end

(1..N).each do |i|
    (0..(N-i)).each do |j|
        d[i][j] = (d[i-1][j+1] - d[i-1][j]) / (x[j+i] - x[j])
    end
end

# ｎ階差分商
a = Array.new(N)
(0..N).each do |j|
    a[j] = d[j][0]
end

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

    # 元の関数と比較
    printf("%5.2f\t%8.5f\t%8.5f\t%8.5f\n", d1, d2, d3, d2 - d3)
end

Z:\>ruby Z:\0703.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 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.

Z:\>fpc -v0 -l- Pas0703.pp

Z:\>Pas0703
-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 Ada0703 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;
    a : Long_Float_Array;
    d : array (0..N, 0..N) of Long_Float;
    d1, d2, d3 : 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;

    -- Newton (ニュートン) 補間
    function newton(d:Long_Float; x:Long_Float_Array; a:Long_Float_Array) return Long_Float is
        sum, prod :Long_Float;
    begin
        sum := a(0);
        for i in 1 .. x'Last loop
            prod := a(i);
            for j in x'First .. i-1 loop
                prod := prod * (d - x(j));
            end loop;
            sum := sum + prod;
        end loop;
        return sum;
    end;
begin
    -- 1.5刻みで -4.5～4.5 まで, ７点だけ値をセット
    for i in x'Range loop
        d1   := Long_Float(i) * 1.5 - 4.5;
        x(i) := d1;
        y(i) := f(d1);
    end loop;

    -- 差分商の表を作る
    for j in y'Range loop
        d(0,j) := y(j);
    end loop;

    for i in 1 .. x'Last loop
        for j in x'First .. x'Last-i loop
            d(i,j) := (d(i-1,j+1) - d(i-1,j)) / (x(j+i) - x(j));
        end loop;
    end loop;

    -- ｎ階差分商
    for j in a'Range loop
        a(j) := d(j,0);
    end loop;

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

        -- 元の関数と比較
        Put(d1,      Fore=>2, Aft=>2, Exp=>0);
        Put(Ascii.HT);
        Put(d2,      Fore=>3, Aft=>5, Exp=>0);
        Put(Ascii.HT);
        Put(d3,      Fore=>3, Aft=>5, Exp=>0);
        Put(Ascii.HT);
        Put(d2 - d3, Fore=>3, Aft=>5, Exp=>0);
        New_Line;
    end loop;
end Ada0703;

xxxxxx@yyyyyy /Z
$ gnatmake Ada0703.adb

xxxxxx@yyyyyy /Z
$ Ada0703
-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 VB0703
    'データ点の数 - 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 まで, ７点だけ値をセット
        For i As Integer = 0 To N
            Dim d1 As Double = i * 1.5 - 4.5
            x(i) = d1
            y(i) = f(d1)
        Next

        '差分商の表を作る
        Dim d(N, N) As Double
        For j As Integer = 0 To N
            d(0,j) = y(j)
        Next

        For i As Integer = 1 To N
            For j As Integer = 0 To (N - i)
                d(i,j) = (d(i-1,j+1) - d(i-1,j)) / (x(j+i) - x(j))
            Next
        Next

        'ｎ階差分商
        Dim a(N) As Double
        For j As Integer = 0 To N
            a(j) = d(j,0)
        Next

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

            '元の関数と比較
            Console.WriteLine(String.Format("{0,5:F2}{4}{1,8:F5}{4}{2,8:F5}{4}{3,8:F5}", d1, d2, d3, d2 - d3, 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

    'Newton (ニュートン) 補間
    Private Function newton(ByVal d As Double, ByVal x() As Double, ByVal a() As Double) As Double
        Dim sum As Double = a(0)
        For i As Integer = 1 To N
            Dim prod As Double = a(i)
            For j As Integer = 0 To (i - 1)
                If j <> i Then
                    prod *= (d - x(j))
                End If
            Next
            sum += prod
        Next
        Return sum
    End Function
End Module

Z:\>vbc -nologo VB0703.vb

Z:\>VB0703
-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 CS0703
{
    // データ点の数
    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 まで, ７点だけ値をセット
        for (int i = 0; i < N; i++)
        {
            double d1 = i * 1.5 - 4.5;
            x[i] = d1;
            y[i] = f(d1);
        }

        // 差分商の表を作る
        double[,] d = new double[N,N];
        for (int j = 0; j < N; j++)
            d[0,j] = y[j];

        for (int i = 1; i < N; i++)
        {
            for (int j = 0; j < N - i; j++)
                d[i,j] = (d[i-1,j+1] - d[i-1,j]) / (x[j+i] - x[j]);
        }

        // ｎ階差分商
        double[] a = new double[N];
        for (int j = 0; j < N; j++)
            a[j] = d[j,0];

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

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

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

    // Newton (ニュートン) 補間
    private static double newton(double d, double[] x, double[] a)
    {
        double sum = a[0];
        for (int i = 1; i < N; i++)
        {
            double prod = a[i];
            for (int j = 0; j < i; j++)
                prod *= (d - x[j]);
            sum += prod;
        }

        return sum;
    }
}

Z:\>csc -nologo CS0703.cs

Z:\>CS0703
-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 Java0703 {

    // データ点の数
    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 まで, ７点だけ値をセット
        for (int i = 0; i < N; i++) {
            double d1 = i * 1.5 - 4.5;
            x[i] = d1;
            y[i] = f(d1);
        }

        // 差分商の表を作る
        double[][] d = new double[N][N];
        for (int j = 0; j < N; j++)
            d[0][j] = y[j];

        for (int i = 1; i < N; i++) {
            for (int j = 0; j < N - i; j++)
                d[i][j] = (d[i-1][j+1] - d[i-1][j]) / (x[j+i] - x[j]);
        }

        // ｎ階差分商
        double[] a = new double[N];
        for (int j = 0; j < N; j++)
            a[j] = d[j][0];

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

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

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

    // Newton (ニュートン) 補間
    private static double newton(double d, double[] x, double[] a) {
        double sum = a[0];
        for (int i = 1; i < N; i++) {
            double prod = a[i];
            for (int j = 0; j < i; j++)
                prod *= (d - x[j]);
            sum += prod;
        }

        return sum;
    }
}

Z:\>javac Java0703.java

Z:\>java Java0703
-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>
using namespace std;

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

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

// Newton (ニュートン) 補間
double newton(double d, double x[], double a[]);

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

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

    // 差分商の表を作る
    double d[N][N];
    for (int j = 0; j < N; j++)
        d[0][j] = y[j];

    for (int i = 1; i < N; i++)
    {
        for (int j = 0; j < N - i; j++)
            d[i][j] = (d[i-1][j+1] - d[i-1][j]) / (x[j+i] - x[j]);
    }

    // ｎ階差分商
    double a[N];
    for (int j = 0; j < N; j++)
        a[j] = d[j][0];

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

        // 元の関数と比較
        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);
}

// Newton (ニュートン) 補間
double newton(double d, double x[], double a[])
{
    double sum = a[0];
    for (int i = 1; i < N; i++)
    {
        double prod = a[i];
        for (int j = 0; j < i; j++)
            prod *= (d - x[j]);
        sum += prod;
    }

    return sum;
}

Z:\>bcc32 -q CP0703.cpp
cp0703.cpp:

Z:\>CP0703
-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);

// Newton (ニュートン) 補間
double newton(double d, double x[], double a[]);

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

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

    // 差分商の表を作る
    double d[N][N];
    for (j = 0; j < N; j++)
        d[0][j] = y[j];

    for (i = 1; i < N; i++)
    {
        for (j = 0; j < N - i; j++)
            d[i][j] = (d[i-1][j+1] - d[i-1][j]) / (x[j+i] - x[j]);
    }

    // ｎ階差分商
    double a[N];
    for (j = 0; j < N; j++)
        a[j] = d[j][0];

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

        // 元の関数と比較
        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);
}

// Newton (ニュートン) 補間
double newton(double d, double x[], double a[])
{
    int i, j;
    double sum = a[0];
    for (i = 1; i < N; i++)
    {
        double prod = a[i];
        for (j = 0; j < i; j++)
            prod *= (d - x[j]);
        sum += prod;
    }

    return sum;
}

xxxxxx@yyyyyy /Z
$ gcc -o OC0703 OC0703.m -lobjc -lgnustep-base -I $INCLUDE -L $LIB $CFLAGS

xxxxxx@yyyyyy /Z
$ OC0703
-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 まで, ７点だけ値をセット
    for (int i = 0; i < N; i++)
    {
        double d = i * 1.5 - 4.5;
        x[i] = d;
        y[i] = f(d);
    }

    // 差分商の表を作る
    double d[N][N];
    for (int j = 0; j < N; j++)
        d[0][j] = y[j];

    for (int i = 1; i < N; i++)
    {
        for (int j = 0; j < N - i; j++)
            d[i][j] = (d[i-1][j+1] - d[i-1][j]) / (x[j+i] - x[j]);
    }

    // ｎ階差分商
    double a[N];
    for (int j = 0; j < N; j++)
        a[j] = d[j][0];

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

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

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

// Newton (ニュートン) 補間
double newton(double d, double x[], double a[])
{
    double sum = a[0];
    for (int i = 1; i < N; i++)
    {
        double prod = a[i];
        for (int j = 0; j < i; j++)
            prod *= (d - x[j]);
        sum += prod;
    }

    return sum;
}

Z:\>dmd D0703.d

Z:\>D0703
-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 まで, ７点だけ値をセット
    for i := 0; i < N; i++ {
        var d float64 = float64(i) * 1.5 - 4.5
        x[i] = d
        y[i] = f(d)
    }

    // 差分商の表を作る
    var d[N][N] float64
    for j := 0; j < N; j++ {
        d[0][j] = y[j]
    }
    for i := 1; i < N; i++ {
        for j := 0; j < N - i; j++ {
            d[i][j] = (d[i-1][j+1] - d[i-1][j]) / (x[j+i] - x[j])
        }
    }

    // ｎ階差分商
    var a [N]float64
    for j := 0; j < N; j++ {
        a[j] = d[j][0]
    }

    // 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 = newton(d, x[:], a[:])

        // 元の関数と比較
        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)
}

// Newton (ニュートン) 補間
func newton(d float64, x []float64, a []float64) float64 {
    var sum float64 = a[0]
    for i := 1; i < N; i++ {
        var prod float64 = a[i]
        for j := 0; j < i; j++ {
            prod *= (d - x[j])
        }
        sum += prod
    }
    return sum
}

Z:\>8g GO0703.go

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

Z:\>GO0703
-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 Scala0703 {
    // データ点の数 - 1
    val N = 6

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

        // 差分商の表を作る
        val d = Array.ofDim[Double](N + 1, N + 1)
        for (j <- 0 to N)
            d(0)(j) = y(j)

        for (i <- 1 to N) {
            for (j <- 0 to N - i)
                d(i)(j) = (d(i-1)(j+1) - d(i-1)(j)) / (x(j+i) - x(j))
        }

        // ｎ階差分商
        val a = (0 to N).map(d(_)(0))

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

        (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)
    }

    // Newton (ニュートン) 補間
    def newton(d:Double, x:IndexedSeq[Double], a:IndexedSeq[Double]) = {
        var sum_list = List(a(0))
        for (i <- 1 to N) {
            var prod_list = List(a(i))
            for (j <- 0 to i - 1) {
                prod_list = (d - x(j))::prod_list
            }
            sum_list = (prod_list.product)::sum_list
        }
        sum_list.sum
    }
}

Z:\>scala Scala0703.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 Fs0703

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))

// Newton (ニュートン) 補間
let newton(d:double) (x:double list) (a:double list) =
    let mutable sum_list = [a.[0]]
    for i in [1..N] do
        let mutable prod_list = [a.[i]]
        for j in [0..i-1] do
            prod_list <- (d - x.[j])::prod_list
        sum_list <- (prod_list |> List.reduce(*))::sum_list
    sum_list |> List.sum

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

// 差分商の表を作る
let d = Array2D.zeroCreate<double> (N+1) (N+1)
for j in [0..N] do
    d.[0,j] <- y.[j]

for i in [1..N] do
    for j in [0..N-i] do
        d.[i,j] <- (d.[i-1,j+1] - d.[i-1,j]) / (x.[j+i] - x.[j])

// ｎ階差分商
let a = [0..N] |> List.map(fun i -> d.[i,0])

// 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 -> (newton d x a))

(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 Fs0703.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))))))

; Newton (ニュートン) 補間
(defn newton [d x a]
    (def sum_list (list (nth a 0)))
    (doseq [i (range 1 N)]
        (def prod_list (list (nth a i)))
        (doseq [j (range 0 i)]
            (def prod_list (cons (- d (nth x j)) prod_list)))
        (def w (reduce * prod_list))
        (def sum_list (cons w sum_list)))
    (reduce + sum_list))

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

; 差分商の表を作る
(def d (cons y nil))
(doseq [i (range 1 N)]
    (def w (nth d 0))
    (def t (list))
    (doseq [j (range 0 (- N i))]
        (def t
            (cons
                (/ (- (nth w (+ j 1)) (nth w j))
                   (- (nth x (+ j i)) (nth x j)))
                t
            )
        )
    )
    (def d (cons (reverse t) d))
)
(def d (reverse d))

; ｎ階差分商
(def a (map #(nth (nth d %) 0) (range 0 N)))

; 0.5刻みで 与えられていない値を補間
(def d1 (map #(- (* % 0.5) 4.5) (range 0 19)))
(def d2 (map #(f %) d1))
(def d3 (map #(newton % x a) 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 Clj0703.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))

-- Newton (ニュートン) 補間
newton::Double->[Double]->[Double]->Double
newton d x a =
    let sum_list = map(\i -> do
        let prod_list = map(\j -> do
            d - x!!j
            ) $ [0..(i-1)::Int]
        product $ a!!i : prod_list
        ) [1..n::Int]
    in
        sum $ a!!0 : sum_list

-- 差分商の表を作る
make_table::[Double]->[Double]->[Double]->Int->[Double]
make_table x d a i =
    let
        w = map(\j -> do
            ((d!!(j+1) - d!!j) / (x!!(j+i) - x!!j))
            ) $ [0..(n-i)::Int]
        t = w!!0:a
    in
        -- ｎ階差分商
        if i == n then t
                  else (make_table x w t (i + 1))

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

    -- 差分商の表を作る
    let a = reverse (make_table x y [y!!0] 1)

    -- 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 -> (newton i x a)) 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 Hs0703.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