(printf "%d\n" (+ 3 5))
(printf "%d\n" (- 3 5))
(printf "%d\n" (* 3 5))
(printf "%f\n" (Math/pow 3 5))
(printf "%f\n" (double (/ 5 3)))
(printf "%f\n" (/ 5.0 3))
(printf "%f\n" (/ 5 3.0))
(printf "%d\n" (quot 5 3))
(printf "%d\n" (rem 5 3))

(print (format "%d\n" (* 3 5)))

(println (format "%3d" (* 3 5)))
(println (format "%23.20f" (/ 5.0 3)))
(def i (* 3 5))
(printf "3 * 5 = %d\n" i)
(loop [n 1]
    (if (< n 10)
        (do
            (printf "%d, " n)
            (recur (inc n)))
        (println "")))
(loop [n 1]
    (if (< n 10)
        (if (= (mod n 3) 0)
            (do
                (printf "%d, " n)
                (recur (inc n)))
            (recur (inc n)))
        (println "")))
(defn total [n]
    (if (zero? n) 0
        (if (= (mod n 3) 0)
            (+' n (total (dec n)))
            (total (dec n)))))

(println (total 99))
(range 1 10)
(filter #(= (mod  % 3) 0)
    (range 1 10))
(filter #(= (rem  % 3) 0)
    (range 1 10))
(apply +
    (filter #(= (mod  % 3) 0)
    (range 1 100)))
(reduce +
    (filter #(= (mod  % 3) 0)
    (range 1 100)))
; 初項:a, 公差:a で, 上限:lim の数列の総和を返す関数
(defn sn [a, lim]
    (def n (quot lim a))    ; 項数:n  =  上限:lim / 公差:a
    (def l (* a   n))       ; 末項:l  =  項数:n   * 公差:a
    (quot (* (+ a l) n) 2)) ; 総和:sn = (初項:a   + 末項:l) * 項数:n / 2
; 3 の倍数の合計
(sn 3 999)
(def n 10000)
(quot (* n (+ n 1)) 2)
(def n 5000)
(* n (+ n 1))
(def n 5000)
(int (Math/pow n 2))
(def n 1000)
(quot (* (* n (+ n 1)) (+ (* 2 n) 1)) 6)
(def n 100)
(int (quot (* (Math/pow n 2) (Math/pow (+ n 1) 2)) 4))
(def n 10)
(def a 2)
(def r 3)
(int (quot (* a (- (Math/pow r n) 1)) (- r 1)))
(range 0 10)
(map
    #(+ 5 (* 3 %))
    (range 0 10))
(apply *
    (map
        #(+ 5 (* 3 %))
        (range 0 10)))
; 等差数列の積
(defn product [m d n]
    (if (zero? n)
        1
        (* m (product (+ m d) d (- n 1)))))
; 初項 5, 公差 3, 項数 10 の数列の積
(product 5 3 10)
; 階乗を求める関数
(defn Fact [n]
    (if (zero? n)
        1
        (* n (Fact (- n 1)))))
; 10の階乗
(Fact 10)
; 下降階乗冪
(defn FallingFact [x n]
    (if (<= n 1)
        x
        (* x (FallingFact (- x 1) (- n 1)))))
; 10 から 6 までの 総乗
(FallingFact 10 5)
; 上昇階乗冪
(defn RisingFact [x n]
    (if (<= n 1)
        x
        (* x (RisingFact (+ x 1) (- n 1)))))
; 10 から 14 までの 総乗
(RisingFact 10 5)
; 階乗
(defn Fact [n]
    (if (zero? n)
        1
        (* n (Fact (- n 1)))))
; 順列 (異なる 10 個のものから 5 個取ってできる順列の総数)
(def n 10)
(def r 5)
(quot (Fact n) (Fact (- n r)))
; 下降階乗冪
(defn FallingFact [x n]
    (if (<= n 1)
        x
        (* x (FallingFact (- x 1) (- n 1)))))
; 順列 (異なる 10 個のものから 5 個取ってできる順列の総数)
(def n 10)
(def r 5)
(FallingFact n r)
; 重複順列 (異なる 10 個のものから重複を許して 5 個取ってできる順列の総数)
(def n 10)
(def r 5)
(Math/pow n r)
; 組合せ
(defn Comb [n r]
    (cond (or (= r 0) (= r n)) 1
          (= r 1)              n
          true                 (+ (Comb (- n 1) (- r 1)) (Comb (- n 1) r))))
; 異なる 10 個のものから 5 個取ってできる組合せの総数
(def n 10)
(def r 5)
(Comb n r)
; 組合せ
(defn Comb [n r]
    (cond (or (= r 0) (= r n)) 1
          (= r 1)              n
          true                 (+ (Comb (- n 1) (- r 1)) (Comb (- n 1) r))))
; 異なる 10 個のものから重複を許して 5 個取る組合せの総数
(def n 10)
(def r 5)
(Comb (+ n (- r 1)) r)
;自作の正弦関数
(defn mySin [x n nega numerator denominator y]
    (def m     (* 2 n))
    (def denom (* denominator (* (+ m 1) m)))
    (def nume  (* numerator (* x x)))
    (def a     (/ nume denom))
    ;十分な精度になったら処理を抜ける
    (if (<= a 0.00000000001)
        y
        (+ y (mySin x (+ n 1) (not nega) nume denom (if nega a (- a))))))
(doseq
    [degree (filter #(or (= (mod % 30) 0)
                         (= (mod % 45) 0))
            (map    #(* % 15)
            (range 0 25)))]
    (do
        (def radian (/ (* degree (. Math PI)) 180.0))
        ;自作の正弦関数
        (def d1     (mySin radian 1 false radian 1.0 radian))
        ;標準の正弦関数
        (def d2     (. Math sin radian))
        ;標準関数との差異
        (println (format "%3d : %13.10f - %13.10f = %13.10f" degree d1 d2 (- d1 d2)))))
