module Fs0101

printfn "%d" (3 + 5)
printfn "%d" (3 - 5)
printfn "%d" (3 * 5)
printfn "%f" (System.Math.Pow(float 3, float 5))
printfn "%d" (5 / 3)
printfn "%f" (5.0 / 3.0)
printfn "%d" (5 % 3)

printf  "%d\n" (3 * 5)

printfn "%3d"     (3 * 5)
printfn "%23.20f" (5.0 / 3.0)

exit 0
module Fs0102

let i = 3 * 5
printfn "3 * 5 = %d" i
exit 0
module Fs0103

for i in [1..9] do
    printf "%d, " i
printfn ""

exit 0
module Fs0104

for i in [1..9] do
    if (i % 3 = 0) then
        printf "%d, " i
printfn ""

exit 0
module Fs0105

let mutable sm = 0
for i in [1..99] do
    if (i % 3 = 0) then
        sm <- sm + i

printfn "%d" sm

exit 0
[1..9]
[1..9]
|> List.filter (fun x -> x % 3 = 0)
[1..99]
|> List.filter (fun x -> x % 3 = 0)
|> List.reduce(+)
[1..99]
|> List.filter (fun x -> x % 3 = 0)
|> List.sum
// 初項:a, 公差:a で, 上限:lim の数列の総和を返す関数
let sn a lim =
    let n = lim / a in  // 項数:n  =  上限:lim / 公差:a
    let l = a * n in    // 末項:l  =  項数:n   * 公差:a
    (a + l) * n / 2     // 総和:sn = (初項:a   + 末項:l) * 項数:n / 2
// 3 の倍数の合計
(sn 3 999)
let n = 10000
n * (n + 1) / 2
let n = 5000
n * (n + 1)
let n = 5000
System.Math.Pow(double n, 2.0)
let n = 1000
n * (n + 1) * (2 * n + 1) / 6
let n = 100
System.Math.Pow(double n, 2.0) * System.Math.Pow((double n + 1.0), 2.0) / 4.0
let n = 10
let a = 2
let r = 3
(a * (int (System.Math.Pow(double r, double n)) - 1)) / (r - 1)
[0..9]
[0..9]
|> List.map (fun n -> n * 3 + 5)
[0L..9L]
|> List.map (fun n -> n * 3L + 5L)
|> List.reduce(*)
// 等差数列の積
let rec product (m:int64) (d:int64) (n:int64) :int64 =
    match n with
        | 0L -> 1L
        | _  -> m * (product (m + d) d (n - 1L))
// 初項 5, 公差 3, 項数 10 の数列の積
product 5L 3L 10L
// 階乗を求める関数
let rec Fact = function
    |   0 -> 1
    |   n -> n * Fact(n - 1)
// 10の階乗
Fact 10
// 下降階乗冪
let rec FallingFact (x:int) (n:int):int =
    match n with
        | 1 -> x
        | _ -> x * (FallingFact (x - 1) (n - 1))
// 10 から 6 までの 総乗
FallingFact 10 5
// 上昇階乗冪
let rec RisingFact (x:int) (n:int):int =
    match n with
        | 1 -> x
        | _ -> x * (RisingFact (x + 1) (n - 1))
// 10 から 14 までの 総乗
RisingFact 10 5
// 階乗
let rec Fact = function
    |   0 -> 1
    |   n -> n * Fact(n - 1)
// 順列 (異なる 10 個のものから 5 個取ってできる順列の総数)
let n = 10
let r = 5
Fact n / Fact (n - r)
// 下降階乗冪
let rec FallingFact (x:int) (n:int):int =
    match n with
        | 1 -> x
        | _ -> x * (FallingFact (x - 1) (n - 1))
// 順列 (異なる 10 個のものから 5 個取ってできる順列の総数)
let n = 10
let r = 5
FallingFact n r
// 重複順列 (異なる 10 個のものから重複を許して 5 個取ってできる順列の総数)
let n = 10
let r = 5
System.Math.Pow(float n, float r)
// 組合せ
let rec Comb(n:int) (r:int):int =
    match n, r with
        | _, 0            -> 1
        | _, 1            -> n
        | _, _ when n = r -> 1
        | _, _            -> Comb (n - 1) (r - 1) + Comb(n - 1) r
