import Text.Printf

main = do
    print (3 + 5)
    print (3 - 5)
    print (3 * 5)
    print (3 ^ 5)
    print (5 / 3)
    print (5 `div` 3)

    print (div  5 3)
    print (mod  5 3)
    print (quot 5 3)
    print (rem  5 3)

    print (div  (-5) 3)
    print (mod  (-5) 3)
    print (quot (-5) 3)
    print (rem  (-5) 3)

    putStrLn ((show (3 * 5)))

    printf "%3d\n"  ((3 * 5)::Int)
    printf "%23.20f\n" ((5 / 3)::Double)
import Text.Printf

main = do
    let i = (3 * 5)::Int
    printf "3 * 5 = %d\n" i
    putStrLn ("3 * 5 = " ++ (show i))
import Text.Printf
import Control.Monad

main = do
    forM [1..9::Int] $ \i -> printf "%d, " i
import Text.Printf
import Control.Monad

main = do
    forM [1..9::Int] $ \i -> do
        if (i `mod` 3 == 0)
            then printf "%d, " i
            else printf ""
total 0 = 0
total i =
    if i `mod` 3 == 0 then
        i + total (i - 1)
    else
        total (i - 1)

main = print (total 99)
[1..9]
filter (\x -> mod x 3 == 0) [1..9]
foldl (+) 0 $ filter (\x -> mod x 3 == 0) [1..99]
sum $ filter (\x -> mod x 3 == 0) [1..99]
-- 初項:a, 公差:a で, 上限:lim の数列の総和を返す関数
sn::Int->Int->Int
sn a lim =
    let
        n = lim `div` a      -- 項数:n  =  上限:lim / 公差:a
        l = a * n            -- 末項:l  =  項数:n   * 公差:a
    in
        (a + l) * n `div` 2  -- 総和:sn = (初項:a   + 末項:l) * 項数:n / 2

-- 3 の倍数の合計
main = putStrLn (show (sn 3 999))
let n = 10000
(n * (n + 1)) `div` 2
let n = 5000
n * (n + 1)
let n = 5000
n ^ 2
let n = 1000
n * (n + 1) * (2 * n + 1) `div` 6
let n = 100
(n ^ 2) * ((n + 1) ^ 2) `div` 4
let n = 10
let a = 2
let r = 3
a * (r ^ n - 1) `div` (r -1)
[0..9]
map(\n -> n * 3 + 5)[0..9]
product $ map(\n -> n * 3 + 5)[0..9]
-- 等差数列の積
prod m d 0 = 1
prod m d n =  m * (prod (m + d) d (n - 1))
-- 初項 5, 公差 3, 項数 10 の数列の積
prod 5 3 10
-- 階乗を求める関数
fact 0 = 1
fact n = n * fact (n-1)
-- 10の階乗
fact 10
-- 下降階乗冪
fallingFact x 1 = x
fallingFact x n = x * (fallingFact (x - 1) (n - 1))
-- 10 から 6 までの 総乗
fallingFact 10 5
-- 上昇階乗冪
risingFact x 1 = x
risingFact x n = x * (risingFact (x + 1) (n - 1))
-- 10 から 14 までの 総乗
risingFact 10 5
-- 階乗
fact 0 = 1
fact n = n * fact (n-1)
-- 順列 (異なる 10 個のものから 5 個取ってできる順列の総数)
let n = 10
let r = 5
(fact n) `div` (fact (n - r))
-- 下降階乗冪
fallingFact x 1 = x
fallingFact x n = 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
n ^ r
-- 組合せ
comb n 0 = 1
comb n 1 = n
comb n r
    | n == r    = 1
    | otherwise = (comb (n - 1) (r - 1)) + (comb (n - 1) r)
-- 異なる 10 個のものから 5 個取ってできる組合せの総数
let n = 10
let r = 5
comb n r
-- 組合せ
comb n 0 = 1
comb n 1 = n
comb n r
    | n == r    = 1
    | otherwise = (comb (n - 1) (r - 1)) + (comb (n - 1) r)
-- 異なる 10 個のものから重複を許して 5 個取ってできる組合せの総数
let n = 10
let r = 5
comb (n + r - 1) r
-- 自作の正弦関数
mySin::Double->Int->Bool->Double->Double->Double->Double
mySin x n nega numerator denominator y =
let
    m     = 2 * n
    denom = denominator * (fromIntegral (m + 1)) * (fromIntegral m)
    num   = numerator   * x * x
    a     = num / denom
in
    -- 十分な精度になったら処理を抜ける
    if a <= 0.00000000001 then
        y
    else
        y + (mySin x (n + 1) (not nega) num denom (if nega then a else -a))
import Text.Printf
import Control.Monad

forM_ (
    filter (\n -> (mod n 30 == 0 || mod n 45 == 0)) $
    map (* 15) $
    [0..24::Int]
) $ \degree -> do
    let radian = (fromIntegral degree) * pi / 180.0
    -- 自作の正弦関数
    let d1 = (mySin radian 1 False radian 1.0 radian)
    -- 標準の正弦関数
    let d2 = sin(radian)
    -- 標準関数との差異
    printf "%3d : %13.10f - %13.10f = %13.10f\n" degree d1 d2 (d1- d2)
-- 自作の余弦関数
myCos::Double->Int->Bool->Double->Double->Double->Double
myCos x n nega numerator denominator y =
let
    m     = 2 * n
    denom = denominator * (fromIntegral m) * (fromIntegral (m - 1))
    num   = numerator   * x * x
    a     = num / denom
in
    -- 十分な精度になったら処理を抜ける
    if a <= 0.00000000001 then
        y
    else
        y + (myCos x (n + 1) (not nega) num denom (if nega then a else -a))
import Text.Printf
import Control.Monad

forM_ (
    filter (\n -> (mod n 30 == 0 || mod n 45 == 0)) $
    map (* 15) $
    [0..24::Int]
) $ \degree -> do
    let radian = (fromIntegral degree) * pi / 180.0
    -- 自作の余弦関数
    let d1 = (myCos radian 1 False 1.0 1.0 1.0)
    -- 標準の余弦関数
    let d2 = cos(radian)
    -- 標準関数との差異
    printf "%3d : %13.10f - %13.10f = %13.10f\n" degree d1 d2 (d1- d2)
-- 自作の正接関数
myTan::Double->Double->Int->Double->Double
myTan x x2 n t =
let
    denom = x2 / ((fromIntegral n) - t)
    num   = n - 2
in
    if num <= 1 then
        x / (1.0 - denom)
    else
        myTan x x2 num denom
import Text.Printf
import Control.Monad

forM_ (
    map    (\n -> n - 90) $
    filter (\n -> mod n 180 /= 0) $
    map    (* 15) $
    [0..12::Int]
) $ \degree -> do
    let radian = (fromIntegral degree) * pi / 180.0
    -- 自作の正接関数
    let x2 = radian * radian
    let d1 = (myTan radian x2 15 0.0) -- 15:必要な精度が得られる十分大きな奇数
    -- 標準の正接関数
    let d2 = tan(radian)
    -- 標準関数との差異
    printf "%3d : %13.10f - %13.10f = %13.10f\n" degree d1 d2 (d1- d2)
-- 自作の指数関数
myExp::Double->Int->Double->Double->Double->Double
myExp x n numerator denominator y =
let
    denom = denominator * (fromIntegral n)
    num   = numerator   * x
    a     = num / denom
in
    -- 十分な精度になったら処理を抜ける
    if abs(a) <= 0.00000000001 then
        y
    else
        y + (myExp x (n + 1) num denom a)
import Text.Printf
import Control.Monad

forM_ (
    map (\n -> (fromIntegral (n - 10)) / 4.0) $
    [0..20::Int]
) $ \x -> do
    -- 標準の指数関数
    let d1 = exp(x)
    -- 自作の指数関数
    let 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)
-- 自作の指数関数
myExp::Double->Double->Int->Double->Double;
myExp x x2 n t =
let
    denom = x2 / ((fromIntegral n) + t)
    num   = n - 4
in
    if num < 6 then
        1.0 + ((2.0 * x) / (2.0 - x + denom))
    else
        myExp x x2 num denom
import Text.Printf
import Control.Monad

forM_ (
    map (\n -> (fromIntegral (n - 10)) / 4.0) $
    [0..20::Int]
) $ \x -> do
    -- 標準の指数関数
    let d1 = exp(x)
    -- 自作の指数関数
    let x2 = x * x
    let 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)
-- 自作の対数関数
myLog::Double->Double->Double->Double->Double;
myLog x2 numerator denominator y =
let
    denom = denominator + 2.0
    num   = numerator   * x2 * x2
    a     = num / denom
in
    -- 十分な精度になったら処理を抜ける
    if abs(a) <= 0.00000000001 then
        y
    else
        y + (myLog x2 num denom a)
import Text.Printf
import Control.Monad

forM_ (
    map (\n -> (fromIntegral n) / 5.0) $
    [1..20::Int]
) $ \x -> do
    -- 標準の対数関数
    let d1 = log(x)
    -- 自作の対数関数
    let x2 = (x - 1.0) / (x + 1.0)
    let d2 = 2.0 * (myLog x2 x2 1.0 x2)
    -- 標準関数との差異
    printf "%5.2f : %13.10f - %13.10f = %13.10f\n" x d1 d2 (d1- d2)
-- 自作の対数関数
myLog::Double->Int->Double->Double
myLog x n t =
let
    n2 = if rem n 2 == 0 then 2
                         else n
    n3 = if n > 3        then n2
                         else n
    x2 = if n > 3        then x * (fromIntegral (div n 2))
                         else x
    t2 = x2 / ((fromIntegral n3) + t)
in
    if n <= 2            then x / (1.0 + t2)
                         else myLog x (n - 1) t2
import Text.Printf
import Control.Monad

forM_ (
    map (\n -> (fromIntegral n) / 5.0) $
    [1..20::Int]
) $ \x -> do
    -- 標準の対数関数
    let d1 = log(x)
    -- 自作の対数関数
    let d2 = (myLog (x - 1.0) 27 0.0) -- 27:必要な精度が得られる十分大きな奇数
    -- 標準関数との差異
    printf "%5.2f : %13.10f - %13.10f = %13.10f\n" x d1 d2 (d1- d2)
-- 自作の双曲線正弦関数
mySinh::Double->Int->Double->Double->Double->Double
mySinh x n numerator denominator y =
let
    m     = 2 * n
    denom = denominator * (fromIntegral (m + 1)) * (fromIntegral m)
    num   = numerator   * x * x
    a     = num / denom
