home > さまざまな言語で数値計算 > 数値積分 >

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

Only Do What Only You Can Do

シンプソン則

関数 $ f(x) $ の $ a, b $ 間の $n$ 個の微小区間を, 二次曲線で近似し面積を合計する.

3点 $ (x_0, f(x_0)), (x_1, f(x_1)), (x_2, f(x_2)) $ を通る二次曲線の方程式 $ g(x) $ は,

で表される(ラグランジュ補間).

$ x_1 - x_0 $ を $h$ として, 右辺第1項を積分すると

同様に右辺第2項は

同じく第3項は

よって, 関数 $ g(x) $ の $ x_0, x_2 $ 区間の定積分は

と表すことができる.

まとめると, 関数 $ f(x) $ の $ a, b $ 間の $n$ 個の微小区間を, 二次曲線で近似し面積を合計する式は以下のとおり.

例題として, $ 4/(1+x^2) $ を $0$ から $1$ まで積分し $\pi$ を求める.

VBScript

Option Explicit

Const PI = 3.14159265359
Const a = 0
Const b = 1

'Simpson則で積分
Dim n: n = 2
Dim i, j
For j = 1 To 5
    Dim h:  h  = (b - a) / n
    Dim s2: s2 = 0
    Dim s4: s4 = 0
    Dim x:  x  = a + h
    For i = 1 To n \ 2
        s4 = s4 + f(x)
        x  = x  + h
        s2 = s2 + f(x)
        x  = x  + h
    Next
    s2 = (s2 - f(b)) * 2 + f(a) + f(b)
    s4 = s4 * 4
    Dim s: s = (s2 + s4) * h / 3
    n = n * 2

    '結果を π と比較
    WScript.StdOut.Write Right(Space(2)  & j,                                   2) & " : "
    WScript.StdOut.Write Right(Space(13) & FormatNumber(s,      10, -1, 0, 0), 13) & ", "
    WScript.StdOut.Write Right(Space(13) & FormatNumber(s - PI, 10, -1, 0, 0), 13) & vbNewLine
Next

Private Function f(x)
    f = 4 / (1 + x * x)
End Function
Z:\>cscript //nologo Z:\0603.vbs
 1 :  3.1333333333, -0.0082593203
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415925025, -0.0000001511
 4 :  3.1415926512, -0.0000000024
 5 :  3.1415926536, -0.0000000000

JScript

var a = 0
var b = 1

// Simpson則で積分
var n = 2
for (var j = 1; j <= 5; j++)
{
    var h  = (b - a) / n
    var s2 = 0
    var s4 = 0
    var x  = a + h
    for (var i = 1; i <= n / 2; i++)
    {
        s4 += f(x)
        x  += h
        s2 += f(x)
        x  += h
    }
    s2 = (s2 - f(b)) * 2 + f(a) + f(b)
    s4 *= 4
    var s = (s2 + s4) * h / 3
    n *= 2

    // 結果を π と比較
    WScript.StdOut.Write(("  "            + j                        ).slice( -2) + " : ");
    WScript.StdOut.Write(("             " + s.toFixed(10)            ).slice(-13) + ", ");
    WScript.StdOut.Write(("             " + (s - Math.PI).toFixed(10)).slice(-13) + "\n" );
}

function f(x)
{
    return 4 / (1 + x * x)
}
Z:\>cscript //nologo Z:\0603.js
 1 :  3.1333333333, -0.0082593203
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415925025, -0.0000001511
 4 :  3.1415926512, -0.0000000024
 5 :  3.1415926536, -0.0000000000

PowerShell

function f($x)
{
    4 / (1 + $x * $x)
}

$a = 0
$b = 1