;自作の余弦関数
(defn myCos [x n nega numerator denominator y]
    (def m     (* 2 n))
    (def denom (* denominator (* m (- m 1))))
    (def nume  (* numerator (* x x)))
    (def a     (/ nume denom))
    ;十分な精度になったら処理を抜ける
    (if (<= a 0.00000000001)
        y
        (+ y (myCos x (+ n 1) (not nega) nume denom (if nega a (- a))))))
(doseq
    [degree (filter #(or (= (mod % 30) 0)
                         (= (mod % 45) 0))
            (map    #(* % 15)
            (range 0 25)))]
    (do
        (def radian (/ (* degree (. Math PI)) 180.0))
        ;自作の余弦関数
        (def d1     (myCos radian 1 false 1.0 1.0 1.0))
        ;標準の余弦関数
        (def d2     (. Math cos radian))
        ;標準関数との差異
        (println (format "%3d : %13.10f - %13.10f = %13.10f" degree d1 d2 (- d1 d2)))))
;自作の正接関数
(defn myTan [x x2 n t]
    (def denom (/ x2 (- n t)))
    (def nume  (- n 2))
    (if (<= nume 1)
        (/ x (- 1 denom))
        (myTan x x2 nume denom)))
(doseq
    [degree (map    #(- % 90)
            (filter #(not (zero? (mod % 180)))
            (map    #(* % 15)
            (range 0 12))))]
    (do (def radian (/ (* degree (. Math PI)) 180.0))
        (def x2     (* radian radian))
        ;自作の正接関数
        (def d1     (myTan radian x2 15 0))
        ;標準の正接関数
        (def d2     (. Math tan radian))
        ;標準関数との差異
        (println (format "%3d : %13.10f - %13.10f = %13.10f" degree d1 d2 (- d1 d2)))))
;自作の指数関数
(defn myExp [x n numerator denominator y]
    (def denom (* denominator n))
    (def nume  (* numerator x))
    (def a     (/ nume denom))
    ;十分な精度になったら処理を抜ける
    (if (<= (. Math abs a) 0.00000000001)
        y
        (+ y (myExp x (+ n 1) nume denom a))))
(doseq
    [x  (map #(/ (- % 10) 4.0)
        (range 0 21))]
    (do
        ;標準の指数関数
        (def d1 (. Math exp x))
        ;自作の指数関数
        (def d2 (myExp x 1 1.0 1.0 1.0))
        ;標準関数との差異
        (println (format "%5.2f : %13.10f - %13.10f = %13.10f" x d1 d2 (- d1 d2)))))