in
    -- 十分な精度になったら処理を抜ける
    if abs(a) <= 0.00000000001 then
        y
    else
        y + (mySinh x (n + 1) num denom a)
import Text.Printf
import Control.Monad

forM_ (
    map ((-) 10) $
    [0..20::Int]
) $ \i -> do
    -- 自作の双曲線正弦関数
    let x  = (fromIntegral i)
    let d1 = (mySinh x 1 x 1.0 x)
    -- 標準の双曲線正弦関数
    let d2 = sinh(x)
    -- 標準関数との差異
    printf "%3d : %17.10f - %17.10f = %13.10f\n" i d1 d2 (d1- d2)
-- 自作の双曲線余弦関数
myCosh::Double->Int->Double->Double->Double->Double;
myCosh x n numerator denominator y =
let
    m     = 2 * n
    denom = denominator * (fromIntegral m) * (fromIntegral (m - 1))
    num   = numerator   * x * x
    a     = num / denom
in
    -- 十分な精度になったら処理を抜ける
    if abs(a) <= 0.00000000001 then
        y
    else
        y + (myCosh x (n + 1) num denom a)
import Text.Printf
import Control.Monad

forM_ (
    map ((-) 10) $
    [0..20::Int]
) $ \i -> do
    -- 自作の双曲線余弦関数
    let x  = (fromIntegral i)
    let d1 = (myCosh x 1 1.0 1.0 1.0)
    -- 標準の双曲線余弦関数
    let d2 = cosh(x)
    -- 標準関数との差異
    printf "%3d : %17.10f - %17.10f = %13.10f\n" i d1 d2 (d1- d2)
-- 自作の逆正接関数
myAtan::Double->Double->Int->Double->Double
myAtan x x2 n t =
let
    m     = (fromIntegral (n `div` 2))
    denom = (m * m * x2) / ((fromIntegral n) + t)
    num   = n - 2
in
    if num <= 1 then
        x / (1.0 + denom)
    else
        myAtan x x2 num denom

import Text.Printf
import Control.Monad

forM_ (
    map (\n -> n - 45) $
    map (* 15) $
    [0..6::Int]
) $ \degree -> do
    let radian = (fromIntegral degree) * pi / 180.0
    -- 自作の逆正接関数
    let x2 = radian * radian
    let d1 = (myAtan radian x2 23 0.0) -- 23:必要な精度が得られる十分大きな奇数
    -- 標準の逆正接関数
    let d2 = atan(radian)
    -- 標準関数との差異
    printf "%3d : %13.10f - %13.10f = %13.10f\n" degree d1 d2 (d1- d2)
-- 自作の逆正接関数
myAtan::Double->Double->Int->Double->Double
myAtan x x2 n t =
let
    m     = (fromIntegral (n `div` 2))
    denom = (m * m * x2) / ((fromIntegral n) + t)
    num   = n - 2
in
    if num <= 1 then
        x / (1.0 + denom)
    else
        myAtan x x2 num denom

import Text.Printf
import Control.Monad

:{
forM_ (
    map (\n -> n * 2 + 1) $
    [5..15::Int]
) $ \i -> do -- 11..31
    let radian = 1
    -- 自作の逆正接関数
    let x2 = radian * radian
    let d1 = (myAtan radian x2 i 0.0) -- i:必要な精度が得られる十分大きな奇数
    -- 標準の逆正接関数
    let d2 = atan(radian)
    -- 標準関数との差異
    printf "%2d : %13.10f, %13.10f\n" i (d1 * 4) (d1 * 4 - pi)
:}


printf "%33.30f\n" pi
 3.14159 26535 89793
 3.14159 26535 89793 23846 26433 83279 50288 …

import Text.Printf
import Control.Monad

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

main = do
    forM_ ([1..10::Integer]) $ \j -> do
        let n = 2 ^ j
        let a = 0.0
        let b = 1.0
        let h = (b - a) / (fromIntegral n)
        -- 台形則で積分
        let w1 = sum $ map(\x -> f x) $ map(\i -> (fromIntegral i) * h + a) $ [1..(n - 1)]
        let w2 = ((f a) + (f b)) / 2
        let t1 = h * (w1 + w2)
        -- 結果を π と比較
        let t2 = t1 - pi
        printf "%3d : %13.10f, %13.10f\n" j t1 t2
import Text.Printf
import Control.Monad

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

main = do
    forM_ ([1..10::Integer]) $ \j -> do
        let n = 2 ^ j
        let a = 0.0
        let b = 1.0
        let h = (b - a) / (fromIntegral n)
        -- 中点則で積分
        let a2 = a + (h / 2)
        let w1 = sum $ map(\x -> f x) $ map(\i -> (fromIntegral i) * h + a2) $ [0..(n-1)]
        let t1 = h * w1
        -- 結果を π と比較
        let t2 = t1 - pi
        printf "%3d : %13.10f, %13.10f\n" j t1 t2
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

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

-- 台形則で積分する関数
trapezoidal::Integer->Integer->Double->Double->Double->Double
trapezoidal i n h x s
    | i <= (n - 1) = trapezoidal (i + 1) n h (x + h) (s + (f (x + h)))
    | otherwise    = s

import Control.Monad

-- 台形則で積分して, 結果を保存

t <- forM ([1..6::Integer]) $ \j -> do
    let n = 2 ^ j
    let a = 0
    let b = 1
    let h = (fromIntegral (b - a)) / (fromIntegral n)
    let x = 0::Double
    let s = 0::Double

    let w1 = trapezoidal 1 n h x 0::Double

    let w2 = ((f (fromIntegral a)) + (f (fromIntegral b))) / 2
    let t1 = h * (w1 + w2)
    let t2 = t1 - pi
    return t1


-- Richardsonの補外法
rich_sub n []     ss = ss
rich_sub n (x:[]) ss = ss
rich_sub n (x:xs) ss =
let
    x2 = (head xs)
    s  = x2 + (x2 - x) / (fromIntegral (n - 1))
in
    rich_sub n xs (s:ss)

richardson n ([])   = 0
richardson n (x:[]) = x
richardson n (x:xs) =
let
    s = rich_sub n (x:xs) []
in
    richardson (n * 4) (reverse s)

-- Richardsonの補外法で積分して, 結果を π と比較
import Text.Printf

forM_ ([2..6::Int]) $ \i -> do
    let s = richardson 4 (take i t)
    printf "%3d : %13.10f, %13.10f\n" i s (s - pi)







:{
let{
trapezoidal::Integer->Integer->Double->Double->Double->Double;
trapezoidal i n h x s
    | i <= (n - 1) = trapezoidal (i + 1) n h (x + h) (s + (f (x + h)))
    | otherwise    = s
}
:}

:{
t <- forM ([1..6::Integer]) $ \j -> do
    let n = 2 ^ j
    let a = 0
    let b = 1
    let h = (fromIntegral (b - a)) / (fromIntegral n)
    let x = 0::Double
    let s = 0::Double

    let w1 = trapezoidal 1 n h x 0::Double

    let w2 = ((f (fromIntegral a)) + (f (fromIntegral b))) / 2
    let t1 = h * (w1 + w2)
    let t2 = t1 - pi
    return t1
:}

Control.Monad> t
[3.1,3.1311764705882354,3.138988494491089,3.140941612041389,3.141429893174975,3.141551963485654]


:{
let {
rich_sub n []     ss = ss;
rich_sub n (x:[]) ss = ss;
rich_sub n (x:xs) ss =
let
    x2 = (head xs)
    s  = x2 + (x2 - x) / (fromIntegral (n - 1))
in
    rich_sub n xs (s:ss)
}
:}

:{
let {
richardson n ([])   = 0;
richardson n (x:[]) = x;
richardson n (x:xs) =
let
    s = rich_sub n (x:xs) []
in
    richardson (n * 4) (reverse s)
}
:}


:{
forM_ ([2..6::Int]) $ \i -> do
    let s = richardson 4 (take i t)
    printf "%3d : %13.10f, %13.10f\n" i s (s - pi)
:}


/*
  2 :  3.1415686275, -0.0000240261
  3 :  3.1415940941,  0.0000014405
  4 :  3.1415926384, -0.0000000152
  5 :  3.1415926536,  0.0000000001
  6 :  3.1415926536, -0.0000000000
*/































midpoint



C:\Users\Ifthen>ghci
GHCi, version 7.6.3: http://www.haskell.org/ghc/  :? for help
Loading package ghc-prim ... linking ... done.
Loading package integer-gmp ... linking ... done.
Loading package base ... linking ... done.
Prelude> :{
Prelude| let {
Prelude| f::Double->Double;
Prelude| f x = 4 / (1 + (x * x))
Prelude| }
Prelude| :}
Prelude> :{
Prelude| let {
Prelude| trape::Integer->Integer->Double->Double->Double->Double;
Prelude| trape i n h x s
Prelude|     | i <= (n - 1) = trape (i + 1) n h (x + h) (s + (f (x + h)))
Prelude|     | otherwise    = s
Prelude| }
Prelude| :}
Prelude> import Text.Printf
Prelude Text.Printf> import Control.Monad
Prelude Text.Printf Control.Monad> :{
Prelude Text.Printf Control.Monad| forM_ ([1..10::Integer]) $ \j -> do
Prelude Text.Printf Control.Monad|     let n = 2 ^ j
Prelude Text.Printf Control.Monad|     let a = 0
Prelude Text.Printf Control.Monad|     let b = 1
Prelude Text.Printf Control.Monad|     let h = (fromIntegral (b - a)) / (fromIntegral n)
Prelude Text.Printf Control.Monad|     let x = 0::Double
Prelude Text.Printf Control.Monad|     let s = 0::Double
Prelude Text.Printf Control.Monad|
Prelude Text.Printf Control.Monad|     let w1 = trape 1 n h x 0::Double
Prelude Text.Printf Control.Monad|
Prelude Text.Printf Control.Monad|     let w2 = ((f (fromIntegral a)) + (f (fromIntegral b))) / 2::Double
Prelude Text.Printf Control.Monad|     let t1 =  (realToFrac h) * (w1 + w2)
Prelude Text.Printf Control.Monad|     let t2 = t1 - pi
Prelude Text.Printf Control.Monad|     printf "%3d : %23.20f, %23.20f\n" j t1 t2
Prelude Text.Printf Control.Monad| :}
  1 :  3.10000000000000000000, -0.04159265358979303000
  2 :  3.13117647058823540000, -0.01041618300155766600
  3 :  3.13898849449108900000, -0.00260415909870426180
  4 :  3.14094161204138900000, -0.00065104154840422980
  5 :  3.14142989317497500000, -0.00016276041481821935
  6 :  3.14155196348565400000, -0.00004069010413898511
  7 :  3.14158248106375200000, -0.00001017252604107455
  8 :  3.14159011045828150000, -0.00000254313151160090
  9 :  3.14159201780692140000, -0.00000063578287168298
 10 :  3.14159249464407250000, -0.00000015894572058528