# Simpson則で積分
$n = 2;
foreach ($j in 1..5)
{
    $h  = ($b - $a) / $n
    $s2 = 0
    $s4 = 0
    $x  = $a + $h
    foreach ($i in 1..($n / 2))
    {
        $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

    # 結果を π と比較
    Write-Host ([String]::Format("{0,2:D} : {1,13:F10}, {2,13:F10}", $j, $s, $s - [Math]::PI))
}
Z:\>powershell -file Z:\0603.ps1
 1 :  3.1333333333, -0.0082593203
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415925025, -0.0000001511
 4 :  3.1415926512, -0.0000000024
 5 :  3.1415926536,  0.0000000000

Perl

use Math::Trig 'pi';

my $a = 0;
my $b = 1;

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

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

sub f
{
    my ($x) = @_;
    4 / (1 + $x * $x);
}
Z:\>perl Z:\0603.pl
 1 :  3.1333333333, -0.0082593203
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415925025, -0.0000001511
 4 :  3.1415926512, -0.0000000024
 5 :  3.1415926536, -0.0000000000

PHP

<?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);
}
?>
Z:\>php Z:\0603.php
 1 :  3.1333333333, -0.0082593203
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415925025, -0.0000001511
 4 :  3.1415926512, -0.0000000024
 5 :  3.1415926536, -0.0000000000

Python

# coding: Shift_JIS

import math

def f(x):
    return 4 / (1 + x * x)

a = 0
b = 1

# Simpson則で積分
n = 2
for j in range(1, 6):
    h  = (b - a) / float(n)
    s2 = 0
    s4 = 0
    x  = a + h
    for i in range(1, (n / 2 + 1)):
        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

    # 結果を π と比較
    print "%2d : %13.10f, %13.10f" % (j, s, s - math.pi)
Z:\>python Z:\0603.py
 1 :  3.1333333333, -0.0082593203
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415925025, -0.0000001511
 4 :  3.1415926512, -0.0000000024
 5 :  3.1415926536, -0.0000000000

Ruby

def f(x)
    4 / (1 + x * x)
end

a = 0
b = 1

# Simpson則で積分
n = 2
(1..5).each do |j|
    h  = (b - a) / n.to_f
    s2 = 0
    s4 = 0
    x  = a + h
    (1..(n / 2)).each do |i|
        s4 += f(x)
        x  += h
        s2 += f(x)
        x  += h
    end
    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 - Math::PI)
end
Z:\>ruby Z:\0603.rb
 1 :  3.1333333333, -0.0082593203
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415925025, -0.0000001511
 4 :  3.1415926512, -0.0000000024
 5 :  3.1415926536, -0.0000000000

Groovy

Pascal

program Pas0603(arg);
{$MODE delphi}

uses
    SysUtils, Math;

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

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

        // 結果を π と比較
        writeln(format('%2d : %13.10f, %13.10f', [j, s, s - PI]));
    end;
end.
Z:\>fpc -v0 -l- Pas0603.pp

Z:\>Pas0603
 1 :  3.1333333333, -0.0082593203
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415925025, -0.0000001511
 4 :  3.1415926512, -0.0000000024
 5 :  3.1415926536,  0.0000000000

Ada

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

procedure Ada0603 is
    a : Constant Long_Float := 0.0;
    b : Constant Long_Float := 1.0;

    h, s, x : Long_Float;
    s2, s4  : Long_Float;
    n       : Integer;

    function f(x:Long_Float) return Long_Float is
    begin
        return 4.0 / (1.0 + x * x);
    end f;

begin
    -- Simpson則で積分
    n := 2;
    for j in 1..5 loop
        h  := (b - a) / Long_Float(n);
        s2 := 0.0;
        s4 := 0.0;
        x  := a + h;
        for i in 1..(n / 2) loop
            s4 := s4 + f(x);
            x  := x  + h;
            s2 := s2 + f(x);
            x  := x  + h;
        end loop;
        s2 := (s2 - f(b)) * 2.0 + f(a) + f(b);
        s4 := s4 * 4.0;
        s  := (s2 + s4) * h / 3.0;
        n  := n * 2;

        -- 結果を π と比較
        Put(j,                   Width=> 2);
        Put(" : ");
        Put(s ,                  Fore=>3, Aft=>10, Exp=>0);
        Put(", ");
        Put(s - Ada.Numerics.Pi, Fore=>3, Aft=>10, Exp=>0);
        New_Line;
    end loop;
end Ada0603;
xxxxxx@yyyyyy /Z
$ gnatmake Ada0603.adb

xxxxxx@yyyyyy /Z
$ Ada0603
 1 :   3.1333333333,  -0.0082593203
 2 :   3.1415686275,  -0.0000240261
 3 :   3.1415925025,  -0.0000001511
 4 :   3.1415926512,  -0.0000000024
 5 :   3.1415926536,  -0.0000000000

VB.NET

Module VB0603
    Public Sub Main()
        Const a As Double = 0
        Const b As Double = 1

        'Simpson則で積分
        Dim n As Integer = 2
        For j As Integer = 1 To 5
            Dim h  As Double = (b - a) / n
            Dim s2 As Double = 0
            Dim s4 As Double = 0
            Dim x  As Double = a + h
            For i As Integer = 1 To n \ 2
                s4 += f(x)
                x  += h
                s2 += f(x)
                x  += h
            Next
            s2 = (s2 - f(b)) * 2 + f(a) + f(b)
            s4 *= 4
            Dim s As Double = (s2 + s4) * h / 3
            n *= 2

            '結果を π と比較
            Console.WriteLine(String.Format("{0,2:D} : {1,13:F10}, {2,13:F10}", j, s, s - Math.PI))
        Next
    End Sub

    Private Function f(ByVal x As Double) As Double
        Return 4 / (1 + x * x)
    End Function
End Module
Z:\>vbc -nologo VB0603.vb

Z:\>VB0603
 1 :  3.1333333333, -0.0082593203
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415925025, -0.0000001511
 4 :  3.1415926512, -0.0000000024
 5 :  3.1415926536,  0.0000000000

C#

using System;

public class CS0603
{
    private static double f(double x)
    {
        return 4 / (1 + x * x);
    }

    public static void Main()
    {
        const double a = 0;
        const double b = 1;

        // Simpson則で積分
        int n = 2;
        for (int j = 1; j <= 5; j++)
        {
            double h  = (b - a) / n;
            double s2 = 0;
            double s4 = 0;
            double x  = a + h;
            for (int i = 1; i <= n / 2; i++)
            {
                s4 += f(x);
                x  += h;
                s2 += f(x);
                x  += h;
            }
            s2 = (s2 - f(b)) * 2 + f(a) + f(b);
            s4 *= 4;
            double s = (s2 + s4) * h / 3;
            n *= 2;

            // 結果を π と比較
            Console.WriteLine(string.Format("{0,2:D} : {1,13:F10}, {2,13:F10}", j, s, s - Math.PI));
        }
    }
}
Z:\>csc -nologo CS0603.cs

Z:\>CS0603
 1 :  3.1333333333, -0.0082593203
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415925025, -0.0000001511
 4 :  3.1415926512, -0.0000000024
 5 :  3.1415926536,  0.0000000000

Java

public class Java0603 {

    private static double f(double x) {
        return 4 / (1 + x * x);
    }

    public static void main(String []args) {
        final double a = 0;
        final double b = 1;

        // Simpson則で積分
        int n = 2;
        for (int j = 1; j <= 5; j++) {
            double h  = (b - a) / n;
            double s2 = 0;
            double s4 = 0;
            double x  = a + h;
            for (int i = 1; i <= n / 2; i++) {
                s4 += f(x);
                x  += h;
                s2 += f(x);
                x  += h;
            }
            s2 = (s2 - f(b)) * 2 + f(a) + f(b);
            s4 *= 4;
            double s = (s2 + s4) * h / 3;
            n *= 2;

            // 結果を π と比較
            System.out.println(String.format("%2d : %13.10f, %13.10f", j, s, s - Math.PI));
        }
    }
}
Z:\>javac Java0603.java

Z:\>java Java0603
 1 :  3.1333333333, -0.0082593203
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415925025, -0.0000001511
 4 :  3.1415926512, -0.0000000024
 5 :  3.1415926536, -0.0000000000

C++

#include <iostream>
#include <iomanip>
#include <math>

using namespace std;

double f(double x)
{
    return 4 / (1 + x * x);
}

int main()
{
    const double a = 0;
    const double b = 1;

    // Simpson則で積分
    int n = 2;
    for (int j = 1; j <= 5; j++)
    {
        double h = (b - a) / n;
        double s2 = 0;
        double s4 = 0;
        double x = a + h;
        for (int i = 1; i <= n / 2; i++)
        {
            s4 += f(x);
            x  += h;
            s2 += f(x);
            x  += h;
        }
        s2 = (s2 - f(b)) * 2 + f(a) + f(b);
        s4 *= 4;
        double s = (s2 + s4) * h / 3;
        n *= 2;

        // 結果を π と比較
        cout << setw(2)  << j << ":";
        cout << setw(13) << fixed << setprecision(10) << s        << ", ";
        cout << setw(13) << fixed << setprecision(10) << s - M_PI << endl;
    }

    return 0;
}
Z:\>bcc32 -q CP0603.cpp
cp0603.cpp:

Z:\>CP0603
 1: 3.1333333333, -0.0082593203
 2: 3.1415686275, -0.0000240261
 3: 3.1415925025, -0.0000001511
 4: 3.1415926512, -0.0000000024
 5: 3.1415926536, -0.0000000000

Objective-C

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

double f(double x)
{
    return 4 / (1 + x * x);
}

int main()
{
    const double a = 0;
    const double b = 1;

    // Simpson則で積分
    int n = 2;
    int i, j;
    for (j = 1; j <= 5; j++)
    {
        double h = (b - a) / n;
        double s2 = 0;
        double s4 = 0;
        double x = a + h;
        for (i = 1; i <= n / 2; i++)
        {
            s4 += f(x);
            x  += h;
            s2 += f(x);
            x  += h;
        }
        s2 = (s2 - f(b)) * 2 + f(a) + f(b);
        s4 *= 4;
        double s = (s2 + s4) * h / 3;
        n *= 2;

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

    return 0;
}
xxxxxx@yyyyyy /Z
$ gcc -o OC0603 OC0603.m -lobjc -lgnustep-base -I $INCLUDE -L $LIB $CFLAGS

xxxxxx@yyyyyy /Z
$ OC0603
 1 :  3.1333333333, -0.0082593203
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415925025, -0.0000001511
 4 :  3.1415926512, -0.0000000024
 5 :  3.1415926536, -0.0000000000

D

import std.stdio;
import std.math;

void main(string[] args)
{
    const double a = 0;
    const double b = 1;

    // Simpson則で積分
    int n = 2;
    for (int j = 1; j <= 5; j++)
    {
        double h = (b - a) / n;
        double s2 = 0;
        double s4 = 0;
        double x = a + h;
        for (int i = 1; i <= n / 2; i++)
        {
            s4 += f(x);
            x  += h;
            s2 += f(x);
            x  += h;
        }
        s2 = (s2 - f(b)) * 2 + f(a) + f(b);
        s4 *= 4;
        double s = (s2 + s4) * h / 3;
        n *= 2;

        // 結果を π と比較
        writefln("%2d : %13.10f, %13.10f", j, s, s - PI);
    }
}

double f(double x)
{
    return 4 / (1 + x * x);
}
Z:\>dmd D0603.d

Z:\>D0603
 1 :  3.1333333333, -0.0082593202
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415925025, -0.0000001511
 4 :  3.1415926512, -0.0000000023
 5 :  3.1415926536, -0.0000000000

Go

package main

import "fmt"
import "math"

func main() {
    const a = 0
    const b = 1

    // Simpson則で積分
    var n int = 2
    for j := 1; j <= 5; j++ {
        var h  float64 = (b - a) / float64(n)
        var s2 float64 = 0
        var s4 float64 = 0
        var x  float64 = a + h
        for i := 1; i <= n / 2; i++ {
            s4 += f(x)
            x  += h
            s2 += f(x)
            x  += h
        }
        s2 = (s2 - f(b)) * 2 + f(a) + f(b)
        s4 *= 4
        var s float64 = (s2 + s4) * h / 3
        n *= 2

        // 結果を π と比較
        fmt.Printf("%2d : %13.10f, %13.10f\n", j, s, s - math.Pi)
    }
}

func f(x float64) float64 {
    return (4 / (1 + x * x))
}
Z:\>8g GO0603.go

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

Z:\>GO0603
 1 :  3.1333333333, -0.0082593203
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415925025, -0.0000001511
 4 :  3.1415926512, -0.0000000024
 5 :  3.1415926536, -0.0000000000

Scala

object Scala0603 {
    def main(args: Array[String]) {
        val a:Double = 0
        val b:Double = 1

        // Simpson則で積分
        var n:Int = 2
        for (j <- 1 to 5) {
            var h:Double  = (b - a) / n
            var s2:Double = 0
            var s4:Double = 0
            var x:Double  = a + h
            for (i <- 1 to (n / 2)) {
                s4 += f(x)
                x  += h
                s2 += f(x)
                x  += h
            }
            s2 = (s2 - f(b)) * 2 + f(a) + f(b)
            s4 *= 4
            var s:Double = (s2 + s4) * h / 3
            n *= 2

            // 結果を π と比較
            println("%3d : %13.10f, %13.10f".format(j, s, s - Math.PI))
        }
    }

    def f(x:Double):Double = {
        4 / (1 + x * x)
    }
}
Z:\>scala Scala0603.scala
  1 :  3.1333333333, -0.0082593203
  2 :  3.1415686275, -0.0000240261
  3 :  3.1415925025, -0.0000001511
  4 :  3.1415926512, -0.0000000024
  5 :  3.1415926536, -0.0000000000

F#

module Fs0603

open System

let f(x:double):double =
    4.0 / (1.0 + x * x)

let a:double = 0.0
let b:double = 1.0

// Simpson則で積分
let mutable n = 2
for j in [1..5] do
    let         h:double  = (b - a) / (double n)
    let mutable s2:double = 0.0
    let mutable s4:double = 0.0
    let mutable x:double  = a + h

    for i in [1..(n / 2)] do
        s4 <- s4 + (f x)
        x  <- x  + h
        s2 <- s2 + (f x)
        x  <- x  + h
    s2 <- (s2 - (f b)) * 2.0 + (f a) + (f b)
    s4 <- s4 * 4.0
    let s = (s2 + s4) * h / 3.0
    n <- n * 2

    // 結果を π と比較
    printfn "%3d : %13.10f, %13.10f" j s (s - System.Math.PI)

exit 0
Z:\>fsi  --nologo --quiet Fs0603.fs
  1 :  3.1333333333, -0.0082593203
  2 :  3.1415686275, -0.0000240261
  3 :  3.1415925025, -0.0000001511
  4 :  3.1415926512, -0.0000000024
  5 :  3.1415926536,  0.0000000000

Clojure

(defn f[x]
    (/ 4 (+ 1 (* x x))))

(def a 0)
(def b 1)

; Simpson則で積分
(doseq [j (range 1 6)]
    (def n (Math/pow 2 j))
    (def h (/ (- b a) n))
    (def x (+ a h))
    (def s2 0)
    (def s4 0)
    (doseq [i (range 1 (+ (/ n 2) 1))]
        (def s4 (+ s4 (f x)))
        (def x  (+ x  h))
        (def s2 (+ s2 (f x)))
        (def x  (+ x  h))
    )
    (def w2 (double (+ (* (- s2 (f b)) 2) (f a) (f b))))
    (def w4 (* s4 4))
    (def t1 (double (/ (* (+ w2 w4) h) 3)))
    (def t2 (- t1 (. Math PI)))
    ; 結果を π と比較
    (println (format "%2d : %13.10f, %13.10f" j t1 t2)))
Z:\>java -cp C:\ProgramFiles\clojure-1.5.1\clojure-1.5.1.jar clojure.main Clj0603.clj
 1 :  3.1333333333, -0.0082593203
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415925025, -0.0000001511
 4 :  3.1415926512, -0.0000000024
 5 :  3.1415926536, -0.0000000000

Haskell

import Text.Printf
import Control.Monad

f::Double->Double
f x = 4 / (1 + (x * x))

main = do
    forM_ ([1..5::Integer]) $ \j -> do
        let n  = 2 ^ j
        let a  = 0.0
        let b  = 1.0
        let h  = (b - a) / (fromIntegral n)
        -- シンプソン則で積分
        let w4 = sum $ map(\x -> (f ((fromIntegral x) * h))) $ map(\i -> i * 2 - 1) $ [1..(div n 2)]
        let w2 = sum $ map(\x -> (f ((fromIntegral x) * h))) $ map(\i -> i * 2)     $ [1..(div n 2)]
        let t1 = (w2 - (f b)) * 2 + (f a) + (f b)
        let t2 = w4 * 4
        let t3 = (t1 + t2) * h / 3
        -- 結果を π と比較
        let t4 = t3 - pi
        printf "%3d : %13.10f, %13.10f\n" j t3 t4
Z:\>runghc Hs0603.hs
  1 :  3.1333333333, -0.0082593203
  2 :  3.1415686275, -0.0000240261
  3 :  3.1415925025, -0.0000001511
  4 :  3.1415926512, -0.0000000024
  5 :  3.1415926536, -0.0000000000
inserted by FC2 system