スプライン補間

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

スプライン補間

複数の点が与えられているとき, これらを１つの多項式でつなげるのではなく, ２点ずつの区間に分けて, それぞれの区間ごとに別々の３次式を設定する.

まず, 連立1次方程式をたてる.

その連立1次方程式を解いて $ s_i $ を求める.

区間 $ x_i $ ～ $ x_{i+1} $ の３次式は, 得られた $ s_i $ を使って, 次のように表すことができる.

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

例題として,

を近似する.

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 a(): ReDim a(N)
Dim b(): ReDim b(N)
Dim c(): ReDim c(N)
Dim d(): ReDim d(N)
For i = 1 To N - 1
    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)))
Next
'３項方程式を解く (ト－マス法)
Dim g(): ReDim g(N)
Dim s(): ReDim s(N)
g(1) = b(1)
s(1) = d(1)
For i = 2 To N - 1
    g(i) = b(i) - a(i) * c(i-1) / g(i-1)
    s(i) = d(i) - a(i) * s(i-1) / g(i-1)
Next
Dim z(): ReDim z(N)
z(0)   = 0
z(N)   = 0
z(N-1) = s(N-1) / g(N-1)
For i = N - 2 To 1 Step -1
    z(i) = (s(i) - c(i) * z(i+1)) / g(i)
Next

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

    '元の関数と比較
    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

'Spline (スプライン) 補間
Private Function spline(ByVal d, ByVal x(), ByVal y(), ByVal z())
    '補間関数値がどの区間にあるか
    Dim k: k = -1
    Dim i
    For i = 1 To N
        If d <= x(i) Then
            k = i - 1
            Exit For
        End If
    Next
    If k < 0 Then k = N

    Dim d1: d1 =  x(k+1) - d
    Dim d2: d2 =  d      - x(k)
    Dim d3: d3 =  x(k+1) - x(k)
    spline     = (z(k) * (d1 ^ 3) + z(k+1) * (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 Function

Z:\>cscript //nologo Z:\0705.vbs
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -2.90573     1.03906
-3.50   -0.73099    -1.41849     0.68750
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.39714    -0.31250
-2.00   -0.93333    -0.70677    -0.22656
-1.50   -1.00078    -1.00078    -0.00000
-1.00   -0.84167    -0.92760     0.08594
-0.50   -0.47943    -0.54193     0.06250
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.54193    -0.06250
 1.00    0.84167     0.92760    -0.08594
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.70677     0.22656
 2.50    0.70964     0.39714     0.31250
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     1.41849    -0.68750
 4.00    1.86667     2.90573    -1.03906
 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 d1 = i * 1.5 - 4.5
    x[i]   = d1
    y[i]   = f(d1)
}

// ３項方程式の係数の表を作る
var a = []
var b = []
var c = []
var d = []
for (var i = 1; i < N - 1; i++)
{
    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]))
}
// ３項方程式を解く (ト－マス法)
var g = []
var s = []
g[1] = b[1]
s[1] = d[1]
for (var i = 2; i < N - 1; i++)
{
    g[i] = b[i] - a[i] * c[i-1] / g[i-1]
    s[i] = d[i] - a[i] * s[i-1] / g[i-1]
}
var z  = []
z[0]   = 0
z[N-1] = 0
z[N-2] = s[N-2] / g[N-2]
for (var i = N - 3; i >= 1; i--)
    z[i] = (s[i] - c[i] * z[i+1]) / g[i]

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

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

// Spline (スプライン) 補間
function spline(d, x, y, z)
{
    // 補間関数値がどの区間にあるか
    var k = -1
    for (var i = 1; i < N; i++)
    {
        if (d <= x[i])
        {
            k = i - 1
            break
        }
    }
    if (k < 0) k = N - 1

    var d1 = x[k+1] - d
    var d2 = d      - x[k]
    var d3 = x[k+1] - x[k]
    return  (z[k] * Math.pow(d1,3) + z[k+1] * Math.pow(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
}

Z:\>cscript //nologo Z:\0705.js
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -2.90573     1.03906
-3.50   -0.73099    -1.41849     0.68750
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.39714    -0.31250
-2.00   -0.93333    -0.70677    -0.22656
-1.50   -1.00078    -1.00078    -0.00000
-1.00   -0.84167    -0.92760     0.08594
-0.50   -0.47943    -0.54193     0.06250
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.54193    -0.06250
 1.00    0.84167     0.92760    -0.08594
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.70677     0.22656
 2.50    0.70964     0.39714     0.31250
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     1.41849    -0.68750
 4.00    1.86667     2.90573    -1.03906
 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)
}

