cauchy-clojure.clj

(ns cauchy-clojure.core
  (:require [clojure.core.match :as cm]
            [clojure.math.numeric-tower :refer [abs]]))

(defn rat [x]
  [:rat x])

(defn rat-val [x]
  (second x))

(defn rat? [x]
  (and (vector? x)
       (= :rat (first x))))

(defn lim [x]
  [:lim x])

(defn lim? [x]
  (and (vector? x)
       (= :lim (first x))))

(defn real? [x]
  (or (rat? x) (lim? x)))

(defn approx [x eps]
  (cm/match
   [x]
   [[:rat a]] a
   [[:lim f]] (f eps)))

(defn fapprox [x eps]
  (float (approx x eps)))

(defn dapprox [x eps]
  (double (approx x eps)))

(defn lift-rat-op [op]
  (fn [a b]
    (cm/match
     [a b]
     [[:rat x] [:rat y]] (rat (op x y))
     [[:lim f] [:rat y]] (lim (fn [eps] (op (f eps) y)))
     [[:rat x] [:lim g]] (lim (fn [eps] (op x (g eps))))
     [[:lim f] [:lim g]] (lim (fn [eps] (op (f eps) (g eps)))))))

(defn times [a b]
  (cm/match
   [a b]
   [[:rat x] [:rat y]] (rat (* x y))
   [[:lim f] [:rat y]] (lim (fn [eps] (* (f (/ eps y)) y)))
   [[:rat x] [:lim g]] (lim (fn [eps] (* x (g (/ eps x)))))
   [[:lim f] [:lim g]] (lim (fn [eps]
                              ;; This doesn't seem quite right
                              (let [a (f eps)]
                                (* a (g (/ eps a))))))))

(def plus (lift-rat-op +))
(def minus (lift-rat-op -))

;; An exponentiation fn that works on rationals
(def leibniz-pi-4
  (lim (fn [eps]
         ;; (println "Effective eps" eps)
         (loop [approx 1
                denom 1]
           (let [t1 (/ -1 (+ denom 2))
                 t2 (/ 1 (+ denom 4))
                 new-approx (+ approx t1 t2)
                 del (abs (- t1 t2))]
             ;; (println "del: " del)
             (if (< del eps )
               new-approx
               (recur new-approx (+ denom 4))))))))

(def leibniz-pi (times leibniz-pi-4 (rat 4)))

(comment
  (fapprox leibniz-pi 1/1000)


  )

;; taylor series expansion of arctan
(defn arctan [x]
  (lim (fn [eps]
         (println "effective eps" eps)
         (loop [approx x
                numer x
                denom 1]
           (let [n1 (* numer x x -1)
                 d1 (+ denom 2)
                 t1 (/ n1 d1)
                 n2 (* n1 x x -1)
                 d2 (+ d1 2)
                 t2 (/ n2 d2)
                 new-approx (+ approx t1 t2)
                 del (abs (- t1 t2))]
             (println "del: " del)
             (if (< del eps )
               new-approx
               (recur new-approx n2 d2)))))))

(def machin-pi-4
  (minus (times (rat 4) (arctan 1/5))
         (arctan 1/239)))


(def machin-pi (times machin-pi-4 (rat 4)))

(comment
  (dapprox machin-pi 1/1000)

  )