;自作の指数関数
(defn myExp [x x2 n t]
    (def denom (/ x2 (+ n t)))
    (def nume  (- n 4))
    (if (< nume 6)
        (+ 1.0 (/ (* 2.0 x) (+ (- 2.0 x) denom)))
        (myExp x x2 nume denom)))
(doseq
    [x  (map #(/ (- % 10) 4.0)
        (range 0 20))]
    (do
        ;標準の指数関数
        (def d1 (. Math exp x))
        ;自作の指数関数
        (def x2 (* x x))
        ;標準関数との差異
        (def d2 (myExp x x2 30 0.0)) ;30:必要な精度が得られるよう, 6から始めて4ずつ増加させる
        (println (format "%5.2f : %13.10f - %13.10f = %13.10f" x d1 d2 (- d1 d2)))))
;自作の対数関数
(defn myLog [x numerator denominator y]
    (def denom (+ denominator 2.0))
    (def nume  (* numerator (* x2 x2)))
    (def a     (/ nume denom))
    ;十分な精度になったら処理を抜ける
    (if (<= (. Math abs a) 0.00000000001)
        y
        (+ y (myLog x nume denom a))))
(doseq
    [x  (map #(/ % 5.0)
        (range 1 21))]
    (do
        ;標準の対数関数
        (def d1 (. Math log x))
        ;自作の対数関数
        (def x2 (/ (- x 1.0) (+ x 1.0)))
        ;標準関数との差異
        (def d2 (* 2.0 (myLog x2 x2 1.0 x2)))
        (println (format "%5.2f : %13.10f - %13.10f = %13.10f" x d1 d2 (- d1 d2)))))
;自作の対数関数
(defn myLog [x n t]
    (do
        (def n2 (if (zero? (rem n 2)) 2
                                      n))
        (def n3 (if (> n 3)           n2
                                      n))
        (def x2 (if (> n 3)           (* x (quot n 2))
                                      x))
        (def t2 (/ x2 (+ n3 t)))

        (if (<= n 2) (/ x (+ 1.0 t2))
                     (myLog x (- n 1) t2))))
(doseq
    [x  (map #(/ % 5.0)
        (range 1 21))]
    (do
        ;標準の対数関数
        (def d1 (. Math log x))
        ;自作の対数関数
        (def d2 (myLog (- x 1.0) 27 0.0)) ;27:必要な精度が得られる十分大きな奇数
        ;標準関数との差異
        (println (format "%5.2f : %13.10f - %13.10f = %13.10f" x d1 d2 (- d1 d2)))))
;自作の双曲線正弦関数
(defn mySinh [x n numerator denominator y]
    (def m     (* 2 n))
    (def denom (* denominator (* (+ m 1) m)))
    (def nume  (* numerator (* x x)))
    (def a     (/ nume denom))
    ;十分な精度になったら処理を抜ける
    (if (<= (. Math abs a) 0.00000000001)
        y
        (+ y (mySinh x (+ n 1) nume denom a))))
(doseq
    [i  (map #(- % 10)
        (range 0 21))]
    (do
        (def x  (double i))
        ;自作の双曲線正弦関数
        (def d1 (mySinh x 1 x 1.0 x))
        ;標準の双曲線正弦関数
        (def d2 (. Math sinh x))
        ;標準関数との差異
        (println (format "%3d : %17.10f - %17.10f = %13.10f" i d1 d2 (- d1 d2)))))
;自作の双曲線余弦関数
(defn myCosh [x n numerator denominator y]
    (def m     (* 2 n))
    (def denom (* denominator (* m (- m 1))))
    (def nume  (* numerator (* x x)))
    (def a     (/ nume denom))
    ;十分な精度になったら処理を抜ける
    (if (<= (. Math abs a) 0.00000000001)
        y
        (+ y (myCosh x (+ n 1) nume denom a))))
