rlm@2: #+TITLE: Building a Classical Mechanics Library in Clojure rlm@2: #+author: Dylan Holmes rlm@2: #+EMAIL: rlm@mit.edu rlm@2: #+MATHJAX: align:"left" mathml:t path:"../MathJax/MathJax.js" rlm@2: #+STYLE: rlm@2: #+OPTIONS: H:3 num:t toc:t \n:nil @:t ::t |:t ^:t -:t f:t *:t <:t rlm@2: #+SETUPFILE: ../templates/level-0.org rlm@2: #+INCLUDE: ../templates/level-0.org rlm@2: #+BABEL: :noweb yes :results silent rlm@2: rlm@2: * Useful data types rlm@2: rlm@2: ** Complex numbers rlm@2: rlm@2: ** Power series rlm@2: rlm@2: ** Tuples and tensors rlm@2: rlm@2: *** Tuples are \ldquo{}sequences with spin\rdquo{} rlm@2: rlm@2: #+srcname: tuples rlm@2: #+begin_src clojure rlm@2: (in-ns 'sicm.utils) rlm@2: rlm@2: ;; Let some objects have spin rlm@2: rlm@2: (defprotocol Spinning rlm@2: (up? [this]) rlm@2: (down? [this])) rlm@2: rlm@2: (defn spin rlm@2: "Returns the spin of the Spinning s, either :up or :down" rlm@2: [#^Spinning s] rlm@2: (cond (up? s) :up (down? s) :down)) rlm@2: rlm@2: rlm@2: ;; DEFINITION: A tuple is a sequence with spin rlm@2: rlm@2: (deftype Tuple rlm@2: [spin coll] rlm@2: rlm@2: clojure.lang.Seqable rlm@2: (seq [this] (seq (.coll this))) rlm@2: rlm@2: clojure.lang.Counted rlm@2: (count [this] (count (.coll this))) rlm@2: rlm@2: Spinning rlm@2: (up? [this] (= ::up (.spin this))) rlm@2: (down? [this] (= ::down (.spin this)))) rlm@2: rlm@2: (defmethod print-method Tuple rlm@2: [o w] rlm@2: (print-simple rlm@2: (if (up? o) rlm@2: (str "u" (.coll o)) rlm@2: (str "d" (vec(.coll o)))) rlm@2: w)) rlm@2: rlm@2: (def tuple? #(= (type %) Tuple)) rlm@2: rlm@2: ;; CONSTRUCTORS rlm@2: rlm@2: (defn up rlm@2: "Create a new up-tuple containing the contents of coll." rlm@2: [coll] rlm@2: (Tuple. ::up coll)) rlm@2: rlm@2: (defn down rlm@2: "Create a new down-tuple containing the contents of coll." rlm@2: [coll] rlm@2: (Tuple. ::down coll)) rlm@2: rlm@2: (defn same-spin rlm@2: "Creates a tuple which has the same spin as tuple and which contains rlm@2: the contents of coll." rlm@2: [tuple coll] rlm@2: (if (up? tuple) rlm@2: (up coll) rlm@2: (down coll))) rlm@2: rlm@2: (defn opposite-spin rlm@2: "Create a tuple which has opposite spin to tuple and which contains rlm@2: the contents of coll." rlm@2: [tuple coll] rlm@2: (if (up? tuple) rlm@2: (down coll) rlm@2: (up coll))) rlm@2: #+end_src rlm@2: rlm@2: rlm@2: rlm@2: *** Matrices rlm@2: #+srcname:matrices rlm@2: #+begin_src clojure rlm@2: (in-ns 'sicm.utils) rlm@2: (require 'incanter.core) ;; use incanter's fast matrices rlm@2: rlm@2: rlm@2: (defn all-equal? [coll] rlm@2: (if (empty? (rest coll)) true rlm@2: (and (= (first coll) (second coll)) rlm@2: (recur (rest coll))))) rlm@2: rlm@2: rlm@2: (defprotocol Matrix rlm@2: (rows [matrix]) rlm@2: (cols [matrix]) rlm@2: (diagonal [matrix]) rlm@2: (trace [matrix]) rlm@2: (determinant [matrix]) rlm@2: (transpose [matrix]) rlm@2: (conjugate [matrix]) rlm@2: ) rlm@2: rlm@2: (extend-protocol Matrix incanter.Matrix rlm@2: (rows [rs] (map down (apply map vector (apply map vector rs)))) rlm@2: (cols [rs] (map up (apply map vector rs))) rlm@2: (diagonal [matrix] (incanter.core/diag matrix) ) rlm@2: (determinant [matrix] (incanter.core/det matrix)) rlm@2: (trace [matrix] (incanter.core/trace matrix)) rlm@2: (transpose [matrix] (incanter.core/trans matrix))) rlm@2: rlm@2: (defn count-rows [matrix] rlm@2: ((comp count rows) matrix)) rlm@2: rlm@2: (defn count-cols [matrix] rlm@2: ((comp count cols) matrix)) rlm@2: rlm@2: (defn square? [matrix] rlm@2: (= (count-rows matrix) (count-cols matrix))) rlm@2: rlm@2: (defn identity-matrix rlm@2: "Define a square matrix of size n-by-n with 1s along the diagonal and rlm@2: 0s everywhere else." rlm@2: [n] rlm@2: (incanter.core/identity-matrix n)) rlm@2: rlm@2: rlm@2: (defn matrix-by-rows rlm@2: "Define a matrix by giving its rows." rlm@2: [& rows] rlm@2: (if rlm@2: (not (all-equal? (map count rows))) rlm@2: (throw (Exception. "All rows in a matrix must have the same number of elements.")) rlm@2: (incanter.core/matrix (vec rows)))) rlm@2: rlm@2: (defn matrix-by-cols rlm@2: "Define a matrix by giving its columns" rlm@2: [& cols] rlm@2: (if (not (all-equal? (map count cols))) rlm@2: (throw (Exception. "All columns in a matrix must have the same number of elements.")) rlm@2: (incanter.core/matrix (vec (apply map vector cols))))) rlm@2: rlm@2: (defn identity-matrix rlm@2: "Define a square matrix of size n-by-n with 1s along the diagonal and rlm@2: 0s everywhere else." rlm@2: [n] rlm@2: (incanter.core/identity-matrix n)) rlm@2: rlm@2: #+end_src rlm@2: rlm@2: ** Generic Operations rlm@2: #+srcname:arith-tuple rlm@2: #+begin_src clojure rlm@2: (in-ns 'sicm.utils) rlm@2: (use 'clojure.contrib.generic.arithmetic rlm@2: 'clojure.contrib.generic.collection rlm@2: 'clojure.contrib.generic.functor rlm@2: 'clojure.contrib.generic.math-functions) rlm@2: rlm@2: (defn numbers? rlm@2: "Returns true if all arguments are numbers, else false." rlm@2: [& xs] rlm@2: (every? number? xs)) rlm@2: rlm@2: (defn tuple-surgery rlm@2: "Applies the function f to the items of tuple and the additional rlm@2: arguments, if any. Returns a Tuple of the same type as tuple." rlm@2: [tuple f & xs] rlm@2: ((if (up? tuple) up down) rlm@2: (apply f (seq tuple) xs))) rlm@2: rlm@2: rlm@2: rlm@2: ;;; CONTRACTION collapses two compatible tuples into a number. rlm@2: rlm@2: (defn contractible? rlm@2: "Returns true if the tuples a and b are compatible for contraction, rlm@2: else false. Tuples are compatible if they have the same number of rlm@2: components, they have opposite spins, and their elements are rlm@2: pairwise-compatible." rlm@2: [a b] rlm@2: (and rlm@2: (isa? (type a) Tuple) rlm@2: (isa? (type b) Tuple) rlm@2: (= (count a) (count b)) rlm@2: (not= (spin a) (spin b)) rlm@2: rlm@2: (not-any? false? rlm@2: (map #(or rlm@2: (numbers? %1 %2) rlm@2: (contractible? %1 %2)) rlm@2: a b)))) rlm@2: rlm@2: (defn contract rlm@2: "Contracts two tuples, returning the sum of the rlm@2: products of the corresponding items. Contraction is recursive on rlm@2: nested tuples." rlm@2: [a b] rlm@2: (if (not (contractible? a b)) rlm@2: (throw rlm@2: (Exception. "Not compatible for contraction.")) rlm@2: (reduce + rlm@2: (map rlm@2: (fn [x y] rlm@2: (if (numbers? x y) rlm@2: (* x y) rlm@2: (contract x y))) rlm@2: a b)))) rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: (defmethod conj Tuple rlm@2: [tuple & xs] rlm@2: (tuple-surgery tuple #(apply conj % xs))) rlm@2: rlm@2: (defmethod fmap Tuple rlm@2: [f tuple] rlm@2: (tuple-surgery tuple (partial map f))) rlm@2: rlm@2: rlm@2: rlm@2: ;; TODO: define Scalar, and add it to the hierarchy above Number and Complex rlm@2: rlm@2: rlm@2: (defmethod * [Tuple Tuple] ; tuple*tuple rlm@2: [a b] rlm@2: (if (contractible? a b) rlm@2: (contract a b) rlm@2: (map (partial * a) b))) rlm@2: rlm@2: rlm@2: (defmethod * [java.lang.Number Tuple] ;; scalar * tuple rlm@2: [a x] (fmap (partial * a) x)) rlm@2: rlm@2: (defmethod * [Tuple java.lang.Number] rlm@2: [x a] (* a x)) rlm@2: rlm@2: (defmethod * [java.lang.Number incanter.Matrix] ;; scalar * matrix rlm@2: [x M] (incanter.core/mult x M)) rlm@2: rlm@2: (defmethod * [incanter.Matrix java.lang.Number] rlm@2: [M x] (* x M)) rlm@2: rlm@2: (defmethod * [incanter.Matrix incanter.Matrix] ;; matrix * matrix rlm@2: [M1 M2] rlm@2: (incanter.core/mmult M1 M2)) rlm@2: rlm@2: (defmethod * [incanter.Matrix Tuple] ;; matrix * tuple rlm@2: [M v] rlm@2: (if (and (apply numbers? v) (up? v)) rlm@2: (* M (matrix-by-cols v)) rlm@2: (throw (Exception. "Currently, you can only multiply a matrix by a tuple of *numbers*")) rlm@2: )) rlm@2: rlm@2: (defmethod * [Tuple incanter.Matrix] ;; tuple * Matrix rlm@2: [v M] rlm@2: (if (and (apply numbers? v) (down? v)) rlm@2: (* (matrix-by-rows v) M) rlm@2: (throw (Exception. "Currently, you can only multiply a matrix by a tuple of *numbers*")) rlm@2: )) rlm@2: rlm@2: rlm@2: (defmethod exp incanter.Matrix rlm@2: [M] rlm@2: (incanter.core/exp M)) rlm@2: rlm@2: #+end_src rlm@2: rlm@2: * Operators and Differentiation rlm@2: ** Operators rlm@2: #+scrname: operators rlm@2: #+begin_src clojure rlm@2: (in-ns 'sicm.utils) rlm@2: (use 'clojure.contrib.seq rlm@2: 'clojure.contrib.generic.arithmetic rlm@2: 'clojure.contrib.generic.collection rlm@2: 'clojure.contrib.generic.math-functions) rlm@2: rlm@2: (defmethod + [clojure.lang.IFn clojure.lang.IFn] rlm@2: [f g] rlm@2: (fn [& args] rlm@2: (+ (apply f args) (apply g args)))) rlm@2: rlm@2: (defmethod * [clojure.lang.IFn clojure.lang.IFn] rlm@2: [f g] rlm@2: (fn [& args] rlm@2: (* (apply f args) (apply g args)))) rlm@2: rlm@2: (defmethod / [clojure.lang.IFn java.lang.Number] rlm@2: [f x] rlm@2: (fn [& args] rlm@2: (/ (apply f args) x))) rlm@2: rlm@2: rlm@2: (defmethod - [clojure.lang.IFn] rlm@2: [f] rlm@2: (fn [& args] rlm@2: (- (apply f args)))) rlm@2: rlm@2: (defmethod - [clojure.lang.IFn clojure.lang.IFn] rlm@2: [f g] rlm@2: (fn [& args] rlm@2: (- (apply f args) (apply g args)))) rlm@2: rlm@2: (defmethod pow [clojure.lang.IFn java.lang.Number] rlm@2: [f x] rlm@2: (fn [& args] rlm@2: (pow (apply f args) x))) rlm@2: rlm@2: rlm@2: (defmethod + [java.lang.Number clojure.lang.IFn] rlm@2: [x f] rlm@2: (fn [& args] rlm@2: (+ x (apply f args)))) rlm@2: rlm@2: (defmethod * [java.lang.Number clojure.lang.IFn] rlm@2: [x f] rlm@2: (fn [& args] rlm@2: (* x (apply f args)))) rlm@2: rlm@2: (defmethod * [clojure.lang.IFn java.lang.Number] rlm@2: [f x] rlm@2: (* x f)) rlm@2: (defmethod + [clojure.lang.IFn java.lang.Number] rlm@2: [f x] rlm@2: (+ x f)) rlm@2: rlm@2: #+end_src rlm@2: rlm@2: ** Differential Terms and Sequences rlm@2: #+srcname: differential rlm@2: #+begin_src clojure rlm@2: (in-ns 'sicm.utils) rlm@2: (use 'clojure.contrib.seq rlm@2: 'clojure.contrib.generic.arithmetic rlm@2: 'clojure.contrib.generic.collection rlm@2: 'clojure.contrib.generic.math-functions) rlm@2: rlm@2: ;;∂ rlm@2: rlm@2: ;; DEFINITION : Differential Term rlm@2: rlm@2: ;; A quantity with infinitesimal components, e.g. x, dxdy, 4dydz. The rlm@2: ;; coefficient of the quantity is returned by the 'coefficient' method, rlm@2: ;; while the sequence of differential parameters is returned by the rlm@2: ;; method 'partials'. rlm@2: rlm@2: ;; Instead of using (potentially ambiguous) letters to denote rlm@2: ;; differential parameters (dx,dy,dz), we use integers. So, dxdz becomes [0 2]. rlm@2: rlm@2: ;; The coefficient can be any arithmetic object; the rlm@2: ;; partials must be a nonrepeating sorted sequence of nonnegative rlm@2: ;; integers. rlm@2: rlm@2: ;; (deftype DifferentialTerm [coefficient partials]) rlm@2: rlm@2: ;; (defn differential-term rlm@2: ;; "Make a differential term from a coefficient and list of partials." rlm@2: ;; [coefficient partials] rlm@2: ;; (if (and (coll? partials) (every? #(and (integer? %) (not(neg? %))) partials)) rlm@2: ;; (DifferentialTerm. coefficient (set partials)) rlm@2: ;; (throw (java.lang.IllegalArgumentException. "Partials must be a collection of integers.")))) rlm@2: rlm@2: rlm@2: ;; DEFINITION : Differential Sequence rlm@2: ;; A differential sequence is a sequence of differential terms, all with different partials. rlm@2: ;; Internally, it is a map from the partials of each term to their coefficients. rlm@2: rlm@2: (deftype DifferentialSeq rlm@2: [terms] rlm@2: ;;clojure.lang.IPersistentMap rlm@2: clojure.lang.Associative rlm@2: (assoc [this key val] rlm@2: (DifferentialSeq. rlm@2: (cons (differential-term val key) terms))) rlm@2: (cons [this x] rlm@2: (DifferentialSeq. (cons x terms))) rlm@2: (containsKey [this key] rlm@2: (not(nil? (find-first #(= (.partials %) key) terms)))) rlm@2: (count [this] (count (.terms this))) rlm@2: (empty [this] (DifferentialSeq. [])) rlm@2: (entryAt [this key] rlm@2: ((juxt #(.partials %) #(.coefficient %)) rlm@2: (find-first #(= (.partials %) key) terms))) rlm@2: (seq [this] (seq (.terms this)))) rlm@2: rlm@2: (def differential? #(= (type %) DifferentialSeq)) rlm@2: rlm@2: (defn zeroth-order? rlm@2: "Returns true if the differential sequence has at most a constant term." rlm@2: [dseq] rlm@2: (and rlm@2: (differential? dseq) rlm@2: (every? rlm@2: #(= #{} %) rlm@2: (keys (.terms dseq))))) rlm@2: rlm@2: (defmethod fmap DifferentialSeq rlm@2: [f dseq] rlm@2: (DifferentialSeq. rlm@2: (fmap f (.terms dseq)))) rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: ;; BUILDING DIFFERENTIAL OBJECTS rlm@2: rlm@2: (defn differential-seq rlm@2: "Define a differential sequence by specifying an alternating rlm@2: sequence of coefficients and lists of partials." rlm@2: ([coefficient partials] rlm@2: (DifferentialSeq. {(set partials) coefficient})) rlm@2: ([coefficient partials & cps] rlm@2: (if (odd? (count cps)) rlm@2: (throw (Exception. "differential-seq requires an even number of terms.")) rlm@2: (DifferentialSeq. rlm@2: (reduce rlm@2: #(assoc %1 (set (second %2)) (first %2)) rlm@2: {(set partials) coefficient} rlm@2: (partition 2 cps)))))) rlm@2: rlm@2: rlm@2: rlm@2: (defn big-part rlm@2: "Returns the part of the differential sequence that is finite, rlm@2: i.e. not infinitely small. If the sequence is zeroth-order, returns rlm@2: the coefficient of the zeroth-order term instead. " rlm@2: [dseq] rlm@2: (if (zeroth-order? dseq) (get (.terms dseq) #{}) rlm@2: (let [m (.terms dseq) rlm@2: keys (sort-by count (keys m)) rlm@2: smallest-var (last (last keys))] rlm@2: (DifferentialSeq. rlm@2: (reduce rlm@2: #(assoc %1 (first %2) (second %2)) rlm@2: {} rlm@2: (remove #((first %) smallest-var) m)))))) rlm@2: rlm@2: rlm@2: (defn small-part rlm@2: "Returns the part of the differential sequence that infinitely rlm@2: small. If the sequence is zeroth-order, returns zero." rlm@2: [dseq] rlm@2: (if (zeroth-order? dseq) 0 rlm@2: (let [m (.terms dseq) rlm@2: keys (sort-by count (keys m)) rlm@2: smallest-var (last (last keys))] rlm@2: (DifferentialSeq. rlm@2: (reduce rlm@2: #(assoc %1 (first %2) (second %2)) {} rlm@2: (filter #((first %) smallest-var) m)))))) rlm@2: rlm@2: rlm@2: rlm@2: (defn cartesian-product [set1 set2] rlm@2: (reduce concat rlm@2: (for [x set1] rlm@2: (for [y set2] rlm@2: [x y])))) rlm@2: rlm@2: (defn nth-subset [n] rlm@2: (if (zero? n) [] rlm@2: (let [lg2 #(/ (log %) (log 2)) rlm@2: k (int(java.lang.Math/floor (lg2 n))) rlm@2: ] rlm@2: (cons k rlm@2: (nth-subset (- n (pow 2 k))))))) rlm@2: rlm@2: (def all-partials rlm@2: (lazy-seq (map nth-subset (range)))) rlm@2: rlm@2: rlm@2: (defn differential-multiply rlm@2: "Multiply two differential sequences. The square of any differential rlm@2: variable is zero since differential variables are infinitesimally rlm@2: small." rlm@2: [dseq1 dseq2] rlm@2: (DifferentialSeq. rlm@2: (reduce rlm@2: (fn [m [[vars1 coeff1] [vars2 coeff2]]] rlm@2: (if (not (empty? (clojure.set/intersection vars1 vars2))) rlm@2: m rlm@2: (assoc m (clojure.set/union vars1 vars2) (* coeff1 coeff2)))) rlm@2: {} rlm@2: (cartesian-product (.terms dseq1) (.terms dseq2))))) rlm@2: rlm@2: rlm@2: rlm@2: (defmethod * [DifferentialSeq DifferentialSeq] rlm@2: [dseq1 dseq2] rlm@2: (differential-multiply dseq1 dseq2)) rlm@2: rlm@2: (defmethod + [DifferentialSeq DifferentialSeq] rlm@2: [dseq1 dseq2] rlm@2: (DifferentialSeq. rlm@2: (merge-with + (.terms dseq1) (.terms dseq2)))) rlm@2: rlm@2: (defmethod * [java.lang.Number DifferentialSeq] rlm@2: [x dseq] rlm@2: (fmap (partial * x) dseq)) rlm@2: rlm@2: (defmethod * [DifferentialSeq java.lang.Number] rlm@2: [dseq x] rlm@2: (fmap (partial * x) dseq)) rlm@2: rlm@2: (defmethod + [java.lang.Number DifferentialSeq] rlm@2: [x dseq] rlm@2: (+ (differential-seq x []) dseq)) rlm@2: (defmethod + [DifferentialSeq java.lang.Number] rlm@2: [dseq x] rlm@2: (+ dseq (differential-seq x []))) rlm@2: rlm@2: (defmethod - DifferentialSeq rlm@2: [x] rlm@2: (fmap - x)) rlm@2: rlm@2: rlm@2: ;; DERIVATIVES rlm@2: rlm@2: rlm@2: rlm@2: (defn linear-approximator rlm@2: "Returns an operator that linearly approximates the given function." rlm@2: ([f df|dx] rlm@2: (fn [x] rlm@2: (let [big-part (big-part x) rlm@2: small-part (small-part x)] rlm@2: ;; f(x+dx) ~= f(x) + f'(x)dx rlm@2: (+ (f big-part) rlm@2: (* (df|dx big-part) small-part) rlm@2: )))) rlm@2: rlm@2: ([f df|dx df|dy] rlm@2: (fn [x y] rlm@2: (let [X (big-part x) rlm@2: Y (big-part y) rlm@2: DX (small-part x) rlm@2: DY (small-part y)] rlm@2: (+ (f X Y) rlm@2: (* DX (f df|dx X Y)) rlm@2: (* DY (f df|dy X Y))))))) rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: (defn D[f] rlm@2: (fn[x] (f (+ x (differential-seq 1 [0] 1 [1] 1 [2]))))) rlm@2: rlm@2: (defn d[partials f] rlm@2: (fn [x] rlm@2: (get rlm@2: (.terms ((D f)x)) rlm@2: (set partials) rlm@2: 0 rlm@2: ))) rlm@2: rlm@2: (defmethod exp DifferentialSeq [x] rlm@2: ((linear-approximator exp exp) x)) rlm@2: rlm@2: (defmethod sin DifferentialSeq rlm@2: [x] rlm@2: ((linear-approximator sin cos) x)) rlm@2: rlm@2: (defmethod cos DifferentialSeq rlm@2: [x] rlm@2: ((linear-approximator cos #(- (sin %))) x)) rlm@2: rlm@2: (defmethod log DifferentialSeq rlm@2: [x] rlm@2: ((linear-approximator log (fn [x] (/ x)) ) x)) rlm@2: rlm@2: (defmethod / [DifferentialSeq DifferentialSeq] rlm@2: [x y] rlm@2: ((linear-approximator / rlm@2: (fn [x y] (/ 1 y)) rlm@2: (fn [x y] (- (/ x (* y y))))) rlm@2: x y)) rlm@2: rlm@2: #+end_src rlm@2: rlm@2: rlm@2: rlm@2: * Derivatives Revisited rlm@2: #+begin_src clojure rlm@2: (in-ns 'sicm.utils) rlm@2: (use 'clojure.contrib.seq rlm@2: 'clojure.contrib.generic.arithmetic rlm@2: 'clojure.contrib.generic.collection rlm@2: 'clojure.contrib.generic.math-functions) rlm@2: rlm@2: (defn replace-in rlm@2: "Replaces the nth item in coll with the given item. If n is larger rlm@2: than the size of coll, adds n to the end of the collection." rlm@2: [coll n item] rlm@2: (concat rlm@2: (take n coll) rlm@2: [item] rlm@2: (drop (inc n) coll))) rlm@2: rlm@2: (defn euclidean-structure [f partials] rlm@2: (fn sd[g v] rlm@2: (cond rlm@2: (tuple? v) rlm@2: (opposite-spin rlm@2: v rlm@2: (map rlm@2: (fn [n] rlm@2: (sd (fn [xn] rlm@2: (g rlm@2: (same-spin v rlm@2: (replace-in v n xn)))) rlm@2: (nth v n))) rlm@2: (range (count v))))))) rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: #+end_src rlm@2: rlm@2: rlm@2: * Symbolic Quantities rlm@2: rlm@2: #+srcname: symbolic rlm@2: #+begin_src clojure rlm@2: (in-ns 'sicm.utils) rlm@2: (use 'clojure.contrib.generic.arithmetic rlm@2: 'clojure.contrib.generic.math-functions) rlm@2: rlm@2: (deftype Symbolic rlm@2: [type rlm@2: expression]) rlm@2: rlm@2: (defn print-expression [s] rlm@2: (print (.expression s))) rlm@2: rlm@2: (defn symbolic-number rlm@2: [symbol] rlm@2: (Symbolic. java.lang.Number symbol)) rlm@2: rlm@2: (defn simplify-expression [s] rlm@2: (let [e (.expression s)] rlm@2: (cond rlm@2: (not(list? e)) e rlm@2: (= (first e) '+) rlm@2: ) rlm@2: rlm@2: rlm@2: rlm@2: (defmethod + [Symbolic Symbolic] rlm@2: [x y] rlm@2: (Symbolic. (.type x) (list '+ (.expression x) (.expression y)))) rlm@2: rlm@2: (defmethod + [Symbolic java.lang.Number] rlm@2: [s x] rlm@2: (if (zero? x) rlm@2: s rlm@2: (Symbolic. (.type s) (list '+ (.expression s) x)))) rlm@2: rlm@2: (defmethod sin Symbolic rlm@2: [s] rlm@2: (Symbolic. (.type s) (list 'sin (.expression s)))) rlm@2: rlm@2: #+end_src rlm@2: rlm@2: * COMMENT To-be-tangled Source rlm@2: #+begin_src clojure :tangle utils.clj rlm@2: (ns sicm.utils) rlm@2: rlm@2: <> rlm@2: <> rlm@2: <> rlm@2: rlm@2: <> rlm@2: #+end_src rlm@2: