<?php
print (3 + 5)."\n";
print (3 + 5)."\n";
print (3 - 5)."\n";
print (3 * 5)."\n";
print pow(3, 5)."\n";
print (5 / 3)."\n";
print (int)(5 / 3)."\n";
print (5 % 3)."\n";

echo  3 * 5, "\n";

printf("%3d\n",     3 * 5);
printf("%23.20f\n", 5 / 3);
?>
<?php
$i = 3 * 5;
print "3 * 5 = $i\n";
print "3 * 5 = ".$i."\n";

printf("3 * 5 = %d\n", $i);

echo "3 * 5 = $i\n";
echo "3 * 5 = ".$i."\n";
echo "3 * 5 = ", $i, "\n";
?>
<?php
foreach (range(1, 9) as $i)
{
    echo "$i, ";
}
echo "\n";
?>
<?php
foreach (range(1, 9) as $i)
{
    if ($i % 3 == 0)
    {
        echo "$i, ";
    }
}
echo "\n";
?>
<?php
error_reporting(E_NOTICE);

$sum = 0;
foreach (range(1, 99) as $i)
{
    if ($i % 3 == 0)
    {
        $sum += $i;
    }
}
echo "$sum\n";
?>
<?php
array_map(function ($n) {echo $n, ",";}, range(1,9))
?>
<?php
$ary = array_filter(range(1,9), function ($n) {return ($n % 3 == 0);});
array_map(function ($n) {echo $n, ",";}, $ary);
?>
<?php
$ary = array_filter(range(1,9), function ($n) {return ($n % 3 == 0);});
echo array_reduce($ary, function($x, $y) {return $x + $y;}, 0);
?>
<?php
# 初項:a, 公差:a で 上限:lim の数列の総和を返す関数
function sn($a, $lim)
{
    $n = (int)($lim / $a);     # 項数:n  =  上限:lim / 公差:a
    $l = $n * $a;              # 末項:l  =  項数:n   * 公差:a
    return ($a + $l) * $n / 2; # 総和:sn = (初項:a   + 末項:l) * 項数:n / 2
}

# 3 の倍数の合計
$sum = sn(3, 999);
echo $sum, "\n";
?>
<?php
# 10000 までの 自然数の和
# 項目数 n = 10000
$n = 10000;
echo ($n * ($n + 1) / 2), "\n";
?>
<?php
# 10000 までの 偶数の和
# 項目数 n = 5000
$n = 10000 / 2;
echo ($n * ($n + 1)), "\n";
?>
<?php
# 10000 までの 奇数の和
# 項目数 n = 5000
$n = 10000 / 2;
echo pow($n, 2), "\n";
?>
<?php
# 1000 までの 自然数の２乗の和
$n = 1000;
echo ($n * ($n + 1) * (2 * $n + 1) / 6), "\n";
?>
<?php
# 100 までの 自然数の３乗の和
$n = 100;
echo (pow($n, 2) * pow(($n + 1), 2) / 4), "\n";
?>
<?php
# 初項 2, 公比 3, 項数 10 の等比数列の和
$n = 10;
$a = 2;
$r = 3;
echo (($a * (pow($r,  $n) - 1)) / ($r - 1)), "\n";
?>
<?php
$a = 5;  # 初項 5
$d = 3;  # 公差 3
$n = 10; # 項数 10
$p = 1;  # 積

foreach (range(1, $n) as $i)
{
    $m = $a + ($d * ($i - 1));
    $p *= $m;
}
echo $p;
?>
<?php
# 初項 5, 公差 3, 項数 10 の数列の積を表示する
echo product(5, 3, 10);

function product($m, $d, $n)
{
    if ($n == 0)
        return 1;
    else
        return $m * product($m + $d, $d, $n - 1);
}
?>
<?php
# 階乗を求める関数
function Fact($n)
{
    if ($n <= 1)
        return 1;
    else
        return $n * Fact($n - 1);
}

# 10の階乗
echo Fact(10), "\n";
echo 10 * 9 * 8 * 7 * 6 * 5 * 4 * 3 * 2 * 1, "\n";
?>
<?php
# 下降階乗冪
function FallingFact($x, $n)
{
    if ($n <= 1)
        return $x;
    else
        return $x * FallingFact($x - 1, $n - 1);
}

# 10 から 6 までの 総乗
echo FallingFact(10, 5), "\n";
echo 10 * 9 * 8 * 7 * 6, "\n";
?>
# 上昇階乗冪
function RisingFact($x, $n)
{
    if ($n <= 1)
        return $x;
    else
        return $x * RisingFact($x + 1, $n - 1);
}

# 10 から 14 までの 総乗
echo RisingFact(10, 5), "\n";
echo 10 * 11 * 12 * 13 * 14, "\n";
<?php
# 階乗
function Fact($n)
{
    if ($n <= 1)
        return 1;
    else
        return $n * Fact($n - 1);
}

# 下降階乗冪
function FallingFact($x, $n)
{
    if ($n <= 1)
        return $x;
    else
        return $x * FallingFact($x - 1, $n - 1);
}

# 順列 (異なる 10 個のものから 5 個取ってできる順列の総数)
$n = 10;
$r = 5;
echo Fact($n) / Fact($n - $r), "\n";
echo FallingFact($n, $r), "\n";
?>
<?php
# 重複順列 (異なる 10 個のものから重複を許して 5 個取ってできる順列の総数)
$n = 10;
$r = 5;
echo pow($n, $r), "\n";
?>
<?php
# 組合せ
function Comb($n, $r)
{
    if ($r == 0 || $r == $n)
        return 1;
    elseif ($r == 1)
        return $n;
    else
        return Comb($n - 1, $r - 1) + Comb($n - 1, $r);
}

# 組合せ (異なる 10 個のものから 5 個取ってできる組合せの総数)
$n = 10;
$r = 5;
echo Comb($n, $r), "\n";
?>
<?php
# 組合せ
function Comb($n, $r)
{
    if ($r == 0 || $r == $n)
        return 1;
    elseif ($r == 1)
        return $n;
    else
        return Comb($n - 1, $r - 1) + Comb($n - 1, $r);
}