# Spline (スプライン) 補間
function spline($d, $x, $y, $z)
{
    # 補間関数値がどの区間にあるか
    $k = -1
    foreach ($i in 1..($N - 1))
    {
        if ($d -le $x[$i])
        {
            $k = $i - 1
            break
        }
    }
    if ($k -lt 0)
    {
        $k = $N - 1
    }

    $d1 = $x[($k+1)] - $d
    $d2 = $d         - $x[$k]
    $d3 = $x[($k+1)] - $x[$k]
    ($z[$k] * [Math]::Pow($d1,3) + $z[($k+1)] * [Math]::Pow($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
}

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

# ３項方程式の係数の表を作る
$a = New-Object double[] $N
$b = New-Object double[] $N
$c = New-Object double[] $N
$d = New-Object double[] $N
foreach ($i in 1..($N - 2))
{
    $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)]))
}
# ３項方程式を解く (ト－マス法)
$g = New-Object double[] $N
$s = New-Object double[] $N
$g[1] = $b[1]
$s[1] = $d[1]
foreach ($i in 2..($N - 2))
{
    $g[$i] = $b[$i] - $a[$i] * $c[($i-1)] / $g[($i-1)]
    $s[$i] = $d[$i] - $a[$i] * $s[($i-1)] / $g[($i-1)]
}
$z = New-Object double[] $N
$z[0] = 0
$z[($N-1)] = 0
$z[($N-2)] = $s[($N-2)] / $g[($N-2)]
for ($i = $N - 3; $i -ge 1; $i--)
{
    $z[$i] = ($s[$i] - $c[$i] * $z[($i+1)]) / $g[$i]
}

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

    # 元の関数と比較
    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:\0705.ps1
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -2.90573     1.03906
-3.50   -0.73099    -1.41849     0.68750
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.39714    -0.31250
-2.00   -0.93333    -0.70677    -0.22656
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.92760     0.08594
-0.50   -0.47943    -0.54193     0.06250
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.54193    -0.06250
 1.00    0.84167     0.92760    -0.08594
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.70677     0.22656
 2.50    0.70964     0.39714     0.31250
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     1.41849    -0.68750
 4.00    1.86667     2.90573    -1.03906
 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 @a = ();
my @b = ();
my @c = ();
my @d = ();
for $i (1..(N - 1))
{
    $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]));
}
# ３項方程式を解く (ト－マス法)
my @g = ();
my @s = ();
$g[1] = $b[1];
$s[1] = $d[1];
for $i (2..(N - 1))
{
    $g[$i] = $b[$i] - $a[$i] * $c[$i-1] / $g[$i-1];
    $s[$i] = $d[$i] - $a[$i] * $s[$i-1] / $g[$i-1];
}
my @z = ();
$z[0] = 0;
$z[N] = 0;
$z[N-1] = $s[N-1] / $g[N-1];
for ($i = N - 2; $i >= 1; $i--)
{
    $z[$i] = ($s[$i] - $c[$i] * $z[$i+1]) / $g[$i];
}

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

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

# Spline (スプライン) 補間
sub spline
{
    my ($d, $x, $y, $z) = @_;

    # 補間関数値がどの区間にあるか
    my $k = -1;
    for $i (1..N)
    {
        if ($d <= $$x[$i])
        {
            $k = $i - 1;
            last;
        }
    }
    if ($k < 0)
    {
        $k = N;
    }

    $d1 = $$x[$k+1] - $d;
    $d2 = $d        - $$x[$k];
    $d3 = $$x[$k+1] - $$x[$k];
    ($$z[$k] * ($d1 ** 3) + $$z[$k+1] * ($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;
}

Z:\>perl Z:\0705.pl
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -2.90573     1.03906
-3.50   -0.73099    -1.41849     0.68750
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.39714    -0.31250
-2.00   -0.93333    -0.70677    -0.22656
-1.50   -1.00078    -1.00078    -0.00000
-1.00   -0.84167    -0.92760     0.08594
-0.50   -0.47943    -0.54193     0.06250
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.54193    -0.06250
 1.00    0.84167     0.92760    -0.08594
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.70677     0.22656
 2.50    0.70964     0.39714     0.31250
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     1.41849    -0.68750
 4.00    1.86667     2.90573    -1.03906
 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);
}

# ３項方程式の係数の表を作る
$a = array();
$b = array();
$c = array();
$d = array();
foreach (range(1, N - 1) as $i)
{
    $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]));
}
# ３項方程式を解く (ト－マス法)
$g = array();
$s = array();
$g[1] = $b[1];
$s[1] = $d[1];
foreach (range(2, N - 1) as $i)
{
    $g[$i] = $b[$i] - $a[$i] * $c[$i-1] / $g[$i-1];
    $s[$i] = $d[$i] - $a[$i] * $s[$i-1] / $g[$i-1];
}
$z = array();
$z[0] = 0;
$z[N] = 0;
$z[N-1] = $s[N-1] / $g[N-1];
for ($i = N - 2; $i >= 1; $i--)
{
    $z[$i] = ($s[$i] - $c[$i] * $z[$i+1]) / $g[$i];
}

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

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

# Spline (スプライン) 補間
function spline($d, $x, $y, $z)
{
    # 補間関数値がどの区間にあるか
    $k = -1;
    foreach (range(1, N) as $i)
    {
        if ($d <= $x[$i])
        {
            $k = $i - 1;
            break;
        }
    }
    if ($k < 0) $k = N;

    $d1 = $x[$k+1] - $d;
    $d2 = $d       - $x[$k];
    $d3 = $x[$k+1] - $x[$k];
    return ($z[$k] * pow($d1, 3) + $z[$k+1] * pow($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;
}
?>

Z:\>php Z:\0705.php
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -2.90573     1.03906
-3.50   -0.73099    -1.41849     0.68750
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.39714    -0.31250
-2.00   -0.93333    -0.70677    -0.22656
-1.50   -1.00078    -1.00078    -0.00000
-1.00   -0.84167    -0.92760     0.08594
-0.50   -0.47943    -0.54193     0.06250
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.54193    -0.06250
 1.00    0.84167     0.92760    -0.08594
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.70677     0.22656
 2.50    0.70964     0.39714     0.31250
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     1.41849    -0.68750
 4.00    1.86667     2.90573    -1.03906
 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)

# Spline (スプライン) 補間
def spline(d, x, y, z):
    # 補間関数値がどの区間にあるか
    k = -1
    for i in range(1, N):
        if d <= x[i]:
            k = i - 1
            break
    if (k < 0):
        k = N - 1

    d1 = x[k+1] - d
    d2 = d      - x[k]
    d3 = x[k+1] - x[k]
    return ( (z[k] * (d1 ** 3) + z[k+1] * (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)

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)

# ３項方程式の係数の表を作る
a = [0 for i in range(N)]
b = [0 for i in range(N)]
c = [0 for i in range(N)]
d = [0 for i in range(N)]
for i in range(1, N - 1):
    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]))