// 異なる 10 個のものから 5 個取ってできる組合せの総数
let n = 10
let r = 5
Comb n r
// 組合せ
let rec Comb(n:int) (r:int):int =
    match n, r with
        | _, 0            -> 1
        | _, 1            -> n
        | _, _ when n = r -> 1
        | _, _            -> Comb (n - 1) (r - 1) + Comb(n - 1) r
// 異なる 10 個のものから重複を許して 5 個取る組合せの総数
let n = 10
let r = 5
Comb (n + r - 1) r
// 自作の正弦関数
let rec mySin (x:double) (n:int) (nega:bool) (numerator:double) (denominator:double) (y:double):double =
    let m     = 2 * n
    let denom = denominator * (double (m + 1)) * (double m)
    let num   = numerator   * x * x
    let a     = num / denom
    // 十分な精度になったら処理を抜ける
    if a <= 0.00000000001 then
        y
    else
        y + (mySin x (n + 1) (not nega) num denom (if nega then a else -a) )
[0..24]
|> List.map    (fun n -> n * 15)
|> List.filter (fun n -> (n % 30 = 0) || (n % 45 = 0))
|> List.iter
    (fun degree ->
        let radian = (double degree) * System.Math.PI / 180.0
        // 自作の正弦関数
        let d1 = (mySin radian 1 false radian 1.0 radian)
        // 標準の正弦関数
        let d2 = System.Math.Sin(radian)
        // 標準関数との差異
        printfn "%3d : %13.10f - %13.10f = %13.10f" degree d1 d2 (d1 - d2)
    )
// 自作の余弦関数
let rec myCos (x:double) (n:int) (nega:bool) (numerator:double) (denominator:double) (y:double):double =
    let m     = 2 * n
    let denom = denominator * (double m) * (double (m - 1))
    let num   = numerator   * x * x
    let a     = num / denom
    // 十分な精度になったら処理を抜ける
    if a <= 0.00000000001 then
        y
    else
        y + (myCos x (n + 1) (not nega) num denom (if nega then a else -a) )
[0..24]
|> List.map    ((*) 15)
|> List.filter (fun n -> (n % 30 = 0) || (n % 45 = 0))
|> List.iter
    (fun degree ->
        let radian = (double degree) * System.Math.PI / 180.0
        // 自作の余弦関数
        let d1 = (myCos radian 1 false 1.0 1.0 1.0)
        // 標準の余弦関数
        let d2 = System.Math.Cos(radian)
        // 標準関数との差異
        printfn "%3d : %13.10f - %13.10f = %13.10f" degree d1 d2 (d1 - d2)
    )
// 自作の正接関数
let rec myTan (x:double) (x2:double) (n:int) (t:double):double =
    let denom = x2 / ((double n) - t)
    let num   = n - 2
    if num <= 1 then
        x / (1.0 - denom)
    else
        myTan x x2 num denom
[0..12]
|> List.map    ((*) 15)
|> List.filter (fun n -> n % 180 <> 0)
|> List.map    ((-) 90)
|> List.iter
    (fun degree ->
        let radian = (double degree) * System.Math.PI / 180.0
        let x2     = radian * radian
        // 自作の正接関数
        let d1 = (myTan radian x2 15 0.0) // 15:必要な精度が得られる十分大きな奇数
        // 標準の正接関数
        let d2 = System.Math.Tan(radian)
        // 標準関数との差異
        printfn "%3d : %13.10f - %13.10f = %13.10f" degree d1 d2 (d1 - d2)
    )
// 自作の指数関数
let rec myExp (x:double) (n:int) (numerator:double) (denominator:double) (y:double):double =
    let denom = denominator * (double n)
    let num   = numerator   * x
    let a     = num / denom
    // 十分な精度になったら処理を抜ける
    if abs(a) <= 0.00000000001 then
        y
    else
        y + (myExp x (n + 1) num denom a)
[0..20]
|> List.map (fun n -> (double (n - 10)) / 4.0)
|> List.iter
    (fun x ->
        // 標準の指数関数
        let d1 = System.Math.Exp(x)
        // 自作の指数関数
        let d2 = (myExp x 1 1.0 1.0 1.0)
        // 標準関数との差異
        printfn "%5.2f : %13.10f - %13.10f = %13.10f" x d1 d2 (d1 - d2)
    )
// 自作の指数関数
let rec myExp (x:double) (x2:double) (n:int) (t:double):double =
    let denom = x2 / ((double n) + t)
    let num   = n - 4;

    if num < 6 then
        1.0 + ((2.0 * x) / (2.0 - x + denom))
    else
        myExp x x2 num denom