Prelude Text.Printf Control.Monad> :{
Prelude Text.Printf Control.Monad| t <- forM ([1..10::Integer]) $ \j -> do
Prelude Text.Printf Control.Monad|     let n = 2 ^ j
Prelude Text.Printf Control.Monad|     let a = 0
Prelude Text.Printf Control.Monad|     let b = 1
Prelude Text.Printf Control.Monad|     let h = (fromIntegral (b - a)) / (fromIntegral n)
Prelude Text.Printf Control.Monad|     let x = 0::Double
Prelude Text.Printf Control.Monad|     let s = 0::Double
Prelude Text.Printf Control.Monad|
Prelude Text.Printf Control.Monad|     let w1 = trape 1 n h x 0::Double
Prelude Text.Printf Control.Monad|
Prelude Text.Printf Control.Monad|     let w2 = ((f (fromIntegral a)) + (f (fromIntegral b))) / 2
Prelude Text.Printf Control.Monad|     let t1 = h * (w1 + w2)
Prelude Text.Printf Control.Monad|     let t2 = t1 - pi
Prelude Text.Printf Control.Monad|     return t1
Prelude Text.Printf Control.Monad| :}
Prelude Text.Printf Control.Monad> mapM_ print t
3.1
3.1311764705882354
3.138988494491089
3.140941612041389
3.141429893174975
3.141551963485654
3.141582481063752
3.1415901104582815
3.1415920178069214
3.1415924946440725
Prelude Text.Printf Control.Monad> forM_ t print
3.1
3.1311764705882354
3.138988494491089
3.140941612041389
3.141429893174975
3.141551963485654
3.141582481063752
3.1415901104582815
3.1415920178069214
3.1415924946440725
Prelude Text.Printf Control.Monad> :{
Prelude Text.Printf Control.Monad| let{
Prelude Text.Printf Control.Monad| rich n []     ss = ss;
Prelude Text.Printf Control.Monad| rich n (x:[]) ss = ss;
Prelude Text.Printf Control.Monad| rich n (x:xs) ss =
Prelude Text.Printf Control.Monad| let
Prelude Text.Printf Control.Monad|     x2 = (head xs)
Prelude Text.Printf Control.Monad|     s  = x2 + (x2- x) / (fromIntegral (n - 1))
Prelude Text.Printf Control.Monad| in
Prelude Text.Printf Control.Monad|     rich n xs (s:ss)
Prelude Text.Printf Control.Monad| }
Prelude Text.Printf Control.Monad| :}
Prelude Text.Printf Control.Monad> :{
Prelude Text.Printf Control.Monad| let{
Prelude Text.Printf Control.Monad| call_rich n ([])   = 0;
Prelude Text.Printf Control.Monad| call_rich n (x:[]) = x;
Prelude Text.Printf Control.Monad| call_rich n (x:xs) =
Prelude Text.Printf Control.Monad| let
Prelude Text.Printf Control.Monad|     s = rich n (x:xs) []
Prelude Text.Printf Control.Monad| in
Prelude Text.Printf Control.Monad|     call_rich (n * 4) (reverse s)
Prelude Text.Printf Control.Monad| }
Prelude Text.Printf Control.Monad| :}
Prelude Text.Printf Control.Monad> :{
Prelude Text.Printf Control.Monad| forM_ ([2..6::Int]) $ \i -> do
Prelude Text.Printf Control.Monad|     let s = call_rich 4 (take i t)
Prelude Text.Printf Control.Monad|     printf "%3d : %13.10f, %13.10f\n" i s (s - pi)
Prelude Text.Printf Control.Monad| :}
  2 :  3.1415686275, -0.0000240261
  3 :  3.1415940941,  0.0000014405
  4 :  3.1415926384, -0.0000000152
  5 :  3.1415926536,  0.0000000001
  6 :  3.1415926536, -0.0000000000
Prelude Text.Printf Control.Monad> :q
Leaving GHCi.

C:\Users\Ifthen>

*/






























import Text.Printf
import Control.Monad

:{
let {
f::Double->Double;
f x = 4 / (1 + (x * x))
}
:}
:{
let {
trape::Integer->Integer->Double->Double->Double->Double;
trape i n h x s
    | i <= (n - 1) = trape (i + 1) n h (x + h) (s + (f (x + h)))
    | otherwise    = s
}
:}
:{
forM_ ([1..6::Integer]) $ \j -> do
    let n = 2 ^ j
    let a = 0
    let b = 1
    let h = (fromIntegral (b - a)) / (fromIntegral n)
    let x = 0::Double
    let s = 0::Double

    let w1 = trape 1 n h x 0::Double

    let w2 = ((f (fromIntegral a)) + (f (fromIntegral b))) / 2::Double
    let t1 =  (realToFrac h) * (w1 + w2)
    let t2 = t1 - pi
    printf "%3d : %23.20f, %23.20f\n" j t1 t2
:}
  1 :  3.1000000000, -0.0415926536
  2 :  3.1311764706, -0.0104161830
  3 :  3.1389884945, -0.0026041591
  4 :  3.1409416120, -0.0006510415
  5 :  3.1414298932, -0.0001627604
  6 :  3.1415519635, -0.0000406901
  7 :  3.1415824811, -0.0000101725
  8 :  3.1415901105, -0.0000025431
  9 :  3.1415920178, -0.0000006358
 10 :  3.1415924946, -0.0000001589

:{
t <- forM ([1..10::Integer]) $ \j -> do
    let n = 2 ^ j
    let a = 0
    let b = 1
    let h = (fromIntegral (b - a)) / (fromIntegral n)
    let x = 0::Double
    let s = 0::Double

    let w1 = trape 1 n h x 0::Double

    let w2 = ((f (fromIntegral a)) + (f (fromIntegral b))) / 2
    let t1 = h * (w1 + w2)
    let t2 = t1 - pi
    return t1
:}


mapM_ print t
forM_ t print

:{
let{
rich n []     ss = ss;
rich n (x:[]) ss = ss;
rich n (x:xs) ss =
let
    x2 = (head xs)
    s  = x2 + (x2- x) / (fromIntegral (n - 1))
in
    rich n xs (s:ss)
}
:}
let s = test1 4 t []
mapM_ print s

let s = test1 4 [3.1,3.13] []
mapM_ print s

let s = test1 4 [3.1] []
mapM_ print s

let s = test1 4 [] []
mapM_ print s

:{
let{
call_rich n ([])   = 0;
call_rich n (x:[]) = x;
call_rich n (x:xs) =
let
    s = rich n (x:xs) []
in
    call_rich (n * 4) (reverse s)
}
:}
let s = test2 4 []
s

let s = test2 4 [3.1]
s

let s = test2 4 [3.1,3.13]
s

let s = test2 4 [3.1,3.13,3.1311]
s

let s = test2 4 t
s


:{
forM_ ([2..6::Int]) $ \i -> do
    let s = call_rich 4 (take i t)
    printf "%3d : %13.10f, %13.10f\n" i s (s - pi)
:}


let x1 = (head t)
let x2 = (head (tail t))
let t2_2  = x2 + ((x2- x1) / (fromIntegral (4 - 1)))
printf "%13.10f\n" t2_2

let x1 = (head (tail t))
let x2 = (head (tail (tail t)))
let t2_3  = x2 + ((x2- x1) / (fromIntegral (4 - 1)))
printf "%13.10f\n" t2_3

let x1 = (head (tail (tail t)))
let x2 = (head (tail (tail (tail t))))
let t2_4  = x2 + ((x2- x1) / (fromIntegral (4 - 1)))
printf "%13.10f\n" t2_4

let x1 = (head (tail (tail (tail t))))
let x2 = (head (tail (tail (tail (tail t)))))
let t2_5  = x2 + ((x2- x1) / (fromIntegral (4 - 1)))
printf "%13.10f\n" t2_5

let x1 = (head (tail (tail (tail (tail t)))))
let x2 = (head (tail (tail (tail (tail (tail t))))))
let t2_6  = x2 + ((x2- x1) / (fromIntegral (4 - 1)))
printf "%13.10f\n" t2_6

let t3_3  = t2_3 + ((t2_3- t2_2) / (fromIntegral (16 - 1)))
printf "%13.10f\n" t3_3


let s = test1 4 t []
mapM_ print s


in
    test1 n xs (s:ss)





:{
forM_ ([t_list]) $ \t -> do
    print  (show t)
:}

:{
forM_ ([1..10::Integer]) $ \j -> do
    n   <- return  (2 ^ j)
    a  <- return  0
    b  <- return  1
    h  <- return  ((fromIntegral (b - a)) / (fromIntegral n))
    x  <- return  (0::Double)
    s  <- return  (0::Double)
    w1  <- return  (trape_sub 1 n h x 0::Double)
    w2  <- return  (((f (fromIntegral a)) + (f (fromIntegral b))) / 2)
    t1 <- return  (h * (w1 + w2))
    t2 <- return  (t1 - pi)
    printf "%3d : %13.10f, %13.10f\n" j t1 t2
:}


:{
forM_ (
    [1..10::Integer]
) $ \j -> do
    let n = 2 ^ j
    let a = 0
    let b = 1
    let h = (b - a) / n
    let w = trape_sub  (toInteger n) (fromIntegral h) 1 (fromIntegral a)
    -- printf "%3d : %13.10f\n" j t1
    printf "%3d\n" j
:}






let {
trape::Double->Double->Integer->Double;
trape a b j =
let
    n = 2 ^ j
    h = (b - a) / n
in
    if j <= 10 then (trape_sub  (toInteger n) (fromIntegral h) 1 (fromIntegral a) 0) then
        trape a b (j + 1)
    else
        0
}
:}



import Text.Printf
import Control.Monad

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

-- Richardsonの補外法
rich_sub n []     ss = ss
rich_sub n (x:[]) ss = ss
rich_sub n (x:xs) ss =
    let
        x2 = (head xs)
        s  = x2 + (x2 - x) / (fromIntegral (n - 1))
    in
        rich_sub n xs (s:ss)