# 重複組合せ (異なる 10 個のものから重複を許して 5 個とる組合せの総数)
$n = 10;
$r = 5;
echo Comb($n + $r - 1, $r), "\n";
?>
<?php
# 自作の正弦関数
function mySin($x, $n, $nega, $numerator, $denominator, $y)
{
    $m           = 2 * $n;
    $denominator = $denominator * ($m + 1) * $m;
    $numerator   = $numerator * $x * $x;
    $a           = $numerator / $denominator;
    # 十分な精度になったら処理を抜ける
    if ($a <= 0.00000000001)
        return $y;
    else
        return $y + mySin($x, ++$n, !$nega, $numerator, $denominator, $nega ? $a : -$a);
}

foreach (range(0, 24) as $i)
{
    $degree = $i * 15;
    if ($degree % 30 == 0 || $degree % 45 == 0)
    {
        $radian = deg2rad($degree);
        # 自作の正弦関数
        $d1     = mySin($radian, 1, false, $radian, 1.0, $radian);
        # 標準の正弦関数
        $d2     = sin($radian);
        # 標準関数との差異
        printf("%3d : %13.10f - %13.10f = %13.10f\n", $degree, $d1, $d2, $d1 - $d2);
    }
}
?>
<?php
# 自作の余弦関数
function myCos($x, $n, $nega, $numerator, $denominator, $y)
{
    $m           = 2 * $n;
    $denominator = $denominator * $m * ($m - 1);
    $numerator   = $numerator * $x * $x;
    $a           = $numerator / $denominator;
    # 十分な精度になったら処理を抜ける
    if ($a <= 0.00000000001)
        return $y;
    else
        return $y + myCos($x, ++$n, !$nega, $numerator, $denominator, $nega ? $a : -$a);
}

foreach (range(0, 24) as $i)
{
    $degree = $i * 15;
    if ($degree % 30 == 0 || $degree % 45 == 0)
    {
        $radian = deg2rad($degree);
        # 自作の余弦関数
        $d1     = myCos($radian, 1, false, 1.0, 1.0, 1.0);
        # 標準の余弦関数
        $d2     = cos($radian);
        # 標準関数との差異
        printf("%3d : %13.10f - %13.10f = %13.10f\n", $degree, $d1, $d2, $d1 - $d2);
    }
}
?>
<?php
# 自作の正接関数
function myTan($x, $x2, $n, $t)
{
    $t = $x2 / ($n - $t);
    $n -= 2;
    if ($n <= 1)
        return $x / (1 - $t);
    else
        return myTan($x, $x2, $n, $t);
}
foreach (range(0, 12) as $i)
{
    if (($i * 15) % 180 != 0)
    {
        $degree = $i * 15 - 90;
        $radian = deg2rad($degree);
        $x2     = $radian * $radian;
        # 自作の正接関数
        $d1     = myTan($radian, $x2, 15, 0.0); # 15:必要な精度が得られる十分大きな奇数
        # 標準の正接関数
        $d2     = tan($radian);
        # 標準関数との差異
        printf("%3d : %13.10f - %13.10f = %13.10f\n", $degree, $d1, $d2, $d1 - $d2);
    }
}
?>
<?php
# 自作の指数関数
function myExp($x, $n, $numerator, $denominator, $y)
{
    $denominator = $denominator * $n;
    $numerator   = $numerator   * $x;
    $a           = $numerator / $denominator;
    # 十分な精度になったら処理を抜ける
    if (abs($a) <= 0.00000000001)
        return $y;
    else
        return $y + myExp($x, ++$n, $numerator, $denominator, $a);
}

foreach (range(0, 20) as $i)
{
    $x  = ($i - 10) / 4.0;
    # 標準の指数関数
    $d1 = exp($x);
    # 自作の指数関数
    $d2 = myExp($x, 1, 1.0, 1.0, 1.0);
    # 標準関数との差異
    printf("%5.2f : %13.10f - %13.10f = %13.10f\n", $x, $d1, $d2, $d1 - $d2);
}
?>
<?php
# 自作の指数関数
function myExp($x, $x2, $n, $t)
{
    $t = $x2 / ($n + $t);
    $n -= 4;

    if ($n < 6)
        return 1 + ((2 * $x) / (2 - $x + $t));
    else
        return myExp($x, $x2, $n, $t);
}

foreach (range(0, 20) as $i)
{
    $x  = ($i - 10) / 4.0;
    # 標準の指数関数
    $d1 = exp($x);
    # 自作の指数関数
    $x2 = $x * $x;
    $d2 = myExp($x, $x2, 30, 0.0); # 30:必要な精度が得られるよう, 6から始めて4ずつ増加させる
    # 標準関数との差異
    printf("%5.2f : %13.10f - %13.10f = %13.10f\n", $x, $d1, $d2, $d1 - $d2);
}
?>
<?php
# 自作の対数関数
function myLog($x2, $numerator, $denominator, $y)
{
    $denominator = $denominator + 2;
    $numerator   = $numerator   * $x2 * $x2;
    $a           = $numerator / $denominator;
    # 十分な精度になったら処理を抜ける
    if (abs($a) <= 0.00000000001)
        return $y;
    else
        return $y + myLog($x2, $numerator, $denominator, $a);
}

foreach (range(1, 20) as $i)
{
    $x  = $i / 5.0;
    # 標準の対数関数
    $d1 = log($x);
    # 自作の対数関数
    $x2 = ($x - 1) / ($x + 1);
    $d2 = 2 * myLog($x2, $x2, 1.0, $x2);
    # 標準関数との差異
    printf("%5.2f : %13.10f - %13.10f = %13.10f\n", $x, $d1, $d2, $d1 - $d2);
}
?>
<?php
# 自作の対数関数
function myLog($x, $n, $t)
{
    $n2 = $n;
    $x2 = $x;
    if ($n > 3)
    {
        if ($n % 2 == 0)
            $n2 = 2;
        $x2 = $x * (int)($n / 2);
    }
    $t = $x2 / ($n2 + $t);

    if ($n <= 2)
        return $x / (1 + $t);
    else
        return myLog($x, --$n, $t);
}