(doseq
    [i  (map #(- % 10)
        (range 0 21))]
    (do
        (def x  (double i))
        ;自作の双曲線余弦関数
        (def d1 (myCosh x 1 1.0 1.0 1.0))
        ;標準の双曲線余弦関数
        (def d2 (. Math cosh x))
        ;標準関数との差異
        (println (format "%3d : %17.10f - %17.10f = %13.10f" i d1 d2 (- d1 d2)))))
;自作の逆正接関数
(defn myAtan [x x2 n t]
    (def m     (quot n 2))
    (def denom (/ (* (* m m) x2) (+ n t)))
    (def nume  (- n 2))
    (if (<= nume 1)
        (/ x (+ 1 denom))
        (myAtan x x2 nume denom)))

(doseq
    [degree (map #(- (* % 15) 45)
            (range 0 7))]
    (do (def radian (/ (* degree (. Math PI)) 180.0))
        (def x2     (* radian radian))
        ;自作の逆正接関数
        (def d1     (myAtan radian x2 23 0))
        ;標準の逆正接関数
        (def d2     (. Math atan radian))
        ;標準関数との差異
        (println (format "%3d : %13.10f - %13.10f = %13.10f" degree d1 d2 (- d1 d2)))))

-45 : -0.6657737500 - -0.6657737500 =  0.0000000000
-30 : -0.4823479071 - -0.4823479071 =  0.0000000000
-15 : -0.2560527700 - -0.2560527700 =  0.0000000000
  0 :  0.0000000000 -  0.0000000000 =  0.0000000000
 15 :  0.2560527700 -  0.2560527700 =  0.0000000000
 30 :  0.4823479071 -  0.4823479071 = -0.0000000000
 45 :  0.6657737500 -  0.6657737500 = -0.0000000000
nil