richardson n ([])   = 0
richardson n (x:[]) = x
richardson n (x:xs) =
    let
        s = rich_sub n (x:xs) []
    in
        richardson (n * 4) (reverse s)

main = do
    let a = 0.0
    let b = 1.0
    t <- forM ([1..6::Integer]) $ \j -> do
        let n = 2 ^ j
        let h = (b - a) / (fromIntegral n)
        -- 台形則で積分
        let w1 = sum $ map(\x -> f x) $ map(\i -> (fromIntegral i) * h + a) $ [1..(n - 1)]
        let w2 = ((f a) + (f b)) / 2
        let t1 = h * (w1 + w2)
        return t1

    -- Richardsonの補外法で積分して, 結果を π と比較
    forM_ ([2..6::Int]) $ \i -> do
        let s = richardson 4 (take i t)
        printf "%3d : %13.10f, %13.10f\n" i s (s - pi)
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))

-- Lagrange (ラグランジュ) 補間
lagrange::Double->[Double]->[Double]->Double
lagrange d x y =
    let sum_list = map(\i -> do
        let prod_list = map(\j -> do
            (d - x!!(fromIntegral j)) /
                (x!!(fromIntegral i)  -
                 x!!(fromIntegral j)))
                $ filter (\j -> j /= i) [0..n]
        product $ y!!(fromIntegral i) : prod_list
        ) [0..n]
    in
        sum $ sum_list

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

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

-- Lagrange (ラグランジュ) 補間
lagrange::Double->[Double]->[Double]->Double
lagrange d x y =
    let sum_list = map(\i -> do
        let prod_list = map(\j -> do
            (d - x!!j) / (x!!i - x!!j)
            ) $ filter (\j -> j /= i) [0..n::Int]
        product $ y!!i : prod_list
        ) [0..n::Int]
    in
        sum $ sum_list

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

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

-- Neville (ネヴィル) 補間
neville::Double->[Double]->[Double]->Int->Double
neville d x w j =
    let
        t = map(\i -> do
                (w!!(i+1) + (w!!(i+1) - w!!i) * (d - x!!(i+j)) / (x!!(i+j) - x!!i))
            ) [0..(n-j)::Int]
    in
        if j == n then t!!0
                  else (neville d x t (j + 1))

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

    -- 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 -> (neville i x y 1)) 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)
import Text.Printf
import Control.Monad

-- データ点の数 - 1
n = 6 :: Int

-- 元の関数
f::Double->Double
f x = x - (x ^ 3) / (fromIntegral (3 * 2)) + (x ^ 5) / (fromIntegral (5 * 4 * 3 * 2))

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

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

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

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

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

    forM_ (zip (zip d1 d2) d3) $ \((d1, d2), d3) -> do
        printf "%5.2f\t%8.5f\t%8.5f\t%8.5f\n" d1 d2 d3 (d2 - d3)
import Text.Printf
import Control.Monad

-- データ点の数 - 1
n   =  6 :: Int
nx2 = 13 :: Int

-- 元の関数
f::Double->Double
f x = x - (x ^ 3) / (fromIntegral (3 * 2)) + (x ^ 5) / (fromIntegral (5 * 4 * 3 * 2))

-- 導関数
fd::Double->Double
fd x = 1.0 - (x ^ 2) / (fromIntegral 2) + (x ^ 4) / (fromIntegral (4 * 3 * 2))

-- Hermite (エルミート) 補間
hermite::Double->[Double]->[Double]->Double
hermite d z a =
    let sum_list = map(\i -> do
        let prod_list = map(\j -> do
            d - z!!j
            ) $ [0..(i-1)::Int]
        product $ a!!i : prod_list
        ) [1..nx2::Int]
    in
        sum $ a!!0 : sum_list

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

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

    -- 差分商の表を作る
    let z = map(\i -> x!!(div i 2)) [0..nx2]
    let w = map(\i -> y!!(div i 2)) [0..nx2]
    let a = reverse (make_table yd z w [w!!0] 1)

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

    forM_ (zip (zip d1 d2) d3) $ \((d1, d2), d3) -> do
        printf "%5.2f\t%8.5f\t%8.5f\t%8.5f\n" d1 d2 d3 (d2 - d3)
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))

-- 3項方程式を解く (ト−マス法)
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 まで, 7点だけ値をセット
    let x = map(\i -> (fromIntegral i) * 1.5 - 4.5) [0..n]
    let y = map(\i -> f(i)) x

    -- 3項方程式の係数の表を作る
    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]

    -- 3項方程式を解く (ト−マス法)
    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)
import Text.Printf

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

-- 空気抵抗による水平方向の減速分
fx::Double->Double->Double
fx vx vy =
    let
        v = sqrt(vx * vx + vy * vy)
    in
        k * v * vx

-- 空気抵抗による鉛直方向の減速分
fy::Double->Double->Double
fy vx vy =
    let
        v = sqrt(vx * vx + vy * vy)
    in
        g + (k * v * vy)

-- Euler法
euler::Integer->Double->Double->Double->Double->IO ()
euler i vx vy x y =
    let
        t = (fromIntegral i) * h
        wx = x + h * vx
        wy = y + h * vy
        wvx = vx + h * (fx vx vy)
        wvy = vy + h * (fy vx vy)
    in
        if y >= 0.0
            then do
                printf "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f\n" t vx vy wx wy
                euler (i+1) wvx wvy wx wy
            else
                return ()

main = do
    -- 角度
    let degree = 45.0 :: Double
    let radian = degree * pi / 180.0
    -- 初速 250 km/h -> 秒速に変換
    let v = (fromIntegral (250 * 1000 `div` 3600))

    -- 水平方向の速度
    let vx = v * cos(radian)
    -- 鉛直方向の速度
    let vy = v * sin(radian)
    -- 位置
    let x = 0.0
    let y = 0.0

    -- Euler法
    euler 1 vx vy x y
import Text.Printf

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

-- 空気抵抗による水平方向の減速分
fx::Double->Double->Double
fx vx vy =
    let
        v = sqrt(vx * vx + vy * vy)
    in
        k * v * vx

-- 空気抵抗による鉛直方向の減速分
fy::Double->Double->Double
fy vx vy =
    let
        v = sqrt(vx * vx + vy * vy)
    in
        g + (k * v * vy)

-- Heun法
heun::Integer->Double->Double->Double->Double->IO ()
heun i vx vy x y =
    let
        t = (fromIntegral i) * h

        wx2  =  x + h *     vx
        wy2  =  y + h *     vy
        wvx2 = vx + h * (fx vx vy)
        wvy2 = vy + h * (fy vx vy)

        wx   =  x + h * (    vx     +     wvx2      ) / 2.0
        wy   =  y + h * (    vy     +     wvy2      ) / 2.0
        wvx  = vx + h * ((fx vx vy) + (fx wvx2 wvy2)) / 2.0
        wvy  = vy + h * ((fy vx vy) + (fy wvx2 wvy2)) / 2.0
    in
        if y >= 0.0
            then do
                printf "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f\n" t wvx wvy wx wy
                heun (i+1) wvx wvy wx wy
            else
                return ()

main = do
    -- 角度
    let degree = 45.0 :: Double
    let radian = degree * pi / 180.0
    -- 初速 250 km/h -> 秒速に変換
    let v = (fromIntegral (250 * 1000 `div` 3600))

    -- 水平方向の速度
    let vx = v * cos(radian)
    -- 鉛直方向の速度
    let vy = v * sin(radian)
    -- 位置
    let x = 0.0
    let y = 0.0

    -- Heun法
    heun 1 vx vy x y
import Text.Printf

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

-- 空気抵抗による水平方向の減速分
fx::Double->Double->Double
fx vx vy =
    let
        v = sqrt(vx * vx + vy * vy)
    in
        k * v * vx

-- 空気抵抗による鉛直方向の減速分
fy::Double->Double->Double
fy vx vy =
    let
        v = sqrt(vx * vx + vy * vy)
    in
        g + (k * v * vy)

-- 中点法
midpoint::Integer->Double->Double->Double->Double->IO ()
midpoint i vx vy x y =
    let
        -- 経過秒数
        t = (fromIntegral i) * h

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

        wvx = vx + h * (fx wvx2 wvy2)
        wvy = vy + h * (fy wvx2 wvy2)
        wx  = x  + h *     wvx2
        wy  = y  + h *     wvy2
    in
        if y >= 0.0
            then do
                printf "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f\n" t wvx wvy wx wy
                midpoint (i+1) wvx wvy wx wy
            else
                return ()

main = do
    -- 角度
    let degree = 45.0 :: Double
    let radian = degree * pi / 180.0
    -- 初速 250 km/h -> 秒速に変換
    let v = (fromIntegral (250 * 1000 `div` 3600))

    -- 水平方向の速度
    let vx = v * cos(radian)
    -- 鉛直方向の速度
    let vy = v * sin(radian)
    -- 位置
    let x = 0.0
    let y = 0.0

    -- 中点法
    midpoint 1 vx vy x y
import Text.Printf

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

-- 空気抵抗による水平方向の減速分
fx::Double->Double->Double
fx vx vy =
    let
        v = sqrt(vx * vx + vy * vy)
    in
        k * v * vx

-- 空気抵抗による鉛直方向の減速分
fy::Double->Double->Double
fy vx vy =
    let
        v = sqrt(vx * vx + vy * vy)
    in
        g + (k * v * vy)

-- Runge-Kutta法
rungekutta::Integer->Double->Double->Double->Double->IO ()
rungekutta i vx vy x y =
    let
        -- 経過秒数
        t = (fromIntegral i) * h

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

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

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

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

        wx  = x  + ( wx1 +  wx2 * 2.0 +  wx3 * 2.0 +  wx4) / 6.0
        wy  = y  + ( wy1 +  wy2 * 2.0 +  wy3 * 2.0 +  wy4) / 6.0
        wvx = vx + (wvx1 + wvx2 * 2.0 + wvx3 * 2.0 + wvx4) / 6.0
        wvy = vy + (wvy1 + wvy2 * 2.0 + wvy3 * 2.0 + wvy4) / 6.0
    in
        if y >= 0.0
            then do
                printf "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f\n" t wvx wvy wx wy
                rungekutta (i+1) wvx wvy wx wy
            else
                return ()