foreach (range(1, 20) as $i)
{
    $x  = $i / 5.0;
    # 標準の対数関数
    $d1 = log($x);
    # 自作の対数関数
    $d2 = myLog($x - 1, 27, 0.0); # 27:必要な精度が得られる十分大きな奇数
    # 標準関数との差異
    printf("%5.2f : %13.10f - %13.10f = %13.10f\n", $x, $d1, $d2, $d1 - $d2);
}
?>
<?php
# 自作の双曲線正弦関数
function mySinh($x, $n, $numerator, $denominator, $y)
{
    $m           = 2 * $n;
    $denominator = $denominator * ($m + 1) * $m;
    $numerator   = $numerator * $x * $x;
    $a           = $numerator / $denominator;

    # 十分な精度になったら処理を抜ける
    if (abs($a) <= 0.00000000001)
        return $y;
    else
        return $y + mySinh($x, ++$n, $numerator, $denominator, $a);
}

foreach (range(0, 20) as $i)
{
    $x  = $i - 10;
    # 自作の双曲線正弦関数
    $d1 = mySinh($x, 1, $x, 1.0, $x);
    # 標準の双曲線正弦関数
    $d2 = sinh($x);
    # 標準関数との差異
    printf("%3d : %17.10f - %17.10f = %13.10f\n", $x, $d1, $d2, $d1 - $d2);
}
?>
<?php
# 自作の双曲線余弦関数
function myCosh($x, $n, $numerator, $denominator, $y)
{
    $m           = 2 * $n;
    $denominator = $denominator * $m * ($m - 1);
    $numerator   = $numerator * $x * $x;
    $a           = $numerator / $denominator;

    # 十分な精度になったら処理を抜ける
    if (abs($a) <= 0.00000000001)
        return $y;
    else
        return $y + myCosh($x, ++$n, $numerator, $denominator, $a);
}

foreach (range(0, 20) as $i)
{
    $x  = $i - 10;
    # 自作の双曲線余弦関数
    $d1 = myCosh($x, 1, 1.0, 1.0, 1.0);
    # 標準の双曲線余弦関数
    $d2 = cosh($x);
    # 標準関数との差異
    printf("%3d : %17.10f - %17.10f = %13.10f\n", $x, $d1, $d2, $d1 - $d2);
}
?>
<?php
$a = 0;
$b = 1;

# 台形則で積分
$n = 2;
foreach (range(1, 10) as $j)
{
    $h = ($b - $a) / $n;
    $s = 0;
    $x = $a;
    foreach (range(1, ($n - 1)) as $i)
    {
        $x += $h;
        $s += f($x);
    }
    $s = $h * ((f($a) + f($b)) / 2 + $s);
    $n *= 2;

    # 結果を π と比較
    printf("%2d : %13.10f, %13.10f\n", $j, $s, $s - M_PI);
}

function f($x)
{
    return 4 / (1 + $x * $x);
}
?>
<?php
$a = 0;
$b = 1;

# 中点則で積分
$n = 2;
foreach (range(1, 10) as $j)
{
    $h = ($b - $a) / $n;
    $s = 0;
    $x = $a + ($h / 2);
    foreach (range(1, $n) as $i)
    {
        $s += f($x);
        $x += $h;
    }
    $s *= $h;
    $n *= 2;

    # 結果を π と比較
    printf("%2d : %13.10f, %13.10f\n", $j, $s, $s - M_PI);
}

function f($x)
{
    return 4 / (1 + $x * $x);
}
?>
<?php
$a = 0;
$b = 1;

# Simpson則で積分
$n = 2;
foreach (range(1, 5) as $j)
{
    $h  = ($b - $a) / $n;
    $s2 = 0;
    $s4 = 0;
    $x  = $a + $h;
    foreach (range(1, ($n / 2)) as $i)
    {
        $s4 += f($x);
        $x  += $h;
        $s2 += f($x);
        $x  += $h;
    }
    $s2 =  ($s2 - f($b)) * 2 + f($a) + f($b);
    $s4 *= 4;
    $s  =  ($s2 + $s4) * $h / 3;
    $n  *= 2;

    # 結果を π と比較
    printf("%2d : %13.10f, %13.10f\n", $j, $s, $s - M_PI);
}

function f($x)
{
    return 4 / (1 + $x * $x);
}
?>
<?php
$a = 0;
$b = 1;

$t = array();

# 台形則で積分
$n = 2;
foreach (range(1, 6) as $i)
{
    $h = ($b - $a) / $n;
    $s = 0;
    $x = $a;
    foreach (range(1, ($n - 1)) as $j)
    {
        $x += $h;
        $s += f($x);
    }
    # 結果を保存
    $t[$i][1] = $h * ((f($a) + f($b)) / 2 + $s);
    $n *= 2;
}

# Richardsonの補外法
$n = 4;
foreach (range(2, 6) as $j)
{
    foreach (range($j, 6) as $i)
    {
        $t[$i][$j] = $t[$i][$j - 1] + ($t[$i][$j - 1] - $t[$i - 1][$j - 1]) / ($n - 1);
        if ($i == $j)
        {
            # 結果を π と比較
            printf("%2d : %13.10f, %13.10f\n", $j, $t[$i][$j], $t[$i][$j] - M_PI);
        }
    }
    $n *= 4;
}

function f($x)
{
    return 4 / (1 + $x * $x);
}
?>
<?php
# データ点の数 - 1
define("N", 6);

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

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

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

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

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