(doseq
    [i (map #(+ (* % 2) 1)
       (range 5 16))]
    (do (def radian 1)
        (def x2 (* radian radian))
        ;自作の逆正接関数
        (def d1 (myAtan radian x2 i 0))
        ;標準関数との差異
        (def pi (* d1 4.0))
        (println (format "%2d : %13.10f, %13.10f" i pi (- (pi (. Math PI))))))



11 :  3.1414634146, -0.0001292390
13 :  3.1416149068,  0.0000222532
15 :  3.1415888251, -0.0000038285
17 :  3.1415933119,  0.0000006583
19 :  3.1415925404, -0.0000001131
21 :  3.1415926730,  0.0000000194
23 :  3.1415926503, -0.0000000033
25 :  3.1415926542,  0.0000000006
27 :  3.1415926535, -0.0000000001
29 :  3.1415926536,  0.0000000000
31 :  3.1415926536, -0.0000000000
nil


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

(def a 0)
(def b 1)

; 台形則で積分
(doseq [j (range 1 11)]
    (def n (Math/pow 2 j))
    (def h (/ (- b a) n))
    (def x a)
    (def s 0)
    (doseq [i (range 1 n)]
        (def x (+ x h))
        (def s (+ s (f x))))
    (def t1 (double (* h (+ (/ (+ (f a) (f b)) 2) s))))
    (def t2 (- t1 (. Math PI)))
    ; 結果を π と比較
    (println (format "%2d : %13.10f, %13.10f" j t1 t2)))
(defn f[x]
    (/ 4 (+ 1 (* x x))))

(def a 0)
(def b 1)

; 中点則で積分
(doseq [j (range 1 11)]
    (def n (Math/pow 2 j))
    (def h (/ (- b a) n))
    (def x (+ a (/ h 2)))
    (def s 0)
    (doseq [i (range 1 (+ n 1))]
        (def s (+ s (f x)))
        (def x (+ x h)))
    (def t1 (double (* h s)))
    (def t2 (- t1 (. Math PI)))
    ; 結果を π と比較
    (println (format "%2d : %13.10f, %13.10f" j t1 t2)))
(defn f[x]
    (/ 4 (+ 1 (* x x))))

(def a 0)
(def b 1)

; Simpson則で積分
(doseq [j (range 1 6)]
    (def n (Math/pow 2 j))
    (def h (/ (- b a) n))
    (def x (+ a h))
    (def s2 0)
    (def s4 0)
    (doseq [i (range 1 (+ (/ n 2) 1))]
        (def s4 (+ s4 (f x)))
        (def x  (+ x  h))
        (def s2 (+ s2 (f x)))
        (def x  (+ x  h))
    )
    (def w2 (double (+ (* (- s2 (f b)) 2) (f a) (f b))))
    (def w4 (* s4 4))
    (def t1 (double (/ (* (+ w2 w4) h) 3)))
    (def t2 (- t1 (. Math PI)))
    ; 結果を π と比較
    (println (format "%2d : %13.10f, %13.10f" j t1 t2)))
(defn f[x]
    (/ 4 (+ 1 (* x x))))

(def a 0)
(def b 1)

(def t (list))
; 台形則で積分
(doseq [j (range 1 11)]
    (def n (Math/pow 2 j))
    (def h (/ (- b a) n))
    (def x a)
    (def s 0)
    (doseq [i (range 1 n)]
        (def x (+ x h))
        (def s (+ s (f x))))
    (def t1 (double (* h (+ (/ (+ (f a) (f b)) 2) s))))
    (def t2 (- t1 (. Math PI)))
    ; 結果を保存
    (def t (cons t1 t)))

; Richardsonの補外法
(def t (reverse t))
(doseq [j (range 2 7)]
    (def n (Math/pow 4 (- j 1)))
    (def s (list))
    (doseq [i (range j 7)]
        (def t1 (+ (second t) (/ (- (second t) (first t)) (- n 1))))
        (def t2 (- t1 (. Math PI)))
        ; 結果を π と比較
        (if (= i j)
            (println (format "%2d : %13.10f, %13.10f" j t1 t2)))
        (def t (rest t))
        (def s (cons t1 s)))
    (def t (reverse s)))
(def N 7)

; 元の関数
(defn f[x]
    (+ (- x (/ (Math/pow x 3.0) (* 3 2))) (/ (Math/pow x 5.0) (* 5 (* 4 (* 3 2))))))

; Lagrange (ラグランジュ) 補間
(defn lagrange [d x y]
    (def sum_list (list 0.0))
    (doseq [i (range 0 N)]
        (def prod_list (list (nth y i)))
        (doseq [j (range 0 N)]
            (if (= i j)
                (def prod_list (cons 1 prod_list))
                (do (def w1 (/ (- d (nth x j)) (- (nth x i) (nth x j))))
                    (def prod_list (cons w1 prod_list)))))
        (def w2 (reduce * prod_list))
        (def sum_list (cons w2 sum_list)))
    (reduce + sum_list))

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

; 0.5刻みで 与えられていない値を補間
(def d1 (map #(- (* % 0.5) 4.5) (range 0 19)))
(def d2 (map #(f %) d1))
(def d3 (map #(lagrange % x y) d1))

(doseq [d (map list d1 d2 d3)]
    (def d1 (nth d 0))
    (def d2 (nth d 1))
    (def d3 (nth d 2))
    (println (format "%5.2f\t%8.5f\t%8.5f\t%8.5f" d1 d2 d3 (- d2 d3))))
(def N 7)

; 元の関数
(defn f[x]
    (+ (- x (/ (Math/pow x 3.0) (* 3 2))) (/ (Math/pow x 5.0) (* 5 (* 4 (* 3 2))))))

; Neville (ネヴィル) 補間
(defn neville [d x y]
    (def w (cons y nil))
    (doseq [j (range 1 N)]
        (def t_list (list))
        (doseq [i (range 0 (- N j))]
            (def t (nth w 0))
            (def t_list
                (cons
                    (+ (nth t (+ i 1))
                       (/ (* (- (nth t (+ i 1)) (nth t i))
                             (- d (nth x (+ i j))))
                          (- (nth x (+ i j)) (nth x i))))
                    t_list
                )
            )
        )
        (def w (cons (reverse t_list) w))
    )
    (nth (nth w 0) 0)
)

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

; 0.5刻みで 与えられていない値を補間
(def d1 (map #(- (* % 0.5) 4.5) (range 0 19)))
(def d2 (map #(f %) d1))
(def d3 (map #(neville % x y) d1))

(doseq [d (map list d1 d2 d3)]
    (def d1 (nth d 0))
    (def d2 (nth d 1))
    (def d3 (nth d 2))
    (println (format "%5.2f\t%8.5f\t%8.5f\t%8.5f" d1 d2 d3 (- d2 d3))))
(def N 7)

; 元の関数
(defn f[x]
    (+ (- x (/ (Math/pow x 3.0) (* 3 2))) (/ (Math/pow x 5.0) (* 5 (* 4 (* 3 2))))))

; Newton (ニュートン) 補間
(defn newton [d x a]
    (def sum_list (list (nth a 0)))
    (doseq [i (range 1 N)]
        (def prod_list (list (nth a i)))
        (doseq [j (range 0 i)]
            (def prod_list (cons (- d (nth x j)) prod_list)))
        (def w (reduce * prod_list))
        (def sum_list (cons w sum_list)))
    (reduce + sum_list))

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

; 差分商の表を作る
(def d (cons y nil))
(doseq [i (range 1 N)]
    (def w (nth d 0))
    (def t (list))
    (doseq [j (range 0 (- N i))]
        (def t
            (cons
                (/ (- (nth w (+ j 1)) (nth w j))
                   (- (nth x (+ j i)) (nth x j)))
                t
            )
        )
    )
    (def d (cons (reverse t) d))
)
(def d (reverse d))

; n階差分商
(def a (map #(nth (nth d %) 0) (range 0 N)))

; 0.5刻みで 与えられていない値を補間
(def d1 (map #(- (* % 0.5) 4.5) (range 0 19)))
(def d2 (map #(f %) d1))
(def d3 (map #(newton % x a) d1))

(doseq [d (map list d1 d2 d3)]
    (def d1 (nth d 0))
    (def d2 (nth d 1))
    (def d3 (nth d 2))
    (println (format "%5.2f\t%8.5f\t%8.5f\t%8.5f" d1 d2 d3 (- d2 d3))))
(def N    7)
(def Nx2 14)

; 元の関数
(defn f[x]
    (+ (- x (/ (Math/pow x 3.0) (* 3 2))) (/ (Math/pow x 5.0) (* 5 (* 4 (* 3 2))))))

; 導関数
(defn fd[x]
    (+ (- 1.0 (/ (Math/pow x 2.0) 2)) (/ (Math/pow x 4.0) (* 4 (* 3 2)))))

; Hermite (エルミート) 補間
(defn hermite [d z a]
    (def sum_list (list (nth a 0)))
    (doseq [i (range 1 Nx2)]
        (def prod_list (list (nth a i)))
        (doseq [j (range 0 i)]
            (def prod_list (cons (- d (nth z j)) prod_list)))
        (def w (reduce * prod_list))
        (def sum_list (cons w sum_list)))
    (reduce + sum_list))

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

; 差分商の表を作る
(def z (map #(nth x %) (map #(quot % 2) (range 0 Nx2))))
(def w (map #(nth y %) (map #(quot % 2) (range 0 Nx2))))
(def d (cons w nil))
(doseq [i (range 1 Nx2)]
    (def w (nth d 0))
    (def t (list))
    (doseq [j (range 0 (- Nx2 i))]
        (def t
            (cons
                (if (and (= i 1) (= (rem j 2) 0))
                    (nth yd (quot j 2))
                    (/ (- (nth w (+ j 1)) (nth w j))
                       (- (nth z (+ j i)) (nth z j)))
                )
                t
            )
        )
    )
    (def d (cons (reverse t) d))
)
(def d (reverse d))

; n階差分商
(def a (map #(nth (nth d %) 0) (range 0 Nx2)))

; 0.5刻みで 与えられていない値を補間
(def d1 (map #(- (* % 0.5) 4.5) (range 0 19)))
(def d2 (map #(f %) d1))
(def d3 (map #(hermite % z a) d1))

(doseq [d (map list d1 d2 d3)]
    (def d1 (nth d 0))
    (def d2 (nth d 1))
    (def d3 (nth d 2))
    (println (format "%5.2f\t%8.5f\t%8.5f\t%8.5f" d1 d2 d3 (- d2 d3))))
(def N 7)

; 元の関数
(defn f[x]
    (+ (- x (/ (Math/pow x 3.0) (* 3 2))) (/ (Math/pow x 5.0) (* 5 (* 4 (* 3 2))))))

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

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

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

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

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

    (+ (+ a3 b4) c4)
)

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

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

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

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

; 0.5刻みで 与えられていない値を補間
(def d1 (map #(- (* % 0.5) 4.5) (range 0 19)))
(def d2 (map #(f %) d1))
(def d3 (map #(spline % x y z) d1))

(doseq [d (map list d1 d2 d3)]
    (def d1 (nth d 0))
    (def d2 (nth d 1))
    (def d3 (nth d 2))
    (println (format "%5.2f\t%8.5f\t%8.5f\t%8.5f" d1 d2 d3 (- d2 d3))))
; 重力加速度
(def g -9.8)
; 空気抵抗係数
(def k -0.01)
; 時間間隔(秒)
(def h 0.01)

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

; 空気抵抗による水平方向の減速分
(defn fx[vx vy]
    (* k (* (. Math sqrt (+ (* vx vx) (* vy vy))) vx)))
; 重力と空気抵抗による鉛直方向の減速分
(defn fy[vx vy]
    (+ g (* k (* (. Math sqrt (+ (* vx vx) (* vy vy))) vy))))

; Euler法
(defn euler[i vx vy x y]
    ; 経過秒数
    (def t (* i h))

    ; 位置
    (def wx (+ x (* h vx)))
    (def wy (+ y (* h vy)))

    (println (format "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f" t vx vy wx wy))

    ; 速度
    (def wvx (+ vx (* h (fx vx vy))))
    (def wvy (+ vy (* h (fy vx vy))))

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

(euler 1 vx vy x y)
; 重力加速度
(def g -9.8)
; 空気抵抗係数
(def k -0.01)
; 時間間隔(秒)
(def h 0.01)

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

; 空気抵抗による水平方向の減速分
(defn fx[vx vy]
    (* k (* (. Math sqrt (+ (* vx vx) (* vy vy))) vx)))
; 重力と空気抵抗による鉛直方向の減速分
(defn fy[vx vy]
    (+ g (* k (* (. Math sqrt (+ (* vx vx) (* vy vy))) vy))))

; Heun法
(defn heun[i vx vy x y]
    ; 経過秒数
    (def t (* i h))

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

    (def wx   (+  x (* h (/ (+     vx         wvx2      ) 2))))
    (def wy   (+  y (* h (/ (+     vy         wvy2      ) 2))))
    (def wvx  (+ vx (* h (/ (+ (fx vx vy) (fx wvx2 wvy2)) 2))))
    (def wvy  (+ vy (* h (/ (+ (fy vx vy) (fy wvx2 wvy2)) 2))))

    (println (format "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%9.5f" t wvx wvy wx wy))

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

(heun 1 vx vy x y)
; 重力加速度
(def g -9.8)
; 空気抵抗係数
(def k -0.01)
; 時間間隔(秒)
(def h 0.01)

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

; 空気抵抗による水平方向の減速分
(defn fx[vx vy]
    (* k (* (. Math sqrt (+ (* vx vx) (* vy vy))) vx)))
; 重力と空気抵抗による鉛直方向の減速分
(defn fy[vx vy]
    (+ g (* k (* (. Math sqrt (+ (* vx vx) (* vy vy))) vy))))

;中点法
(defn midpoint[i vx vy x y]
    ; 経過秒数
    (def t (* i h))

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

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

    (println (format "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%9.5f" t wvx wvy wx wy))

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

(midpoint 1 vx vy x y)
; 重力加速度
(def g -9.8)
; 空気抵抗係数
(def k -0.01)
; 時間間隔(秒)
(def h 0.01)

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

; 空気抵抗による水平方向の減速分
(defn fx[vx vy]
    (* k (* (. Math sqrt (+ (* vx vx) (* vy vy))) vx)))
; 重力と空気抵抗による鉛直方向の減速分
(defn fy[vx vy]
    (+ g (* k (* (. Math sqrt (+ (* vx vx) (* vy vy))) vy))))

; Runge-Kutta法
(defn rungekutta[i vx vy x y]
    ; 経過秒数
    (def t (* i h))

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

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

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

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

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

    (println (format "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%9.5f" t wvx wvy wx wy))

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

(rungekutta 1 vx vy x y)
; 重力加速度
(def g -9.8)
; 空気抵抗係数
(def k -0.01)
; 時間間隔(秒)
(def h 0.01)

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

; 空気抵抗による水平方向の減速分
(defn fx[vx vy]
    (* k (* (. Math sqrt (+ (* vx vx) (* vy vy))) vx)))
; 重力と空気抵抗による鉛直方向の減速分
(defn fy[vx vy]
    (+ g (* k (* (. Math sqrt (+ (* vx vx) (* vy vy))) vy))))

; Runge-Kutta-Gill法
(defn rungekuttagill[i vx vy x y]
    ; 経過秒数
    (def t (* i h))

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

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

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

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

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

    (println (format "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%9.5f" t wvx wvy wx wy))

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

(rungekuttagill 1 vx vy x y)
(defn f[x]
    (- (* x x) 2.0))

(defn bisection [a b]
    ; 区間 (a, b) の中点 c = (a + b) / 2
    (def c (/ (+ a b) 2.0))
    (println (format "%12.10f\t%13.10f" c (- c (Math/sqrt 2))))

    (def fc (f c))
    (if (< (. Math abs fc) 0.00000000001)
        c
        (if (< fc 0)
            (bisection c b)
            (bisection a c))))

(def a 1.0)
(def b 2.0)
(println (format "%12.10f" (bisection a b)))
(defn f[x]
    (- (* x x) 2.0))

(defn falseposition [a b]
    ; 点 (a,f(a)) と 点 (b,f(b)) を結ぶ直線と x軸の交点
    (def c (/ (- (* a (f b)) (* b (f a))) (- (f b) (f a))))
    (println (format "%12.10f\t%13.10f" c (- c (Math/sqrt 2))))

    (def fc (f c))
    (if (< (. Math abs fc) 0.00000000001)
        c
        (if (< fc 0)
            (falseposition c b)
            (falseposition a c))))

(def a 1.0)
(def b 2.0)
(println (format "%12.10f" (falseposition a b)))
(defn g[x]
    (+ (/ x 2.0) (/ 1.0 x)))

(defn iterative [x0]
    (def x1 (g x0))
    (println (format "%12.10f\t%13.10f" x1 (- x1 (Math/sqrt 2))))

    (if (< (. Math abs (- x1 x0)) 0.00000000001)
        x1
        (iterative x1)))

(def x 1.0)
(println (format "%12.10f" (iterative x)))
(defn f0[x]
    (- (* x x) 2.0))

(defn f1[x]
    (* 2.0 x))

(defn newton [x0]
    (def x1 (- x0 (/ (f0 x0) (f1 x0))))
    (println (format "%12.10f\t%13.10f" x1 (- x1 (Math/sqrt 2))))

    (if (< (. Math abs (- x1 x0)) 0.00000000001)
        x1
        (newton x1)))

(def x 2.0)
(println (format "%12.10f" (newton x)))
(defn f0[x]
    (- (* x x) 2.0))

(defn f1[x]
    (* 2.0 x))

(defn f2[x]
    2.0)

(defn bailey [x0]
    (def x1 (- x0
               (/ (f0 x0)
                  (- (f1 x0)
                     (/ (* (f0 x0) (f2 x0))
                        (* 2.0 (f1 x0)))))))
    (println (format "%12.10f\t%13.10f" x1 (- x1 (Math/sqrt 2))))

    (if (< (. Math abs (- x1 x0)) 0.00000000001)
        x1
        (bailey x1)))

(def x 2.0)
(println (format "%12.10f" (bailey x)))
(defn f[x]
    (- (* x x) 2.0))

(defn secant [x0 x1]
    (def x2 (- x1
               (/ (* (f x1) (- x1 x0))
                  (- (f x1) (f x0)))))
    (println (format "%12.10f\t%13.10f" x2 (- x2 (Math/sqrt 2))))

    (if (< (. Math abs (- x2 x1)) 0.00000000001)
        x2
        (secant x1 x2)))

(def x0 1.0)
(def x1 2.0)
(println (format "%12.10f" (secant x0 x1)))
inserted by FC2 system