main = do
    -- 角度
    let degree = 45.0 :: Double
    let radian = degree * pi / 180.0
    -- 初速 250 km/h -> 秒速に変換
    let v = (fromIntegral (250 * 1000 `div` 3600))

    -- 水平方向の速度
    let vx = v * cos(radian)
    -- 鉛直方向の速度
    let vy = v * sin(radian)
    -- 位置
    let x = 0.0
    let y = 0.0

    -- Runge-Kutta法
    rungekutta 1 vx vy x y
import Text.Printf

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

-- 空気抵抗による水平方向の減速分
fx::Double->Double->Double
fx vx vy =
    let
        v = sqrt(vx * vx + vy * vy)
    in
        k * v * vx

-- 空気抵抗による鉛直方向の減速分
fy::Double->Double->Double
fy vx vy =
    let
        v = sqrt(vx * vx + vy * vy)
    in
        g + (k * v * vy)

-- Runge-Kutta-Gill法
rungekuttagill::Integer->Double->Double->Double->Double->IO ()
rungekuttagill i vx vy x y =
    let
        -- 経過秒数
        t = (fromIntegral i) * h

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

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

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

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

        wx  = x  + ( wx1 +  wx2 * (2.0 - sqrt(2.0)) +  wx3 * (2.0 + sqrt(2.0)) +  wx4) / 6.0
        wy  = y  + ( wy1 +  wy2 * (2.0 - sqrt(2.0)) +  wy3 * (2.0 + sqrt(2.0)) +  wy4) / 6.0
        wvx = vx + (wvx1 + wvx2 * (2.0 - sqrt(2.0)) + wvx3 * (2.0 + sqrt(2.0)) + wvx4) / 6.0
        wvy = vy + (wvy1 + wvy2 * (2.0 - sqrt(2.0)) + wvy3 * (2.0 + sqrt(2.0)) + wvy4) / 6.0
    in
        if y >= 0.0
            then do
                printf "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f\n" t wvx wvy wx wy
                rungekuttagill (i+1) wvx wvy wx wy
            else
                return ()

main = do
    -- 角度
    let degree = 45.0 :: Double
    let radian = degree * pi / 180.0
    -- 初速 250 km/h -> 秒速に変換
    let v = (fromIntegral (250 * 1000 `div` 3600))

    -- 水平方向の速度
    let vx = v * cos(radian)
    -- 鉛直方向の速度
    let vy = v * sin(radian)
    -- 位置
    let x = 0.0
    let y = 0.0

    -- Runge-Kutta-Gill法
    rungekuttagill 1 vx vy x y
import Text.Printf
import Debug.Trace

f::Double->Double
f x = x * x - 2.0

bisection::Double->Double->Double
bisection a b =
    let
        -- 区間 (a, b) の中点 c = (a + b) / 2
        c = (a + b) / 2.0
        fc = (f c) :: Double
    in
        if abs(fc) < 0.0000000001 then
            c
        else
            trace (printf "%12.10f\t%13.10f" c (c - ((sqrt 2.0)::Double))) $
            if fc < 0.0
                then
                    -- f(c) < 0 であれば, 解は区間 (c, b) の中に存在
                    bisection c b
                else
                    -- f(c) > 0 であれば, 解は区間 (a, c) の中に存在
                    bisection a c

main = do
    let a = 1.0
    let b = 2.0
    printf "%12.10f\n" (bisection a b)
import Text.Printf
import Debug.Trace

f::Double->Double
f x = x * x - 2.0

falseposition::Double->Double->Double
falseposition a b =
    let
        -- 点 (a,f(a)) と 点 (b,f(b)) を結ぶ直線と x軸の交点
        c = (a * (f b) - b * (f a)) / ((f b) - (f a))
        fc = (f c) :: Double
    in
        if abs(fc) < 0.0000000001 then
            c
        else
            trace (printf "%12.10f\t%13.10f" c (c - ((sqrt 2.0)::Double))) $
            if fc < 0.0
                then
                    -- f(c) < 0 であれば, 解は区間 (c, b) の中に存在
                    falseposition c b
                else
                    -- f(c) > 0 であれば, 解は区間 (a, c) の中に存在
                    falseposition a c

main = do
    let a = 1.0
    let b = 2.0
    printf "%12.10f\n" (falseposition a b)
import Text.Printf
import Debug.Trace

g::Double->Double
g x = (x / 2.0) + (1.0 / x)

iterative::Double->Double
iterative x0 =
    let
        x1 = (g x0)
    in
        if abs(x1 - x0) < 0.0000000001 then
            x1
        else
            trace (printf "%12.10f\t%13.10f" x1 (x1 - ((sqrt 2.0)::Double)))
            iterative x1

main = do
    let x = 1.0
    printf "%12.10f\n" (iterative x)
import Text.Printf
import Debug.Trace

f0::Double->Double
f0 x = x * x - 2.0

f1::Double->Double
f1 x = 2.0 * x

newton::Double->Double
newton x0 =
    let
        x1 = x0 - ((f0 x0) / (f1 x0))
    in
        if abs(x1 - x0) < 0.0000000001 then
            x1
        else
            trace (printf "%12.10f\t%13.10f" x1 (x1 - ((sqrt 2.0)::Double)))
            newton x1

main = do
    let x = 2.0
    printf "%12.10f\n" (newton x)
import Text.Printf
import Debug.Trace

f0::Double->Double
f0 x = x * x - 2.0

f1::Double->Double
f1 x = 2.0 * x

f2::Double->Double
f2 x = 2.0

bailey::Double->Double
bailey x0 =
    let
        x1 = x0 - ((f0 x0) / ((f1 x0) - ((f0 x0) * (f2 x0) / (2.0 * (f1 x0)))))
    in
        if abs(x1 - x0) < 0.0000000001 then
            x1
        else
            trace (printf "%12.10f\t%13.10f" x1 (x1 - ((sqrt 2.0)::Double)))
            bailey x1

main = do
    let x = 2.0
    printf "%12.10f\n" (bailey x)
import Text.Printf
import Debug.Trace

f::Double->Double
f x = x * x - 2.0

secant::Double->Double->Double
secant x0 x1 =
    let
        x2 = x1 - (f x1) * (x1 - x0) / ((f x1) - (f x0))
    in
        if abs(x2 - x1) < 0.0000000001 then
            x2
        else
            trace (printf "%12.10f\t%13.10f" x2 (x2 - ((sqrt 2.0)::Double)))
            secant x1 x2

main = do
    let x0 = 1.0
    let x1 = 2.0
    printf "%12.10f\n" (secant x0 x1)
import Data.Array
import Text.Printf
import Debug.Trace
import Control.Monad

n = 3 :: Int

jacobi::Bool->[[Double]]->[Double]->[Double]->IO()
jacobi finish a b x0 = do
--    let
--        x1 = array (0,3) [(0,0.0),(1,0.0),(2,0.0),(3,0.0)]
--    in
        forM_ [0..3] $ \row -> do
            let a1 = zip (a!!row) x0
            let a2 = take row a1 ++ drop (row+1) a1
            let a3 = map (\(x, y) -> x * y) $ a2
            let s  = sum a3
            printf "%f\t" s
            let x1 = ((b!!row) - s) / a!!row!!row
            printf "%f\n" x1

            --x1 // [((fromIntegral row), ((b!!row) - s) / a!!row!!row)]
        -- printf "%d" (length $ filter (< 0) $ map (\(x,y) -> x - y) $ zip x0 (elems x1))
        --putStrLn ""

        if finish
            then putStrLn ""
            else jacobi True a b x0


disp_progress::[[Double]]->[Double]->IO()
disp_progress a b = do
    forM ([0..n]) $ \i -> do
        forM ([0..n]) $ \j -> do
            printf "%14.10f\t" (a!!i!!j)
        printf "%14.10f\t" (b!!i)
        putStrLn ""
    putStrLn ""

disp_result::[Double]->IO()
disp_result a = do
    forM ([0..n]) $ \i -> do
        printf "%14.10f\t" (a!!i)
    putStrLn ""

main = do
    let a  = [[9,2,1,1],[2,8,-2,1],[-1,-2,7,-2],[1,-1,-2,6::Double]]
    let b  = [20,16,8,17::Double]
    let x0 = [0,0,0,0::Double]
    disp_progress a b
    putStrLn ""

    jacobi False a b x0
    disp_progress a b
    putStrLn ""


--  let a1 = map (\pivot -> map(\x -> x        / (a!!pivot!!pivot)) (a!!pivot)) [0..n]
--  let b1 = map (\pivot ->           b!!pivot / (a!!pivot!!pivot)            ) [0..n]
    --disp_progress a  b1
    --disp_result b1



    putStrLn ""


--      while (true)
--      {
--          double[] x1 = new double[N];
--          bool finish = true;
--          for (int i = 0; i < N; i++)
--          {
--              x1[i] = 0;
--              for (int j = 0; j < N; j++)
--                  if (j != i)
--                    map (\x -> filter(/= x)[0..3])[0..3]
--                      x1[i] += a[i,j] * x0[j];
--
--              x1[i] = (b[i] - x1[i]) / a[i,i];
--              if (Math.Abs(x1[i] - x0[i]) > 0.0000000001) finish = false;
--          }
--
--          for (int i = 0; i < N; i++)
--          {
--              x0[i] = x1[i];
--              Console.Write(string.Format("{0,14:F10}\t", x0[i]));
--          }
--          Console.WriteLine();
--          if (finish) return;
--      }

--      while (true)
--      {
--          double[] x1 = new double[N];
--          bool finish = true;
--          for (int i = 0; i < N; i++)
--          {
--              x1[i] = 0;
--              for (int j = 0; j < N; j++)
--                  if (j != i)
--                      x1[i] += a[i,j] * x0[j];
--
--              x1[i] = (b[i] - x1[i]) / a[i,i];
--              if (Math.Abs(x1[i] - x0[i]) > 0.0000000001) finish = false;
--          }
--
--          for (int i = 0; i < N; i++)
--          {
--              x0[i] = x1[i];
--              Console.Write(string.Format("{0,14:F10}\t", x0[i]));
--          }
--          Console.WriteLine();
--          if (finish) return;
--      }
import Data.Array
import Text.Printf
import Debug.Trace
import Control.Monad

n = 3 :: Int

disp_progress::[[Double]]->[Double]->IO()
disp_progress a b = do
    forM ([0..n]) $ \i -> do
        forM ([0..n]) $ \j -> do
            printf "%14.10f\t" (a!!i!!j)
        printf "%14.10f\t" (b!!i)
        putStrLn ""
    putStrLn ""

disp_result::[Double]->IO()
disp_result a = do
    forM ([0..n]) $ \i -> do
        printf "%14.10f\t" (a!!i)
    putStrLn ""