[0..20]
|> List.map    (fun n -> (double (n - 10)) / 4.0)
|> List.iter
    (fun x ->
        // 標準の指数関数
        let d1 = System.Math.Exp(x)
        // 自作の指数関数
        let x2 = x * x
        let d2 = (myExp x x2 30 0.0) // 30:必要な精度が得られるよう, 6から始めて4ずつ増加させる
        // 標準関数との差異
        printfn "%5.2f : %13.10f - %13.10f = %13.10f" x d1 d2 (d1 - d2)
    )
// 自作の対数関数
let rec myLog (x2:double) (numerator:double) (denominator:double) (y:double):double =
    let denom = denominator + 2.0
    let num   = numerator   * x2 * x2
    let a     = num / denom
    // 十分な精度になったら処理を抜ける
    if abs(a) <= 0.00000000001 then
        y
    else
        y + (myLog x2 num denom a)
[1..20]
|> List.map  (fun n -> (double n) / 5.0)
|> List.iter
    (fun x ->
        // 標準の対数関数
        let d1 = System.Math.Log(x)
        // 自作の対数関数
        let x2 = (x - 1.0) / (x + 1.0)
        let d2 = 2.0 * (myLog x2 x2 1.0 x2)
        // 標準関数との差異
        printfn "%5.2f : %13.10f - %13.10f = %13.10f" x d1 d2 (d1 - d2)
    )
// 自作の対数関数
let rec myLog (x:double) (n:int) (t:double):double =
    let mutable n2 = n
    let mutable x2 = x
    if n > 3 then
        if n % 2 = 0 then
            n2 <- 2
        x2 <- x * (double (n / 2))

    let t2 = x2 / ((double n2) + t)

    if n <= 2 then
        x / (1.0 + t2)
    else
        myLog x (n - 1) t2
[1..20]
|> List.map (fun n -> (double n) / 5.0)
|> List.iter
    (fun x ->
        // 標準の対数関数
        let d1 = System.Math.Log(x)
        // 自作の対数関数
        let d2 = (myLog (x - 1.0) 27 0.0) // 27:必要な精度が得られる十分大きな奇数
        // 標準関数との差異
        printfn "%5.2f : %13.10f - %13.10f = %13.10f" x d1 d2 (d1 - d2)
    )
// 自作の双曲線正弦関数
let rec mySinh (x:double) (n:int) (numerator:double) (denominator:double) (y:double):double =
    let m     = 2 * n
    let denom = denominator * (double (m + 1)) * (double m)
    let num   = numerator   * x * x
    let a     = num / denom
    // 十分な精度になったら処理を抜ける
    if abs(a) <= 0.00000000001 then
        y
    else
        y + (mySinh x (n + 1) num denom a)
[0..20]
|> List.map((-) 10)
|> List.iter
    (fun i ->
        // 自作の双曲線正弦関数
        let x  = (double i)
        let d1 = (mySinh x 1 x 1.0 x)
        // 標準の双曲線正弦関数
        let d2 = System.Math.Sinh(x)
        // 標準関数との差異
        printfn "%3d : %17.10f - %17.10f = %13.10f" i d1 d2 (d1 - d2)
    )
// 自作の双曲線余弦関数
let rec myCosh (x:double) (n:int) (numerator:double) (denominator:double) (y:double):double =
    let m     = 2 * n
    let denom = denominator * (double m) * (double (m - 1))
    let num   = numerator   * x * x
    let a     = num / denom
    // 十分な精度になったら処理を抜ける
    if abs(a) <= 0.00000000001 then
        y
    else
        y + (myCosh x (n + 1) num denom a)
[0..20]
|> List.map((-) 10)
|> List.iter
    (fun i ->
        let x  = (double i)
        // 自作の双曲線余弦関数
        let d1 = (myCosh x 1 1.0 1.0 1.0)
        // 標準の双曲線余弦関数
        let d2 = System.Math.Cosh(x)
        // 標準関数との差異
        printfn "%3d : %17.10f - %17.10f = %13.10f" i d1 d2 (d1 - d2)
    )
module Fs0601

open System

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

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

// 台形則で積分
let mutable n = 2
for j in [1..10] do
    let         h:double = (b - a) / (double n)
    let mutable s:double = 0.0
    let mutable x:double = a

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

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

exit 0
module Fs0602

open System

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

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

// 中点則で積分
let mutable n = 2
for j in [1..10] do
    let         h:double = (b - a) / (double n)
    let mutable s:double = 0.0
    let mutable x:double = a + (h / 2.0)

    for i in [1..n] do
        s <- s + (f x)
        x <- x + h
    s <- s * h
    n <- n * 2

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

exit 0
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
module Fs0604

open System

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

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

let t = Array2D.zeroCreate<double> 7 7

// 台形則で積分
let mutable n = 2
for i in [1..6] do
    let         h:double = (b - a) / (double n)
    let mutable s:double = 0.0
    let mutable x:double = a

    for j in [1..(n - 1)] do
        x <- x + h
        s <- s + (f x)

    // 結果を保存
    t.[i,1] <- h * (((f a) + (f b)) / 2.0 + s)
    n <- n * 2

// Richardsonの補外法
n <- 4
for j in [2..6] do
    for i in [j..6] do
        t.[i,j] <- t.[i, j - 1] + (t.[i, j - 1] - t.[i - 1, j - 1]) / (double (n - 1))
        if i = j then
            // 結果を π と比較
            printfn "%3d : %13.10f, %13.10f" j t.[i,j] (t.[i,j] - Math.PI)
    n <- n * 4

exit 0
module Fs0701

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

// Lagrange (ラグランジュ) 補間
let lagrange(d:double) (x:double list) (y:double list) =
    let mutable sum_list = [0.0]
    for i in [0..N] do
        let mutable prod_list = [y.[i]]
        for j in [0..N] do
            if i <> j then
                prod_list <- ((d - x.[j]) / (x.[i] - x.[j]))::prod_list
        sum_list <- (prod_list |> List.reduce(*))::sum_list
    sum_list |> List.sum

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

// 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 -> (lagrange d x y))

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

exit 0
module Fs0702

open System

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

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

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

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

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

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

exit 0
module Fs0703

open System

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

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

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

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

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

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

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

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

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

exit 0
module Fs0704

open System

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

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

// Hermite (エルミート) 補間
let hermite(d:double) (z:double list) (a:double list) =
    let mutable sum_list = [a.[0]]
    for i in [1..Nx2] do
        let mutable prod_list = [a.[i]]
        for j in [0..i-1] do
            prod_list <- (d - z.[j])::prod_list
        sum_list <- (prod_list |> List.reduce(*))::sum_list
    sum_list |> List.sum

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

// 差分商の表を作る
let z =  [0..Nx2] |> List.map(fun i -> i / 2) |> List.map(fun j -> x.[j])
let d = Array2D.zeroCreate<double> (Nx2+1) (Nx2+1)
for i in [0..Nx2] do
    d.[0,i] <- y.[i / 2]

for i in [1..Nx2] do
    for j in [0..Nx2-i] do
        if (i = 1 && j % 2 = 0) then
            d.[i,j] <- yd.[j / 2]
        else
            d.[i,j] <- (d.[i-1,j+1] - d.[i-1,j]) / (z.[j+i] - z.[j])

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

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

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

exit 0
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
module Fs0801

open System

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

// 角度
let degree = 45.0
let radian = degree * Math.PI / 180.0
// 初速 250 km/h -> 秒速に変換
let v = float(250 * 1000 / 3600)
// 水平方向の速度
let vx = v * Math.Cos(radian)
// 鉛直方向の速度
let vy = v * Math.Sin(radian)
// 位置
let x = 0.0
let y = 0.0

// 空気抵抗による水平方向の減速分
let fx(vx:Double) (vy:Double) =
    k * Math.Sqrt(vx * vx + vy * vy) * vx
// 重力と空気抵抗による鉛直方向の減速分
let fy(vx:Double) (vy:Double) =
    g + (k * Math.Sqrt(vx * vx + vy * vy) * vy)

// Euler法
let rec euler(i:int) (vx:double) (vy:double) (x:double) (y:double):unit =
    // 経過秒数
    let t = float(i) * h

    // 位置
    let wx = x + h * vx
    let wy = y + h * vy

    printfn "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f" t vx vy wx wy

    // 速度
    let wvx = vx + h * (fx vx vy)
    let wvy = vy + h * (fy vx vy)

    if wy >= 0.0 then
        (euler (i+1) wvx wvy wx wy)
    else
        ()

// Euler法
(euler 1 vx vy x y)