# ３項方程式を解く (ト－マス法)
g = [0 for i in range(N)]
s = [0 for i in range(N)]
g[1] = b[1]
s[1] = d[1]
for i in range(2, N - 1):
    g[i] = b[i] - a[i] * c[i-1] / g[i-1]
    s[i] = d[i] - a[i] * s[i-1] / g[i-1]

z = [0 for i in range(N)]
z[0]   = 0
z[N-1] = 0
z[N-2] = s[N-2] / g[N-2]
i = N - 3
while i >= 1:
    z[i] = (s[i] - c[i] * z[i+1]) / g[i]
    i -= 1

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

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

Z:\>python Z:\0705.py
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -2.90573     1.03906
-3.50   -0.73099    -1.41849     0.68750
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.39714    -0.31250
-2.00   -0.93333    -0.70677    -0.22656
-1.50   -1.00078    -1.00078    -0.00000
-1.00   -0.84167    -0.92760     0.08594
-0.50   -0.47943    -0.54193     0.06250
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.54193    -0.06250
 1.00    0.84167     0.92760    -0.08594
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.70677     0.22656
 2.50    0.70964     0.39714     0.31250
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     1.41849    -0.68750
 4.00    1.86667     2.90573    -1.03906
 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

# Spline (スプライン) 補間
def spline(d, x, y, z)
    # 補間関数値がどの区間にあるか
    k = -1
    (1..N).each do |i|
        if d <= x[i]
            k = i - 1
            break
        end
    end
    if k < 0
        k = N
    end

    d1 =  x[k+1] - d
    d2 =  d      - x[k]
    d3 =  x[k+1] - x[k]
    (z[k] * (d1 ** 3) + z[k+1] * (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

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

# ３項方程式の係数の表を作る
a = Array.new(N)
b = Array.new(N)
c = Array.new(N)
d = Array.new(N)
(1..(N - 1)).each do |i|
    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 = Array.new(N)
s = Array.new(N)
g[1] = b[1]
s[1] = d[1]
(2..(N - 1)).each do |i|
    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 = Array.new(N)
z[0]   = 0
z[N]   = 0
z[N-1] = s[N-1] / g[N-1]
i = N - 2
while i > 0 do
    z[i] = (s[i] - c[i] * z[i+1]) / g[i]
    i -= 1
end

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

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

Z:\>ruby Z:\0705.rb
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -2.90573     1.03906
-3.50   -0.73099    -1.41849     0.68750
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.39714    -0.31250
-2.00   -0.93333    -0.70677    -0.22656
-1.50   -1.00078    -1.00078    -0.00000
-1.00   -0.84167    -0.92760     0.08594
-0.50   -0.47943    -0.54193     0.06250
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.54193    -0.06250
 1.00    0.84167     0.92760    -0.08594
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.70677     0.22656
 2.50    0.70964     0.39714     0.31250
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     1.41849    -0.68750
 4.00    1.86667     2.90573    -1.03906
 4.50    4.68984     4.68984     0.00000

Groovy

Pascal

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.

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

Z:\>Pas0705
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -2.90573     1.03906
-3.50   -0.73099    -1.41849     0.68750
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.39714    -0.31250
-2.00   -0.93333    -0.70677    -0.22656
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.92760     0.08594
-0.50   -0.47943    -0.54193     0.06250
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.54193    -0.06250
 1.00    0.84167     0.92760    -0.08594
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.70677     0.22656
 2.50    0.70964     0.39714     0.31250
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     1.41849    -0.68750
 4.00    1.86667     2.90573    -1.03906
 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 Ada0705 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;
    b : Long_Float_Array;
    c : Long_Float_Array;
    d : Long_Float_Array;
    g : Long_Float_Array;
    s : Long_Float_Array;
    z : Long_Float_Array;
    d1, d2, d3 : Long_Float;
    k : Integer;

    -- 元の関数
    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;

    -- Spline (スプライン) 補間
    function spline(d:Long_Float; x:Long_Float_Array; y:Long_Float_Array; z:Long_Float_Array) return Long_Float is
        d1, d2, d3 : Long_Float;
        k : Integer;
    begin
        -- 補間関数値がどの区間にあるか
        k := -1;
        for i in 1 .. x'Last loop
            if d <= x(i) then
                k := i - 1;
                exit;
            end if;
        end loop;
        if k < 0 then
            k := N;
        end if;

        d1   := x(k+1) - d;
        d2   := d      - x(k);
        d3   := x(k+1) - x(k);
        return (z(k) * (d1 ** 3) + z(k+1) * (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;
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 i in 1..x'Last-1 loop
        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 loop;
    -- ３項方程式を解く (ト－マス法)
    g(1) := b(1);
    s(1) := d(1);
    for i in 2..a'Last-1 loop
        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 loop;
    z(0)   := 0.0;
    z(N)   := 0.0;
    z(N-1) := s(N-1) / g(N-1);
    k := N - 2;
    while k > 0 loop
        z(k) := (s(k) - c(k) * z(k+1)) / g(k);
        k := k - 1;
    end loop;

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

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

xxxxxx@yyyyyy /Z
$ gnatmake Ada0705.adb

xxxxxx@yyyyyy /Z
$ Ada0705
-4.50    -4.68984    -4.68984     0.00000
-4.00    -1.86667    -2.90573     1.03906
-3.50    -0.73099    -1.41849     0.68750
-3.00    -0.52500    -0.52500     0.00000
-2.50    -0.70964    -0.39714    -0.31250
-2.00    -0.93333    -0.70677    -0.22656
-1.50    -1.00078    -1.00078     0.00000
-1.00    -0.84167    -0.92760     0.08594
-0.50    -0.47943    -0.54193     0.06250
 0.00     0.00000     0.00000     0.00000
 0.50     0.47943     0.54193    -0.06250
 1.00     0.84167     0.92760    -0.08594
 1.50     1.00078     1.00078     0.00000
 2.00     0.93333     0.70677     0.22656
 2.50     0.70964     0.39714     0.31250
 3.00     0.52500     0.52500     0.00000
 3.50     0.73099     1.41849    -0.68750
 4.00     1.86667     2.90573    -1.03906
 4.50     4.68984     4.68984     0.00000

VB.NET

Module VB0705
    'データ点の数 - 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 a(N) As Double
        Dim b(N) As Double
        Dim c(N) As Double
        Dim d(N) As Double
        For i As Integer = 1 To N - 1
            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)))
        Next
        '３項方程式を解く (ト－マス法)
        Dim g(N) As Double
        Dim s(N) As Double
        g(1) = b(1)
        s(1) = d(1)
        For i As Integer = 2 To N - 1
            g(i) = b(i) - a(i) * c(i-1) / g(i-1)
            s(i) = d(i) - a(i) * s(i-1) / g(i-1)
        Next
        Dim z(N) As Double
        z(0)   = 0
        z(N)   = 0
        z(N-1) = s(N-1) / g(N-1)
        For i As Integer = N - 2 To 1 Step -1
            z(i) = (s(i) - c(i) * z(i+1)) / g(i)
        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 = spline(d1, x, y, z)

            '元の関数と比較
            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

    'Spline (スプライン) 補間
    Private Function spline(ByVal d As Double, ByVal x() As Double, ByVal y() As Double, ByVal z() As Double) As Double
        '補間関数値がどの区間にあるか
        Dim k As Integer = -1
        For i As Integer = 1 To N
            If d <= x(i) Then
                k = i - 1
                Exit For
            End If
        Next
        If (k < 0) Then k = N

        Dim d1 As Double = x(k+1) - d
        Dim d2 As Double = d      - x(k)
        Dim d3 As Double = x(k+1) - x(k)
        Return            (z(k) * (d1 ^ 3) + z(k+1) * (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 Function
End Module

Z:\>vbc -nologo VB0705.vb

Z:\>VB0705
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -2.90573     1.03906
-3.50   -0.73099    -1.41849     0.68750
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.39714    -0.31250
-2.00   -0.93333    -0.70677    -0.22656
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.92760     0.08594
-0.50   -0.47943    -0.54193     0.06250
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.54193    -0.06250
 1.00    0.84167     0.92760    -0.08594
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.70677     0.22656
 2.50    0.70964     0.39714     0.31250
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     1.41849    -0.68750
 4.00    1.86667     2.90573    -1.03906
 4.50    4.68984     4.68984     0.00000

C#

using System;

public class CS0705
{
    // データ点の数
    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[] a = new double[N];
        double[] b = new double[N];
        double[] c = new double[N];
        double[] d = new double[N];
        for (int i = 1; i < N - 1; i++)
        {
            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]));
        }
        // ３項方程式を解く (ト－マス法)
        double[] g = new double[N];
        double[] s = new double[N];
        g[1] = b[1];
        s[1] = d[1];
        for (int i = 2; i < N - 1; i++)
        {
            g[i] = b[i] - a[i] * c[i-1] / g[i-1];
            s[i] = d[i] - a[i] * s[i-1] / g[i-1];
        }
        double[] z = new double[N];
        z[0]   = 0;
        z[N-1] = 0;
        z[N-2] = s[N-2] / g[N-2];
        for (int i = N - 3; i >= 1; i--)
        {
            z[i] = (s[i] - c[i] * z[i+1]) / g[i];
        }

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

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

    // Spline (スプライン) 補間
    private static double spline(double d, double[] x, double[] y, double[] z)
    {
        // 補間関数値がどの区間にあるか
        int k = -1;
        for (int i = 1; i < N; i++)
        {
            if (d <= x[i])
            {
                k = i - 1;
                break;
            }
        }
        if (k < 0) k = N - 1;

        double d1 = x[k+1] - d;
        double d2 = d      - x[k];
        double d3 = x[k+1] - x[k];
        return      (z[k] * Math.Pow(d1,3) + z[k+1] * Math.Pow(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;
    }
}

Z:\>csc -nologo CS0705.cs

Z:\>CS0705
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -2.90573     1.03906
-3.50   -0.73099    -1.41849     0.68750
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.39714    -0.31250
-2.00   -0.93333    -0.70677    -0.22656
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.92760     0.08594
-0.50   -0.47943    -0.54193     0.06250
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.54193    -0.06250
 1.00    0.84167     0.92760    -0.08594
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.70677     0.22656
 2.50    0.70964     0.39714     0.31250
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     1.41849    -0.68750
 4.00    1.86667     2.90573    -1.03906
 4.50    4.68984     4.68984     0.00000

Java

public class Java0705 {

    // データ点の数
    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[] a = new double[N];
        double[] b = new double[N];
        double[] c = new double[N];
        double[] d = new double[N];
        for (int i = 1; i < N - 1; i++) {
            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]));
        }
        // ３項方程式を解く (ト－マス法)
        double[] g = new double[N];
        double[] s = new double[N];
        g[1] = b[1];
        s[1] = d[1];
        for (int i = 2; i < N - 1; i++) {
            g[i] = b[i] - a[i] * c[i-1] / g[i-1];
            s[i] = d[i] - a[i] * s[i-1] / g[i-1];
        }
        double[] z = new double[N];
        z[0]   = 0;
        z[N-1] = 0;
        z[N-2] = s[N-2] / g[N-2];
        for (int i = N - 3; i >= 1; i--)
            z[i] = (s[i] - c[i] * z[i+1]) / g[i];

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

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

    // Spline (スプライン) 補間
    private static double spline(double d, double[] x, double[] y, double[] z) {
        // 補間関数値がどの区間にあるか
        int k = -1;
        for (int i = 1; i < N; i++) {
            if (d <= x[i]) {
                k = i - 1;
                break;
            }
        }
        if (k < 0) k = N - 1;

        double d1 = x[k+1] - d;
        double d2 = d      - x[k];
        double d3 = x[k+1] - x[k];
        return      (z[k] * Math.pow(d1,3) + z[k+1] * Math.pow(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;
    }
}

Z:\>javac Java0705.java

Z:\>java Java0705
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -2.90573     1.03906
-3.50   -0.73099    -1.41849     0.68750
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.39714    -0.31250
-2.00   -0.93333    -0.70677    -0.22656
-1.50   -1.00078    -1.00078    -0.00000
-1.00   -0.84167    -0.92760     0.08594
-0.50   -0.47943    -0.54193     0.06250
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.54193    -0.06250
 1.00    0.84167     0.92760    -0.08594
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.70677     0.22656
 2.50    0.70964     0.39714     0.31250
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     1.41849    -0.68750
 4.00    1.86667     2.90573    -1.03906
 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);