# Lagrange (ラグランジュ) 補間
function lagrange($d, $x, $y)
{
    $sum = 0;
    foreach (range(0, N) as $i)
    {
        $prod = $y[$i];
        foreach (range(0, N) as $j)
        {
            if ($j != $i)
            {
                $prod *= ($d - $x[$j]) / ($x[$i] - $x[$j]);
            }
        }
        $sum += $prod;
    }
    return $sum;
}
?>
<?php
# データ点の数 - 1
define("N", 6);

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

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

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

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

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

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

    foreach (range(1, N) as $j)
    {
        foreach (range(0, N - $j) as $i)
        {
            $w[$j][$i] = $w[$j-1][$i+1] + ($w[$j-1][$i+1] - $w[$j-1][$i]) * ($d - $x[$i+$j]) / ($x[$i+$j] - $x[$i]);
        }
    }
    return $w[N][0];
}
?>
<?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;
}
?>
<?php
# データ点の数 - 1
define("N",    6);
define("Nx2", 13);

$x  = array();
$y  = array();
$yd = 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);
    $yd[$i] = fd($d1);
}

# 差分商の表を作る
$z = array();
$d = array();
foreach (range(0, Nx2) as $i)
{
    $j        = (int)($i / 2);
    $z[$i]    = $x[$j];
    $d[0][$i] = $y[$j];
}
foreach (range(1, Nx2) as $i)
{
    foreach (range(0, Nx2 - $i) as $j)
    {
        if ($i == 1 && $j % 2 == 0)
            $d[$i][$j] = $yd[(int)($j / 2)];
        else
            $d[$i][$j] = ($d[$i-1][$j+1] - $d[$i-1][$j]) / ($z[$j+$i] - $z[$j]);
    }
}
# ｎ階差分商
$a = array();
foreach (range(0, Nx2) 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 = hermite($d1, $z, $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);
}
# 導関数
function fd($x)
{
    return 1 - pow($x, 2) / 2 + pow($x, 4) / (4 * 3 * 2);
}

# Hermite (エルミート) 補間
function hermite($d, $z, $a)
{
    $sum = $a[0];
    foreach (range(1, Nx2) as $i)
    {
        $prod = $a[$i];
        foreach (range(0, $i - 1) as $j)
        {
            $prod *= ($d - $z[$j]);
        }
        $sum += $prod;
    }
    return $sum;
}
?>
<?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;
}
?>
<?php

# 重力加速度
$g = -9.8;
# 空気抵抗係数
$k = -0.01;
# 時間間隔(秒)
$h = 0.01;

# 角度
$degree = 45;
$radian = $degree * M_PI / 180.0;
# 初速 250 km/h -> 秒速に変換
$v = (int)(250 * 1000 / 3600);
# 水平方向の速度
$vx = array();
$vx[0] = $v * cos($radian);
# 鉛直方向の速度
$vy = array();
$vy[0] = $v * sin($radian);
# 経過秒数
$t = 0.0;
# 位置
$x = 0.0;
$y = 0.0;

# Euler法
for ($i = 1; $y >= 0.0; $i++)
{
    # 経過秒数
    $t = $i * $h;

    # 位置
    $x += $h * $vx[0];
    $y += $h * $vy[0];

    printf("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f\n", $t, $vx[0], $vy[0], $x, $y);

    # 速度
    $vx[1] = $vx[0] + $h * fx($vx[0], $vy[0]);
    $vy[1] = $vy[0] + $h * fy($vx[0], $vy[0]);
    $vx[0] = $vx[1];
    $vy[0] = $vy[1];
}

# 空気抵抗による水平方向の減速分
function fx($vx, $vy)
{
    global $k;
    return $k * sqrt($vx * $vx + $vy * $vy) * $vx;
}
# 重力と空気抵抗による鉛直方向の減速分
function fy($vx, $vy)
{
    global $g, $k;
    return $g + ($k * sqrt($vx * $vx + $vy * $vy) * $vy);
}
?>
<?php

# 重力加速度
$g = -9.8;
# 空気抵抗係数
$k = -0.01;
# 時間間隔(秒)
$h = 0.01;

# 角度
$degree = 45;
$radian = $degree * M_PI / 180.0;
# 初速 250 km/h -> 秒速に変換
$v = (int)(250 * 1000 / 3600);
# 水平方向の速度
$vx = array();
$vx[0] = $v * cos($radian);
# 鉛直方向の速度
$vy = array();
$vy[0] = $v * sin($radian);
# 経過秒数
$t = 0.0;
# 位置
$x = array();
$x[0] = 0.0;
$y = array();
$y[0] = 0.0;

# Heun法
for ($i = 1; $y[0] >= 0.0; $i++)
{
    # 経過秒数
    $t = $i * $h;

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

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

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

    printf("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f\n", $t, $vx[0], $vy[0], $x[0], $y[0]);
}

# 空気抵抗による水平方向の減速分
function fx($vx, $vy)
{
    global $k;
    return $k * sqrt($vx * $vx + $vy * $vy) * $vx;
}
# 重力と空気抵抗による鉛直方向の減速分
function fy($vx, $vy)
{
    global $g, $k;
    return $g + ($k * sqrt($vx * $vx + $vy * $vy) * $vy);
}
?>
<?php

# 重力加速度
$g = -9.8;
# 空気抵抗係数
$k = -0.01;
# 時間間隔(秒)
$h = 0.01;

# 角度
$degree = 45;
$radian = $degree * M_PI / 180.0;
# 初速 250 km/h -> 秒速に変換
$v = (int)(250 * 1000 / 3600);
# 水平方向の速度
$vx = array();
$vx[0] = $v * cos($radian);
# 鉛直方向の速度
$vy = array();
$vy[0] = $v * sin($radian);
# 経過秒数
$t = 0.0;
# 位置
$x = array();
$x[0] = 0.0;
$y = array();
$y[0] = 0.0;

# 中点法
for ($i = 1; $y[0] >= 0.0; $i++)
{
    # 経過秒数
    $t = $i * $h;

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

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

    printf("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%9.5f\n", $t, $vx[0], $vy[0], $x[0], $y[0]);
}

# 空気抵抗による水平方向の減速分
function fx($vx, $vy)
{
    global $k;
    return $k * sqrt($vx * $vx + $vy * $vy) * $vx;
}
# 重力と空気抵抗による鉛直方向の減速分
function fy($vx, $vy)
{
    global $g, $k;
    return $g + ($k * sqrt($vx * $vx + $vy * $vy) * $vy);
}
?>
<?php

