sicmutils · sritchie · Mar 27, 2022 · Mar 23, 2022 · Mar 23, 2022 · Mar 23, 2022
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -2,8 +2,61 @@
 
 ## unreleased
 
+- #503:
+
+  - Changes the default implementation of `square` and `cube` for differentials
+    to use `(expt <> 2)` etc instead of `(* <> <>)`.
+
+    This is a _big deal!_ For certain expressions there's a huge blowup when you
+    square a big symbolic term, and taking the derivative of it TWICE is very
+    messy.
+
+    With this change, differentials use the chain rule to calculuate the
+    derivative of $x^2$ as $2x*x'$, instead of using the product rule and
+    achieving a SECOND differentiatation of the same form and another
+    multiplication: $xx' + x'x$.
+
+    Before a judicious `simplify` call I added, this change dropped the runtime
+    of the `sicmutils.sicm.ch3-test` suite down by 6x. After the simplify change
+    in `sicmutils.examples.top` the tests were still 40% faster in that
+    namespace.
+
+  - Fixes a bug where the `RationalFunction` cube implementation actually called
+    `square`.
+
+  - adds `sicmutils.mechanics.lagrange/Lagrangian` for building function
+    signatures of Lagrangians.
+
+  - adds the `sicmutils.mechanics.time-evolution` namespace
+
+  - adds `sicmutils.mechanics.lagrange/L-axisymmetric-top`, more efficient than
+    the version in `sicmutils.examples.top`
+
+  - Fleshes out `sicmutils.mechanics.hamilton`:
+
+    - New functions: `H-state?`, `compatible-H-state?`, `state->p`, `momenta`,
+      `P`, `literal-Hamiltonian-state`, `L-state->H-state`, `H-state->L-state`,
+      `H-state->matrix`, `matrix->H-state`, `make-Hamiltonian`, `D-phase-space`,
+      `Hamiltonian->Lagrangian-procedure`, `Hamiltonian->Lagrangian`,
+      `flow-derivative`, `flow-transform`, `standard-map-inverse`, `F->K`,
+      `J-func`, `T-func`, `canonical-H?`, `canonical-K?`,
+      `linear-function->multiplier`, `Phi`, `Phi*`, `qp-canonical?`,
+      `polar-canonical-inverse`, `two-particle-center-of-mass` ,
+      `two-particle-center-of-mass-canonical`, `transpose-function`,
+      `multiplicative-transpose`, `symplectic-two-form`, `canonical-transform?`,
+      `J-matrix`, `symplectic?`
+
+    - `F->CH` moves to `F->CT` (`F->CT` is now an alias)
+
+    - `Legendre-transform-fn` becomes `Legendre-transform-procedure` and gains
+      more correctness tests, toggled on and off by the
+      `*validate-Legendre-transform?*` dynamic variable.
+
 - #508 adds `sicmutils.mechanics.noether` namespace, with `Noether-integral`.
 
+- #506 tidies up the build by removing unneeded reader conditionals and
+  replacing renames like `core-=` with a proper require of `clojure.core`.
+
 - #502 begins the port of the remaining items in the scmutils `mechanics`
   package over the Clojure. This PR focuses on `sicmutils.mechanics.lagrange`,
   which contains functions from many files in the original `mechanics` folder.

diff --git a/src/sicmutils/calculus/covariant.cljc b/src/sicmutils/calculus/covariant.cljc
@@ -19,35 +19,34 @@
 ;; This comes from `Lie.scm`.
 
 (defn- vector-field-Lie-derivative [X]
-  (o/make-operator
-   (fn [Y]
-     (cond (f/function? Y)      (X Y)
-           (vf/vector-field? Y) (o/commutator X Y)
-
-           (ff/form-field? Y)
-           (let [k (ff/get-rank Y)
-                 op (fn [& vectors]
-                      (let [vectors (into [] vectors)]
-                        (assert (= k (count vectors))
-                                `(~'≠ ~k ~(count vectors)
-                                  ~@vectors
-                                  ~@(map meta vectors)))
-                        (g/- ((g/Lie-derivative X) (apply Y vectors))
-                             (ua/generic-sum
-                              (fn [i]
-                                (let [xs (update vectors i (g/Lie-derivative X))]
-                                  (apply Y xs)))
-                              0 k))))
-                 name `((~'Lie-derivative ~(v/freeze X))
-                        ~(v/freeze Y))]
-             (ff/procedure->nform-field op k name))
-
-           (s/structure? Y)
-           (s/mapr (vector-field-Lie-derivative X) Y)
-
-           :else (u/unsupported "Bad argument: Lie Derivative")))
-   `(~'Lie-derivative
-     ~(v/freeze X))))
+  (let [freeze-X (v/freeze X)
+        op-name `(~'Lie-derivative ~freeze-X)]
+    (-> (fn rec [Y]
+          (cond (f/function? Y)      (X Y)
+                (vf/vector-field? Y) (o/commutator X Y)
+
+                (ff/form-field? Y)
+                (let [k (ff/get-rank Y)
+                      op (fn [& vectors]
+                           (let [vectors (into [] vectors)]
+                             (assert (= k (count vectors))
+                                     `(~'≠ ~k ~(count vectors)
+                                       ~@vectors
+                                       ~@(map meta vectors)))
+                             (g/- ((g/Lie-derivative X) (apply Y vectors))
+                                  (ua/generic-sum
+                                   (fn [i]
+                                     (let [xs (update vectors i (g/Lie-derivative X))]
+                                       (apply Y xs)))
+                                   0 k))))
+                      name `((~'Lie-derivative ~freeze-X) ~(v/freeze Y))]
+                  (ff/procedure->nform-field op k name))
+
+                (s/structure? Y)
+                (s/mapr (vector-field-Lie-derivative X) Y)
+
+                :else (u/unsupported "Bad argument: Lie Derivative")))
+        (o/make-operator op-name))))
 
 (defmethod g/Lie-derivative [::vf/vector-field] [V]
   (vector-field-Lie-derivative V))

diff --git a/src/sicmutils/differential.cljc b/src/sicmutils/differential.cljc
@@ -1214,11 +1214,18 @@
 (defunary g/negate (lift-1 g/negate))
 (defbinary g/sub (lift-2 g/sub))
 
-(let [mul (lift-2 g/mul)]
+(let [mul  (lift-2 g/mul)
+      expt (lift-2 g/expt)]
   (defbinary g/mul mul)
-  (defunary g/square (fn [x] (mul x x)))
-  (defunary g/cube (fn [x] (mul x (mul x x))))
-  (defbinary g/dot-product mul))
+  (defbinary g/dot-product mul)
+  (defbinary g/expt expt)
+
+  ;; NOTE that it's important that these stay exponents vs repeated
+  ;; multiplication. It is cheaper to compute the derivative just once using the
+  ;; chain rule, vs performing a full multiplication and then taking the
+  ;; derivative of the result.
+  (defunary g/square (fn [x] (expt x 2)))
+  (defunary g/cube (fn [x] (expt x 3))))
 
 (defunary g/invert (lift-1 g/invert))
 (defbinary g/div (lift-2 g/div))
@@ -1264,7 +1271,6 @@
   (defbinary g/solve-linear-right div))
 
 (defunary g/sqrt (lift-1 g/sqrt))
-(defbinary g/expt (lift-2 g/expt))
 (defunary g/log (lift-1 g/log))
 (defunary g/exp (lift-1 g/exp))
 

diff --git a/src/sicmutils/examples/top.cljc b/src/sicmutils/examples/top.cljc
@@ -7,7 +7,7 @@
 
 (defn L
   [A B C gMR]
-  (let [T (rigid/T-rigid-body A B C)
+  (let [T (comp e/simplify (rigid/T-rigid-body A B C))
         V (fn [[_ [theta _ _]]]
             (* gMR (cos theta)))]
     (- T V)))

diff --git a/src/sicmutils/generic.cljc b/src/sicmutils/generic.cljc
@@ -850,7 +850,7 @@ defaults to `ln((1 + sqrt(1+x^2)) / x)`."
            (if (v/zero? x)
              x
              (sub (div (cos x) x)
-                  (div (sin x) (* x x)))))})
+                  (div (sin x) (square x)))))})
 
 (defmethod sinc :default [x]
   (if (v/zero? x)
@@ -886,7 +886,7 @@ defaults to `ln((1 + sqrt(1+x^2)) / x)`."
              x
              (let [sx (sec x)]
                (sub (div (* sx sx) x)
-                    (div (tan x) (* x x))))))})
+                    (div (tan x) (square x))))))})
 
 (defmethod tanc :default [x]
   (if (v/zero? x)
@@ -907,7 +907,7 @@ defaults to `ln((1 + sqrt(1+x^2)) / x)`."
            (if (v/zero? x)
              x
              (sub (div (cosh x) x)
-                  (div (sinh x) (* x x)))))})
+                  (div (sinh x) (square x)))))})
 
 (defmethod sinhc :default [x]
   (if (v/zero? x)
@@ -927,7 +927,7 @@ defaults to `ln((1 + sqrt(1+x^2)) / x)`."
              x
              (let [sx (sech x)]
                (sub (div (* sx sx) x)
-                    (div (tanh x) (* x x))))))})
+                    (div (tanh x) (square x))))))})
 
 (defmethod tanhc :default [x]
   (if (v/zero? x)