// Spline (スプライン) 補間
double spline(double d, double x[], double y[], double z[]);

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 a[N], b[N], c[N], d[N];
    for (int i = 1; i < N - 1; i++)
    {
        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]));
    }
    // ３項方程式を解く (ト－マス法)
    double g[N], s[N];
    g[1] = b[1];
    s[1] = d[1];
    for (int i = 2; i < N - 1; i++)
    {
        g[i] = b[i] - a[i] * c[i-1] / g[i-1];
        s[i] = d[i] - a[i] * s[i-1] / g[i-1];
    }
    double z[N];
    z[0]   = 0;
    z[N-1] = 0;
    z[N-2] = s[N-2] / g[N-2];
    for (int i = N - 3; i >= 1; i--)
    {
        z[i] = (s[i] - c[i] * z[i+1]) / g[i];
    }

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

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

// Spline (スプライン) 補間
double spline(double d, double x[], double y[], double z[])
{
    // 補間関数値がどの区間にあるか
    int k = -1;
    for (int i = 1; i < N; i++)
    {
        if (d <= x[i])
        {
            k = i - 1;
            break;
        }
    }
    if (k < 0) k = N - 1;

    double d1 = x[k+1] - d;
    double d2 = d      - x[k];
    double d3 = x[k+1] - x[k];
    return      (z[k] * pow(d1,3) + z[k+1] * pow(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;
}

Z:\>bcc32 -q CP0705.cpp
cp0705.cpp:

Z:\>CP0705
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -2.90573     1.03906
-3.50   -0.73099    -1.41849     0.68750
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.39714    -0.31250
-2.00   -0.93333    -0.70677    -0.22656
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.92760     0.08594
-0.50   -0.47943    -0.54193     0.06250
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.54193    -0.06250
 1.00    0.84167     0.92760    -0.08594
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.70677     0.22656
 2.50    0.70964     0.39714     0.31250
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     1.41849    -0.68750
 4.00    1.86667     2.90573    -1.03906
 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);

// Spline (スプライン) 補間
double spline(double d, double x[], double y[], double z[]);

int main()
{
    int i;
    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 a[N], b[N], c[N], d[N];
    for (i = 1; i < N - 1; i++)
    {
        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]));
    }
    // ３項方程式を解く (ト－マス法)
    double g[N], s[N];
    g[1] = b[1];
    s[1] = d[1];
    for (i = 2; i < N - 1; i++)
    {
        g[i] = b[i] - a[i] * c[i-1] / g[i-1];
        s[i] = d[i] - a[i] * s[i-1] / g[i-1];
    }
    double z[N];
    z[0]   = 0;
    z[N-1] = 0;
    z[N-2] = s[N-2] / g[N-2];
    for (i = N - 3; i >= 1; i--)
    {
        z[i] = (s[i] - c[i] * z[i+1]) / g[i];
    }

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

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