# 重力加速度
$g = -9.8;
# 空気抵抗係数
$k = -0.01;
# 時間間隔(秒)
$h = 0.01;

# 角度
$degree = 45;
$radian = $degree * M_PI / 180.0;
# 初速 250 km/h -> 秒速に変換
$v = (int)(250 * 1000 / 3600);
# 水平方向の速度
$vx = array();
$vx[0] = $v * cos($radian);
# 鉛直方向の速度
$vy = array();
$vy[0] = $v * sin($radian);
# 経過秒数
$t = 0.0;
# 位置
$x = array();
$x[0] = 0.0;
$y = array();
$y[0] = 0.0;

# Runge-Kutta法
for ($i = 1; $y[0] >= 0.0; $i++)
{
    # 経過秒数
    $t = $i * $h;

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

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

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

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

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

    printf("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%9.5f\n", $t, $vx[0], $vy[0], $x[0], $y[0]);
}

# 空気抵抗による水平方向の減速分
function fx($vx, $vy)
{
    global $k;
    return $k * sqrt($vx * $vx + $vy * $vy) * $vx;
}
# 重力と空気抵抗による鉛直方向の減速分
function fy($vx, $vy)
{
    global $g, $k;
    return $g + ($k * sqrt($vx * $vx + $vy * $vy) * $vy);
}
?>
<?php

# 重力加速度
$g = -9.8;
# 空気抵抗係数
$k = -0.01;
# 時間間隔(秒)
$h = 0.01;

# 角度
$degree = 45;
$radian = $degree * M_PI / 180.0;
# 初速 250 km/h -> 秒速に変換
$v = (int)(250 * 1000 / 3600);
# 水平方向の速度
$vx = array();
$vx[0] = $v * cos($radian);
# 鉛直方向の速度
$vy = array();
$vy[0] = $v * sin($radian);
# 経過秒数
$t = 0.0;
# 位置
$x = array();
$x[0] = 0.0;
$y = array();
$y[0] = 0.0;

# Runge-Kutta-Gill法
for ($i = 1; $y[0] >= 0.0; $i++)
{
    # 経過秒数
    $t = $i * $h;

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

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

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

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

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

    printf("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%9.5f\n", $t, $vx[0], $vy[0], $x[0], $y[0]);
}

# 空気抵抗による水平方向の減速分
function fx($vx, $vy)
{
    global $k;
    return $k * sqrt($vx * $vx + $vy * $vy) * $vx;
}
# 重力と空気抵抗による鉛直方向の減速分
function fy($vx, $vy)
{
    global $g, $k;
    return $g + ($k * sqrt($vx * $vx + $vy * $vy) * $vy);
}
?>
<?php
$a = 1;
$b = 2;
printf("%12.10f\n", bisection($a, $b));

