Mercurial > dylan

     1 <?xml version="1.0" encoding="utf-8"?>

     2 <!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN"

     3                "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">

     4 <html xmlns="http://www.w3.org/1999/xhtml"

     5 lang="en" xml:lang="en">

     6 <head>

     7 <title>Building a Classical Mechanics Library in Clojure</title>

     8 <meta http-equiv="Content-Type" content="text/html;charset=utf-8"/>

     9 <meta name="generator" content="Org-mode"/>

    10 <meta name="generated" content="2011-08-09 18:41:37 EDT"/>

    11 <meta name="author" content="Robert McIntyre & Dylan Holmes"/>

    12 <meta name="description" content=""/>

    13 <meta name="keywords" content=""/>

    14 <style type="text/css">

    15  <!--/*--><![CDATA[/*><!--*/

    16   html { font-family: Times, serif; font-size: 12pt; }

    17   .title  { text-align: center; }

    18   .todo   { color: red; }

    19   .done   { color: green; }

    20   .tag    { background-color: #add8e6; font-weight:normal }

    21   .target { }

    22   .timestamp { color: #bebebe; }

    23   .timestamp-kwd { color: #5f9ea0; }

    24   .right  {margin-left:auto; margin-right:0px;  text-align:right;}

    25   .left   {margin-left:0px;  margin-right:auto; text-align:left;}

    26   .center {margin-left:auto; margin-right:auto; text-align:center;}

    27   p.verse { margin-left: 3% }

    28   pre {

    29 	border: 1pt solid #AEBDCC;

    30 	background-color: #F3F5F7;

    31 	padding: 5pt;

    32 	font-family: courier, monospace;

    33         font-size: 90%;

    34         overflow:auto;

    35   }

    36   table { border-collapse: collapse; }

    37   td, th { vertical-align: top;  }

    38   th.right  { text-align:center;  }

    39   th.left   { text-align:center;   }

    40   th.center { text-align:center; }

    41   td.right  { text-align:right;  }

    42   td.left   { text-align:left;   }

    43   td.center { text-align:center; }

    44   dt { font-weight: bold; }

    45   div.figure { padding: 0.5em; }

    46   div.figure p { text-align: center; }

    47   textarea { overflow-x: auto; }

    48   .linenr { font-size:smaller }

    49   .code-highlighted {background-color:#ffff00;}

    50   .org-info-js_info-navigation { border-style:none; }

    51   #org-info-js_console-label { font-size:10px; font-weight:bold;

    52                                white-space:nowrap; }

    53   .org-info-js_search-highlight {background-color:#ffff00; color:#000000;

    54                                  font-weight:bold; }

    55   /*]]>*/-->

    56 </style>

    57 <link rel="stylesheet" type="text/css" href="../css/aurellem.css" />

    58 <script type="text/javascript">

    59 <!--/*--><![CDATA[/*><!--*/

    60  function CodeHighlightOn(elem, id)

    61  {

    62    var target = document.getElementById(id);

    63    if(null != target) {

    64      elem.cacheClassElem = elem.className;

    65      elem.cacheClassTarget = target.className;

    66      target.className = "code-highlighted";

    67      elem.className   = "code-highlighted";

    68    }

    69  }

    70  function CodeHighlightOff(elem, id)

    71  {

    72    var target = document.getElementById(id);

    73    if(elem.cacheClassElem)

    74      elem.className = elem.cacheClassElem;

    75    if(elem.cacheClassTarget)

    76      target.className = elem.cacheClassTarget;

    77  }

    78 /*]]>*///-->

    79 </script>

    80 <script type="text/javascript" src="../MathJax/MathJax.js">

    81 <!--/*--><![CDATA[/*><!--*/

    82     MathJax.Hub.Config({

    83         // Only one of the two following lines, depending on user settings

    84         // First allows browser-native MathML display, second forces HTML/CSS

    85             config: ["MMLorHTML.js"], jax: ["input/TeX"],

    86         //  jax: ["input/TeX", "output/HTML-CSS"],

    87         extensions: ["tex2jax.js","TeX/AMSmath.js","TeX/AMSsymbols.js",

    88                      "TeX/noUndefined.js"],

    89         tex2jax: {

    90             inlineMath: [ ["\\(","\\)"] ],

    91             displayMath: [ ['$$','$$'], ["\\[","\\]"], ["\\begin{displaymath}","\\end{displaymath}"] ],

    92             skipTags: ["script","noscript","style","textarea","pre","code"],

    93             ignoreClass: "tex2jax_ignore",

    94             processEscapes: false,

    95             processEnvironments: true,

    96             preview: "TeX"

    97         },

    98         showProcessingMessages: true,

    99         displayAlign: "left",

   100         displayIndent: "2em",

   101 

   102         "HTML-CSS": {

   103              scale: 100,

   104              availableFonts: ["STIX","TeX"],

   105              preferredFont: "TeX",

   106              webFont: "TeX",

   107              imageFont: "TeX",

   108              showMathMenu: true,

   109         },

   110         MMLorHTML: {

   111              prefer: {

   112                  MSIE:    "MML",

   113                  Firefox: "MML",

   114                  Opera:   "HTML",

   115                  other:   "HTML"

   116              }

   117         }

   118     });

   119 /*]]>*///-->

   120 </script>

   121 </head>

   122 <body>

   123 

   124 <div id="content">

   125 

   126 

   127 

   128 <div class="header">

   129   <div class="float-right">	

   130     <!-- 

   131     <form>

   132       <input type="text"/><input type="submit" value="search the blog &raquo;"/> 

   133     </form>

   134     -->

   135   </div>

   136 

   137   <h1>aurellem <em>&#x2609;</em></h1>

   138   <ul class="nav">

   139     <li><a href="/">read the blog &raquo;</a></li>

   140     <!-- li><a href="#">learn about us &raquo;</a></li-->

   141   </ul>

   142 </div>

   143 

   144 <h1 class="title">Building a Classical Mechanics Library in Clojure</h1>

   145 

   146 

   147 

   148 

   149 

   150 

   151 

   152 <div id="table-of-contents">

   153 <h2>Table of Contents</h2>

   154 <div id="text-table-of-contents">

   155 <ul>

   156 <li><a href="#sec-1">1 Generic Arithmetic </a></li>

   157 <li><a href="#sec-2">2 Useful Data Types </a>

   158 <ul>

   159 <li><a href="#sec-2-1">2.1 Complex Numbers </a></li>

   160 <li><a href="#sec-2-2">2.2 Tuples and Tensors </a>

   161 <ul>

   162 <li><a href="#sec-2-2-1">2.2.1 Contraction </a></li>

   163 <li><a href="#sec-2-2-2">2.2.2 Matrices </a></li>

   164 </ul>

   165 </li>

   166 <li><a href="#sec-2-3">2.3 Power Series </a></li>

   167 </ul>

   168 </li>

   169 <li><a href="#sec-3">3 Basic Utilities </a>

   170 <ul>

   171 <li><a href="#sec-3-1">3.1 Sequence manipulation </a></li>

   172 <li><a href="#sec-3-2">3.2 Ranges, Airity and Function Composition </a></li>

   173 </ul>

   174 </li>

   175 <li><a href="#sec-4">4 Numerical Methods </a></li>

   176 <li><a href="#sec-5">5 Differentiation </a></li>

   177 <li><a href="#sec-6">6 Symbolic Manipulation </a></li>

   178 </ul>

   179 </div>

   180 </div>

   181 

   182 <div id="outline-container-1" class="outline-2">

   183 <h2 id="sec-1"><span class="section-number-2">1</span> Generic Arithmetic </h2>

   184 <div class="outline-text-2" id="text-1">

   185 

   186 

   187 

   188 

   189 

   190 <pre class="src src-clojure">(<span style="color: #f0dfaf; font-weight: bold;">ns</span> sicm.utils)

   191 (<span style="color: #f0dfaf; font-weight: bold;">in-ns</span> 'sicm.utils)

   192 

   193 

   194 

   195 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">all-equal?</span> [coll]

   196   (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #8cd0d3;">empty?</span> (<span style="color: #8cd0d3;">rest</span> coll)) true

   197       (<span style="color: #f0dfaf; font-weight: bold;">and</span> (<span style="color: #8cd0d3;">=</span> (<span style="color: #8cd0d3;">first</span> coll) (<span style="color: #8cd0d3;">second</span> coll))

   198            (<span style="color: #f0dfaf; font-weight: bold;">recur</span> (<span style="color: #8cd0d3;">rest</span> coll)))))

   199 

   200 

   201 (<span style="color: #f0dfaf; font-weight: bold;">defprotocol</span> <span style="color: #f0dfaf;">Arithmetic</span>

   202   (zero [this])

   203   (one [this])

   204   (negate [this])

   205   (invert [this])

   206   (conjugate [this]))

   207 

   208 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">zero?</span> [x]

   209   (<span style="color: #8cd0d3;">=</span> x (zero x)))

   210 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">one?</span> [x]

   211   (<span style="color: #8cd0d3;">=</span> x (one x)))

   212 

   213 

   214 

   215 (<span style="color: #8cd0d3;">extend-protocol</span> Arithmetic

   216   java.lang.Number

   217   (zero [x] 0)

   218   (one [x] 1)

   219   (negate [x] (<span style="color: #8cd0d3;">-</span> x))

   220   (invert [x] (<span style="color: #8cd0d3;">/</span> x))

   221   )

   222 

   223 (<span style="color: #8cd0d3;">extend-protocol</span> Arithmetic

   224   clojure.lang.IFn

   225   (one [f] identity)

   226   (negate [f] (<span style="color: #8cd0d3;">comp</span> negate f))

   227   (invert [f] (<span style="color: #8cd0d3;">comp</span> invert f)))

   228 

   229     

   230 (<span style="color: #8cd0d3;">extend-protocol</span> Arithmetic

   231   clojure.lang.Seqable

   232   (zero [this] (<span style="color: #8cd0d3;">map</span> zero this))

   233   (one [this] (<span style="color: #8cd0d3;">map</span> one this))

   234   (invert [this] (<span style="color: #8cd0d3;">map</span> invert this))

   235   (negate [this] (<span style="color: #8cd0d3;">map</span> negate this)))

   236 

   237 

   238 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">ordered-like</span>

   239   <span style="color: #8fb28f;">"Create a comparator using the sorted collection as an</span>

   240 <span style="color: #8fb28f;">  example. Elements not in the sorted collection are sorted to the</span>

   241 <span style="color: #8fb28f;">  end."</span>

   242   [sorted-coll]

   243   (<span style="color: #f0dfaf; font-weight: bold;">let</span> [

   244         sorted-coll? (<span style="color: #8cd0d3;">set</span> sorted-coll) 

   245         ascending-pair?

   246         (<span style="color: #8cd0d3;">set</span>(<span style="color: #8cd0d3;">reduce</span> concat

   247                     (map-indexed

   248                      (<span style="color: #8cd0d3;">fn</span> [n x]

   249                        (<span style="color: #8cd0d3;">map</span> #(<span style="color: #8cd0d3;">vector</span> x %) (<span style="color: #8cd0d3;">nthnext</span> sorted-coll n)))

   250                      sorted-coll)))]

   251     (<span style="color: #8cd0d3;">fn</span> [x y]

   252       (<span style="color: #f0dfaf; font-weight: bold;">cond</span>

   253        (<span style="color: #8cd0d3;">=</span> x y) 0

   254        (ascending-pair? [x y]) -1

   255        (ascending-pair? [y x]) 1

   256        (sorted-coll? x) -1

   257        (sorted-coll? y) 1))))

   258   

   259 

   260 

   261 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">type-precedence</span>

   262   (ordered-like [incanter.Matrix]))

   263 

   264 (<span style="color: #f0dfaf; font-weight: bold;">defmulti</span> <span style="color: #f0dfaf;">add</span>

   265     (<span style="color: #8cd0d3;">fn</span> [x y]

   266       (<span style="color: #8cd0d3;">sort</span> type-precedence [(<span style="color: #8cd0d3;">type</span> x)(<span style="color: #8cd0d3;">type</span> y)]))) 

   267 

   268 (<span style="color: #f0dfaf; font-weight: bold;">defmulti</span> <span style="color: #f0dfaf;">multiply</span>

   269   (<span style="color: #8cd0d3;">fn</span> [x y]

   270     (<span style="color: #8cd0d3;">sort</span> type-precedence [(<span style="color: #8cd0d3;">type</span> x) (<span style="color: #8cd0d3;">type</span> y)])))

   271 

   272 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">add</span> [java.lang.Number java.lang.Number] [x y] (<span style="color: #8cd0d3;">+</span> x y))

   273 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">multiply</span> [java.lang.Number java.lang.Number] [x y] (<span style="color: #8cd0d3;">*</span> x y))

   274 

   275 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">multiply</span> [incanter.Matrix java.lang.Integer] [x y]

   276            (<span style="color: #f0dfaf; font-weight: bold;">let</span> [args (<span style="color: #8cd0d3;">sort</span> #(type-precedence (<span style="color: #8cd0d3;">type</span> %1)(<span style="color: #8cd0d3;">type</span> %2)) [x y])

   277                  matrix (<span style="color: #8cd0d3;">first</span> args)

   278                  scalar (<span style="color: #8cd0d3;">second</span> args)]

   279              (incanter.core/matrix (<span style="color: #8cd0d3;">map</span> (<span style="color: #8cd0d3;">partial</span> map (<span style="color: #8cd0d3;">partial</span> multiply scalar)) matrix))))

   280            

   281 </pre>

   282 

   283 

   284 

   285 

   286 

   287 </div>

   288 

   289 </div>

   290 

   291 <div id="outline-container-2" class="outline-2">

   292 <h2 id="sec-2"><span class="section-number-2">2</span> Useful Data Types </h2>

   293 <div class="outline-text-2" id="text-2">

   294 

   295 

   296 

   297 </div>

   298 

   299 <div id="outline-container-2-1" class="outline-3">

   300 <h3 id="sec-2-1"><span class="section-number-3">2.1</span> Complex Numbers </h3>

   301 <div class="outline-text-3" id="text-2-1">

   302 

   303 

   304 

   305 

   306 <pre class="src src-clojure">(<span style="color: #f0dfaf; font-weight: bold;">in-ns</span> 'sicm.utils)

   307 

   308 (<span style="color: #f0dfaf; font-weight: bold;">defprotocol</span> <span style="color: #f0dfaf;">Complex</span>

   309   (real-part [z])

   310   (imaginary-part [z])

   311   (magnitude-squared [z])

   312   (angle [z])

   313   (conjugate [z])

   314   (norm [z]))

   315 

   316 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">complex-rectangular</span>

   317   <span style="color: #8fb28f;">"Define a complex number with the given real and imaginary</span>

   318 <span style="color: #8fb28f;">  components."</span>

   319   [re im]

   320   (<span style="color: #8cd0d3;">reify</span> Complex

   321          (real-part [z] re)

   322          (imaginary-part [z] im)

   323          (magnitude-squared [z] (<span style="color: #8cd0d3;">+</span> (<span style="color: #8cd0d3;">*</span> re re) (<span style="color: #8cd0d3;">*</span> im im)))

   324          (angle [z] (java.lang.Math/atan2 im re))

   325          (conjugate [z] (complex-rectangular re (<span style="color: #8cd0d3;">-</span> im)))

   326 

   327          Arithmetic

   328          (zero [z] (complex-rectangular 0 0))

   329          (one [z] (complex-rectangular 1 0))

   330          (negate [z] (complex-rectangular (<span style="color: #8cd0d3;">-</span> re) (<span style="color: #8cd0d3;">-</span> im)))

   331          (invert [z] (complex-rectangular

   332                       (<span style="color: #8cd0d3;">/</span> re (magnitude-squared z))

   333                       (<span style="color: #8cd0d3;">/</span> (<span style="color: #8cd0d3;">-</span> im) (magnitude-squared z))))

   334 

   335          Object

   336          (toString [_]

   337                    (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #f0dfaf; font-weight: bold;">and</span> (zero? re) (zero? im)) (<span style="color: #8cd0d3;">str</span> 0)

   338                        (<span style="color: #8cd0d3;">str</span>

   339                         (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #8cd0d3;">not</span>(zero? re))

   340                           re)

   341                         (<span style="color: #f0dfaf; font-weight: bold;">if</span> ((<span style="color: #8cd0d3;">comp</span> not zero?) im)

   342                           (<span style="color: #8cd0d3;">str</span>

   343                            (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #8cd0d3;">neg?</span> im) <span style="color: #cc9393;">"-"</span> <span style="color: #cc9393;">"+"</span>)

   344                            (<span style="color: #f0dfaf; font-weight: bold;">if</span> ((<span style="color: #8cd0d3;">comp</span> not one?) (java.lang.Math/abs im))

   345                            (java.lang.Math/abs im))

   346                            <span style="color: #cc9393;">"i"</span>)))))))

   347 

   348 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">complex-polar</span>

   349   <span style="color: #8fb28f;">"Define a complex number with the given magnitude and angle."</span>

   350   [mag ang]

   351   (<span style="color: #8cd0d3;">reify</span> Complex

   352          (magnitude-squared [z] (<span style="color: #8cd0d3;">*</span> mag mag))

   353          (angle [z] angle)

   354          (real-part [z] (<span style="color: #8cd0d3;">*</span> mag (java.lang.Math/cos ang)))

   355          (imaginary-part [z] (<span style="color: #8cd0d3;">*</span> mag (java.lang.Math/sin ang)))

   356          (conjugate [z] (complex-polar mag (<span style="color: #8cd0d3;">-</span> ang)))

   357 

   358          Arithmetic

   359          (zero [z] (complex-polar 0 0))

   360          (one [z] (complex-polar 1 0))

   361          (negate [z] (complex-polar (<span style="color: #8cd0d3;">-</span> mag) ang))

   362          (invert [z] (complex-polar (<span style="color: #8cd0d3;">/</span> mag) (<span style="color: #8cd0d3;">-</span> ang)))

   363 

   364          Object

   365          (toString [_] (<span style="color: #8cd0d3;">str</span> mag <span style="color: #cc9393;">" * e^(i"</span> ang<span style="color: #cc9393;">")"</span>))

   366          ))

   367 

   368 

   369 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">Numbers are complex quantities</span>

   370 

   371 (<span style="color: #8cd0d3;">extend-protocol</span> Complex

   372   java.lang.Number

   373   (real-part [x] x)

   374   (imaginary-part [x] 0)

   375   (magnitude [x] x)

   376   (angle [x] 0)

   377   (conjugate [x] x))

   378 

   379 

   380 </pre>

   381 

   382 

   383 

   384 

   385 

   386 

   387 </div>

   388 

   389 </div>

   390 

   391 <div id="outline-container-2-2" class="outline-3">

   392 <h3 id="sec-2-2"><span class="section-number-3">2.2</span> Tuples and Tensors </h3>

   393 <div class="outline-text-3" id="text-2-2">

   394 

   395 

   396 <p>

   397 A tuple is a vector which is spinable&mdash;it can be either <i>spin up</i> or <i>spin down</i>. (Covariant, contravariant; dual vectors)

   398 </p>

   399 

   400 

   401 <pre class="src src-clojure">(<span style="color: #f0dfaf; font-weight: bold;">in-ns</span> 'sicm.utils)

   402 

   403 (<span style="color: #f0dfaf; font-weight: bold;">defprotocol</span> <span style="color: #f0dfaf;">Spinning</span>

   404   (up? [this])

   405   (down? [this]))

   406 

   407 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">spin</span>

   408   <span style="color: #8fb28f;">"Returns the spin of the Spinning s, either :up or :down"</span>

   409   [<span style="color: #dfdfbf; font-weight: bold;">#^Spinning</span> s]

   410   (<span style="color: #f0dfaf; font-weight: bold;">cond</span> (up? s) <span style="color: #8cd0d3;">:up</span> (down? s) <span style="color: #8cd0d3;">:down</span>))

   411 

   412 

   413 (<span style="color: #f0dfaf; font-weight: bold;">deftype</span> <span style="color: #f0dfaf;">Tuple</span>

   414   [spin coll]

   415   clojure.lang.Seqable

   416   (<span style="color: #8cd0d3;">seq</span> [this] (<span style="color: #8cd0d3;">seq</span> (.coll this)))

   417   clojure.lang.Counted

   418   (<span style="color: #8cd0d3;">count</span> [this] (<span style="color: #8cd0d3;">count</span> (.coll this))))

   419 

   420 (<span style="color: #8cd0d3;">extend-type</span> Tuple

   421   Spinning

   422   (up? [this] (<span style="color: #8cd0d3;">=</span> <span style="color: #8cd0d3;">::up</span> (.spin this)))

   423   (down? [this] (<span style="color: #8cd0d3;">=</span> <span style="color: #8cd0d3;">::down</span> (.spin this))))

   424 

   425 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">print-method</span> Tuple

   426   [o w]

   427   (<span style="color: #8cd0d3;">print-simple</span> (<span style="color: #8cd0d3;">str</span> (<span style="color: #f0dfaf; font-weight: bold;">if</span> (up? o) 'u 'd) (.coll o))  w))

   428 

   429 

   430 

   431 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">up</span>

   432   <span style="color: #8fb28f;">"Create a new up-tuple containing the contents of coll."</span>

   433   [coll]

   434   (Tuple. <span style="color: #8cd0d3;">::up</span> coll))       

   435 

   436 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">down</span>

   437   <span style="color: #8fb28f;">"Create a new down-tuple containing the contents of coll."</span>

   438   [coll]

   439   (Tuple. <span style="color: #8cd0d3;">::down</span> coll))

   440 

   441 

   442 </pre>

   443 

   444 

   445 

   446 

   447 

   448 </div>

   449 

   450 <div id="outline-container-2-2-1" class="outline-4">

   451 <h4 id="sec-2-2-1"><span class="section-number-4">2.2.1</span> Contraction </h4>

   452 <div class="outline-text-4" id="text-2-2-1">

   453 

   454 <p>Contraction is a binary operation that you can apply to compatible

   455  tuples. Tuples are compatible for contraction if they have the same

   456  length and opposite spins, and if the corresponding items in each

   457  tuple are both numbers or both compatible tuples.

   458 </p>

   459 

   460 

   461 

   462 <pre class="src src-clojure">(<span style="color: #f0dfaf; font-weight: bold;">in-ns</span> 'sicm.utils)

   463 

   464 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">numbers?</span>

   465   <span style="color: #8fb28f;">"Returns true if all arguments are numbers, else false."</span>

   466   [&amp; xs]

   467   (<span style="color: #8cd0d3;">every?</span> number? xs))

   468 

   469 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">contractible?</span>

   470   <span style="color: #8fb28f;">"Returns true if the tuples a and b are compatible for contraction,</span>

   471 <span style="color: #8fb28f;">  else false. Tuples are compatible if they have the same number of</span>

   472 <span style="color: #8fb28f;">  components, they have opposite spins, and their elements are</span>

   473 <span style="color: #8fb28f;">  pairwise-compatible."</span>

   474   [a b]

   475   (<span style="color: #f0dfaf; font-weight: bold;">and</span>

   476    (<span style="color: #8cd0d3;">isa?</span> (<span style="color: #8cd0d3;">type</span> a) Tuple)

   477    (<span style="color: #8cd0d3;">isa?</span> (<span style="color: #8cd0d3;">type</span> b) Tuple)

   478    (<span style="color: #8cd0d3;">=</span> (<span style="color: #8cd0d3;">count</span> a) (<span style="color: #8cd0d3;">count</span> b))

   479    (<span style="color: #8cd0d3;">not=</span> (spin a) (spin b))

   480    

   481    (<span style="color: #8cd0d3;">not-any?</span> false?

   482              (<span style="color: #8cd0d3;">map</span> #(<span style="color: #f0dfaf; font-weight: bold;">or</span>

   483                     (numbers? %1 %2)

   484                     (contractible? %1 %2))

   485                   a b))))

   486 

   487 

   488 

   489 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">contract</span>

   490   <span style="color: #8fb28f;">"Contracts two tuples, returning the sum of the</span>

   491 <span style="color: #8fb28f;">  products of the corresponding items. Contraction is recursive on</span>

   492 <span style="color: #8fb28f;">  nested tuples."</span>

   493   [a b]

   494   (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #8cd0d3;">not</span> (contractible? a b))

   495     (<span style="color: #f0dfaf; font-weight: bold;">throw</span>

   496      (Exception. <span style="color: #cc9393;">"Not compatible for contraction."</span>))

   497     (<span style="color: #8cd0d3;">reduce</span> +

   498             (<span style="color: #8cd0d3;">map</span>

   499              (<span style="color: #8cd0d3;">fn</span> [x y]

   500                (<span style="color: #f0dfaf; font-weight: bold;">if</span> (numbers? x y)

   501                  (<span style="color: #8cd0d3;">*</span> x y)

   502                  (contract x y)))

   503              a b))))

   504 

   505 </pre>

   506 

   507 

   508 

   509 

   510 </div>

   511 

   512 </div>

   513 

   514 <div id="outline-container-2-2-2" class="outline-4">

   515 <h4 id="sec-2-2-2"><span class="section-number-4">2.2.2</span> Matrices </h4>

   516 <div class="outline-text-4" id="text-2-2-2">

   517 

   518 

   519 

   520 

   521 <pre class="src src-clojure">(<span style="color: #f0dfaf; font-weight: bold;">in-ns</span> 'sicm.utils)

   522 (<span style="color: #8cd0d3;">require</span> 'incanter.core) <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">use incanter's fast matrices</span>

   523 

   524 (<span style="color: #f0dfaf; font-weight: bold;">defprotocol</span> <span style="color: #f0dfaf;">Matrix</span>

   525   (rows [matrix])

   526   (cols [matrix])

   527   (diagonal [matrix])

   528   (trace [matrix])

   529   (determinant [matrix])

   530   (transpose [matrix])

   531   (conjugate [matrix])

   532 )

   533 

   534 (<span style="color: #8cd0d3;">extend-protocol</span> Matrix

   535   incanter.Matrix

   536   (rows [rs] (<span style="color: #8cd0d3;">map</span> down (<span style="color: #8cd0d3;">apply</span> map vector (<span style="color: #8cd0d3;">apply</span> map vector rs))))

   537   (cols [rs] (<span style="color: #8cd0d3;">map</span> up (<span style="color: #8cd0d3;">apply</span> map vector rs)))

   538   (diagonal [matrix] (incanter.core/diag matrix) )

   539   (determinant [matrix] (incanter.core/det matrix))

   540   (trace [matrix] (incanter.core/trace matrix))

   541   (transpose [matrix] (incanter.core/trans matrix))

   542   )

   543 

   544 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">count-rows</span> [matrix]

   545   ((<span style="color: #8cd0d3;">comp</span> count rows) matrix))

   546 

   547 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">count-cols</span> [matrix]

   548   ((<span style="color: #8cd0d3;">comp</span> count cols) matrix))

   549 

   550 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">square?</span> [matrix]

   551   (<span style="color: #8cd0d3;">=</span> (count-rows matrix) (count-cols matrix)))

   552 

   553 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">identity-matrix</span>

   554   <span style="color: #8fb28f;">"Define a square matrix of size n-by-n with 1s along the diagonal and</span>

   555 <span style="color: #8fb28f;">  0s everywhere else."</span>

   556   [n]

   557   (incanter.core/identity-matrix n))

   558 

   559 

   560 

   561 

   562 

   563 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">matrix-by-rows</span>

   564   <span style="color: #8fb28f;">"Define a matrix by giving its rows."</span>

   565   [&amp; rows]

   566   (<span style="color: #f0dfaf; font-weight: bold;">if</span>

   567    (<span style="color: #8cd0d3;">not</span> (all-equal? (<span style="color: #8cd0d3;">map</span> count rows)))

   568    (<span style="color: #f0dfaf; font-weight: bold;">throw</span> (Exception. <span style="color: #cc9393;">"All rows in a matrix must have the same number of elements."</span>))

   569    (incanter.core/matrix (<span style="color: #8cd0d3;">vec</span> rows))))

   570 

   571 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">matrix-by-cols</span>

   572   <span style="color: #8fb28f;">"Define a matrix by giving its columns"</span>

   573   [&amp; cols]

   574   (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #8cd0d3;">not</span> (all-equal? (<span style="color: #8cd0d3;">map</span> count cols)))

   575    (<span style="color: #f0dfaf; font-weight: bold;">throw</span> (Exception. <span style="color: #cc9393;">"All columns in a matrix must have the same number of elements."</span>))

   576    (incanter.core/matrix (<span style="color: #8cd0d3;">vec</span> (<span style="color: #8cd0d3;">apply</span> map vector cols)))))

   577 

   578 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">identity-matrix</span>

   579   <span style="color: #8fb28f;">"Define a square matrix of size n-by-n with 1s along the diagonal and</span>

   580 <span style="color: #8fb28f;">  0s everywhere else."</span>

   581   [n]

   582   (incanter.core/identity-matrix n))

   583 

   584 

   585 

   586 (<span style="color: #8cd0d3;">extend-protocol</span> Arithmetic

   587   incanter.Matrix

   588   (one [matrix]

   589        (<span style="color: #f0dfaf; font-weight: bold;">if</span> (square? matrix)

   590          (identity-matrix (count-rows matrix))

   591          (<span style="color: #f0dfaf; font-weight: bold;">throw</span> (Exception. <span style="color: #cc9393;">"Non-square matrices have no multiplicative unit."</span>))))

   592   (zero [matrix]

   593         (<span style="color: #8cd0d3;">apply</span> matrix-by-rows (<span style="color: #8cd0d3;">map</span> zero (rows matrix))))

   594   (negate [matrix]

   595           (<span style="color: #8cd0d3;">apply</span> matrix-by-rows (<span style="color: #8cd0d3;">map</span> negate (rows matrix))))

   596   (invert [matrix]

   597           (incanter.core/solve matrix)))

   598 

   599 

   600 

   601 (<span style="color: #f0dfaf; font-weight: bold;">defmulti</span> <span style="color: #f0dfaf;">coerce-to-matrix</span>

   602  <span style="color: #8fb28f;">"Converts x into a matrix, if possible."</span>

   603  type)

   604 

   605 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">coerce-to-matrix</span> incanter.Matrix [x] x)

   606 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">coerce-to-matrix</span> Tuple [x]

   607            (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #8cd0d3;">apply</span> numbers? (<span style="color: #8cd0d3;">seq</span> x))

   608              (<span style="color: #f0dfaf; font-weight: bold;">if</span> (up? x)

   609              (matrix-by-cols (<span style="color: #8cd0d3;">seq</span> x))

   610              (matrix-by-rows (<span style="color: #8cd0d3;">seq</span> x)))

   611              (<span style="color: #f0dfaf; font-weight: bold;">throw</span> (Exception. <span style="color: #cc9393;">"Non-numerical tuple cannot be converted into a matrix."</span>)))) 

   612              

   613   

   614 

   615 

   616 

   617 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">(defn matrix-by-cols</span>

   618 <span style="color: #708070;">;;   </span><span style="color: #7f9f7f;">"Define a matrix by giving its columns."</span>

   619 <span style="color: #708070;">;;   </span><span style="color: #7f9f7f;">[&amp; cols]</span>

   620 <span style="color: #708070;">;;   </span><span style="color: #7f9f7f;">(cond</span>

   621 <span style="color: #708070;">;;    </span><span style="color: #7f9f7f;">(not (all-equal? (map count cols)))</span>

   622 <span style="color: #708070;">;;    </span><span style="color: #7f9f7f;">(throw (Exception. "All columns in a matrix must have the same number of elements."))</span>

   623 <span style="color: #708070;">;;    </span><span style="color: #7f9f7f;">:else</span>

   624 <span style="color: #708070;">;;    </span><span style="color: #7f9f7f;">(reify Matrix</span>

   625 <span style="color: #708070;">;;        </span><span style="color: #7f9f7f;">(cols [this] (map up cols))</span>

   626 <span style="color: #708070;">;;        </span><span style="color: #7f9f7f;">(rows [this] (map down (apply map vector cols)))</span>

   627 <span style="color: #708070;">;;        </span><span style="color: #7f9f7f;">(diagonal [this] (map-indexed (fn [i col] (nth col i) cols)))</span>

   628 <span style="color: #708070;">;;        </span><span style="color: #7f9f7f;">(trace [this]</span>

   629 <span style="color: #708070;">;;               </span><span style="color: #7f9f7f;">(if (not= (count-cols this) (count-rows this))</span>

   630 <span style="color: #708070;">;;                 </span><span style="color: #7f9f7f;">(throw (Exception.</span>

   631 <span style="color: #708070;">;;                         </span><span style="color: #7f9f7f;">"Cannot take the trace of a non-square matrix."))</span>

   632 <span style="color: #708070;">;;                 </span><span style="color: #7f9f7f;">(reduce + (diagonal this))))</span>

   633           

   634 <span style="color: #708070;">;;        </span><span style="color: #7f9f7f;">(determinant [this]</span>

   635 <span style="color: #708070;">;;                     </span><span style="color: #7f9f7f;">(if (not= (count-cols this) (count-rows this))</span>

   636 <span style="color: #708070;">;;                       </span><span style="color: #7f9f7f;">(throw (Exception.</span>

   637 <span style="color: #708070;">;;                               </span><span style="color: #7f9f7f;">"Cannot take the determinant of a non-square matrix."))</span>

   638 <span style="color: #708070;">;;                       </span><span style="color: #7f9f7f;">(reduce * (map-indexed (fn [i col] (nth col i)) cols))))</span>

   639 <span style="color: #708070;">;;        </span><span style="color: #7f9f7f;">)))</span>

   640 

   641 (<span style="color: #8cd0d3;">extend-protocol</span> Matrix  Tuple

   642                  (rows [this] (<span style="color: #f0dfaf; font-weight: bold;">if</span> (down? this)

   643                                 (<span style="color: #8cd0d3;">list</span> this)

   644                                 (<span style="color: #8cd0d3;">map</span> (<span style="color: #8cd0d3;">comp</span> up vector) this)))

   645 

   646                  (cols [this] (<span style="color: #f0dfaf; font-weight: bold;">if</span> (up? this)

   647                                 (<span style="color: #8cd0d3;">list</span> this)

   648                                 (<span style="color: #8cd0d3;">map</span> (<span style="color: #8cd0d3;">comp</span> down vector) this))

   649                  ))

   650   

   651 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">matrix-multiply</span>

   652   <span style="color: #8fb28f;">"Returns the matrix resulting from the matrix multiplication of the given arguments."</span>

   653   ([A] (coerce-to-matrix A))

   654   ([A B] (incanter.core/mmult (coerce-to-matrix A) (coerce-to-matrix B)))

   655   ([M1 M2 &amp; Ms] (<span style="color: #8cd0d3;">reduce</span> matrix-multiply (matrix-multiply M1 M2) Ms)))

   656 

   657 </pre>

   658 

   659 

   660 

   661 

   662 

   663 </div>

   664 </div>

   665 

   666 </div>

   667 

   668 <div id="outline-container-2-3" class="outline-3">

   669 <h3 id="sec-2-3"><span class="section-number-3">2.3</span> Power Series </h3>

   670 <div class="outline-text-3" id="text-2-3">

   671 

   672 

   673 

   674 

   675 <pre class="src src-clojure">(<span style="color: #f0dfaf; font-weight: bold;">in-ns</span> 'sicm.utils)

   676 (<span style="color: #8cd0d3;">use</span> 'clojure.contrib.def)

   677 

   678 

   679 

   680 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">series-fn</span>

   681   <span style="color: #8fb28f;">"The function corresponding to the given power series."</span>

   682   [series]

   683   (<span style="color: #8cd0d3;">fn</span> [x]

   684     (<span style="color: #8cd0d3;">reduce</span> +

   685             (map-indexed  (<span style="color: #8cd0d3;">fn</span>[n x] (<span style="color: #8cd0d3;">*</span> (<span style="color: #8cd0d3;">float</span> (<span style="color: #8cd0d3;">nth</span> series n)) (<span style="color: #8cd0d3;">float</span>(java.lang.Math/pow (<span style="color: #8cd0d3;">float</span> x) n)) ))

   686                           (<span style="color: #8cd0d3;">range</span> 20)))))

   687 

   688 (<span style="color: #f0dfaf; font-weight: bold;">deftype</span> <span style="color: #f0dfaf;">PowerSeries</span>

   689   [coll]

   690   clojure.lang.Seqable

   691   (<span style="color: #8cd0d3;">seq</span> [this] (<span style="color: #8cd0d3;">seq</span> (.coll this)))

   692   

   693   clojure.lang.Indexed

   694   (<span style="color: #8cd0d3;">nth</span> [this n] (<span style="color: #8cd0d3;">nth</span> (.coll this) n 0))

   695   (<span style="color: #8cd0d3;">nth</span> [this n not-found] (<span style="color: #8cd0d3;">nth</span> (.coll this) n not-found))

   696 

   697   <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">clojure.lang.IFn</span>

   698   <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">(call [this] (throw(Exception.)))</span>

   699   <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">(invoke [this &amp; args] args</span>

   700   <span style="color: #708070;">;;      </span><span style="color: #7f9f7f;">(let [f </span>

   701   <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">)</span>

   702   <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">(run [this] (throw(Exception.)))</span>

   703   )

   704 

   705 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">power-series</span>

   706   <span style="color: #8fb28f;">"Returns a power series with the items of the coll as its</span>

   707 <span style="color: #8fb28f;">  coefficients. Trailing zeros are added to the end of coll."</span>

   708   [coeffs]

   709   (PowerSeries. coeffs))

   710 

   711 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">power-series-indexed</span>

   712   <span style="color: #8fb28f;">"Returns a power series consisting of the result of mapping f to the non-negative integers."</span>

   713   [f]

   714   (PowerSeries. (<span style="color: #8cd0d3;">map</span> f (<span style="color: #8cd0d3;">range</span>))))

   715 

   716 

   717 (<span style="color: #f0dfaf; font-weight: bold;">defn-memo</span> <span style="color: #f0dfaf;">nth-partial-sum</span>

   718   ([series n]

   719      (<span style="color: #f0dfaf; font-weight: bold;">if</span> (zero? n) (<span style="color: #8cd0d3;">first</span> series)

   720          (<span style="color: #8cd0d3;">+</span> (<span style="color: #8cd0d3;">nth</span> series n)

   721             (nth-partial-sum series (<span style="color: #8cd0d3;">dec</span> n))))))

   722 

   723 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">partial-sums</span> [series]

   724   (<span style="color: #8cd0d3;">lazy-seq</span> (<span style="color: #8cd0d3;">map</span> nth-partial-sum (<span style="color: #8cd0d3;">range</span>))))

   725 

   726 

   727 

   728 

   729 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">cos-series</span>

   730      (power-series-indexed

   731       (<span style="color: #8cd0d3;">fn</span>[n]

   732         (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #8cd0d3;">odd?</span> n) 0

   733             (<span style="color: #8cd0d3;">/</span>

   734              (<span style="color: #8cd0d3;">reduce</span> *

   735                      (<span style="color: #8cd0d3;">reduce</span> * (<span style="color: #8cd0d3;">repeat</span> (<span style="color: #8cd0d3;">/</span> n 2) -1))

   736                      (<span style="color: #8cd0d3;">range</span> 1 (<span style="color: #8cd0d3;">inc</span> n)))

   737              )))))

   738 

   739 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">sin-series</span>

   740      (power-series-indexed

   741       (<span style="color: #8cd0d3;">fn</span>[n]

   742         (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #8cd0d3;">even?</span> n) 0

   743             (<span style="color: #8cd0d3;">/</span>

   744              (<span style="color: #8cd0d3;">reduce</span> *

   745                      (<span style="color: #8cd0d3;">reduce</span> * (<span style="color: #8cd0d3;">repeat</span> (<span style="color: #8cd0d3;">/</span> (<span style="color: #8cd0d3;">dec</span> n) 2) -1))

   746                      (<span style="color: #8cd0d3;">range</span> 1 (<span style="color: #8cd0d3;">inc</span> n)))

   747              )))))

   748 

   749 </pre>

   750 

   751 

   752 

   753 

   754 

   755 </div>

   756 </div>

   757 

   758 </div>

   759 

   760 <div id="outline-container-3" class="outline-2">

   761 <h2 id="sec-3"><span class="section-number-2">3</span> Basic Utilities </h2>

   762 <div class="outline-text-2" id="text-3">

   763 

   764 

   765 

   766 </div>

   767 

   768 <div id="outline-container-3-1" class="outline-3">

   769 <h3 id="sec-3-1"><span class="section-number-3">3.1</span> Sequence manipulation </h3>

   770 <div class="outline-text-3" id="text-3-1">

   771 

   772 

   773 

   774 

   775 

   776 <pre class="src src-clojure">(<span style="color: #f0dfaf; font-weight: bold;">ns</span> sicm.utils)

   777 

   778 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">do-up</span>

   779   <span style="color: #8fb28f;">"Apply f to each number from low to high, presumably for</span>

   780 <span style="color: #8fb28f;">  side-effects."</span>

   781   [f low high]

   782   (<span style="color: #f0dfaf; font-weight: bold;">doseq</span> [i (<span style="color: #8cd0d3;">range</span> low high)] (f i)))

   783    

   784 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">do-down</span>

   785   <span style="color: #8fb28f;">"Apply f to each number from high to low, presumably for</span>

   786 <span style="color: #8fb28f;">   side-effects."</span>

   787   [f high low]

   788   (<span style="color: #f0dfaf; font-weight: bold;">doseq</span> [i (<span style="color: #8cd0d3;">range</span> high low -1)] (f i)))

   789 

   790 

   791 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">all-equal?</span> [coll]

   792   (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #8cd0d3;">empty?</span> (<span style="color: #8cd0d3;">rest</span> coll)) true

   793       (<span style="color: #f0dfaf; font-weight: bold;">and</span> (<span style="color: #8cd0d3;">=</span> (<span style="color: #8cd0d3;">first</span> coll) (<span style="color: #8cd0d3;">second</span> coll))

   794            (<span style="color: #f0dfaf; font-weight: bold;">recur</span> (<span style="color: #8cd0d3;">rest</span> coll))))))

   795 

   796 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">multiplier</span>

   797   <span style="color: #8fb28f;">"Returns a function that 'multiplies' the members of a collection,</span>

   798 <span style="color: #8fb28f;">returning unit if called on an empty collection."</span>

   799   [multiply unit]

   800   (<span style="color: #8cd0d3;">fn</span> [coll] ((<span style="color: #8cd0d3;">partial</span> reduce multiply unit) coll)))

   801 

   802 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">divider</span>

   803   <span style="color: #8fb28f;">"Returns a function that 'divides' the first element of a collection</span>

   804 <span style="color: #8fb28f;">by the 'product' of the rest of the collection."</span>

   805   [divide multiply invert unit]

   806   (<span style="color: #8cd0d3;">fn</span> [coll]

   807     (<span style="color: #8cd0d3;">apply</span>

   808      (<span style="color: #8cd0d3;">fn</span>

   809        ([] unit)

   810        ([x] (invert x))

   811        ([x y] (divide x y))

   812        ([x y &amp; zs] (divide x (<span style="color: #8cd0d3;">reduce</span> multiply y zs))))

   813      coll)))

   814 

   815 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">left-circular-shift</span>

   816   <span style="color: #8fb28f;">"Remove the first element of coll, adding it to the end of coll."</span>

   817   [coll]

   818   (<span style="color: #8cd0d3;">concat</span> (<span style="color: #8cd0d3;">rest</span> coll) (<span style="color: #8cd0d3;">take</span> 1 coll)))

   819 

   820 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">right-circular-shift</span>

   821   <span style="color: #8fb28f;">"Remove the last element of coll, adding it to the front of coll."</span>

   822   [coll]

   823   (<span style="color: #8cd0d3;">cons</span> (<span style="color: #8cd0d3;">last</span> coll) (<span style="color: #8cd0d3;">butlast</span> coll)))

   824 </pre>

   825 

   826 

   827 

   828 

   829 

   830 

   831 

   832 </div>

   833 

   834 </div>

   835 

   836 <div id="outline-container-3-2" class="outline-3">

   837 <h3 id="sec-3-2"><span class="section-number-3">3.2</span> Ranges, Airity and Function Composition </h3>

   838 <div class="outline-text-3" id="text-3-2">

   839 

   840 

   841 

   842 

   843 <pre class="src src-clojure">(<span style="color: #f0dfaf; font-weight: bold;">in-ns</span> 'sicm.utils)

   844 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">infinity</span> Double/POSITIVE_INFINITY)

   845 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">infinite?</span> [x] (Double/isInfinite x))

   846 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">finite?</span> (<span style="color: #8cd0d3;">comp</span> not infinite?)) 

   847 

   848 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">arity-min</span>

   849   <span style="color: #8fb28f;">"Returns the smallest number of arguments f can take."</span>

   850   [f]

   851   (<span style="color: #8cd0d3;">apply</span>

   852    min

   853    (<span style="color: #8cd0d3;">map</span> (<span style="color: #8cd0d3;">comp</span> alength #(.getParameterTypes %))

   854         (<span style="color: #8cd0d3;">filter</span> (<span style="color: #8cd0d3;">comp</span> (<span style="color: #8cd0d3;">partial</span> = <span style="color: #cc9393;">"invoke"</span>) #(.getName %))

   855                 (.getDeclaredMethods (<span style="color: #8cd0d3;">class</span> f))))))

   856   

   857 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">arity-max</span>

   858   <span style="color: #8fb28f;">"Returns the largest number of arguments f can take, possibly</span>

   859 <span style="color: #8fb28f;">  Infinity."</span>

   860   [f]

   861   (<span style="color: #f0dfaf; font-weight: bold;">let</span> [methods (.getDeclaredMethods (<span style="color: #8cd0d3;">class</span> f))]

   862     (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #8cd0d3;">not-any?</span> (<span style="color: #8cd0d3;">partial</span> = <span style="color: #cc9393;">"doInvoke"</span>) (<span style="color: #8cd0d3;">map</span> #(.getName %) methods))

   863       (<span style="color: #8cd0d3;">apply</span> max

   864              (<span style="color: #8cd0d3;">map</span> (<span style="color: #8cd0d3;">comp</span> alength #(.getParameterTypes %))

   865                   (<span style="color: #8cd0d3;">filter</span> (<span style="color: #8cd0d3;">comp</span> (<span style="color: #8cd0d3;">partial</span> = <span style="color: #cc9393;">"invoke"</span>) #(.getName %))  methods)))

   866       infinity)))

   867 

   868 

   869 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">^</span>{<span style="color: #8cd0d3;">:arglists</span> '([f])

   870        <span style="color: #8cd0d3;">:doc</span> <span style="color: #cc9393;">"Returns a two-element list containing the minimum and</span>

   871 <span style="color: #cc9393;">     maximum number of args that f can take."</span>}

   872      arity-interval

   873      (<span style="color: #8cd0d3;">juxt</span> arity-min arity-max))

   874 

   875 

   876 

   877 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">--- intervals</span>

   878 

   879 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">intersect</span>

   880   <span style="color: #8fb28f;">"Returns the interval of overlap between interval-1 and interval-2"</span>

   881   [interval-1 interval-2]

   882   (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #f0dfaf; font-weight: bold;">or</span> (<span style="color: #8cd0d3;">empty?</span> interval-1) (<span style="color: #8cd0d3;">empty?</span> interval-2)) []

   883       (<span style="color: #f0dfaf; font-weight: bold;">let</span> [left (<span style="color: #8cd0d3;">max</span> (<span style="color: #8cd0d3;">first</span> interval-1) (<span style="color: #8cd0d3;">first</span> interval-2))

   884             right (<span style="color: #8cd0d3;">min</span> (<span style="color: #8cd0d3;">second</span> interval-1) (<span style="color: #8cd0d3;">second</span> interval-2))]

   885         (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #8cd0d3;">&gt;</span> left right) []

   886             [left right]))))

   887 

   888 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">same-endpoints?</span>

   889   <span style="color: #8fb28f;">"Returns true if the left endpoint is the same as the right</span>

   890 <span style="color: #8fb28f;">  endpoint."</span>

   891   [interval]

   892   (<span style="color: #8cd0d3;">=</span> (<span style="color: #8cd0d3;">first</span> interval) (<span style="color: #8cd0d3;">second</span> interval)))

   893 

   894 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">naturals?</span>

   895   <span style="color: #8fb28f;">"Returns true if the left endpoint is 0 and the right endpoint is</span>

   896 <span style="color: #8fb28f;">infinite."</span>

   897   [interval]

   898   (<span style="color: #f0dfaf; font-weight: bold;">and</span> (zero? (<span style="color: #8cd0d3;">first</span> interval))

   899        (infinite? (<span style="color: #8cd0d3;">second</span> interval))))

   900  

   901 

   902 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">fan-in</span>

   903   <span style="color: #8fb28f;">"Returns a function that pipes its input to each of the gs, then</span>

   904 <span style="color: #8fb28f;">  applies f to the list of results. Consequently, f must be able to</span>

   905 <span style="color: #8fb28f;">  take a number of arguments equal to the number of gs."</span>

   906   [f &amp; gs]

   907   (<span style="color: #8cd0d3;">fn</span> [&amp; args]

   908     (<span style="color: #8cd0d3;">apply</span> f (<span style="color: #8cd0d3;">apply</span> (<span style="color: #8cd0d3;">apply</span> juxt gs) args))))

   909 

   910 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">fan-in</span>

   911   <span style="color: #8fb28f;">"Returns a function that pipes its input to each of the gs, then applies f to the list of results. The resulting function takes any number of arguments, but will fail if given arguments that are incompatible with any of the gs."</span>

   912   [f &amp; gs]

   913   (<span style="color: #8cd0d3;">comp</span> (<span style="color: #8cd0d3;">partial</span> apply f) (<span style="color: #8cd0d3;">apply</span> juxt gs)))

   914 

   915 

   916 

   917 (<span style="color: #f0dfaf; font-weight: bold;">defmacro</span> <span style="color: #f0dfaf;">airty-blah-sad</span> [f n more?]

   918   (<span style="color: #f0dfaf; font-weight: bold;">let</span> [syms (<span style="color: #8cd0d3;">vec</span> (<span style="color: #8cd0d3;">map</span> (<span style="color: #8cd0d3;">comp</span> gensym (<span style="color: #8cd0d3;">partial</span> str <span style="color: #cc9393;">"x"</span>)) (<span style="color: #8cd0d3;">range</span> n)))

   919          optional (<span style="color: #8cd0d3;">gensym</span> <span style="color: #cc9393;">"xs"</span>)]

   920     (<span style="color: #f0dfaf; font-weight: bold;">if</span> more?

   921       `(<span style="color: #8cd0d3;">fn</span> <span style="color: #f0dfaf;">~</span>(<span style="color: #8cd0d3;">conj</span> syms '&amp; optional)

   922          (<span style="color: #8cd0d3;">apply</span> ~f ~@syms ~optional))

   923       `(<span style="color: #8cd0d3;">fn</span> <span style="color: #f0dfaf;">~syms</span> (~f ~@syms)))))

   924 

   925 (<span style="color: #f0dfaf; font-weight: bold;">defmacro</span> <span style="color: #f0dfaf;">airt-whaa*</span> [f n more?]

   926   `(airty-blah-sad ~f ~n ~more?))

   927 

   928 

   929 

   930 

   931 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">fan-in*</span>

   932   <span style="color: #8fb28f;">"Returns a function that pipes its input to each of the gs, then</span>

   933 <span style="color: #8fb28f;">  applies f to the list of results. Unlike fan-in, fan-in* strictly</span>

   934 <span style="color: #8fb28f;">  enforces arity: it will fail if the gs do not have compatible</span>

   935 <span style="color: #8fb28f;">  arities."</span>

   936   [f &amp; gs]

   937   (<span style="color: #f0dfaf; font-weight: bold;">let</span> [arity-in (<span style="color: #8cd0d3;">reduce</span> intersect (<span style="color: #8cd0d3;">map</span> arity-interval gs))

   938         left (<span style="color: #8cd0d3;">first</span> arity-in)

   939         right (<span style="color: #8cd0d3;">second</span> arity-in)

   940         composite (fan-in f gs)

   941         ]

   942     (<span style="color: #f0dfaf; font-weight: bold;">cond</span>

   943      (<span style="color: #8cd0d3;">empty?</span> arity-in)

   944      (<span style="color: #f0dfaf; font-weight: bold;">throw</span> (Exception. <span style="color: #cc9393;">"Cannot compose functions with incompatible arities."</span>))

   945 

   946      (<span style="color: #8cd0d3;">not</span>

   947       (<span style="color: #f0dfaf; font-weight: bold;">or</span> (<span style="color: #8cd0d3;">=</span> left right)

   948           (<span style="color: #f0dfaf; font-weight: bold;">and</span> (finite? left)

   949                (<span style="color: #8cd0d3;">=</span> right infinity))))

   950          

   951      (<span style="color: #f0dfaf; font-weight: bold;">throw</span> (Exception.

   952              <span style="color: #cc9393;">"Compose can only handle arities of the form [n n] or [n infinity]"</span>))

   953      <span style="color: #8cd0d3;">:else</span>

   954      (airty-blah-sad composite left (<span style="color: #8cd0d3;">=</span> right infinity)))))

   955     

   956 

   957 

   958 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">compose-n</span> <span style="color: #8fb28f;">"Compose any number of functions together."</span>

   959   ([] identity)

   960   ([f] f)

   961   ([f &amp; fs]

   962   (<span style="color: #f0dfaf; font-weight: bold;">let</span> [fns (<span style="color: #8cd0d3;">cons</span> f fs)]

   963     (compose-bin (<span style="color: #8cd0d3;">reduce</span> fan-in (<span style="color: #8cd0d3;">butlast</span> fs)) (<span style="color: #8cd0d3;">last</span> fs))))

   964 )

   965 

   966 

   967 

   968 

   969            

   970 

   971 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">iterated</span>

   972   ([f n id] (<span style="color: #8cd0d3;">reduce</span> comp id (<span style="color: #8cd0d3;">repeat</span> n f)))

   973   ([f n] (<span style="color: #8cd0d3;">reduce</span> comp identity (<span style="color: #8cd0d3;">repeat</span> n f))))

   974 

   975 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">iterate-until-stable</span>

   976   <span style="color: #8fb28f;">"Repeatedly applies f to x, returning the first result that is close</span>

   977 <span style="color: #8fb28f;">enough to its predecessor."</span>

   978   [f close-enough? x]

   979   (<span style="color: #8cd0d3;">second</span> (swank.util/find-first

   980            (<span style="color: #8cd0d3;">partial</span> apply close-enough?)

   981            (<span style="color: #8cd0d3;">partition</span> 2 1 (<span style="color: #8cd0d3;">iterate</span> f x)))))

   982 

   983 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">lexical&lt;</span> [x y]

   984   (<span style="color: #8cd0d3;">neg?</span> (<span style="color: #8cd0d3;">compare</span> (<span style="color: #8cd0d3;">str</span> x) (<span style="color: #8cd0d3;">str</span> y))))

   985 

   986 

   987 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">do-up</span>

   988 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">do-down</span>

   989 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">make-pairwise-test</span> comparator)

   990 <span style="color: #708070;">;;</span><span style="color: #7f9f7f;">all-equal?</span>

   991 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">accumulation</span> multiplier)

   992 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">inverse-accumulation</span> divider)

   993 <span style="color: #708070;">;;</span><span style="color: #7f9f7f;">left-circular-shift</span>

   994 <span style="color: #708070;">;;</span><span style="color: #7f9f7f;">right-circular-shift</span>

   995 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">exactly-n?</span> same-endpoints?)

   996 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">any-number?</span> naturals?)

   997 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">TODO compose</span>

   998 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">TODO compose-n</span>

   999 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">identity</span>

  1000 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">compose-2</span> fan-in)

  1001 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">compose-bin</span> fan-in*)

  1002 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">any?</span> (<span style="color: #8cd0d3;">constantly</span> true))

  1003 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">none?</span> (<span style="color: #8cd0d3;">constantly</span> false))

  1004 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">constant</span> constantly)

  1005 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">joint-arity</span> intersect)

  1006 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">a-reduce</span> reduce)

  1007 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">filter</span>

  1008 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">make-map</span> (<span style="color: #8cd0d3;">partial</span> partial map) )

  1009 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">bracket</span> juxt)

  1010 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">TODO apply-to-all</span>

  1011 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">TODO nary-combine</span>

  1012 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">TODO binary-combine</span>

  1013 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">TODO unary-combine</span>

  1014 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">iterated</span>

  1015 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">iterate-until-stable</span>

  1016 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">make-function-of-vector</span> (<span style="color: #8cd0d3;">partial</span> partial map))

  1017 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">make-function-of-arguments</span> (<span style="color: #8cd0d3;">fn</span> [f] (<span style="color: #8cd0d3;">fn</span> [&amp; args] (f args))))

  1018 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">alphaless</span> lexical&lt;)

  1019 

  1020 </pre>

  1021 

  1022 

  1023 

  1024 

  1025 

  1026 

  1027 

  1028 

  1029 

  1030 

  1031 

  1032 

  1033 </div>

  1034 </div>

  1035 

  1036 </div>

  1037 

  1038 <div id="outline-container-4" class="outline-2">

  1039 <h2 id="sec-4"><span class="section-number-2">4</span> Numerical Methods </h2>

  1040 <div class="outline-text-2" id="text-4">

  1041 

  1042 

  1043 

  1044 

  1045 <pre class="src src-clojure">(<span style="color: #f0dfaf; font-weight: bold;">in-ns</span> 'sicm.utils)

  1046 (<span style="color: #f0dfaf; font-weight: bold;">import</span> java.lang.Math)

  1047 (<span style="color: #8cd0d3;">use</span> 'clojure.contrib.def)

  1048 

  1049 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">---- USEFUL CONSTANTS</span>

  1050 

  1051 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">machine-epsilon</span>

  1052   <span style="color: #8fb28f;">"Smallest value representable on your machine, as determined by</span>

  1053 <span style="color: #8fb28f;">successively dividing a number in half until consecutive results are</span>

  1054 <span style="color: #8fb28f;">indistinguishable."</span>

  1055   []

  1056   (<span style="color: #8cd0d3;">ffirst</span>

  1057    (<span style="color: #8cd0d3;">drop-while</span>

  1058     (<span style="color: #8cd0d3;">comp</span> not zero? second)

  1059     (<span style="color: #8cd0d3;">partition</span> 2 1

  1060                (<span style="color: #8cd0d3;">iterate</span> (<span style="color: #8cd0d3;">partial</span> * 0.5) 1)))))

  1061 

  1062 

  1063 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">pi</span> (Math/PI))

  1064 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">two-pi</span> (<span style="color: #8cd0d3;">*</span> 2 pi))

  1065 

  1066 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">eulers-gamma</span> 0.5772156649015328606065)

  1067 

  1068 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">phi</span> (<span style="color: #8cd0d3;">/</span> (<span style="color: #8cd0d3;">inc</span> (Math/sqrt 5)) 2))

  1069 

  1070 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">ln2</span> (Math/log 2))

  1071 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">ln10</span> (Math/log 10))

  1072 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">exp10</span> #(Math/pow 10 %))

  1073 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">exp2</span> #(Math/pow 2 %))

  1074 

  1075 

  1076 <span style="color: #708070;">;;</span>

  1077 

  1078 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">---- ANGLES AND TRIGONOMETRY</span>

  1079 

  1080 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">angle-restrictor</span>

  1081   <span style="color: #8fb28f;">"Returns a function that ensures that angles lie in the specified interval of length two-pi."</span>

  1082   [max-angle]

  1083   (<span style="color: #f0dfaf; font-weight: bold;">let</span> [min-angle (<span style="color: #8cd0d3;">-</span> max-angle two-pi)]

  1084     (<span style="color: #8cd0d3;">fn</span> [x]

  1085       (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #f0dfaf; font-weight: bold;">and</span>

  1086            (<span style="color: #8cd0d3;">&lt;=</span> min-angle x)

  1087            (<span style="color: #8cd0d3;">&lt;</span> x max-angle))

  1088         x

  1089         (<span style="color: #f0dfaf; font-weight: bold;">let</span> [corrected-x (<span style="color: #8cd0d3;">-</span> x (<span style="color: #8cd0d3;">*</span> two-pi (Math/floor (<span style="color: #8cd0d3;">/</span> x two-pi))))]

  1090           (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #8cd0d3;">&lt;</span> corrected-x max-angle)

  1091             corrected-x

  1092             (<span style="color: #8cd0d3;">-</span> corrected-x two-pi)))))))

  1093 

  1094 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">angle-restrict-pi</span>

  1095   <span style="color: #8fb28f;">"Coerces angles to lie in the interval from -pi to pi."</span>

  1096   [angle]

  1097   ((angle-restrictor pi) angle))

  1098 

  1099 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">angle-restrict-two-pi</span>

  1100   <span style="color: #8fb28f;">"Coerces angles to lie in the interval from zero to two-pi"</span>

  1101   [angle]

  1102   ((angle-restrictor two-pi) angle))

  1103 

  1104 

  1105 

  1106 

  1107 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">invert</span> [x] (<span style="color: #8cd0d3;">/</span> x))

  1108 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">negate</span> [x] (<span style="color: #8cd0d3;">-</span> x))

  1109 

  1110 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">exp</span> [x] (Math/exp x))

  1111 

  1112 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">sin</span> [x] (Math/sin x))

  1113 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">cos</span> [x] (Math/cos x))

  1114 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">tan</span> [x] (Math/tan x))

  1115 

  1116 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">sec</span> (<span style="color: #8cd0d3;">comp</span> invert cos))

  1117 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">csc</span> (<span style="color: #8cd0d3;">comp</span> invert sin))

  1118 

  1119 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">sinh</span> [x] (Math/sinh x))

  1120 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">cosh</span> [x] (Math/cosh x))

  1121 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">tanh</span> [x] (Math/tanh x))

  1122 

  1123 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">sech</span> (<span style="color: #8cd0d3;">comp</span> invert cosh))

  1124 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">csch</span> (<span style="color: #8cd0d3;">comp</span> invert sinh))

  1125 

  1126 

  1127 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">------------</span>

  1128 

  1129 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">factorial</span>

  1130   <span style="color: #8fb28f;">"Computes the factorial of the nonnegative integer n."</span>

  1131   [n]

  1132   (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #8cd0d3;">neg?</span> n)

  1133     (<span style="color: #f0dfaf; font-weight: bold;">throw</span> (Exception. <span style="color: #cc9393;">"Cannot compute the factorial of a negative number."</span>))

  1134     (<span style="color: #8cd0d3;">reduce</span> * 1 (<span style="color: #8cd0d3;">range</span> 1 (<span style="color: #8cd0d3;">inc</span> n)))))

  1135 

  1136 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">exact-quotient</span> [n d] (<span style="color: #8cd0d3;">/</span> n d))

  1137 

  1138 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">binomial-coefficient</span>

  1139   <span style="color: #8fb28f;">"Computes the number of different ways to choose m elements from n."</span>

  1140   [n m]

  1141   (<span style="color: #8cd0d3;">assert</span> (<span style="color: #8cd0d3;">&lt;=</span> 0 m n))

  1142   (<span style="color: #f0dfaf; font-weight: bold;">let</span> [difference (<span style="color: #8cd0d3;">-</span> n m)]

  1143     (exact-quotient

  1144      (<span style="color: #8cd0d3;">reduce</span> * (<span style="color: #8cd0d3;">range</span> n (<span style="color: #8cd0d3;">max</span> difference m) -1 ))

  1145      (factorial (<span style="color: #8cd0d3;">min</span> difference m)))))

  1146 

  1147 (<span style="color: #f0dfaf; font-weight: bold;">defn-memo</span> <span style="color: #f0dfaf;">stirling-1</span>

  1148   <span style="color: #cc9393;">"Stirling numbers of the first kind: the number of permutations of n</span>

  1149 <span style="color: #cc9393;">  elements with exactly m permutation cycles. "</span>

  1150   [n k]

  1151   <span style="color: #708070;">;</span><span style="color: #7f9f7f;">(assert (&lt;= 1 k n))</span>

  1152   (<span style="color: #f0dfaf; font-weight: bold;">if</span> (zero? n)

  1153     (<span style="color: #f0dfaf; font-weight: bold;">if</span> (zero? k) 1 0)

  1154     (<span style="color: #8cd0d3;">+</span> (stirling-1 (<span style="color: #8cd0d3;">dec</span> n) (<span style="color: #8cd0d3;">dec</span> k))

  1155        (<span style="color: #8cd0d3;">*</span> (<span style="color: #8cd0d3;">dec</span> n) (stirling-1 (<span style="color: #8cd0d3;">dec</span> n) k)))))

  1156 

  1157 (<span style="color: #f0dfaf; font-weight: bold;">defn-memo</span> <span style="color: #f0dfaf;">stirling-2</span> <span style="color: #708070;">;;</span><span style="color: #7f9f7f;">count-partitions</span>

  1158   <span style="color: #cc9393;">"Stirling numbers of the second kind: the number of ways to partition a set of n elements into k subsets."</span>

  1159   [n k]

  1160   (<span style="color: #f0dfaf; font-weight: bold;">cond</span>

  1161    (<span style="color: #8cd0d3;">=</span> k 1) 1

  1162    (<span style="color: #8cd0d3;">=</span> k n) 1

  1163    <span style="color: #8cd0d3;">:else</span> (<span style="color: #8cd0d3;">+</span> (stirling-2 (<span style="color: #8cd0d3;">dec</span> n) (<span style="color: #8cd0d3;">dec</span> k))

  1164             (<span style="color: #8cd0d3;">*</span> k (stirling-2 (<span style="color: #8cd0d3;">dec</span> n) k)))))

  1165 

  1166 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">harmonic-number</span> [n]

  1167   (<span style="color: #8cd0d3;">/</span> (stirling-1 (<span style="color: #8cd0d3;">inc</span> n) 2)

  1168      (factorial n)))

  1169 

  1170 

  1171 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">sum</span>

  1172   [f low high]

  1173   (<span style="color: #8cd0d3;">reduce</span> + (<span style="color: #8cd0d3;">map</span> f (<span style="color: #8cd0d3;">range</span> low (<span style="color: #8cd0d3;">inc</span> high)))))

  1174 

  1175 </pre>

  1176 

  1177 

  1178 

  1179 

  1180 

  1181 

  1182 

  1183 

  1184 

  1185 

  1186 

  1187 

  1188 

  1189 

  1190 

  1191 </div>

  1192 

  1193 </div>

  1194 

  1195 <div id="outline-container-5" class="outline-2">

  1196 <h2 id="sec-5"><span class="section-number-2">5</span> Differentiation </h2>

  1197 <div class="outline-text-2" id="text-5">

  1198 

  1199 

  1200 <p>

  1201 We compute derivatives by passing special <b>differential objects</b> \([x,

  1202 dx]\) through functions. Roughly speaking, applying a function \(f\) to a

  1203 differential object \([x, dx]\) should produce a new differential

  1204 object \([f(x),\,Df(x)\cdot dx]\).

  1205 </p>

  1206 

  1207 

  1208 \([x,\,dx]\xrightarrow{\quad f \quad}[f(x),\,Df(x)\cdot dx]\)

  1209 <p>

  1210 Notice that you can obtain the derivative of \(f\) from this

  1211 differential object, as it is the coefficient of the \(dx\) term. Also,

  1212 as you apply successive functions using this rule, you get the

  1213 chain-rule answer you expect:

  1214 </p>

  1215 

  1216 

  1217 \([f(x),\,Df(x)\cdot dx]\xrightarrow{\quad g\quad} [gf(x),\,

  1218 Dgf(x)\cdot Df(x) \cdot dx ]\)

  1219 

  1220 <p>

  1221 In order to generalize to multiple variables and multiple derivatives,

  1222 we use a <b>power series of differentials</b>, a sortred infinite sequence which

  1223 contains all terms like \(dx\cdot dy\), \(dx^2\cdot dy\), etc.

  1224 </p>

  1225 

  1226 

  1227 

  1228 

  1229 <pre class="src src-clojure">(<span style="color: #f0dfaf; font-weight: bold;">in-ns</span> 'sicm.utils)

  1230 

  1231 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">differential-seq</span>

  1232   <span style="color: #8fb28f;">"Constructs a sequence of differential terms from a numerical</span>

  1233 <span style="color: #8fb28f;">coefficient and a list of keys for variables. If no coefficient is supplied, uses 1."</span>

  1234   ([variables] (differential-seq* 1 variables))

  1235   ([coefficient variables &amp; cvs]

  1236      (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #8cd0d3;">number?</span> coefficient)

  1237        (<span style="color: #8cd0d3;">conj</span> (<span style="color: #8cd0d3;">assoc</span> {} (<span style="color: #8cd0d3;">apply</span> sorted-set variables) coefficient)

  1238              (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #8cd0d3;">empty?</span> cvs)

  1239                nil

  1240                (<span style="color: #8cd0d3;">apply</span> differential-seq* cvs)))

  1241        (<span style="color: #8cd0d3;">apply</span> differential-seq* 1 coefficient 1 variables cvs)  

  1242        )))

  1243 

  1244 

  1245 

  1246 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">differential-add</span>

  1247   <span style="color: #8fb28f;">"Add two differential sequences by combining like terms."</span>

  1248   [dseq1 dseq2]

  1249   (<span style="color: #8cd0d3;">merge-with</span> + dseq1 dseq2))

  1250 

  1251 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">differential-multiply</span>

  1252   <span style="color: #8fb28f;">"Multiply two differential sequences. The square of any differential variable is zero since differential variables are infinitesimally small."</span>

  1253   [dseq1 dseq2]

  1254   (<span style="color: #8cd0d3;">reduce</span>

  1255    (<span style="color: #8cd0d3;">fn</span> [m [[vars1 coeff1] [vars2 coeff2]]]

  1256      (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #8cd0d3;">empty?</span> (clojure.set/<span style="color: #dfdfbf; font-weight: bold;">intersection</span> vars1 vars2))

  1257        (<span style="color: #8cd0d3;">assoc</span> m (clojure.set/<span style="color: #dfdfbf; font-weight: bold;">union</span> vars1 vars2) (<span style="color: #8cd0d3;">*</span> coeff1 coeff2))

  1258        m))

  1259    {}

  1260   (clojure.contrib.combinatorics/cartesian-product

  1261    dseq1

  1262    dseq2)))

  1263 

  1264 

  1265 

  1266 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">big-part</span>

  1267   <span style="color: #8fb28f;">"Returns the 'finite' part of the differential sequence."</span>

  1268   [dseq]

  1269   (<span style="color: #f0dfaf; font-weight: bold;">let</span>

  1270       [keys (<span style="color: #8cd0d3;">sort-by</span> count (<span style="color: #8cd0d3;">keys</span> dseq))

  1271        pivot-key (<span style="color: #8cd0d3;">first</span> (<span style="color: #8cd0d3;">last</span> keys))]

  1272 

  1273     (<span style="color: #8cd0d3;">apply</span> hash-map

  1274     (<span style="color: #8cd0d3;">reduce</span> concat

  1275             (<span style="color: #8cd0d3;">filter</span> (<span style="color: #8cd0d3;">comp</span> pivot-key first) dseq)

  1276             ))))

  1277 

  1278 

  1279 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">small-part</span>

  1280   <span style="color: #8fb28f;">"Returns the 'infinitesimal' part of the differential sequence."</span>

  1281   [dseq]

  1282     (<span style="color: #f0dfaf; font-weight: bold;">let</span>

  1283       [keys (<span style="color: #8cd0d3;">sort-by</span> count (<span style="color: #8cd0d3;">keys</span> dseq))

  1284        pivot-key (<span style="color: #8cd0d3;">first</span> (<span style="color: #8cd0d3;">last</span> keys))]

  1285 

  1286     (<span style="color: #8cd0d3;">apply</span> hash-map

  1287     (<span style="color: #8cd0d3;">reduce</span> concat

  1288             (<span style="color: #8cd0d3;">remove</span> (<span style="color: #8cd0d3;">comp</span> pivot-key first) dseq)

  1289             ))))

  1290 

  1291 

  1292 

  1293 

  1294 

  1295 

  1296 

  1297 

  1298 

  1299 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">;; A differential term consists of a numerical coefficient and a</span>

  1300 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">;; sorted  </span>

  1301 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">(defrecord DifferentialTerm [coefficient variables])</span>

  1302 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">(defmethod print-method DifferentialTerm</span>

  1303 <span style="color: #708070;">;;    </span><span style="color: #7f9f7f;">[o w]</span>

  1304 <span style="color: #708070;">;;     </span><span style="color: #7f9f7f;">(print-simple</span>

  1305 <span style="color: #708070;">;;      </span><span style="color: #7f9f7f;">(apply str (.coefficient o)(map (comp (partial str "d") name) (.variables o)))</span>

  1306 <span style="color: #708070;">;;      </span><span style="color: #7f9f7f;">w))</span>

  1307 

  1308 

  1309 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">(defn differential-seq</span>

  1310 <span style="color: #708070;">;;   </span><span style="color: #7f9f7f;">"Constructs a sequence of differential terms from a numerical</span>

  1311 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">coefficient and a list of keywords for variables. If no coefficient is</span>

  1312 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">supplied, uses 1."</span>

  1313 <span style="color: #708070;">;;   </span><span style="color: #7f9f7f;">([variables] (differential-seq 1 variables))</span>

  1314 <span style="color: #708070;">;;   </span><span style="color: #7f9f7f;">([coefficient variables]</span>

  1315 <span style="color: #708070;">;;      </span><span style="color: #7f9f7f;">(list</span>

  1316 <span style="color: #708070;">;;       </span><span style="color: #7f9f7f;">(DifferentialTerm. coefficient (apply sorted-set variables))))</span>

  1317 <span style="color: #708070;">;;   </span><span style="color: #7f9f7f;">([coefficient variables &amp; cvs]</span>

  1318 <span style="color: #708070;">;;      </span><span style="color: #7f9f7f;">(sort-by</span>

  1319 <span style="color: #708070;">;;       </span><span style="color: #7f9f7f;">#(vec(.variables %))</span>

  1320 <span style="color: #708070;">;;       </span><span style="color: #7f9f7f;">(concat (differential-seq coefficient variables) (apply differential-seq cvs)))))</span>

  1321 

  1322 </pre>

  1323 

  1324 

  1325 

  1326 

  1327 

  1328 

  1329 

  1330 

  1331 

  1332 

  1333 

  1334 

  1335 

  1336 

  1337 

  1338 

  1339 </div>

  1340 

  1341 </div>

  1342 

  1343 <div id="outline-container-6" class="outline-2">

  1344 <h2 id="sec-6"><span class="section-number-2">6</span> Symbolic Manipulation </h2>

  1345 <div class="outline-text-2" id="text-6">

  1346 

  1347 

  1348 

  1349 

  1350 

  1351 <pre class="src src-clojure">(<span style="color: #f0dfaf; font-weight: bold;">in-ns</span> 'sicm.utils)

  1352 

  1353 (<span style="color: #f0dfaf; font-weight: bold;">deftype</span> <span style="color: #f0dfaf;">Symbolic</span> [type expression]

  1354   Object

  1355   (equals [this that]

  1356           (<span style="color: #f0dfaf; font-weight: bold;">cond</span>

  1357            (<span style="color: #8cd0d3;">=</span> (.expression this) (.expression that)) true

  1358            <span style="color: #8cd0d3;">:else</span>

  1359            (Symbolic.

  1360             java.lang.Boolean

  1361             (<span style="color: #8cd0d3;">list</span> '= (.expression this) (.expression that)))

  1362           )))

  1363 

  1364 

  1365 

  1366 

  1367 (<span style="color: #f0dfaf; font-weight: bold;">deftype</span> <span style="color: #f0dfaf;">AbstractSet</span> [glyph membership-test])

  1368 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">member?</span> [abstract-set x]

  1369   ((.membership-test abstract-set) x))

  1370 

  1371 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">------------ Some important AbstractSets</span>

  1372 

  1373 

  1374 (<span style="color: #f0dfaf; font-weight: bold;">def</span> <span style="color: #f0dfaf;">Real</span>

  1375      (AbstractSet.

  1376       'R

  1377       (<span style="color: #8cd0d3;">fn</span>[x](<span style="color: #8cd0d3;">number?</span> x))))

  1378 

  1379 

  1380 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">------------ Create new AbstractSets from existing ones</span>

  1381 

  1382 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">abstract-product</span>

  1383   <span style="color: #8fb28f;">"Gives the cartesian product of abstract sets."</span>

  1384   ([sets]

  1385        (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #8cd0d3;">=</span> (<span style="color: #8cd0d3;">count</span> sets) 1) (<span style="color: #8cd0d3;">first</span> sets)

  1386            (AbstractSet.

  1387             (<span style="color: #8cd0d3;">symbol</span>

  1388              (<span style="color: #8cd0d3;">apply</span> str

  1389                     (<span style="color: #8cd0d3;">interpose</span> 'x (<span style="color: #8cd0d3;">map</span> #(.glyph %) sets))))

  1390            (<span style="color: #8cd0d3;">fn</span> [x]

  1391              (<span style="color: #f0dfaf; font-weight: bold;">and</span>

  1392               (<span style="color: #8cd0d3;">coll?</span> x)

  1393               (<span style="color: #8cd0d3;">=</span> (<span style="color: #8cd0d3;">count</span> sets) (<span style="color: #8cd0d3;">count</span> x))

  1394               (<span style="color: #8cd0d3;">reduce</span> #(<span style="color: #f0dfaf; font-weight: bold;">and</span> %1 %2)

  1395                       true

  1396                       (<span style="color: #8cd0d3;">map</span> #(member? %1 %2) sets x)))))))

  1397   ([abstract-set n]

  1398      (abstract-product (<span style="color: #8cd0d3;">repeat</span> n abstract-set))))

  1399 

  1400 

  1401               

  1402 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">abstract-up</span>

  1403   <span style="color: #8fb28f;">"Returns the abstract set of all up-tuples whose items belong to the</span>

  1404 <span style="color: #8fb28f;">  corresponding abstract sets in coll."</span>

  1405   ([coll]

  1406      (AbstractSet.

  1407       (<span style="color: #8cd0d3;">symbol</span> (<span style="color: #8cd0d3;">str</span> <span style="color: #cc9393;">"u["</span>

  1408                    (<span style="color: #8cd0d3;">apply</span> str

  1409                           (<span style="color: #8cd0d3;">interpose</span> <span style="color: #cc9393;">" "</span>

  1410                                      (<span style="color: #8cd0d3;">map</span> #(.glyph %) coll)))

  1411                    <span style="color: #cc9393;">"]"</span>))

  1412       (<span style="color: #8cd0d3;">fn</span> [x]

  1413         (<span style="color: #f0dfaf; font-weight: bold;">and</span>

  1414          (<span style="color: #8cd0d3;">satisfies?</span> Spinning x)

  1415          (up? x)

  1416          (<span style="color: #8cd0d3;">=</span> (<span style="color: #8cd0d3;">count</span> coll) (<span style="color: #8cd0d3;">count</span> x))

  1417          (<span style="color: #8cd0d3;">reduce</span>

  1418           #(<span style="color: #f0dfaf; font-weight: bold;">and</span> %1 %2)

  1419           true

  1420           (<span style="color: #8cd0d3;">map</span> #(member? %1 %2) coll x))))))

  1421   ([abstract-set n]

  1422      (abstract-up (<span style="color: #8cd0d3;">repeat</span> n abstract-set))))

  1423 

  1424 

  1425 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">abstract-down</span>

  1426   <span style="color: #8fb28f;">"Returns the abstract set of all down-tuples whose items belong to the</span>

  1427 <span style="color: #8fb28f;">  corresponding abstract sets in coll."</span>

  1428   ([coll]

  1429      (AbstractSet.

  1430       (<span style="color: #8cd0d3;">symbol</span> (<span style="color: #8cd0d3;">str</span> <span style="color: #cc9393;">"d["</span>

  1431                    (<span style="color: #8cd0d3;">apply</span> str

  1432                           (<span style="color: #8cd0d3;">interpose</span> <span style="color: #cc9393;">" "</span>

  1433                                      (<span style="color: #8cd0d3;">map</span> #(.glyph %) coll)))

  1434                    <span style="color: #cc9393;">"]"</span>))

  1435       (<span style="color: #8cd0d3;">fn</span> [x]

  1436         (<span style="color: #f0dfaf; font-weight: bold;">and</span>

  1437          (<span style="color: #8cd0d3;">satisfies?</span> Spinning x)

  1438          (down? x)

  1439          (<span style="color: #8cd0d3;">=</span> (<span style="color: #8cd0d3;">count</span> coll) (<span style="color: #8cd0d3;">count</span> x))

  1440          (<span style="color: #8cd0d3;">reduce</span>

  1441           #(<span style="color: #f0dfaf; font-weight: bold;">and</span> %1 %2)

  1442           true

  1443           (<span style="color: #8cd0d3;">map</span> #(member? %1 %2) coll x))))))

  1444   ([abstract-set n]

  1445      (abstract-down (<span style="color: #8cd0d3;">repeat</span> n abstract-set))))

  1446 

  1447               

  1448 

  1449 

  1450 

  1451                                         <span style="color: #708070;">;</span><span style="color: #7f9f7f;">-------ABSTRACT FUNCTIONS</span>

  1452 (<span style="color: #f0dfaf; font-weight: bold;">defrecord</span> <span style="color: #f0dfaf;">AbstractFn</span>

  1453   [<span style="color: #dfdfbf; font-weight: bold;">#^AbstractSet</span> domain <span style="color: #dfdfbf; font-weight: bold;">#^AbstractSet</span> codomain])

  1454 

  1455 

  1456 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">print-method</span> AbstractFn

  1457   [o w]

  1458   (<span style="color: #8cd0d3;">print-simple</span>

  1459    (<span style="color: #8cd0d3;">str</span>

  1460     <span style="color: #cc9393;">"f:"</span>

  1461    (.glyph (<span style="color: #8cd0d3;">:domain</span> o))

  1462    <span style="color: #cc9393;">"--&gt;"</span>

  1463    (.glyph (<span style="color: #8cd0d3;">:codomain</span> o))) w))

  1464 </pre>

  1465 

  1466 

  1467 

  1468 

  1469 

  1470 

  1471 

  1472 </div>

  1473 </div>

  1474 <div id="postamble">

  1475 <p class="date">Date: 2011-08-09 18:41:37 EDT</p>

  1476 <p class="author">Author: Robert McIntyre & Dylan Holmes</p>

  1477 <p class="creator">Org version 7.6 with Emacs version 23</p>

  1478 <a href="http://validator.w3.org/check?uri=referer">Validate XHTML 1.0</a>

  1479 </div>

  1480 </div>

  1481 </body>

  1482 </html>
author	Robert McIntyre <rlm@mit.edu>
date	Fri, 28 Oct 2011 00:03:05 -0700
parents
children