// Spline (スプライン) 補間
double spline(double d, double x[], double y[], double z[])
{
    // 補間関数値がどの区間にあるか
    int i, k = -1;
    for (i = 1; i < N; i++)
    {
        if (d <= x[i])
        {
            k = i - 1;
            break;
        }
    }
    if (k < 0) k = N - 1;

    double d1 = x[k+1] - d;
    double d2 = d      - x[k];
    double d3 = x[k+1] - x[k];
    return      (z[k] * pow(d1,3) + z[k+1] * pow(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;
}

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

xxxxxx@yyyyyy /Z
$ OC0705
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -2.90573     1.03906
-3.50   -0.73099    -1.41849     0.68750
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.39714    -0.31250
-2.00   -0.93333    -0.70677    -0.22656
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.92760     0.08594
-0.50   -0.47943    -0.54193     0.06250
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.54193    -0.06250
 1.00    0.84167     0.92760    -0.08594
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.70677     0.22656
 2.50    0.70964     0.39714     0.31250
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     1.41849    -0.68750
 4.00    1.86667     2.90573    -1.03906
 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 a[N];
    double b[N];
    double c[N];
    double d[N];
    for (int i = 1; i < N - 1; i++)
    {
        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]));
    }
    // ３項方程式を解く (ト－マス法)
    double g[N];
    double s[N];
    g[1] = b[1];
    s[1] = d[1];
    for (int i = 2; i < N - 1; i++)
    {
        g[i] = b[i] - a[i] * c[i-1] / g[i-1];
        s[i] = d[i] - a[i] * s[i-1] / g[i-1];
    }
    double z[N];
    z[0]   = 0;
    z[N-1] = 0;
    z[N-2] = s[N-2] / g[N-2];
    for (int i = N - 3; i >= 1; i--)
    {
        z[i] = (s[i] - c[i] * z[i+1]) / g[i];
    }

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

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