exit 0
module Fs0802

open System

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

// 角度
let degree = 45.0
let radian = degree * Math.PI / 180.0
// 初速 250 km/h -> 秒速に変換
let v = float(250 * 1000 / 3600)
// 水平方向の速度
let vx = v * Math.Cos(radian)
// 鉛直方向の速度
let vy = v * Math.Sin(radian)
// 位置
let x = 0.0
let y = 0.0

// 空気抵抗による水平方向の減速分
let fx(vx:Double) (vy:Double) =
    k * Math.Sqrt(vx * vx + vy * vy) * vx
// 重力と空気抵抗による鉛直方向の減速分
let fy(vx:Double) (vy:Double) =
    g + (k * Math.Sqrt(vx * vx + vy * vy) * vy)

// Heun法
let rec heun(i:int) (vx:double) (vy:double) (x:double) (y:double):unit =
    // 経過秒数
    let t = float(i) * h

    // 位置・速度
    let wx2  =  x + h *     vx
    let wy2  =  y + h *     vy
    let wvx2 = vx + h * (fx vx vy)
    let wvy2 = vy + h * (fy vx vy)

    let wx   =  x + h * (    vx     +     wvx2      ) / 2.0
    let wy   =  y + h * (    vy     +     wvy2      ) / 2.0
    let wvx  = vx + h * ((fx vx vy) + (fx wvx2 wvy2)) / 2.0
    let wvy  = vy + h * ((fy vx vy) + (fy wvx2 wvy2)) / 2.0

    printfn "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f" t wvx wvy wx wy

    if wy >= 0.0 then
        (heun (i+1) wvx wvy wx wy)
    else
        ()

// Heun法
(heun 1 vx vy x y)

exit 0
module Fs0803

open System

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

// 角度
let degree = 45.0
let radian = degree * Math.PI / 180.0
// 初速 250 km/h -> 秒速に変換
let v = float(250 * 1000 / 3600)
// 水平方向の速度
let vx = v * Math.Cos(radian)
// 鉛直方向の速度
let vy = v * Math.Sin(radian)
// 位置
let x = 0.0
let y = 0.0

// 空気抵抗による水平方向の減速分
let fx(vx:Double) (vy:Double) =
    k * Math.Sqrt(vx * vx + vy * vy) * vx
// 重力と空気抵抗による鉛直方向の減速分
let fy(vx:Double) (vy:Double) =
    g + (k * Math.Sqrt(vx * vx + vy * vy) * vy)

// 中点法
let rec midpoint(i:int) (vx:double) (vy:double) (x:double) (y:double):unit =
    // 経過秒数
    let t = float(i) * h

    // 位置・速度
    let wvx1 = h * (fx vx vy)
    let wvy1 = h * (fy vx vy)
    let wvx2 = vx + wvx1 / 2.0
    let wvy2 = vy + wvy1 / 2.0

    let wvx = vx + h * (fx wvx2 wvy2)
    let wvy = vy + h * (fy wvx2 wvy2)
    let wx  = x  + h *     wvx2
    let wy  = y  + h *     wvy2

    printfn "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f" t wvx wvy wx wy

    if wy >= 0.0 then
        (midpoint (i+1) wvx wvy wx wy)
    else
        ()

// 中点法
(midpoint 1 vx vy x y)

exit 0
module Fs0804

open System

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

// 角度
let degree = 45.0
let radian = degree * Math.PI / 180.0
// 初速 250 km/h -> 秒速に変換
let v = float(250 * 1000 / 3600)
// 水平方向の速度
let vx = v * Math.Cos(radian)
// 鉛直方向の速度
let vy = v * Math.Sin(radian)
// 位置
let x = 0.0
let y = 0.0

// 空気抵抗による水平方向の減速分
let fx(vx:Double) (vy:Double) =
    k * Math.Sqrt(vx * vx + vy * vy) * vx
// 重力と空気抵抗による鉛直方向の減速分
let fy(vx:Double) (vy:Double) =
    g + (k * Math.Sqrt(vx * vx + vy * vy) * vy)

