sicmutils · sritchie · Apr 7, 2021 · Apr 7, 2021 · Apr 7, 2021 · Apr 7, 2021
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -2,6 +2,28 @@
 
 ## [unreleased]
 
+- From #342:
+
+  - Added `sicmutils.calculus.derivative/D-as-matrix` and
+    `sicmutils.matrix/as-matrix`, ported from scmutils.
+
+  - converted `sicmutils.modint.ModInt` to a `deftype`; this allows `ModInt`
+    instances to be `=` to non-`ModInt` numbers on the right, if the right side
+    is equal to the residue plus any integer multiple of the modulus. `v/=`
+    gives us this behavior with numbers on the LEFT too, and `ModInt` on the
+    right.
+
+    - This change means that `:i` and `:m` won't return the residue and modulus
+      anymore. `sicmutils.modint` gains new `residue` and `modulus` functions to
+      access these attributes.
+
+  - The JVM version of sicmutils gains more efficient `gcd` implementations
+    for `Integer` and `Long` (in addition to the existing native `BigInteger`
+    `gcd`), thanks to our existing Apache Commons-Math dependency.
+
+  - `sicmutils.structure/dual-zero` aliases `compatible-zero` to match the
+    scmutils interface. Both are now aliased into `sicmutils.env`.
+
 - From #339:
 
   - `Structure` instances can now hold metadata.

diff --git a/src/sicmutils/calculus/derivative.cljc b/src/sicmutils/calculus/derivative.cljc
@@ -456,6 +456,12 @@
   (o/make-operator #(g/partial-derivative % [])
                    g/derivative-symbol))
 