function bisection($a, $b)
{
    while (TRUE)
    {
        # 区間 (a, b) の中点 c = (a + b) / 2
        $c = ($a + $b) / 2;
        printf("%12.10f\t%13.10f\n", $c, $c - sqrt(2));

        $fc = f($c);
        if (abs($fc) < 0.0000000001) break;

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

function f($x)
{
    return $x * $x - 2;
}
?>
<?php
$a = 1;
$b = 2;
printf("%12.10f\n", falseposition($a, $b));

function falseposition($a, $b)
{
    while (TRUE)
    {
        # 点 (a,f(a)) と 点 (b,f(b)) を結ぶ直線と x軸の交点
        $c = ($a * f($b) - $b * f($a)) / (f($b) - f($a));
        printf("%12.10f\t%13.10f\n", $c, $c - sqrt(2));

        $fc = f($c);
        if (abs($fc) < 0.0000000001) break;

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

function f($x)
{
    return $x * $x - 2;
}
?>
<?php
$x = 1;
printf("%12.10f\n", iterative($x));

function iterative($x0)
{
    while (TRUE)
    {
        $x1 = g($x0);
        printf("%12.10f\t%13.10f\n", $x1, $x1 - sqrt(2));

        if (abs($x1 - $x0) < 0.0000000001) break;

        $x0 = $x1;
    }
    return $x1;
}

function g($x)
{
    return ($x / 2) + (1 / $x);
}
?>
<?php
$x = 2;
printf("%12.10f\n", newton($x));

function newton($x0)
{
    while (TRUE)
    {
        $x1 = $x0 - (f0($x0) / f1($x0));
        printf("%12.10f\t%13.10f\n", $x1, $x1 - sqrt(2));

        if (abs($x1 - $x0) < 0.0000000001) break;

        $x0 = $x1;
    }
    return $x1;
}

function f0($x)
{
    return $x * $x - 2;
}

function f1($x)
{
    return 2 * $x;
}
?>
<?php
$x = 2;
printf("%12.10f\n", bailey($x));

function bailey($x0)
{
    while (TRUE)
    {
        $x1 = $x0 - (f0($x0) / (f1($x0) - (f0($x0) * f2($x0) / (2 * f1($x0)))));
        printf("%12.10f\t%13.10f\n", $x1, $x1 - sqrt(2));

        if (abs($x1 - $x0) < 0.0000000001) break;

        $x0 = $x1;
    }
    return $x1;
}

function f0($x)
{
    return $x * $x - 2;
}

function f1($x)
{
    return 2 * $x;
}

function f2($x)
{
    return 2;
}
?>
<?php
$x0 = 1;
$x1 = 2;
printf("%12.10f\n", secant($x0, $x1));

function secant($x0, $x1)
{
    while (TRUE)
    {
        $x2 = $x1 - f($x1) * ($x1 - $x0) / (f($x1) - f($x0));
        printf("%12.10f\t%13.10f\n", $x2, $x2 - sqrt(2));

        if (abs($x2 - $x1) < 0.0000000001) break;

        $x0 = $x1;
        $x1 = $x2;
    }
    return $x2;
}

function f($x)
{
    return $x * $x - 2;
}
?>
<?php
define("N", 3);

$a = [[9,2,1,1],[2,8,-2,1],[-1,-2,7,-2],[1,-1,-2,6]];
$b = [20,16,8,17];
$c = [0,0,0,0];

# ヤコビの反復法
jacobi($a, $b, $c);

print "解\n";
foreach (range(0, N) as $i)
    printf("%14.10f\t", $c[$i]);
print "\n";

# ヤコビの反復法
function jacobi($a, $b, &$x0)
{
    while (true)
    {
        $x1 = array();
        $finish = true;
        foreach (range(0, N) as $i)
        {
            $x1[$i] = 0;
            foreach (range(0, N) as $j)
            {
                if ($j != $i)
                    $x1[$i] += $a[$i][$j] * $x0[$j];
            }
            $x1[$i] = ($b[$i] - $x1[$i]) / $a[$i][$i];
            if (abs($x1[$i] - $x0[$i]) > 0.0000000001) $finish = false;
        }

        foreach (range(0, N) as $i)
        {
            $x0[$i] = $x1[$i];
            printf("%14.10f\t", $x0[$i]);
        }
        print "\n";
        if ($finish) return;
    }
}
?>
<?php
define("N", 3);

$a = [[9,2,1,1],[2,8,-2,1],[-1,-2,7,-2],[1,-1,-2,6]];
$b = [20,16,8,17];
$c = [0,0,0,0];

# ガウス・ザイデル法
gauss($a, $b, $c);

print "解\n";
foreach (range(0, N) as $i)
    printf("%14.10f\t", $c[$i]);
print "\n";

# ガウス・ザイデル法
function gauss($a, $b, &$x0)
{
    while (true)
    {
        $finish = true;
        foreach (range(0, N) as $i)
        {
            $x1 = 0;
            foreach (range(0, N) as $j)
            {
                if ($j != $i)
                    $x1 += $a[$i][$j] * $x0[$j];
            }
            $x1 = ($b[$i] - $x1) / $a[$i][$i];
            if (abs($x1 - $x0[$i]) > 0.0000000001) $finish = false;
            $x0[$i] = $x1;
        }

        foreach (range(0, N) as $i)
        {
            printf("%14.10f\t", $x0[$i]);
        }
        print "\n";
        if ($finish) return;
    }
}
?>
<?php
define("N", 3);

$a = [[-1,-2,7,-2],[1,-1,-2,6],[9,2,1,1],[2,8,-2,1]];
$b = [8,17,20,16];

# ピボット選択
pivoting($a,$b);
print "ピボット選択後\n";
disp_progress($a,$b);

# 前進消去
forward_elimination($a,$b);
print "前進消去後\n";
disp_progress($a,$b);

# 後退代入
backward_substitution($a,$b);
print "解\n";
foreach (range(0, N) as $i)
    printf("%14.10f\t", $b[$i]);
print "\n";

# 前進消去
function forward_elimination(&$a, &$b)
{
    foreach (range(0, N - 1) as $pivot)
    {
        foreach (range($pivot + 1, N) as $row)
        {
            $s = $a[$row][$pivot] / $a[$pivot][$pivot];
            foreach (range($pivot, N) as $col)
                    $a[$row][$col] -= $a[$pivot][$col]    * $s;
            $b[$row]               -= $b[$pivot]          * $s;
        }
    }
}
# 後退代入
function backward_substitution(&$a, &$b)
{
     for ($row = N; $row >= 0; $row--)
    {
        for ($col = N; $col > $row; $col--)
            $b[$row] -= $a[$row][$col] * $b[$col];
        $b[$row] /= $a[$row][$row];
    }
}
# ピボット選択
function pivoting(&$a, &$b)
{
    foreach (range(0, N) as $pivot)
    {
        # 各列で 一番値が大きい行を 探す
        $max_row =   $pivot;
        $max_val =   0;
        foreach (range($pivot, N) as $row)
        {
            if (abs($a[$row][$pivot]) > $max_val)
            {
                # 一番値が大きい行
                $max_val =   abs($a[$row][$pivot]);
                $max_row =   $row;
            }
        }

        # 一番値が大きい行と入れ替え
        if ($max_row != $pivot)
        {
            foreach (range(0, N) as $col)
            {
                $tmp               =   $a[$max_row][$col];
                $a[$max_row][$col] =   $a[$pivot][$col];
                $a[$pivot][$col]   =   $tmp;
            }
            $tmp         =   $b[$max_row];
            $b[$max_row] =   $b[$pivot];
            $b[$pivot]   =   $tmp;
        }
    }
}
# 状態を確認
function disp_progress($a, $b)
{
  foreach (range(0, N) as $i)
 {
    foreach (range(0, N) as $j)
            printf("%14.10f\t", $a[$i][$j]);
    printf("%14.10f\n", $b[$i]);
  }
  print "\n";
}
?>
<?php
define("N", 3);

$a = [[-1,-2,7,-2],[1,-1,-2,6],[9,2,1,1],[2,8,-2,1]];
$b = [8,17,20,16];

# ピボット選択
pivoting($a,$b);
print "ピボット選択後\n";
disp_progress($a,$b);

# 前進消去
forward_elimination($a,$b);
print "前進消去後\n";
disp_progress($a,$b);

# 後退代入
backward_substitution($a,$b);
print "解\n";
foreach (range(0, N) as $i)
    printf("%14.10f\t", $b[$i]);
print "\n";

# 前進消去
function forward_elimination(&$a, &$b)
{
    # 対角線上の係数を 1 にする
    foreach (range(0, N) as $pivot)
    {
        $s = $a[$pivot][$pivot];
        foreach (range(0, N) as $col)
            $a[$pivot][$col] /= $s;
        $b[$pivot]           /= $s;
    }
    print "対角線上の係数を 1 にする\n";
    disp_progress($a,$b);

    # 対角行列にする
    foreach (range(0, N) as $pivot)
    {
        foreach (range(0, N) as $row)
        {
            if ($row == $pivot) continue;

            $s = $a[$row][$pivot] / $a[$pivot][$pivot];
        foreach (range($pivot, N) as $col)
                $a[$row][$col] -= $a[$pivot][$col] * $s;
            $b[$row] -= $b[$pivot] * $s;
        }
    }
}
# 後退代入
function backward_substitution(&$a, &$b)
{
    foreach (range(0, N) as $pivot)
        $b[$pivot]  /= $a[$pivot][$pivot];
}
# ピボット選択
function pivoting(&$a, &$b)
{
    foreach (range(0, N) as $pivot)
    {
        # 各列で 一番値が大きい行を 探す
        $max_row =   $pivot;
        $max_val =   0;
        foreach (range($pivot, N) as $row)
        {
            if (abs($a[$row][$pivot]) > $max_val)
            {
                # 一番値が大きい行
                $max_val =   abs($a[$row][$pivot]);
                $max_row =   $row;
            }
        }

        # 一番値が大きい行と入れ替え
        if ($max_row != $pivot)
        {
            foreach (range(0, N) as $col)
            {
                $tmp               =   $a[$max_row][$col];
                $a[$max_row][$col] =   $a[$pivot][$col];
                $a[$pivot][$col]   =   $tmp;
            }
            $tmp         =   $b[$max_row];
            $b[$max_row] =   $b[$pivot];
            $b[$pivot]   =   $tmp;
        }
    }
}
# 状態を確認
function disp_progress($a, $b)
{
  foreach (range(0, N) as $i)
 {
    foreach (range(0, N) as $j)
    {
            printf("%14.10f\t", $a[$i][$j]);
    }
    printf("%14.10f\n", $b[$i]);
  }
  print "\n";
}
?>
<?php
define("N", 3);

$a = [[-1,-2,7,-2],[1,-1,-2,6],[9,2,1,1],[2,8,-2,1]];
$b = [8,17,20,16];

# ピボット選択
pivoting($a,$b);
print "ピボット選択後\n";
disp_progress($a,$b);

# LU分解
$x = [0,0,0,0];
decomp($a,$b,$x);
print "X\n";
foreach (range(0, N) as $i)
    printf("%14.10f\t", $x[$i]);
print "\n";

# LU分解
function decomp(&$a, &$b, &$x)
{
    # 前進消去
    forward_elimination($a,$b);
    # 下三角行列を確認
    print "L\n";
    foreach (range(0, N) as $i)
    {
        foreach (range(0, N) as $j)
        {
            if ($i > $j)
            {
                printf("%14.10f\t", $a[$i][$j]);
            }
            elseif ($i == $j)
            {
                printf("%14.10f\t", 1.0);
            }
            else
            {
                printf("%14.10f\t", 0.0);
            }
        }
        print "\n";
    }
    print "\n";
    # 上三角行列を確認
    print "U\n";
    foreach (range(0, N) as $i)
    {
        foreach (range(0, N) as $j)
        {
            if ($i <= $j)
            {
                printf("%14.10f\t", $a[$i][$j]);
            }
            else
            {
                printf("%14.10f\t", 0.0);
            }
        }
        print "\n";
    }
    print "\n";

    # Ly=b から y を求める (前進代入)
    $y = [0,0,0,0];
    forward_substitution($a,$b,$y);
    # y の 値 を確認
    print "Y\n";
    foreach (range(0, N) as $i)
    {
        printf("%14.10f\t", $y[$i]);
    }
    print "\n";

    # Ux=y から x を求める (後退代入)
    backward_substitution($a,$y,$x);
}

# 前進消去
function forward_elimination(&$a, $b)
{
    foreach (range(0, N - 1) as $pivot)
    {
        foreach (range($pivot + 1, N) as $row)
        {
            $s = $a[$row][$pivot] / $a[$pivot][$pivot];
            foreach (range($pivot, N) as $col)
                    $a[$row][$col] -= $a[$pivot][$col]    * $s; # これが 上三角行列
            $a[$row][$pivot] = $s;                              # これが 下三角行列
            # $b[$row]               -= $b[$pivot]          * $s;# この値は変更しない
        }
    }
}
# Ly=b から y を求める (前進代入)
function forward_substitution($a, $b, &$y)
{
    foreach (range(0, N) as $row)
    {
        foreach (range(0, $row) as $col)
        {
            $b[$row] -= $a[$row][$col] * $y[$col];
        }
        $y[$row] = $b[$row];
    }
}
# Ux=y から x を求める (後退代入)
function backward_substitution($a, $y, &$x)
{
    for ($row = N; $row >= 0; $row--)
    {
        for ($col = N; $col > $row; $col--)
        {
            $y[$row] -= $a[$row][$col] * $x[$col];
        }
        $x[$row] = $y[$row] / $a[$row][$row];
    }
}
# ピボット選択
function pivoting(&$a, &$b)
{
    foreach (range(0, N) as $pivot)
    {
        # 各列で 一番値が大きい行を 探す
        $max_row =   $pivot;
        $max_val =   0;
        foreach (range($pivot, N) as $row)
        {
            if (abs($a[$row][$pivot]) > $max_val)
            {
                # 一番値が大きい行
                $max_val =   abs($a[$row][$pivot]);
                $max_row =   $row;
            }
        }

        # 一番値が大きい行と入れ替え
        if ($max_row != $pivot)
        {
            foreach (range(0, N) as $col)
            {
                $tmp               =   $a[$max_row][$col];
                $a[$max_row][$col] =   $a[$pivot][$col];
                $a[$pivot][$col]   =   $tmp;
            }
            $tmp         =   $b[$max_row];
            $b[$max_row] =   $b[$pivot];
            $b[$pivot]   =   $tmp;
        }
    }
}
# 状態を確認
function disp_progress($a, $b)
{
  foreach (range(0, N) as $i)
 {
    foreach (range(0, N) as $j)
            printf("%14.10f\t", $a[$i][$j]);
    printf("%14.10f\n", $b[$i]);
  }
  print "\n";
}
?>
<?php
define("N", 3);

$a = [[5,2,3,4],[2,10,6,7],[3,6,15,9],[4,7,9,20]];
$b = [34,68,96,125];
print "A\n";
disp_progress($a);

# LL^T分解
$x = [0,0,0,0];
decomp($a,$b,$x);
print "X\n";
foreach (range(0, N) as $i)
  printf("%14.10f\t", $x[$i]);
print "\n";

# LL^T分解
function decomp(&$a, &$b, &$x)
{
  # 前進消去
    forward_elimination($a,$b);
    print "L と L^T\n";
    foreach (range(0, N) as $i)
    {
        foreach (range(0, N) as $j)
        {
            if ($j <= $i)
            {
                printf("%14.10f\t", $a[$i][$j]);
            }
            else
            {
                printf("%14.10f\t", $a[$j][$i]);
            }
        }
        print "\n";
    }
    print "\n";

    # Ly=b から y を求める (前進代入)
    $y = [0,0,0,0];
    forward_substitution($a,$b,$y);
    # y の 値 を確認
    print "Y\n";
    foreach (range(0, N) as $i)
    {
        printf("%14.10f\t", $y[$i]);
    }
    print "\n";

    # Ux=y から x を求める (後退代入)
    backward_substitution($a,$y,$x);
}

# 前進消去
function forward_elimination(&$a, $b)
{
    foreach (range(0, N) as $pivot)
    {
        $s = 0;
        for ($col = 0; $col < $pivot; $col++)
            $s += $a[$pivot][$col] * $a[$pivot][$col];
        # ここで根号の中が負の値になると計算できない！
        $a[$pivot][$pivot] = sqrt($a[$pivot][$pivot] - $s);

        for ($row = $pivot + 1; $row <= N; $row++)
        {
            $s = 0;
            for ($col = 0; $col < $pivot; $col++)
                $s += $a[$row][$col] * $a[$pivot][$col];
            $a[$row][$pivot] = ($a[$row][$pivot] - $s) / $a[$pivot][$pivot];
        }
    }
}
# Ly=b から y を求める (前進代入)
function forward_substitution($a, $b, &$y)
{
    foreach (range(0, N) as $row)
    {
        foreach (range(0, $row - 1) as $col)
        {
            $b[$row] -= $a[$row][$col] * $y[$col];
        }
        $y[$row] = $b[$row] / $a[$row][$row];
    }
}
# Ux=y から x を求める (後退代入)
function backward_substitution($a, $y, &$x)
{
    for ($row = N; $row >= 0; $row--)
    {
        for ($col = N; $col > $row; $col--)
        {
            $y[$row] -= $a[$col][$row] * $x[$col];
        }
        $x[$row] = $y[$row] / $a[$row][$row];
    }
}
# 状態を確認
function disp_progress($a)
{
  foreach (range(0, N) as $i)
 {
    foreach (range(0, N) as $j)
        printf("%14.10f\t", $a[$i][$j]);
    print "\n";
  }
  print "\n";
}
?>
<?php
define("N", 3);

$a = [[5,2,3,4],[2,10,6,7],[3,6,15,9],[4,7,9,20]];
$b = [34,68,96,125];

print "A\n";
disp_progress($a);

# LDL^T分解
$x = [0,0,0,0];
decomp($a,$b,$x);
print "X\n";
foreach (range(0, N) as $i)
  printf("%14.10f\t", $x[$i]);
print "\n";

# LDL^T分解
function decomp(&$a, &$b, &$x)
{
  # 前進消去
  forward_elimination($a,$b);
    print "L と D\n";
    foreach (range(0, N) as $i)
    {
        foreach (range(0, N) as $j)
        {
            if ($j <= $i)
            {
                printf("%14.10f\t", $a[$i][$j]);
            }
            else
            {
                printf("%14.10f\t", $a[$j][$i]);
            }
        }
        print "\n";
    }
    print "\n";

    # Ly=b から y を求める (前進代入)
    $y = [0,0,0,0];
    forward_substitution($a,$b,$y);
    # y の 値 を確認
    print "Y\n";
    foreach (range(0, N) as $i)
    {
        printf("%14.10f\t", $y[$i]);
    }
    print "\n";

    # DL^Tx=y から x を求める (後退代入)
    backward_substitution($a,$y,$x);
}

# 前進消去
function forward_elimination(&$a, $b)
{
    foreach (range(0, N) as $pivot)
    {
        # pivot < k の場合
        $s = 0;
        for ($col = 0; $col < $pivot; $col++)
        {
            $s = $a[$pivot][$col];
            for ($k = 0; $k < $col; $k++)
            {
                $s -= $a[$pivot][$k] * $a[$col][$k] * $a[$k][$k];
            }
            $a[$pivot][$col] = $s / $a[$col][$col];
        }

        # pivot == k の場合
        $s = $a[$pivot][$pivot];
        for ($k = 0; $k < $pivot; $k++)
        {
            $s -= $a[$pivot][$k] * $a[$pivot][$k] * $a[$k][$k];
        }
        $a[$pivot][$pivot] = $s;
    }
}
# Ly=b から y を求める (前進代入)
function forward_substitution($a, $b, &$y)
{
    foreach (range(0, N) as $row)
    {
        foreach (range(0, $row - 1) as $col)
        {
            $b[$row] -= $a[$row][$col] * $y[$col];
        }
        $y[$row] = $b[$row];
    }
}
# Ux=y から x を求める (後退代入)
function backward_substitution($a, $y, &$x)
{
    for ($row = N; $row >= 0; $row--)
    {
        for ($col = N; $col > $row; $col--)
        {
            $y[$row] -= $a[$col][$row] * $a[$row][$row] * $x[$col];
        }
        $x[$row] = $y[$row] / $a[$row][$row];
    }
}
# 状態を確認
function disp_progress($a)
{
  foreach (range(0, N) as $i)
 {
    foreach (range(0, N) as $j)
        printf("%14.10f\t", $a[$i][$j]);
    print "\n";
  }
  print "\n";
}
?>