// Spline (スプライン) 補間
double spline(double d, double x[], double y[], double z[])
{
    // 補間関数値がどの区間にあるか
    int k = -1;
    for (int i = 1; i < N; i++)
    {
        if (d <= x[i])
        {
            k = i - 1;
            break;
        }
    }
    if (k < 0) k = N - 1;

    double d1 = x[k+1] - d;
    double d2 = d      - x[k];
    double d3 = x[k+1] - x[k];
    return      (z[k] * pow(d1,3) + z[k+1] * pow(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;
}

Z:\>dmd D0705.d

Z:\>D0705
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -2.90573     1.03906
-3.50   -0.73099    -1.41849     0.68750
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.39714    -0.31250
-2.00   -0.93333    -0.70677    -0.22656
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.92760     0.08593
-0.50   -0.47943    -0.54193     0.06250
 0.00    0.00000    -0.00000     0.00000
 0.50    0.47943     0.54193    -0.06250
 1.00    0.84167     0.92760    -0.08593
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.70677     0.22656
 2.50    0.70964     0.39714     0.31250
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     1.41849    -0.68750
 4.00    1.86667     2.90573    -1.03906
 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 a [N]float64
    var b [N]float64
    var c [N]float64
    var d [N]float64
    for i := 1; i < N - 1; i++ {
        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]))
    }

    // ３項方程式を解く (ト－マス法)
    var g [N]float64
    var s [N]float64
    g[1] = b[1]
    s[1] = d[1]
    for i := 2; i < N - 1; i++ {
        g[i] = b[i] - a[i] * c[i-1] / g[i-1]
        s[i] = d[i] - a[i] * s[i-1] / g[i-1]
    }
    var z [N]float64
    z[0]   = 0
    z[N-1] = 0
    z[N-2] = s[N-2] / g[N-2]
    for i := N - 3; i >= 1; i-- {
        z[i] = (s[i] - c[i] * z[i+1]) / g[i]
    }

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

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

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