// Runge-Kutta法
let rec rungekutta(i:int) (vx:double) (vy:double) (x:double) (y:double):unit =
    // 経過秒数
    let t = float(i) * h

    // 位置・速度
    let wx1  = h *     vx
    let wy1  = h *     vy
    let wvx1 = h * (fx vx vy)
    let wvy1 = h * (fy vx vy)

    let wvx5 = vx + wvx1 / 2.0
    let wvy5 = vy + wvy1 / 2.0
    let wx2  = h *     wvx5
    let wy2  = h *     wvy5
    let wvx2 = h * (fx wvx5 wvy5)
    let wvy2 = h * (fy wvx5 wvy5)

    let wvx6 = vx + wvx2 / 2.0
    let wvy6 = vy + wvy2 / 2.0
    let wx3  = h *     wvx6
    let wy3  = h *     wvy6
    let wvx3 = h * (fx wvx6 wvy6)
    let wvy3 = h * (fy wvx6 wvy6)

    let wvx7 = vx + wvx3
    let wvy7 = vy + wvy3
    let wx4  = h *     wvx7
    let wy4  = h *     wvy7
    let wvx4 = h * (fx wvx7 wvy7)
    let wvy4 = h * (fy wvx7 wvy7)

    let wx  = x  + ( wx1 +  wx2 * 2.0 +  wx3 * 2.0 +  wx4) / 6.0
    let wy  = y  + ( wy1 +  wy2 * 2.0 +  wy3 * 2.0 +  wy4) / 6.0
    let wvx = vx + (wvx1 + wvx2 * 2.0 + wvx3 * 2.0 + wvx4) / 6.0
    let wvy = vy + (wvy1 + wvy2 * 2.0 + wvy3 * 2.0 + wvy4) / 6.0

    printfn "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f" t wvx wvy wx wy

    if wy >= 0.0 then
        (rungekutta (i+1) wvx wvy wx wy)
    else
        ()

// Runge-Kutta法
(rungekutta 1 vx vy x y)

exit 0
module Fs0805

open System

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

// 角度
let degree = 45.0
let radian = degree * Math.PI / 180.0
// 初速 250 km/h -> 秒速に変換
let v = float(250 * 1000 / 3600)
// 水平方向の速度
let vx = v * Math.Cos(radian)
// 鉛直方向の速度
let vy = v * Math.Sin(radian)
// 位置
let x = 0.0
let y = 0.0

// 空気抵抗による水平方向の減速分
let fx(vx:Double) (vy:Double) =
    k * Math.Sqrt(vx * vx + vy * vy) * vx
// 重力と空気抵抗による鉛直方向の減速分
let fy(vx:Double) (vy:Double) =
    g + (k * Math.Sqrt(vx * vx + vy * vy) * vy)

// Runge-Kutta-Gill法
let rec rungekuttagill(i:int) (vx:double) (vy:double) (x:double) (y:double):unit =
    // 経過秒数
    let t = float(i) * h

    // 位置・速度
    let wx1  = h *     vx
    let wy1  = h *     vy
    let wvx1 = h * (fx vx vy)
    let wvy1 = h * (fy vx vy)

    let wvx5 = vx + wvx1 / 2.0
    let wvy5 = vy + wvy1 / 2.0
    let wx2  = h *     wvx5
    let wy2  = h *     wvy5
    let wvx2 = h * (fx wvx5 wvy5)
    let wvy2 = h * (fy wvx5 wvy5)

    let wvx6 = vx + wvx1 * ((Math.Sqrt(2.0) - 1.0) / 2.0) + wvx2 * (1.0 - 1.0 / Math.Sqrt(2.0))
    let wvy6 = vy + wvy1 * ((Math.Sqrt(2.0) - 1.0) / 2.0) + wvy2 * (1.0 - 1.0 / Math.Sqrt(2.0))
    let wx3  = h *     wvx6
    let wy3  = h *     wvy6
    let wvx3 = h * (fx wvx6 wvy6)
    let wvy3 = h * (fy wvx6 wvy6)

    let wvx7 = vx - wvx2 / Math.Sqrt(2.0) + wvx3 * (1.0 + 1.0 / Math.Sqrt(2.0))
    let wvy7 = vy - wvy2 / Math.Sqrt(2.0) + wvy3 * (1.0 + 1.0 / Math.Sqrt(2.0))
    let wx4  = h *     wvx7
    let wy4  = h *     wvy7
    let wvx4 = h * (fx wvx7 wvy7)
    let wvy4 = h * (fy wvx7 wvy7)

    let wx  = x  + ( wx1 +  wx2 * (2.0 - Math.Sqrt(2.0)) +  wx3 * (2.0 + Math.Sqrt(2.0)) +  wx4) / 6.0
    let wy  = y  + ( wy1 +  wy2 * (2.0 - Math.Sqrt(2.0)) +  wy3 * (2.0 + Math.Sqrt(2.0)) +  wy4) / 6.0
    let wvx = vx + (wvx1 + wvx2 * (2.0 - Math.Sqrt(2.0)) + wvx3 * (2.0 + Math.Sqrt(2.0)) + wvx4) / 6.0
    let wvy = vy + (wvy1 + wvy2 * (2.0 - Math.Sqrt(2.0)) + wvy3 * (2.0 + Math.Sqrt(2.0)) + wvy4) / 6.0

    printfn "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f" t wvx wvy wx wy

    if wy >= 0.0 then
        (rungekuttagill (i+1) wvx wvy wx wy)
    else
        ()