row_loop::Int->[[Double]]->[Double]->[Double]->[Double]->[Double]
row_loop row a b x0 x1 =
    let
        a1 = zip (a!!row) x0
        a2 = take row a1 ++ drop (row+1) a1
        a3 = map (\(x, y) -> x * y) $ a2
        s  = sum a3
        x = ((b!!row) - s) / a!!row!!row
    in
        if row >= n
            then x:x1
            else (row_loop (row+1) a b x0 (x:x1))


main = do
    let a  = [[9,2,1,1],[2,8,-2,1],[-1,-2,7,-2],[1,-1,-2,6::Double]]
    let b  = [20,16,8,17::Double]
    let x0 = [0,0,0,0::Double]
    disp_progress a b
    putStrLn ""

    let x1 = row_loop 0 a b x0 []
    disp_result x1
    putStrLn ""


--      while (true)
--      {
--          double[] x1 = new double[N];
--          bool finish = true;
--          for (int i = 0; i < N; i++)
--          {
--              x1[i] = 0;
--              for (int j = 0; j < N; j++)
--                  if (j != i)
--                      x1[i] += a[i,j] * x0[j];
--
--              x1[i] = (b[i] - x1[i]) / a[i,i];
--              if (Math.Abs(x1[i] - x0[i]) > 0.0000000001) finish = false;
--          }
--
--          for (int i = 0; i < N; i++)
--          {
--              x0[i] = x1[i];
--              Console.Write(string.Format("{0,14:F10}\t", x0[i]));
--          }
--          Console.WriteLine();
--          if (finish) return;
--      }
import Text.Printf
import Debug.Trace
import Control.Monad

n = 3::Int

disp_progress::[[Double]]->[Double]->IO()
disp_progress a b = do
    forM_ ([0..n]) $ \i -> do
        forM_ ([0..n]) $ \j -> do
            printf "%14.10f\t" (a!!i!!j)
        printf "%14.10f\t" (b!!i)
        putStrLn ""
    putStrLn ""

disp_result::[Double]->IO()
disp_result a = do
    forM_ ([0..n]) $ \i -> do
        printf "%14.10f\t" (a!!i)
    putStrLn ""

row_loop::Int->[[Double]]->[Double]->[Double]->[Double]->[Double]
row_loop row a b x0 x1 =
    let
        a1 = zip (a!!row) x0
        a2 = take row a1 ++ drop (row+1) a1
        a3 = map (\(x, y) -> x * y) $ a2
        s  = sum a3
        x = ((b!!row) - s) / a!!row!!row
    in
        if row <= 0
            then
                x:x1
            else
                (row_loop (row-1) a b x0 (x:x1))

jacobi::[[Double]]->[Double]->[Double]->[[Double]]->([[Double]],[Double])
jacobi a b x0 xs =
    let
        x1  = (row_loop n a b x0 [])
        cnt = length $ filter (>= 0.0000000001) $ map (\(x,y) -> abs(x - y)) $ zip x0 x1
    in
        if cnt < 1
            then
                (reverse xs, x1)
            else
                (jacobi a b x1 (x1:xs))

main = do
    let a  = [[9,2,1,1],[2,8,-2,1],[-1,-2,7,-2],[1,-1,-2,6::Double]]
    let b  = [20,16,8,17::Double]
    let x0 = [0,0,0,0::Double]
    disp_progress a b
    putStrLn ""

    let (xs, x1) = jacobi a b x0 []

    putStrLn "mapM_"
    mapM_ (\x -> (printf "%14.10f\t" x)) $ x1
    putStrLn ""

    putStrLn "forM_1"
    forM_ x1 $ \x -> (printf "%14.10f\t" x)
    putStrLn ""

    putStrLn "forM_2"
    forM_ x1 $ \x -> do
        (printf "%14.10f\t" x)
    putStrLn ""

    putStrLn "forM_3"
    forM_ xs $ \x -> do
        forM_ x $ \x1 -> do
            (printf "%14.10f\t" x1)
        putStrLn ""
    putStrLn ""

    mapM_ (\x -> (printf "%14.10f\t" x)) $ x1
import Text.Printf
import Debug.Trace
import Control.Monad

n = 3 :: Int

disp_progress::[[Double]]->IO()
disp_progress a = do
    forM ([0..n]) $ \i -> do
        forM ([0..n]) $ \j -> do
            printf "%14.10f\t" (a!!i!!j)
        putStrLn ""
    putStrLn ""
    --return ()

jacobi::Bool->[[Double]]->[Double]->([[Double]], [Double])
jacobi finish a b = -- do
    if finish then
        (a, b)
    else
        --trace (printf "%14.10f", 1.0 ) $
        jacobi True a b

jacobi2::Bool->[[Double]]->[Double]->IO()
jacobi2 finish a b = -- do
    if finish then
        putStrLn ""
    else
        jacobi2 True a b

test::(Show a)=>a->a
test x =
        trace (show x) $
        x

main_loop::Bool->[[Double]]->[Double]->IO()
main_loop finish a b = -- do
    if finish then
        putStrLn ""
    else
--      double[] x1 = new double[N];
--      bool finish = true;
--      for (int i = 0; i < N; i++)
--      {
--          x1[i] = 0;
--          for (int j = 0; j < N; j++)
--              if (j != i)
--                  x1[i] += a[i,j] * x0[j];
--
--          x1[i] = (b[i] - x1[i]) / a[i,i];
--          if (Math.Abs(x1[i] - x0[i]) > 0.0000000001) finish = false;
--      }
--
--      for (int i = 0; i < N; i++)
--      {
--          x0[i] = x1[i];
--          Console.Write(string.Format("{0,14:F10}\t", x0[i]));
--      }
        main_loop True a b

--main :: IO ()
main = do
    let a = [[9,2,1,1],[2,8,-2,1],[-1,-2,7,-2],[1,-1,-2,6::Double]]
    let b = [20,16,8,17]
    disp_progress a

    let (x, y) = jacobi False a b
    printf "%14.10f\n" (x!!0!!0)

    main_loop False a b

--      while (true)
--      {
--          double[] x1 = new double[N];
--          bool finish = true;
--          for (int i = 0; i < N; i++)
--          {
--              x1[i] = 0;
--              for (int j = 0; j < N; j++)
--                  if (j != i)
--                      x1[i] += a[i,j] * x0[j];
--
--              x1[i] = (b[i] - x1[i]) / a[i,i];
--              if (Math.Abs(x1[i] - x0[i]) > 0.0000000001) finish = false;
--          }
--
--          for (int i = 0; i < N; i++)
--          {
--              x0[i] = x1[i];
--              Console.Write(string.Format("{0,14:F10}\t", x0[i]));
--          }
--          Console.WriteLine();
--          if (finish) return;
--      }

    forM ([0..n]) $ \i -> do
        forM ([0..n]) $ \j -> do
            printf "%14.10f\t" (test (a!!i!!j))
        putStrLn ""
    putStrLn ""

  --print $ id' (8::Int)
  --print $ sum' ([1,2,3]::[Int])
  --print $ append "abc" "def"
  --print $ log2 10
  --print $ getHeight tree1
  --return ()
import Text.Printf
import Control.Monad

n = 4::Int

-- 2次元配列を表示
disp_matrix::[[Double]]->IO()
disp_matrix matrix = do
    forM_ matrix $ \row -> do
        forM_ row $ \elem -> do
            printf "%14.10f\t" elem
        putStrLn ""

-- 1次元配列を表示
disp_vector::[Double]->IO()
disp_vector vector = do
    forM_ vector $ \elem -> do
        printf "%14.10f\t" elem
    putStrLn ""

-- ヤコビの反復法
row_loop::Int->[[Double]]->[Double]->[Double]->[Double]->[Double]
row_loop row a b x0 x1 =
    let
        a1 = zip (a!!row) x0
        a2 = take row a1 ++ drop (row + 1) a1
        a3 = map (\(x, y) -> x * y) $ a2
        s  = sum a3
        x  = ((b!!row) - s) / a!!row!!row
        xs = x:x1
    in
        if row >= (n - 1)
            then
                reverse xs
            else
                (row_loop (row+1) a b x0 xs)

jacobi::[[Double]]->[Double]->[Double]->[[Double]]->([[Double]],[Double])
jacobi a b x0 xs =
    let
        x1  = (row_loop 0 a b x0 [])
        cnt = length $ filter (>= 0.0000000001) $ map (\(x,y) -> abs(x - y)) $ zip x0 x1
    in
        if cnt < 1
            then
                (reverse xs, x1)
            else
                (jacobi a b x1 (x1:xs))

main = do
    let a  = [[9,2,1,1],[2,8,-2,1],[-1,-2,7,-2],[1,-1,-2,6::Double]]
    let b  = [20,16,8,17::Double]
    let x0 = [0,0,0,0::Double]

    putStrLn "A"
    disp_matrix a
    putStrLn "B"
    disp_vector b
    putStrLn ""

    -- ヤコビの反復法
    let (xs, x1) = jacobi a b x0 []
    disp_matrix xs

    putStrLn "X"
    disp_vector x1
import Text.Printf
import Control.Monad

n = 4::Int

-- 2次元配列を表示
disp_matrix::[[Double]]->IO()
disp_matrix matrix = do
    forM_ matrix $ \row -> do
        forM_ row $ \elem -> do
            printf "%14.10f\t" elem
        putStrLn ""

-- 1次元配列を表示
disp_vector::[Double]->IO()
disp_vector vector = do
    forM_ vector $ \elem -> do
        printf "%14.10f\t" elem
    putStrLn ""

-- ガウス・ザイデル法
row_loop::Int->[[Double]]->[Double]->[Double]->[Double]->[Double]
row_loop row a b x0 x1 =
    let
        a1 = zip (a!!row) x0
        a2 = take row a1 ++ drop (row + 1) a1
        a3 = map (\(x, y) -> x * y) $ a2
        s  = sum a3
        x  = ((b!!row) - s) / a!!row!!row
        xs = x:x1
    in
        if row >= (n - 1)
            then
                reverse xs
            else
                (row_loop (row+1) a b ((take (row+1) (reverse xs)) ++ (drop (row+1) x0)) xs)

gauss::[[Double]]->[Double]->[Double]->[[Double]]->([[Double]],[Double])
gauss a b x0 xs =
    let
        x1  = (row_loop 0 a b x0 [])
        cnt = length $ filter (>= 0.0000000001) $ map (\(x,y) -> abs(x - y)) $ zip x0 x1
    in
        if cnt < 1
            then
                (reverse xs, x1)
            else
                (gauss a b x1 (x1:xs))