// Spline (スプライン) 補間
func spline(d float64, x []float64, y []float64, z []float64) float64 {
    // 補間関数値がどの区間にあるか
    k := -1
    for i := 1; i < N; i++ {
        if d <= x[i] {
            k = i - 1
            break
        }
    }
    if k < 0 {
        k = N - 1
    }

    var d1 float64 = x[k+1] - d
    var d2 float64 = d      - x[k]
    var d3 float64 = x[k+1] - x[k]
    return (z[k] * math.Pow(d1,3) + z[k+1] * math.Pow(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
}

Z:\>8g GO0705.go

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

Z:\>GO0705
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -2.90573     1.03906
-3.50   -0.73099    -1.41849     0.68750
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.39714    -0.31250
-2.00   -0.93333    -0.70677    -0.22656
-1.50   -1.00078    -1.00078    -0.00000
-1.00   -0.84167    -0.92760     0.08594
-0.50   -0.47943    -0.54193     0.06250
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.54193    -0.06250
 1.00    0.84167     0.92760    -0.08594
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.70677     0.22656
 2.50    0.70964     0.39714     0.31250
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     1.41849    -0.68750
 4.00    1.86667     2.90573    -1.03906
 4.50    4.68984     4.68984     0.00000

Scala

object Scala0705 {
    // データ点の数 - 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 a = Array.ofDim[Double](N)
        val b = Array.ofDim[Double](N)
        val c = Array.ofDim[Double](N)
        val d = Array.ofDim[Double](N)
        for (i <- 1 to N - 1) {
            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)))
        }
        // ３項方程式を解く (ト－マス法)
        val g = Array.ofDim[Double](N)
        val s = Array.ofDim[Double](N)
        g(1) = b(1)
        s(1) = d(1)
        for (i <- 2 to N - 1) {
            g(i) = b(i) - a(i) * c(i-1) / g(i-1)
            s(i) = d(i) - a(i) * s(i-1) / g(i-1)
        }
        val z = Array.ofDim[Double](N + 1)
        z(0)   = 0
        z(N)   = 0
        z(N-1) = s(N-1) / g(N-1)
        for (i <- N - 2 to 1 by -1)
            z(i) = (s(i) - c(i) * z(i+1)) / g(i)

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

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

    // Spline (スプライン) 補間
    def spline(d:Double, x:IndexedSeq[Double], y:IndexedSeq[Double], z:IndexedSeq[Double]) = {
        // 補間関数値がどの区間にあるか
        var k = -1
        for (i <- N to 1 by -1) {
            if (d <= x(i)) k = i - 1
        }
        if (k < 0) k = N

        val d1 = x(k+1) - d
        val d2 = d      - x(k)
        val d3 = x(k+1) - x(k)
        (z(k) * Math.pow(d1,3) + z(k+1) * Math.pow(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
    }
}

Z:\>scala Scala0705.scala
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -2.90573     1.03906
-3.50   -0.73099    -1.41849     0.68750
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.39714    -0.31250
-2.00   -0.93333    -0.70677    -0.22656
-1.50   -1.00078    -1.00078    -0.00000
-1.00   -0.84167    -0.92760     0.08594
-0.50   -0.47943    -0.54193     0.06250
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.54193    -0.06250
 1.00    0.84167     0.92760    -0.08594
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.70677     0.22656
 2.50    0.70964     0.39714     0.31250
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     1.41849    -0.68750
 4.00    1.86667     2.90573    -1.03906
 4.50    4.68984     4.68984     0.00000

F#

module Fs0705

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

// Spline (スプライン) 補間
let spline(d:double) (x:double list) (y:double list) (z:double list) =
    // 補間関数値がどの区間にあるか
    let mutable k = -1
    for i in [N..(-1)..1] do
        if d <= x.[i] then k <- i - 1
    if k < 0 then k <- N

    let d1 = x.[k+1] - d
    let d2 = d       - x.[k]
    let d3 = x.[k+1] - x.[k]
    (z.[k] * Math.Pow(d1,3.0) + z.[k+1] * Math.Pow(d2,3.0)) / (6.0 * d3) +
    (y.[k]   / d3 - z.[k]   * d3 / 6.0) * d1 +
    (y.[k+1] / d3 - z.[k+1] * d3 / 6.0) * d2

// 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 a = Array.zeroCreate<double> N
let b = Array.zeroCreate<double> N
let c = Array.zeroCreate<double> N
let d = Array.zeroCreate<double> N
for i in [1..N - 1] do
    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]))
// ３項方程式を解く (ト－マス法)
let g = Array.zeroCreate<double> N
let s = Array.zeroCreate<double> N
g.[1] <- b.[1]
s.[1] <- d.[1]
for i in [2..N - 1] do
    g.[i] <- b.[i] - a.[i] * c.[i-1] / g.[i-1]
    s.[i] <- d.[i] - a.[i] * s.[i-1] / g.[i-1]
let z = Array.zeroCreate<double> (N+1)
z.[0]   <- 0.0
z.[N]   <- 0.0
z.[N-1] <- s.[N-1] / g.[N-1]
for i in [(N-2)..(-1)..1] do
    z.[i] <- (s.[i] - c.[i] * z.[i+1]) / g.[i]

// 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 -> (spline d x y (Array.toList(z))))

(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 Fs0705.fs
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -2.90573     1.03906
-3.50   -0.73099    -1.41849     0.68750
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.39714    -0.31250
-2.00   -0.93333    -0.70677    -0.22656
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.92760     0.08594
-0.50   -0.47943    -0.54193     0.06250
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.54193    -0.06250
 1.00    0.84167     0.92760    -0.08594
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.70677     0.22656
 2.50    0.70964     0.39714     0.31250
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     1.41849    -0.68750
 4.00    1.86667     2.90573    -1.03906
 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))))))

; Spline (スプライン) 補間
(defn spline [d x y z]
    ; 補間関数値がどの区間にあるか
    (def k -1)
    (doseq [i (range 1 N)]
        (if (and (<= d (nth x i)) (< k 0))
            (def k (- i 1))
        )
    )
    (if (< k 0) (def k N))

    (def d1 (- (nth x (+ k 1)) d))
    (def d2 (- d               (nth x k)))
    (def d3 (- (nth x (+ k 1)) (nth x k)))

    (def a1 (*  (nth z k)       (Math/pow d1 3.0)))
    (def a2 (*  (nth z (+ k 1)) (Math/pow d2 3.0)))
    (def a3 (/ (+ a1 a2) (* 6.0 d3)))

    (def b1 (/ (nth y k) d3))
    (def b2 (* (nth z k) d3))
    (def b3 (/ b2 6.0))
    (def b4 (* (- b1 b3) d1))

    (def c1 (/ (nth y (+ k 1)) d3))
    (def c2 (* (nth z (+ k 1)) d3))
    (def c3 (/ c2 6.0))
    (def c4 (* (- c1 c3) d2))

    (+ (+ a3 b4) c4)
)

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