// Runge-Kutta-Gill法
(rungekuttagill 1 vx vy x y)

exit 0
module Fs0901

open System

let f (x:double):double =
    x * x - 2.0

let rec bisection (a:double) (b:double):double =
    // 区間 (a, b) の中点 c = (a + b) / 2
    let c = (a + b) / 2.0
    printfn "%12.10f\t%13.10f" c (c - Math.Sqrt(2.0))

    let fc = f c
    if abs(fc) < 0.0000000001 then
        c
    else
        if fc < 0.0 then
            // f(c) < 0 であれば, 解は区間 (c, b) の中に存在
            bisection c b
        else
            // f(c) > 0 であれば, 解は区間 (a, c) の中に存在
            bisection a c

let a = 1.0
let b = 2.0
printfn "%12.10f" (bisection a b)

exit 0
module Fs0902

open System

let f (x:double):double =
    x * x - 2.0

let rec falseposition (a:double) (b:double):double =
    // 点 (a,f(a)) と 点 (b,f(b)) を結ぶ直線と x軸の交点
    let c = (a * f(b) - b * f(a)) / (f(b) - f(a))
    printfn "%12.10f\t%13.10f" c (c - Math.Sqrt(2.0))

    let fc = f c
    if abs(fc) < 0.0000000001 then
        c
    else
        if fc < 0.0 then
            // f(c) < 0 であれば, 解は区間 (c, b) の中に存在
            falseposition c b
        else
            // f(c) > 0 であれば, 解は区間 (a, c) の中に存在
            falseposition a c

let a = 1.0
let b = 2.0
printfn "%12.10f" (falseposition a b)

exit 0
module Fs0903

open System

let g (x:double):double =
    (x / 2.0) + (1.0 / x)

let rec iterative (x0:double):double =
    let x1 = g(x0)
    printfn "%12.10f\t%13.10f" x1 (x1 - Math.Sqrt(2.0))

    if abs(x1 - x0) < 0.0000000001 then
        x1
    else
        iterative x1

let x = 1.0
printfn "%12.10f" (iterative x)

exit 0
module Fs0904

open System

let f0 (x:double):double =
    x * x - 2.0

let f1 (x:double):double =
    2.0 * x

let rec newton (x0:double):double =
    let x1 = x0 - (f0(x0) / f1(x0))
    printfn "%12.10f\t%13.10f" x1 (x1 - Math.Sqrt(2.0))

    if abs(x1 - x0) < 0.0000000001 then
        x1
    else
        newton x1

let x = 2.0
printfn "%12.10f" (newton x)

exit 0
module Fs0905

open System

let f0 (x:double):double =
    x * x - 2.0

let f1 (x:double):double =
    2.0 * x

let f2 (x:double):double =
    2.0

let rec bailey (x0:double):double =
    let x1 = x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2.0 * f1(x0)))))
    printfn "%12.10f\t%13.10f" x1 (x1 - Math.Sqrt(2.0))

    if abs(x1 - x0) < 0.0000000001 then
        x1
    else
        bailey x1

let x = 2.0
printfn "%12.10f" (bailey x)

exit 0
module Fs0906

open System

let f (x:double):double =
    x * x - 2.0

let rec secant (x0:double)(x1:double):double =
    let x2 = x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0))
    printfn "%12.10f\t%13.10f" x2 (x2 - Math.Sqrt(2.0))

    if abs(x2 - x1) < 0.0000000001 then
        x2
    else
        secant x1 x2

let x0 = 1.0
let x1 = 2.0
printfn "%12.10f" (secant x0 x1)

exit 0