+(defn D-as-matrix [F]
+  (fn [s]
+    (matrix/s->m (s/compatible-shape (F s))
+	               ((D F) s)
+	               s)))
+
 (defn partial
   "Returns an operator that, when applied to a function `f`, produces a function
   that computes the partial derivative of `f` at the (zero-based) slot index
@@ -464,51 +470,6 @@
   (o/make-operator #(g/partial-derivative % selectors)
                    `(~'partial ~@selectors)))
 
-(def ^{:doc "Operator that takes a function `f` and returns a new function that
-  calculates the [Gradient](https://en.wikipedia.org/wiki/Gradient) of `f`.
-
-  The related [[D]] operator returns a function that produces a structure of the
-  opposite orientation as [[Grad]]. Both of these functions use forward-mode
-  automatic differentiation."}
-  Grad
-  (-> (fn [f]
-        (f/compose s/opposite
-                   (g/partial-derivative f [])))
-      (o/make-operator 'Grad)))
-
-(def ^{:doc "Operator that takes a function `f` and returns a function that
-  calculates the [Divergence](https://en.wikipedia.org/wiki/Divergence) of
-  `f` at its input point.
-
- The divergence is a one-level contraction of the gradient."}
-  Div
-  (-> (f/compose g/trace Grad)
-      (o/make-operator 'Div)))
-
-(def ^{:doc "Operator that takes a function `f` and returns a function that
-  calculates the [Curl](https://en.wikipedia.org/wiki/Curl_(mathematics)) of `f`
-  at its input point.
-
-  `f` must be a function from $\\mathbb{R}^3 \\to \\mathbb{R}^3$."}
-  Curl
-  (-> (fn [f-triple]
-        (let [[Dx Dy Dz] (map partial [0 1 2])
-              fx (f/get f-triple 0)
-              fy (f/get f-triple 1)
-              fz (f/get f-triple 2)]
-          (s/up (g/- (Dy fz) (Dz fy))
-                (g/- (Dz fx) (Dx fz))
-                (g/- (Dx fy) (Dy fx)))))
-      (o/make-operator 'Curl)))
-
-(def ^{:doc "Operator that takes a function `f` and returns a function that
-  calculates the [Vector
-  Laplacian](https://en.wikipedia.org/wiki/Laplace_operator#Vector_Laplacian) of
-  `f` at its input point."}
-  Lap
-  (-> (f/compose g/trace (g/* Grad Grad))
-      (o/make-operator 'Lap)))
-
 ;; ## Derivative Utilities
 ;;
 ;; Functions that make use of the differential operators defined above in

diff --git a/src/sicmutils/calculus/form_field.cljc b/src/sicmutils/calculus/form_field.cljc
@@ -450,11 +450,11 @@
                 (assert (= (count args) n)
                         "Wrong number of args to alternation")
                 (g/* (/ 1 (factorial n))
-                     (apply g/+
-                            (map (fn [permutation parity]
-                                   (g/* parity (apply form permutation)))
-                                 (permute/permutation-sequence args)
-                                 (cycle [1 -1])))))]
+                     (ua/generic-sum
+                      (map (fn [permutation parity]
+                             (g/* parity (apply form permutation)))
+                           (permute/permutation-sequence args)
+                           (cycle [1 -1])))))]
         (procedure->nform-field
          alternation n `(~'Alt ~(v/freeze form)))))))
 

diff --git a/src/sicmutils/calculus/vector_calculus.cljc b/src/sicmutils/calculus/vector_calculus.cljc
@@ -25,6 +25,7 @@
   metric and basis."
   (:refer-clojure :exclude [+ - * /])
   (:require [sicmutils.calculus.basis :as b]
+            [sicmutils.calculus.derivative :as d]
             [sicmutils.calculus.covariant :as cov]
             [sicmutils.calculus.form-field :as ff]
             [sicmutils.calculus.hodge-star :as hs]
@@ -90,13 +91,13 @@
   `f` must be a function from $\\mathbb{R}^3 \\to \\mathbb{R}^3$."}
   Curl
   (-> (fn [f-triple]
-        (let [[Dx Dy Dz] (map partial [0 1 2])
+        (let [[Dx Dy Dz] (map d/partial [0 1 2])
               fx (f/get f-triple 0)
               fy (f/get f-triple 1)
               fz (f/get f-triple 2)]
-          (s/up (g/- (Dy fz) (Dz fy))
-                (g/- (Dz fx) (Dx fz))
-                (g/- (Dx fy) (Dy fx)))))
+          (s/up (- (Dy fz) (Dz fy))
+                (- (Dz fx) (Dx fz))
+                (- (Dx fy) (Dy fx)))))
       (o/make-operator 'Curl)))
 
 (defn curl
@@ -114,7 +115,7 @@
   Laplacian](https://en.wikipedia.org/wiki/Laplace_operator#Vector_Laplacian) of
   `f` at its input point."}
   Lap
-  (-> (f/compose g/trace (g/* Grad Grad))
+  (-> (f/compose g/trace (* Grad Grad))
       (o/make-operator 'Lap)))
 
 (defn Laplacian [metric orthonormal-basis]

diff --git a/src/sicmutils/env.cljc b/src/sicmutils/env.cljc
@@ -276,6 +276,7 @@ constant [Pi](https://en.wikipedia.org/wiki/Pi)."}
   simplify]
  [sicmutils.structure
   compatible-shape
+  compatible-zero dual-zero
   down
   mapr sumr
   orientation
@@ -333,7 +334,7 @@ constant [Pi](https://en.wikipedia.org/wiki/Pi)."}
   curvature-components]
 
  [sicmutils.calculus.derivative
-  derivative D taylor-series]
+  derivative D D-as-matrix taylor-series]
 
  [sicmutils.calculus.form-field
   form-field? nform-field? oneform-field?

diff --git a/src/sicmutils/euclid.cljc b/src/sicmutils/euclid.cljc
@@ -57,7 +57,5 @@
 
 ;; multimethod implementation for basic numeric types.
 
-(defmethod g/gcd :default [a b] (gcd a b))
-
-(defmethod g/lcm :default [a b]
-  (g/abs (g/divide (g/* a b) (g/gcd a b))))
+(defmethod g/gcd :default [a b]
+  (gcd a b))
diff --git a/src/sicmutils/generic.cljc b/src/sicmutils/generic.cljc
@@ -178,11 +178,13 @@
   ([] 1)
   ([x] x)
   ([x y]
-   (cond (and (v/numerical? x) (v/zero? x)) (v/zero-like y)
-         (and (v/numerical? y) (v/zero? y)) (v/zero-like x)
-         (v/one? x) y
-         (v/one? y) x
-         :else (mul x y)))
+   (let [numx? (v/numerical? x)
+         numy? (v/numerical? y)]
+     (cond (and numx? (v/zero? x)) (v/zero-like y)
+           (and numy? (v/zero? y)) (v/zero-like x)
+           (and numx? (v/one? x)) y
+           (and numy? (v/one? y)) x
+           :else (mul x y))))
   ([x y & more]
    (reduce * (* x y) more)))
 
@@ -457,6 +459,11 @@
   multiple](https://en.wikipedia.org/wiki/Least_common_multiple) of the two
   inputs `a` and `b`.")
 
+(defmethod lcm :default [a b]
+  (abs
+   (divide (* a b)
+           (gcd a b))))
+
 ;; ### Trigonometric functions
 
 (declare sin)

diff --git a/src/sicmutils/matrix.cljc b/src/sicmutils/matrix.cljc
@@ -27,7 +27,6 @@
                            some core-some}
                   #?@(:cljs [:exclude [get-in some]]))
   (:require [sicmutils.differential :as d]
-            [sicmutils.expression :as x]
             [sicmutils.function :as f]
             [sicmutils.generic :as g]
             [sicmutils.series :as series]
@@ -38,7 +37,7 @@
   #?(:clj
      (:import [clojure.lang Associative AFn IFn Sequential])))
 
-(declare fmap generate I identity-like identity? m:=)
+(declare fmap identity-like identity? m:=)
 
 (derive ::square-matrix ::matrix)
 (derive ::column-matrix ::matrix)
@@ -81,7 +80,7 @@
        Associative
        (assoc [_ k entry] (Matrix. r c (assoc v k entry)))
        (containsKey [_ k] (contains? v k))
-       (entryAt [_ k] (.entryAt v k))
+       (entryAt [_ k] (.entryAt ^Associative v k))
        (count [_] (count v))
        (seq [_] (seq v))
        (valAt [_ key] (get v key))
@@ -540,11 +539,14 @@
   [m u]
   (when (not= (num-cols m) (count u))
     (u/illegal "matrix and tuple incompatible for multiplication"))
-  (s/up* (map (fn [i]
-                (reduce g/+ (for [k (range (num-cols m))]
-                              (g/* (core-get-in m [i k])
-                                   (get u k)))))
-              (range (num-rows m)))))
+  (s/up*
+   (map (fn [i]
+          (ua/generic-sum
+           (fn [k]
+             (g/* (core-get-in m [i k])
+                  (get u k)))
+           0 (num-cols m)))
+        (range (num-rows m)))))
 
 (defn- d*M
   "Multiply a matrix `m` by down tuple `d` on the left. The return value has
@@ -554,10 +556,11 @@
     (u/illegal "matrix and tuple incompatible for multiplication"))
   (s/down*
    (map (fn [i]
-          (reduce g/+ (for [k (range (num-rows m))]
-                        (g/* (get d k)
-                             (core-get-in m [i k])
-                             ))))
+          (ua/generic-sum
+           (fn [k]
+             (g/* (get d k)
+                  (core-get-in m [i k])))
+           0 (num-rows m)))
         (range (num-cols m)))))
 
 (def ^{:dynamic true
@@ -580,6 +583,15 @@
                      (g/* ms (s/unflatten
                               (s/basis-unit nups j) rs)))))))
 
+(defn as-matrix
+  "Any one argument function of a structure can be seen as a matrix. This is only
+  useful if the function has a linear multiplier (e.g. derivative)"
+  [F]
+  (fn [s]
+    (let [v (F s)]
+      (s->m
+       (s/compatible-shape (g/* v s)) v s))))
+
 (defn nth-row
   "Returns the `n`-th row of the supplied matrix `m` as a `down` structure."
   [m n]
@@ -678,16 +690,20 @@
   (assoc m i v))
 
 (defn submatrix
-  "Returns the submatrix of `m` generated by taking
+  "Returns the submatrix of the matrix (or matrix-like structure) `s` generated by
+  taking
 
   - rows from `lowrow` -> `hirow`,
   - columns from `lowcol` -> `hicol`"
-  [m lowrow hirow lowcol hicol]
-  (generate (inc (- hirow lowrow))
-            (inc (- hicol lowcol))
-            (fn [i j]
-              (core-ge
B8FF
t-in m [(+ i lowrow)
-                              (+ j lowcol)]))))
+  [x lowrow hirow lowcol hicol]
+  (let [m (if (s/structure? x)
+            (square-structure-> x (fn [m _] m))
+            x)]
+    (generate (inc (- hirow lowrow))
+              (inc (- hicol lowcol))
+              (fn [i j]
+                (core-get-in m [(+ i lowrow)
+                                (+ j lowcol)])))))
 
 (defn without
   "Returns the matrix formed by deleting the `i`-th row and `j`-th column of the
@@ -733,11 +749,12 @@
     2 (let [[[a b] [c d]] m]
         (g/- (g/* a d)
              (g/* b c)))
-    (reduce g/+ (map g/*
-                     (cycle [1 -1])
-                     (nth m 0)
-                     (for [i (range (num-rows m))]
-                       (determinant (without m 0 i)))))))
+    (ua/generic-sum
+     (map g/*
+          (cycle [1 -1])
+          (nth m 0)
+          (for [i (range (num-rows m))]
+            (determinant (without m 0 i)))))))
 
 (defn cofactors
   "Returns the matrix of cofactors of the supplied square matrix `m`."
@@ -763,7 +780,7 @@
       0 m
       1 (->Matrix 1 1 [[(g/invert (core-get-in m [0 0]))]])
       (let [C (cofactors m)
-            Δ (reduce g/+ (map g/* (nth m 0) (nth C 0)))]
+            Δ (ua/generic-sum (map g/* (nth m 0) (nth C 0)))]
         (fmap #(g/divide % Δ)
               (transpose C))))))