; ３項方程式の係数の表を作る
(def a (list 0.0))
(def b (list 0.0))
(def c (list 0.0))
(def d (list 0.0))
(doseq [i (range 1 (- N 1))]
    (def a (cons          (- (nth x i      ) (nth x (- i 1)))  a))
    (def b (cons  (* 2.0  (- (nth x (+ i 1)) (nth x (- i 1)))) b))
    (def c (cons          (- (nth x (+ i 1)) (nth x i      ))  c))
    (def d (cons  (* 6.0  (- (/ (- (nth y(+ i 1)) (nth y i))
                                (- (nth x(+ i 1)) (nth x i)))
                            (/  (- (nth y i)      (nth y (- i 1)))
                                (- (nth x i)      (nth x (- i 1)))))) d)))
(def a (reverse a))
(def b (reverse b))
(def c (reverse c))
(def d (reverse d))

; ３項方程式を解く (ト－マス法)
(def g (list (nth b 1) 0.0))
(def s (list (nth d 1) 0.0))
(doseq [i (range 2 (- N 1))]
    (def w (nth g 0))
    (def g (cons    (-  (nth b i)
                        (/  (*  (nth a i)
                                (nth c (- i 1)))
                            w))
                    g))
    (def s (cons    (-  (nth d i)
                        (/  (*  (nth a i)
                                (nth s 0))
                            w))
                    s))
)
(def g (reverse g))
(def s (reverse s))

(def z (list 0.0))
(def z (cons (/ (nth s (- N 2)) (nth g (- N 2))) z))
(doseq [i (reverse (range 1 (- N 2)))]
    (def z (cons    (/  (-  (nth s i)
                            (* (nth c i) (nth z 0)))
                        (nth g i))
                    z))
)
(def z (cons 0.0 z))

; 0.5刻みで 与えられていない値を補間
(def d1 (map #(- (* % 0.5) 4.5) (range 0 19)))
(def d2 (map #(f %) d1))
(def d3 (map #(spline % x y z) 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 Clj0705.clj
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -2.90573     1.03906
-3.50   -0.73099    -1.41849     0.68750
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.39714    -0.31250
-2.00   -0.93333    -0.70677    -0.22656
-1.50   -1.00078    -1.00078    -0.00000
-1.00   -0.84167    -0.92760     0.08594
-0.50   -0.47943    -0.54193     0.06250
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.54193    -0.06250
 1.00    0.84167     0.92760    -0.08594
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.70677     0.22656
 2.50    0.70964     0.39714     0.31250
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     1.41849    -0.68750
 4.00    1.86667     2.90573    -1.03906
 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))

-- ３項方程式を解く (ト－マス法)
thomas::Int->[Double]->[Double]->[Double]->[Double]->[Double]->[Double]->([Double], [Double])
thomas i a b c d gs ss =
    let
        g = b!!i - a!!i * c!!(i-1) / gs!!0
        s = d!!i - a!!i * ss!!0    / gs!!0
    in
        if i == (n - 1) then (reverse (g:gs), reverse (s:ss))
                        else (thomas (i + 1) a b c d (g:gs) (s:ss))

thomas2::Int->[Double]->[Double]->[Double]->[Double]->[Double]
thomas2 i g s c zs =
    let
        z = (s!!i - c!!i * zs!!0) / g!!i
    in
        if i == 1 then 0.0 : z : zs
                  else thomas2 (i - 1) g s c (z:zs)

-- Spline (スプライン) 補間
spline::Double->[Double]->[Double]->[Double]->Double
spline d x y z =
    let
        -- 補間関数値がどの区間にあるか
        k = if d <= x!!0 then 0
                        else (length $ filter (\i -> i < d) $ x) - 1

        d1 = x!!(k+1) - d
        d2 = d        - x!!k
        d3 = x!!(k+1) - x!!k
    in
        (z!!k * (d1 ^ 3) + z!!(k+1) * (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


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 = 0.0 : map(\i ->       x!!i      - x!!(i-1))  [1..n - 1]
    let b = 0.0 : map(\i -> 2 *  (x!!(i+1)  - x!!(i-1))) [1..n - 1]
    let c = 0.0 : map(\i ->       x!!(i+1)  - x!!i)      [1..n - 1]
    let d = 0.0 : map(\i -> 6 *  ((y!!(i+1) - y!!i) / (x!!(i+1) - x!!i) - (y!!i - y!!(i-1)) / (x!!i - x!!(i-1)))) [1..n - 1]

    -- ３項方程式を解く (ト－マス法)
    let (g, s) = (thomas 2 a b c d [b!!1, 0.0] [d!!1, 0.0])
    let z      = (thomas2 (n - 2) g s c [s!!(n-1) / g!!(n-1), 0.0])

    -- 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 -> (spline i x y z)) 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 Hs0705.hs
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -2.90573     1.03906
-3.50   -0.73099    -1.41849     0.68750
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.39714    -0.31250
-2.00   -0.93333    -0.70677    -0.22656
-1.50   -1.00078    -1.00078    -0.00000
-1.00   -0.84167    -0.92760     0.08594
-0.50   -0.47943    -0.54193     0.06250
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.54193    -0.06250
 1.00    0.84167     0.92760    -0.08594
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.70677     0.22656
 2.50    0.70964     0.39714     0.31250
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     1.41849    -0.68750
 4.00    1.86667     2.90573    -1.03906
 4.50    4.68984     4.68984     0.00000

VBScript
JScript
PowerShell
Perl
PHP
Python
Ruby
Groovy

Pascal
Ada
VB.NET

C#
Java
C++
Objective-C
D
Go

Scala
F#
Clojure
Haskell