main = do
    let a  = [[9,2,1,1],[2,8,-2,1],[-1,-2,7,-2],[1,-1,-2,6::Double]]
    let b  = [20,16,8,17::Double]
    let x0 = [0,0,0,0::Double]

    putStrLn "A"
    disp_matrix a
    putStrLn "B"
    disp_vector b
    putStrLn ""

    -- ガウス・ザイデル法
    let (xs, x1) = gauss a b x0 []
    disp_matrix xs

    putStrLn "X"
    disp_vector x1
import Text.Printf
import Debug.Trace
import Control.Monad

n = 4::Int

-- 2次元配列を表示
disp_matrix::[[Double]]->IO()
disp_matrix matrix = do
    forM_ matrix $ \row -> do
        forM_ row $ \elem -> do
            printf "%14.10f\t" elem
        putStrLn ""

-- 1次元配列を表示
disp_vector::[Double]->IO()
disp_vector vector = do
    forM_ vector $ \elem -> do
        printf "%14.10f\t" elem
    putStrLn ""

-- 各列で 一番値が大きい行を 探す
get_max_row::Int->Int->[[Double]]->Int->Double->(Int, Double)
get_max_row row col a max_row max_val =
    let
        -- 一番値が大きい行
        (max_row2, max_val2) =
            if max_val < abs(a!!row!!col)
                then (row,     abs(a!!row!!col))
                else (max_row, max_val)
    in
        if row >= length(a) - 1
            then
                (max_row2, max_val2)
            else
                get_max_row (row+1) col a max_row2 max_val2

-- ピボット選択
pivoting::Int->[[Double]]->[Double]->[[Double]]->[Double]->([[Double]],[Double])
pivoting pivot a b a2 b2 =
    let
        (max_row, max_val) = get_max_row 0 pivot a 0 0.0
        a3 = (a!!max_row):a2
        b3 = (b!!max_row):b2
        a4 = (take (max_row) a) ++ (drop (max_row+1) a)
        b4 = (take (max_row) b) ++ (drop (max_row+1) b)
    in
        if pivot >= (n - 1)
            then
                (reverse a3, reverse b3)
            else
                pivoting (pivot+1) a4 b4 a3 b3

-- 前進消去
forward_elim_loop::Int->Int->[[Double]]->[Double]->([[Double]],[Double])
forward_elim_loop pivot row a b =
    let
        s = a!!row!!pivot / a!!pivot!!pivot

        a_zip   = zip (drop pivot (a!!pivot)) (drop pivot (a!!row))
        a_map   = map (\(a_pivot, a_row) -> a_row - a_pivot * s) a_zip
        new_row = (take pivot (a!!row)) ++ a_map
        a2      = (take row   a       ) ++ (new_row:(drop (row+1) a))

        x  = b!!row - b!!pivot * s
        b2 = (take row b) ++ (x:(drop (row+1) b))
    in
        if row < (n - 1)
            then
                forward_elim_loop pivot (row+1) a2 b2
            else
                (a2, b2)

forward_elimination::Int->[[Double]]->[Double]->([[Double]],[Double])
forward_elimination pivot a b =
    let
        (a2, b2) =
            if pivot < (n - 1)
                then forward_elim_loop pivot (pivot+1) a b
                else (a, b)
    in
        if pivot < (n - 1)
            then
                forward_elimination (pivot+1) a2 b2
            else
                (a2, b2)

-- 後退代入
backward_substitution::Int->[[Double]]->[Double]->[Double]
backward_substitution row a b =
    let
        a_zip = zip (drop (row+1) (a!!row)) (drop (row+1) b)
        s  = sum $ map(\(a_col, b_col) -> a_col * b_col) a_zip
        x  = (b!!row - s) / (a!!row!!row)
        b2 = (take row b) ++ (x:(drop (row+1) b))
    in
        if row > 0
            then x:(backward_substitution (row-1) a b2)
            else x:[]

main = do
    let a  = [[-1,-2,7,-2],[1,-1,-2,6],[9,2,1,1],[2,8,-2,1::Double]]
    let b  = [8,17,20,16::Double]

    -- ピボット選択
    let (a1, b1) = pivoting 0 a b [] []

    putStrLn "Pivoting"
    putStrLn "A"
    disp_matrix a1
    putStrLn "B"
    disp_vector b1
    putStrLn ""

    -- 前進消去
    let (a2, b2) = forward_elimination 0 a1 b1

    putStrLn "Forward Elimination"
    putStrLn "A"
    disp_matrix a2
    putStrLn "B"
    disp_vector b2
    putStrLn ""

    -- 後退代入
    let x = backward_substitution (n-1) a2 b2

    putStrLn "X"
    disp_vector (reverse x)
import Text.Printf
import Debug.Trace
import Control.Monad

n = 4::Int

-- 2次元配列を表示
disp_matrix::[[Double]]->IO()
disp_matrix matrix = do
    forM_ matrix $ \row -> do
        forM_ row $ \elem -> do
            printf "%14.10f\t" elem
        putStrLn ""

-- 1次元配列を表示
disp_vector::[Double]->IO()
disp_vector vector = do
    forM_ vector $ \elem -> do
        printf "%14.10f\t" elem
    putStrLn ""

-- 各列で 一番値が大きい行を 探す
get_max_row::Int->Int->[[Double]]->Int->Double->(Int, Double)
get_max_row row col a max_row max_val =
    let
        -- 一番値が大きい行
        (max_row2, max_val2) =
            if max_val < abs(a!!row!!col)
                then (row,     abs(a!!row!!col))
                else (max_row, max_val)
    in
        if row >= length(a) - 1
            then
                (max_row2, max_val2)
            else
                get_max_row (row+1) col a max_row2 max_val2

-- ピボット選択
pivoting::Int->[[Double]]->[Double]->[[Double]]->[Double]->([[Double]],[Double])
pivoting pivot a b a2 b2 =
    let
        (max_row, max_val) = get_max_row 0 pivot a 0 0.0
        a3 = (a!!max_row):a2
        b3 = (b!!max_row):b2
        a4 = (take (max_row) a) ++ (drop (max_row+1) a)
        b4 = (take (max_row) b) ++ (drop (max_row+1) b)
    in
        if pivot >= (n - 1)
            then
                (reverse a3, reverse b3)
            else
                pivoting (pivot+1) a4 b4 a3 b3

-- 前進消去
forward_elim_loop::Int->Int->[[Double]]->[Double]->([[Double]],[Double])
forward_elim_loop pivot row a b =
    let
        s = a!!row!!pivot / a!!pivot!!pivot

        a_zip   = zip (a!!pivot) (a!!row)
        new_row = map (\(a_pivot, a_row) -> a_row - a_pivot * s) a_zip
        a2      = (take row a) ++ (new_row:(drop (row+1) a))

        x  = b!!row - b!!pivot * s
        b2 = (take row b) ++ (x:(drop (row+1) b))

        next_row =
            if row + 1 == pivot
                then row + 2
                else row + 1
    in
        if next_row < n
            then
                forward_elim_loop pivot next_row a2 b2
            else
                (a2, b2)

forward_elimination::Int->[[Double]]->[Double]->([[Double]],[Double])
forward_elimination pivot a b =
    let
        row = if pivot == 0 then 1 else 0
        (a2, b2) =
            if pivot < n
                then forward_elim_loop pivot row a b
                else (a, b)
    in
        if pivot < (n - 1)
            then
                forward_elimination (pivot+1) a2 b2
            else
                (a2, b2)

-- 後退代入
backward_substitution::Int->[[Double]]->[Double]->[Double]
backward_substitution row a b =
    let
        x  = b!!row / (a!!row!!row)
        b2 = (take row b) ++ (x:(drop (row+1) b))
    in
        if row > 0
            then x:(backward_substitution (row-1) a b2)
            else x:[]

main = do
    let a  = [[-1,-2,7,-2],[1,-1,-2,6],[9,2,1,1],[2,8,-2,1::Double]]
    let b  = [8,17,20,16::Double]

    -- ピボット選択
    let (a1, b1) = pivoting 0 a b [] []

    putStrLn "Pivoting"
    putStrLn "A"
    disp_matrix a1
    putStrLn "B"
    disp_vector b1
    putStrLn ""

    -- 前進消去
    let (a2, b2) = forward_elimination 0 a1 b1

    putStrLn "Forward Elimination"
    putStrLn "A"
    disp_matrix a2
    putStrLn "B"
    disp_vector b2
    putStrLn ""

    -- 後退代入
    let x = backward_substitution (n-1) a2 b2

    putStrLn "X"
    disp_vector (reverse x)
import Text.Printf
import Debug.Trace
import Control.Monad

n = 4::Int

-- 2次元配列を表示
disp_matrix::[[Double]]->IO()
disp_matrix matrix = do
    forM_ matrix $ \row -> do
        forM_ row $ \elem -> do
            printf "%14.10f\t" elem
        putStrLn ""

-- 1次元配列を表示
disp_vector::[Double]->IO()
disp_vector vector = do
    forM_ vector $ \elem -> do
        printf "%14.10f\t" elem
    putStrLn ""

-- 各列で 一番値が大きい行を 探す
get_max_row::Int->Int->[[Double]]->Int->Double->(Int, Double)
get_max_row row col a max_row max_val =
    let
        -- 一番値が大きい行
        (max_row2, max_val2) =
            if max_val < abs(a!!row!!col)
                then (row,     abs(a!!row!!col))
                else (max_row, max_val)
    in
        if row >= length(a) - 1
            then
                (max_row2, max_val2)
            else
                get_max_row (row+1) col a max_row2 max_val2

-- ピボット選択
pivoting::Int->[[Double]]->[Double]->[[Double]]->[Double]->([[Double]],[Double])
pivoting pivot a b a2 b2 =
    let
        (max_row, max_val) = get_max_row 0 pivot a 0 0.0
        a3 = (a!!max_row):a2
        b3 = (b!!max_row):b2
        a4 = (take (max_row) a) ++ (drop (max_row+1) a)
        b4 = (take (max_row) b) ++ (drop (max_row+1) b)
    in
        if pivot >= (n - 1)
            then
                (reverse a3, reverse b3)
            else
                pivoting (pivot+1) a4 b4 a3 b3

