sicmutils · sritchie · May 27, 2021 · May 3, 2021 · May 5, 2021 · May 11, 2021
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -2,6 +2,183 @@
 
 ## [unreleased]
 
+- #341 takes on a large rewrite of the rational function and polynomial
+  simplfiers. One goal of this project was to improve the performance of the
+  Bianchi Identities in `sicmutils.fdg.bianchi-test`, and I'm happy to say that
+  they are now a good bit faster than the original scmutils implementation.
+
+  `sicmutils.polynomial` and `sicmutils.rational-function` are now solid data
+  structures of their own, with many operations installed into the generic
+  system. These are now valuable and useful outside of their role in the
+  simplifier.
+
+  This was a large project, and many small improvements and bugfixes snuck in.
+  Here is the full list:
+
+  - `v/kind` now works for `sorted-map` instances.
+
+  - GCD in Clojurescript is now fast and efficient between all combinations of
+    `js/BigInt` and `js/Number`, and in Clojure between all combinations of
+    `clojure.lang.BigInt`, `BigInteger`, `Long` and `Integer`.
+
+  - on the JVM, GCD now works properly with rational numbers. Previously
+    anything non-integral would return `1`; now `(gcd 1/2 1/3)` properly returns
+    `1/6`.
+
+  - `g/exact-divide` now succeeds for all non-exact `::v/scalar` types (symbols,
+    floats, etc) either if the denominator is zero, or if the two arguments are
+    equal. Else, it throws, just like before.
+
+  - A multi-arity call to `sicmutils.generic/*` now stops if it encounters a
+    0, rather than attempting to multiply all remaining items by 0.
+
+  - The default function for `sicmutils.generic/lcm` protects against overflow
+    by dividing only a single one of its arguments `a` and `b` by `(gcd a b)`.
+
+  - `(g/lcm 0 0)` now properly returns 0.
+
+  - New `sicmutils.util.aggregate/{monoid,group}` functions let you build
+    multi-arity aggregations out of binary combination functions, with an option
+    to bail early at "annihilator" values, like 0 for multiplication.
+
+  - New multi-arity `lcm` and `gcd` implementations for symbolic expressions
+    appropriately handle `0` and `1` on either side, as well as the case where
+    both arguments are equal.
+
+  - In the `sicmutils.numsymb` namespace, thanks to `monoid` and `group`, the
+    `'*`, `'/`, `'-`, `'+`, `'or`, `'and`, `'gcd`, `'lcm` and `'=` operations
+    now have efficient multi-arity implementations that stop computing when they
+    receive an annihilator, like `0` for multiplication or `true` for `or`.
+    Access these via `(sicmutils.numsymb/symbolic-operator <symbol>)`.
+
+  - `sicmutils.series/PowerSeries` gains `arg-scale` and `arg-shift` functions;
+    these are identical to `sicmutils.function/arg-{scale,shift}`, but preserve
+    the `PowerSeries` type. (#367 proposes making these functions generic.)
+
+  - New `sicmutils.ratio/IRational` protocol, with `numerator` and `denominator`
+    functions implemented for ratios and for the `RationalFunction` data type.
+    These two are now exposed in `sicmutils.env`.
+
+  - `sicmutils.simplify.rules/*divide-numbers-through-simplify?*` is now `true`
+    by default; numbers in the denominator will now automatically pull up into
+    the numerator. All tests now reflect this setting.
+
+  - Any analyzer generated from `sicmutils.expression.analyze` can now act on
+    both bare, unwrapped expressions (raw lists etc) and on
+    `sicmutils.expression.Literal` instances. This means that you can now call
+    `sicmutils.simplify/{*rf-simplify*,*poly-simplify*}` as functions and
+    canonicalize some form with either simplifier without triggering a full
+    simplification. A small win, but ice.
+
+  - `sicmutils.polynomial.factor` got a major rewrite, and now exposes a few
+    functions like `poly->factored-expression`, `factor-expression` and
+    `factor`.
+
+      - `factor` is _tremendously useful_! Call `factor` (it's aliased into
+        `sicmutils.env`) on any expression to factor out all possible terms.
+        This makes it much easier to see where there is some cancellation
+        lurking, in, say, some expression you know should equal zero (a
+        residual).
+
+  - bugfix: `sicmutils.expression.Literal` instances now compare their contained
+    expression via `sicmutils.value/=`.
+
+  - `sicmutils.rules/constant-elimination` can now eliminate constants from
+    expressions with any arity, not just binary forms.
+
+
+  Now, the three big namespaces... `sicmutils.polynomial`,
+  `sicmutils.rational-function` and `sicmutils.polynomial.gcd` all got a big
+  overhaul.
+
+  - `sicmutils.polynomial` notes:
+
+    - `Polynomial` uses a new sparse representation for its "power product"
+      term; this, plus an arithmetic rewrite, makes the whole system much faster
+      for larger numbers of variables (for all #s, really).
+
+    - `Polynomial` instances implement many more Clojure(script) protocols. They
+      can hold metadata; they can be evaluated as functions of their
+      indeterminates, and `seq` now returns a sequence of terms.
+
+    - `Polynomial` extends `sicmutils.function/IArity` and
+      `differential/IPerturbed`, so you can use `sicmutils.function/arity`, and
+      take derivatives of functions that return polynomials.
+
+    - In their arithmetic, `Polynomial` instances will drop down to bare
+      coefficients whenever some multiplication or addition removes all
+      indeterminates. All binary arithmetic exposed in the namespace can handle
+      non-`Polynomial` instances on either or both sides, so this is fine.
+      Coefficients are treated as constant polynomials.
+
+    - The namespace holds many new functions. Some choice ones are:
+
+      - constructors: `make`, `constant`, `linear`, `c*xn`, `identity`, and
+        `new-variables`
+
+      - accessor functions: `arity`, `degree`, `coefficients`, `leading-term`,
+        `leading-coefficient`, `leading-exponents`, `leading-base-coefficient`,
+        `trailing-coefficient`, `lowest-degree`
+
+      - predicates: `monomial?`, `monic?`, `univariate?`, `multivariate?`,
+        `negative?`
+
+      - functions to generate new polynomials: `map-coefficients`,
+        `map-exponents`, `scale`, `scale-l`, `normalize`, `reciprocal`,
+        `drop-leading-term`, `contract` and `extend` alongside `contractible?`,
+        `lower-arity`, `raise-arity`, `with-lower-arity`, `arg-scale`,
+        `arg-shift`
+
+      - arithmetic: `negate`, `abs`, `add`, `sub`, `mul`, `square`, `cube`,
+        `expt`, `divide` along with `divisible?`, `evenly-divide`,
+        `pseudo-remainder`, and _lots_ of functions installed into the generic
+        arithmetic system.
+
+      - different ways to evaluate polynomials: `evaluate`, `horner-with-error`
+
+      - calculus! `partial-derivative` and `partial-derivatives` are alive and
+        well, and work with the `D` operator.
+
+      - Functions to get in and out of polynomials from other types:
+        `univariate->dense`, `->power-series`, `expression->`, `->expression`
+
+  - `sicmutils.polynomial.gcd` also got a rewrite; it's fairly clear to read
+    now, and prepared for the eventual addition of the sparse multivariate GCD
+    routine that scmutils uses. There are some efficiency gains here too that
+    let us turn a number of tests back on, or demote them from `:long` markers.
+
+  - `sicmutils.rational-function` notes:
+
+    - `RationalFunction` instances implement many more Clojure(script)
+      protocols. They can hold metadata; they can be evaluated as functions of
+      their indeterminates, and `seq` now returns a pair of `numerator`,
+      `denominator`.
+
+    - `RationalFunction` extends `sicmutils.function/IArity` and
+      `sicmutils.ratio/IRational`, so our generic `arity`, `numerator` and
+      `denominator` work on these instances.
+
+    - Here are some new functions from the `RationalFunction` namespace:
+
+      - constructor: `make`, drops to polynomial or coefficient where needed
+        just like `Polynomial` functions
+
+      - functions to generate new rational functions: `arg-scale`, `arg-shift`
+
+      - predicates: `negative?`
+
+      - arithmetic: `negate`, `abs`, `add`, `sub`, `mul`, `square`, `cube`,
+        `expt`, `invert`, `div`, `gcd`, and many functions installed into the
+        generic arithmetic system.
+
+      - evaluation via `evaluate`
+
+      - calculus! `partial-derivative` and `partial-derivatives` are alive and
+        well, and work with the `D` operator.
+
+      - Functions to get in and out of rational functions from symbolic
+        expressions: `expression->`, `->expression`.
+
 - #360 introduces a number of performance improvements to the
   `sicmutils.differential.Differential` implementation, primarily in `terms:+`
   and `terms:*`. thanks again to @ptaoussanis and the

diff --git a/src/sicmutils/abstract/number.cljc b/src/sicmutils/abstract/number.cljc
@@ -28,7 +28,6 @@
             [sicmutils.expression :as x]
             [sicmutils.generic :as g]
             [sicmutils.numsymb :as sym]
-            [sicmutils.simplify :as s]
             [sicmutils.util :as u]
             [sicmutils.value :as v])
   #?(:clj
@@ -217,8 +216,6 @@
 (defunary g/conjugate 'conjugate)
 
 (defbinary g/gcd 'gcd)
+(defbinary g/lcm 'lcm)
 
 (defmethod g
BEA9
/simplify [Symbol] [a] a)
-(defmethod g/simplify [::x/numeric] [a]
-  (literal-number
-   (s/simplify-expression (v/freeze a))))
diff --git a/src/sicmutils/calculus/manifold.cljc b/src/sicmutils/calculus/manifold.cljc
@@ -33,6 +33,7 @@
             [sicmutils.function :as f]
             [sicmutils.generic :as g]
             [sicmutils.matrix :as matrix]
+            [sicmutils.simplify :refer [simplify-numerical-expression]]
             [sicmutils.structure :as s]
             [sicmutils.util :as u]
             [sicmutils.value :as v]
@@ -297,7 +298,7 @@
   [manifold-point coordinate-system thunk]
   (let [reps (:coordinate-representations manifold-point)]
     (or (@reps coordinate-system)
-        (let [rep (g/simplify (thunk))]
+        (let [rep (s/mapr simplify-numerical-expression (thunk))]
           (swap! reps assoc coordinate-system rep)
           rep))))
 

diff --git a/src/sicmutils/collection.cljc b/src/sicmutils/collection.cljc
@@ -186,7 +186,9 @@
    (identity-like [m] (u/unsupported (str "identity-like: " m)))
    (exact? [m] (every? v/exact? (vals m)))
    (freeze [m] (u/map-vals v/freeze m))
-   (kind [m] (:type m (type m)))
+   (kind [m] (if (sorted? m)
+               (type m)
+               (:type m (type m))))
 
    f/IArity
    (arity [_] [:between 1 2])

diff --git a/src/sicmutils/differential.cljc b/src/sicmutils/differential.cljc
@@ -29,6 +29,7 @@
   (:require [clojure.string :refer [join]]
             [sicmutils.generic :as g]
             [sicmutils.util :as u]
+            [sicmutils.util.aggregate :as ua]
             [sicmutils.util.stream :as us]
             [sicmutils.util.vector-set :as uv]
             [sicmutils.value :as v]))
@@ -254,7 +255,7 @@
 ;; al. (2019).
 ;;
 ;; The solution is to introduce a new $\varepsilon$ for every level, and allow
-;; different $\varepsilon$ instances to multiply without annihalating. Each
+;; different $\varepsilon$ instances to multiply without annihilating. Each
 ;; $\varepsilon$ is called a "tag". [[Differential]] (implemented below) is a
 ;; generalized dual number that can track many tags at once, allowing nested
 ;; derivatives like the one described above to work.
@@ -475,48 +476,14 @@
 ;; NOTE: Clojure vectors already implement this ordering properly, so we can
 ;; use [[clojure.core/compare]] to determine an ordering on a tag list.
 
-(defn- terms:+
-  "Returns the sum of the two supplied sequences of differential terms; any terms
+(def ^{:private true
+       :doc "Returns the sum of the two supplied sequences of differential terms; any terms
   in the result with a zero coefficient will be removed.
 
   Each input must be sequence of `[tag-set, coefficient]` pairs, sorted by
-  `tag-set`.
-
-  NOTE that this function recurs on increasing indices internally instead of
-  walking through the lists directly. This method of traversing vectors is more
-  efficient, and this function is called so often that the performance gain is
-  worth it, and reads almost like the explicit sequence traversal."
-  ([] [])
-  ([xs] xs)
-  ([xs ys]
-   (loop [i (long 0)
-          j (long 0)
-          result (transient [])]
-     (let [x (nth xs i nil)
-           y (nth ys j nil)]
-       (cond (not x) (into (persistent! result) (subvec ys j))
-             (not y) (into (persistent! result) (subvec xs i))
-             :else (let [[x-tags x-coef] x
-                         [y-tags y-coef] y
-                         compare-flag (core-compare x-tags y-tags)]
-                     (cond
-                       ;; If the terms have the same tag set, add the coefficients
-                       ;; together. Include the term in the result only if the new
-                       ;; coefficient is non-zero.
-                       (zero? compare-flag)
-                       (let [sum (g/add x-coef y-coef)]
-                         (recur (inc i)
-                                (inc j)
-                                (if (v/zero? sum)
-                                  result
-                                  (conj! result (make-term x-tags sum)))))
-
-                       ;; Else, pass the smaller term on unchanged and proceed.
-                       (neg? compare-flag)
-                       (recur (inc i) j (conj! result x))
-
-                       :else
-                       (recur i (inc j) (conj! result y)))))))))
+  `tag-set`."}
+  terms:+
+  (ua/merge-fn core-compare g/add v/zero? make-term))
 
 ;; Because we've decided to store terms as a vector, we can multiply two vectors
 ;; of terms by:

diff --git a/src/sicmutils/env.cljc b/src/sicmutils/env.cljc
@@ -41,7 +41,9 @@
                            partial core-partial
                            compare core-compare
                            = core=}
-                  :exclude [+ - * / zero? compare divide #?@(:cljs [= partial])])
+                  :exclude [+ - * / zero? compare divide
+                            numerator denominator
+                            #?@(:cljs [= partial])])
   (:require #?(:clj [potemkin :refer [import-def import-vars]])
             #?(:clj [nrepl.middleware.print])
             [sicmutils.abstract.function :as af #?@(:cljs [:include-macros true])]
@@ -53,6 +55,7 @@
             [sicmutils.generic :as g]
             [sicmutils.modint]
             [sicmutils.operator :as o]
+            [sicmutils.polynomial.factor :as pf]
             [sicmutils.ratio]
             [sicmutils.simplify]
             [sicmutils.structure :as structure]
@@ -244,8 +247,8 @@ constant [Pi](https://en.wikipedia.org/wiki/Pi)."}
  [sicmutils.function arity compose arg-shift arg-scale I]
  [sicmutils.modint chinese-remainder]
  [sicmutils.operator commutator anticommutator]
- #?(:cljs [sicmutils.ratio
-           ratio? rationalize numerator denominator])
+ [sicmutils.polynomial.factor factor]
+ [sicmutils.ratio numerator denominator #?@(:cljs [ratio? rationalize])]
  [sicmutils.series binomial-series partial-sums]
  [sicmutils.util bigint? #?@(:cljs [bigint])]
 

diff --git a/src/sicmutils/expression.cljc b/src/sicmutils/expression.cljc
@@ -87,7 +87,7 @@
              (if (instance? Literal b)
                (let [b ^Literal b]
                  (and (= type (.-type b))
-                      (= expression (.-expression b))
+                      (v/= expression (.-expression b))
                       (= m (.-m b))))
                (v/= expression b))))
 
@@ -114,7 +114,7 @@
                (if (instance? Literal b)
                  (let [b ^Literal b]
                    (and (= type (.-type b))
-                        (= expression (.-expression b))
+                        (v/= expression (.-expression b))
                         (= m (.-m b))))
                  (v/= expression b)))
 
@@ -221,7 +221,11 @@
                   (sequential? node)
                   (let [[f-sym & args] node]
                     (if-let [f (sym->f f-sym)]
-                      (apply f (map walk args))
+                      ;; NOTE: I'm not sure why this `doall` is required.
+                      ;; Without it, we were getting heisenbugs in the rational
+                      ;; function simplifier, and `mismatched-arity` notes.
+                      (apply f (doall
+                                (map walk args)))
                       (u/illegal (str "Missing fn for symbol - " f-sym))))
                   :else node))]
     (walk