-- 前進消去
forward_elim_loop::Int->Int->[[Double]]->[Double]->([[Double]],[Double])
forward_elim_loop pivot row a b =
    let
        s = a!!row!!pivot / a!!pivot!!pivot

        a_zip   = zip (drop (pivot+1) (a!!pivot)) (drop (pivot+1) (a!!row))
        a_map   = map (\(a_pivot, a_row) -> a_row - a_pivot * s) a_zip
        new_row = (take pivot (a!!row)) ++ (s:a_map                   )
        a2      = (take row   a       ) ++ (new_row:(drop (row+1)   a))
    in
        if row < (n - 1)
            then
                forward_elim_loop pivot (row+1) a2 b
            else
                (a2, b)

forward_elimination::Int->[[Double]]->[Double]->([[Double]],[Double])
forward_elimination pivot a b =
    let
        (a2, b2) =
            if pivot < (n - 1)
                then forward_elim_loop pivot (pivot+1) a b
                else (a, b)
    in
        if pivot < (n - 1)
            then
                forward_elimination (pivot+1) a2 b2
            else
                (a2, b2)

-- Ly=b から y を求める (前進代入)
forward_substitution::Int->[[Double]]->[Double]->[Double]
forward_substitution row a b =
    let
        a_zip = zip (take row (a!!row)) b
        s  = sum $ map(\(a_col, b_col) -> a_col * b_col) a_zip
        y  = b!!row - s
        b2 = (take row b) ++ (y:(drop (row+1) b))
    in
        if row < (n - 1)
            then y:(forward_substitution (row+1) a b2)
            else y:[]

-- Ux=y から x を求める (後退代入)
backward_substitution::Int->[[Double]]->[Double]->[Double]
backward_substitution row a b =
    let
        a_zip = zip (drop (row+1) (a!!row)) (drop (row+1) b)
        s  = sum $ map(\(a_col, b_col) -> a_col * b_col) a_zip
        x  = (b!!row - s) / (a!!row!!row)
        b2 = (take row b) ++ (x:(drop (row+1) b))
    in
        if row > 0
            then x:(backward_substitution (row-1) a b2)
            else x:[]

main = do
    let a  = [[-1,-2,7,-2],[1,-1,-2,6],[9,2,1,1],[2,8,-2,1::Double]]
    let b  = [8,17,20,16::Double]

    -- ピボット選択
    let (a1, b1) = pivoting 0 a b [] []

    putStrLn "Pivoting"
    putStrLn "A"
    disp_matrix a1
    putStrLn "B"
    disp_vector b1
    putStrLn ""

    -- 前進消去
    let (a2, b2) = forward_elimination 0 a1 b1

    putStrLn "LU"
    disp_matrix a2

    -- Ly=b から y を求める (前進代入)
    let y = forward_substitution 0 a2 b2

    putStrLn "Y"
    disp_vector (reverse y)

    -- Ux=y から x を求める (後退代入)
    let x = backward_substitution (n-1) a2 y

    putStrLn "X"
    disp_vector (reverse x)
import Text.Printf
import Debug.Trace
import Control.Monad

n = 4::Int

-- 2次元配列を表示
disp_matrix::[[Double]]->IO()
disp_matrix matrix = do
    forM_ matrix $ \row -> do
        forM_ row $ \elem -> do
            printf "%14.10f\t" elem
        putStrLn ""

-- 1次元配列を表示
disp_vector::[Double]->IO()
disp_vector vector = do
    forM_ vector $ \elem -> do
        printf "%14.10f\t" elem
    putStrLn ""

-- 前進消去
forward_elim_loop::Int->Int->[[Double]]->Double->[[Double]]
forward_elim_loop pivot row a a_sqrt =
    let
        a_zip = zip (take pivot (a!!row)) (take pivot (a!!pivot))
        s     = sum $ map(\(a_row, a_pivot) -> a_row * a_pivot) a_zip
        a_col = (a!!row!!pivot - s) / a_sqrt

        a_row = (take pivot (a!!row)  ) ++ (a_col:(drop (pivot+1) (a!!row)))
        a_piv = (take row   (a!!pivot)) ++ (a_col:(drop (row+1)   (a!!pivot)))

        a1 = (take pivot a ) ++ (a_piv:(drop (pivot+1) a))
        a2 = (take row   a1) ++ (a_row:(drop (row+1)   a))
    in
        if row < (n - 1)
            then forward_elim_loop pivot (row+1) a2 a_sqrt
            else a2

forward_elimination::Int->[[Double]]->[[Double]]
forward_elimination pivot a =
    let
        s      = sum $ map(\a_col -> a_col * a_col) $ take pivot (a!!pivot)
        -- ここで根号の中が負の値になると計算できない!
        a_sqrt = (sqrt (a!!pivot!!pivot - s))
        a2     = forward_elim_loop pivot pivot a a_sqrt
    in
        if pivot < (n - 1)
            then
                forward_elimination (pivot+1) a2
            else
                a2

-- Ly=b から y を求める (前進代入)
forward_substitution::Int->[[Double]]->[Double]->[Double]
forward_substitution row a b =
    let
        a_zip = zip (take row (a!!row)) b
        s  = sum $ map(\(a_col, b_col) -> a_col * b_col) a_zip
        y  = (b!!row - s) / a!!row!!row
        b2 = (take row b) ++ (y:(drop (row+1) b))
    in
        if row < (n - 1)
            then y:(forward_substitution (row+1) a b2)
            else y:[]

-- Ux=y から x を求める (後退代入)
backward_substitution::Int->[[Double]]->[Double]->[Double]
backward_substitution row a b =
    let
        a_zip = zip (drop (row+1) (a!!row)) (drop (row+1) b)
        s  = sum $ map(\(a_col, b_col) -> a_col * b_col) a_zip
        x  = (b!!row - s) / (a!!row!!row)
        b2 = (take row b) ++ (x:(drop (row+1) b))
    in
        if row > 0
            then x:(backward_substitution (row-1) a b2)
            else x:[]

main = do
    let a  = [[5,2,3,4],[2,10,6,7],[3,6,15,9],[4,7,9,20::Double]]
    let b  = [34,68,96,125::Double]

    putStrLn "A"
    disp_matrix a
    putStrLn "B"
    disp_vector b
    putStrLn ""

    -- 前進消去
    let a2 = forward_elimination 0 a

    putStrLn "LL^T"
    disp_matrix a2

    -- Ly=b から y を求める (前進代入)
    let y = forward_substitution 0 a2 b

    putStrLn "Y"
    disp_vector y

    -- Ux=y から x を求める (後退代入)
    let x = backward_substitution (n-1) a2 y

    putStrLn "X"
    disp_vector (reverse x)
import Text.Printf
import Debug.Trace
import Control.Monad

n = 4::Int

-- 2次元配列を表示
disp_matrix::[[Double]]->IO()
disp_matrix matrix = do
    forM_ matrix $ \row -> do
        forM_ row $ \elem -> do
            printf "%14.10f\t" elem
        putStrLn ""

-- 1次元配列を表示
disp_vector::[Double]->IO()
disp_vector vector = do
    forM_ vector $ \elem -> do
        printf "%14.10f\t" elem
    putStrLn ""

-- 前進消去
forward_elim_loop::Int->Int->[[Double]]->[[Double]]
forward_elim_loop pivot col a =
    let
        a_zip = zip [0..col] $ zip (take col (a!!pivot)) (take col (a!!col))
        s     = sum $ map(\(k, (a_pivot, a_col)) -> a_pivot * a_col * (a!!k!!k)) a_zip
        a_col = (a!!pivot!!col - s) / (a!!col!!col)

        a_row = (take col   (a!!pivot)) ++ (a_col:(drop (col+1)   (a!!pivot)))
        a_piv = (take pivot (a!!col  )) ++ (a_col:(drop (pivot+1) (a!!col)))

        a1 = (take col   a ) ++ (a_piv:(drop (col+1)   a))
        a2 = (take pivot a1) ++ (a_row:(drop (pivot+1) a))

        next_col = if col + 1 == pivot then pivot + 1
                                       else col   + 1
    in
        if next_col < pivot then forward_elim_loop pivot next_col a2
                            else a2

forward_elimination::Int->[[Double]]->[[Double]]
forward_elimination pivot a =
    let
        next_col = if pivot == 0 then 1
                                 else 0
        a1    = forward_elim_loop pivot next_col a
        s     = sum $ map(\(k, a_col) -> a_col * a_col * a1!!k!!k) $ zip [0..pivot] $ take pivot (a1!!pivot)
        a_col = a1!!pivot!!pivot - s

        a_row = (take pivot (a1!!pivot)) ++ (a_col:(drop (pivot+1) (a1!!pivot)))
        a2    = (take pivot a1         ) ++ (a_row:(drop (pivot+1) a1         ))
    in
        if pivot < (n - 1)
            then forward_elimination (pivot+1) a2
            else a2

-- Ly=b から y を求める (前進代入)
forward_substitution::Int->[[Double]]->[Double]->[Double]
forward_substitution row a b =
    let
        a_zip = zip (take row (a!!row)) b
        s  = sum $ map(\(a_col, b_col) -> a_col * b_col) a_zip
        y  = (b!!row - s)
        b2 = (take row b) ++ (y:(drop (row+1) b))
    in
        if row < (n - 1)
            then y:(forward_substitution (row+1) a b2)
            else y:[]

-- Ux=y から x を求める (後退代入)
backward_substitution::Int->[[Double]]->[Double]->[Double]
backward_substitution row a b =
    let
        a_zip = zip (drop (row+1) (a!!row)) (drop (row+1) b)
        s  = sum $ map(\(a_col, b_col) -> a_col * a!!row!!row * b_col) a_zip
        x  = (b!!row - s) / (a!!row!!row)
        b2 = (take row b) ++ (x:(drop (row+1) b))
    in
        if row > 0
            then x:(backward_substitution (row-1) a b2)
            else x:[]

main = do
    let a  = [[5,2,3,4],[2,10,6,7],[3,6,15,9],[4,7,9,20::Double]]
    let b  = [34,68,96,125::Double]

    putStrLn "A"
    disp_matrix a
    putStrLn "B"
    disp_vector b
    putStrLn ""

    -- 前進消去
    let a2 = forward_elimination 0 a

    putStrLn "LDL^T"
    disp_matrix a2

    -- Ly=b から y を求める (前進代入)
    let y = forward_substitution 0 a2 b

    putStrLn "Y"
    disp_vector y

    -- Ux=y から x を求める (後退代入)
    let x = backward_substitution (n-1) a2 y

    putStrLn "X"
    disp_vector (reverse x)
inserted by FC2 system