mentat-collective / emmy Goto Github PK

The Emmy Computer Algebra System.

License: GNU General Public License v3.0

Clojure 99.96% Emacs Lisp 0.02% Shell 0.02%

automatic-differentiation calculus clojure clojurescript computer-algebra differential-geometry explorable-explanations hamiltonian lagrangian-mechanics mathematics physics physics-simulation sussman symbolic-math

emmy's People

Contributors

Stargazers

Watchers

Forkers

plurigrid mk allofphysicsgraph danieltanfh95 sritchie vnredmon davidalphafox mhuebert jamestiotio brucerennie mathisto aakarsh v0idmatr1x eightysteele danielrmeyer alexgian ingydotnet

emmy's Issues

Handle structural args in sicmutils.expression.compile/compile-fn

Right now, compile-fn (link) can handle functions of different arities, and the GUTS of compile-state-fn can handle arguments with any structure.

It would be great to unify these, so compile-state-fn called the guts of compile-fn, and the latter could take an arbitrary argument template instead of just n for its arity.

Some nice design work waiting here!

implement ->graphviz renderer, tune up ->TeX?

I've been thinking about what other renderers might make sense, in addition to "infix" and "TeX"... a couple of thoughts (cc @mk).

If there are significant differences in what katex and mathjax are able to render, we could take a flag that customized the latex renderer to each one. @mk or @littleredcomputer , have either of you hit cases here?
An expression tree => graphviz renderer would be very cool.

Here is some prior art for the latter. Here's an example of sympy doing this trick.

add Bessel function implementation to `sicmutils.numerical`

The goal here is to get the bessel functions J0, Y0, J1 and Y1 into Clojure. Once we do this, we can use the recurrence relations in bessel.scm to implement the binary functions.

I found a nice paper by John Harrison called "Fast and Accurate Bessel Function Computation", with an associated presentation here.

He used a power series approach for part of it, so I spent some time this morning getting that working using our power series library. Too good to pass up, right?

This task is not done, at all, but this is a fast, promising approach that someone should feel free to take on as a challenge. Happy to help stare at the code!

Here's how far I got (also stashed at https://github.com/sicmutils/sicmutils/tree/sritchie/bessel_exploration):

(ns sicmutils.numerical.bessel
  (:require [sicmutils.generic :as g]
            [sicmutils.util :as u]
            [sicmutils.util.stream :as us]
            [sicmutils.series :as s]
            [sicmutils.value :as v]))

(defn coef
  "This inefficiently implements `(n, m)` from section 4 (Asymptotic Expansions)
  of the paper, just to check that we have the right idea."
  [n m]
  (let [num (reduce (fn [acc i]
                      (* acc (- (* 4 (* n n))
                                (g/square
                                 (dec (* 2 (inc i)))))))
                    1
                    (range 0 m))
        den (* (g/expt 2 (* 2 m))
               (g/factorial m))]
    (/ num den)))

(comment
  ;; next, these are the example expansions for $P_0(x)$ and $Q_0(x)$, just to
  ;; verify:
  (let [p0-series (-> (s/power-series*
                       (map (fn [m]
                              (* (g/expt (- 1) m)
                                 (/ (coef 0 (* 2 m))
                                    (g/expt 2 (* 2 m)))))
                            (range)))
                      (s/inflate 2))
        q0-series (-> (s/power-series*
                       (map (fn [m]
                              (* (g/expt (- 1) m)
                                 (/ (coef 0 (inc (* 2 m)))
                                    (g/expt 2 (inc (* 2 m))))))
                            (range)))
                      (s/inflate 2))]

    (= '(1
         0
         (/ -9 128)
         0
         (/ 3675 32768)
         0
         (/ -2401245 4194304)
         0
         (/ 13043905875 2147483648))
       (v/freeze (take 9 p0-series)))

    (= '((/ -1 8)
         0
         (/ 75 1024)
         0
         (/ -59535 262144)
         0
         (/ 57972915 33554432)
         0)
       (v/freeze (take 8 q0-series)))))

(defn updater
  "Returns a function that generates $(n, m+2)$ given a triple of the numerator,
  denominator and `m` from the previous evaluation.

  (Returns a triple of the same form.)"
  [n]
  (let [four*n**2 (* 4 n n)]
    (fn [[num denom m]]
      (let [two-m     (* 2 m)
            num'      (* num
                         (- four*n**2 (g/square (+ two-m 1)))
                         (- four*n**2 (g/square (+ two-m 3))))
            den'      (* denom
                         16
                         (+ m 1)
                         (+ m 2))]
        [num' den' (+ 2 m)]))))

(defn nm-seq
  "Given a value of `n` and an optional initial `m` (defaults to 0, only 0 or 1
  accepted), returns an infinite sequence of:

  [(n, m), (n, m+2), (n, m+4), ...]

  Using the definition of (n, m) from the paper."
  ([n]
   (raw-coef** n 0))
  ([n m]
   {:pre [(#{0 1} m)]}
   (let [num (if (zero? m)
               1
               (- (* 4 n n) (g/square (- (* 2 m) 1))))
         den (g/expt 2 (* 2 m))
         m   (u/bigint m)]
     (iterate (updater n) [num den m]))))

(defn p-series
  "Returns the power series $P_n(x)$ from section 4 of the paper."
  [n]
  (-> (s/power-series*
       (map-indexed
        (fn [i [num den two-m]]
          (* (g/expt (- 1) i)
             (/ num (* den (g/expt 2 two-m)))))
        (nm-seq n 0)))
      (s/inflate 2)))

(defn q-series
  "Returns the power series $Q_n(x)$ from section 4 of the paper.

  TODO this uses a goofy `cons 0` to do an efficient multiplication by `x`...
  bake this into the power series library."
  [n]
  (let [unshifted (-> (s/power-series*
                       (map-indexed
                        (fn [i [num den two-m+1]]
                          (* (g/expt (- 1) i)
                             (/ num (* den (g/expt 2 two-m+1)))))
                        (raw-coef** n 1)))
                      (s/inflate 2))]
    (s/power-series*
     (cons 0 unshifted))))

(defn alpha-series
  "Returns the power series $\\alpha(n)$ defined in the 'modified expansions'
  section of the paper."
  [n]
  (s/compose s/atan-series
             (g// (g/- (q-series n))
                  (p-series n))))

(defn beta-series
  "Returns the power series $\\beta(n)$ defined in the 'modified expansions'
  section of the paper."
  [n]
  (g/sqrt (g/+ (g/square (p-series n))
               (g/square (q-series n)))))

(comment
  ;; Once again, validate that our efficient version calculates the correct
  ;; series from the paper.
  (= '(1
       0
       (/ -9 128)
       0
       (/ 3675 32768)
       0
       (/ -2401245 4194304)
       0
       (/ 13043905875 2147483648))
     (v/freeze (take 9 (p-series 0))))

  (= '((/ -1 8)
       0
       (/ 75 1024)
       0
       (/ -59535 262144)
       0
       (/ 57972915 33554432)
       0)
     (v/freeze
      (take 8 (q-series 0)))))

;; ## Locating the Zeros
;;
;; This doesn't QUITE read like the paper because those power series are in
;; $\frac{1}{x}$, so we need to be a bit careful when trying to multiply both
;; sides by $x$. `identity` here stands for $\frac{1}{x}$.

(def x-factor
  (g/+ 1 (g/* s/identity (g/- (alpha-series 0)))))

(def p-factor
  (g/invert x-factor))

(def one-over-p
  (g/* s/identity p-factor))

(def one-over-x
  (s/revert one-over-p))

(def p-over-x
  (g// one-over-x s/identity))

(def x-over-p
  (g/invert p-over-x))

(defn x-as-fn-of-inv-p
  "This is a little wonky; because the series is in 1/p, we have to evaluate it
  first then add the final p term back on."
  [p]
  (g/+ (s/series 0 p)
       ((g// (g/- x-over-p 1) s/identity)
        (g/invert p))))

(comment
  (= (0 (/ (+ (* 8N (expt p 2)) 1) (* 8N p))
        0
        (/ -31N (* 384N (expt p 3)))
        0
        (/ 3779N (* 15360N (expt p 5)))
        0
        (/ -6277237N (* 3440640N (expt p 7)))
        0
        (/ 2092163573N (* 82575360N (expt p 9))))
     ((g/simplify (take 10 (x-as-fn-of-inv-p 'p))))))

Make a function to generate x↑3 etc

We need to:

Bind (sicmutils.structure/orientation->separator :sicmutils.structure/up) to some easier variable and expose it
Make a function like (raise-sym 'x 3) that will generate x↑3

So users don't have to type that in!

Constructor-simplification versions of numsymb functions

NOTE: Please also

These live in the kernel of scmutils. Find these by looking for enable-constructor-simplifications? and seeing which branch that gives. This is ON by default in scmutils, so we might find some behavior change. Here are some, but please check https://github.com/Tipoca/scmutils/blob/master/src/kernel/numsymb.scm#L34 for more.

;;; The following are more hairy construction simplifications for the
;;; most important functions.

(define (symb1:+ a1 a2)
  (cond ((sum? a1)
	 (cond ((sum? a2)
		(addup-args (append (operands a1) (operands a2)) '()))
	       ((difference? a2)
		(if (null? (cdr (operands a2)))
		    (addup-args (operands a1) (operands a2))
		    (addup-args (append (operands a1)
                                        (list (car (operands a2))))
				(cdr (operands a2)))))
	       (else (addup-args (append (operands a1)
                                         (list a2)) '()))))
	((difference? a1)
	 (if (null? (cdr (operands a1)))
	     (cond ((sum? a2) (addup-args (operands a2) (operands a1)))
		   ((difference? a2)
		    (if (null? (cdr (operands a2)))
			(addup-args '() (append (operands a1) (operands a2)))
			(addup-args (list (car (operands a2)))
				    (append (operands a1) (cdr (operands a2))))))
		   (else (addup-args (list a2) (operands a1))))
	     (cond ((sum? a2)
		    (addup-args (append (list (car (operands a1))) (operands a2))
				(cdr (operands a1))))
		   ((difference? a2)
		    (if (null? (cdr (operands a2)))
			(addup-args (list (car (operands a1)))
				    (append (cdr (operands a1))
                                            (operands a2)))
			(addup-args (list (car (operands a1))
                                          (car (operands a2)))
				    (append (cdr (operands a1))
                                            (cdr (operands a2))))))
		   (else (addup-args (list (car (operands a1)) a2)
				     (cdr (operands a1)))))))
	(else
	 (cond ((sum? a2)
		(addup-args (append (list a1) (operands a2)) '()))
	       ((difference? a2)
		(if (null? (cdr (operands a2)))
		    (addup-args (list a1) (operands a2))
		    (addup-args (append (list a1) (list (car (operands a2))))
				(cdr (operands a2)))))
	       (else (addup-args (list a1 a2) '()))))))

(define (addup-args pos neg)
  (define (make-answer sum pos neg)
    (if (zero? sum)
	(if (null? pos)
	    (if (null? neg)
		:zero
		(if (null? (cdr neg))
		    `(- ,(car neg))
		    `(- (+ ,@neg))))
	    (if (null? neg)
		(if (null? (cdr pos))
		    (car pos)
		    `(+ ,@pos))
		(if (null? (cdr pos))
		    (if (null? (cdr neg))
                        (if (equal? (car pos) (car neg))
                            :zero
                            `(- ,(car pos) ,(car neg)))
			`(- ,(car pos) (+ ,@neg)))
		    (if (null? (cdr neg))
			`(- (+ ,@pos) ,(car neg))
                        (if (equal? pos neg)
                            :zero
                            `(- (+ ,@pos) (+ ,@neg)))))))
	(if (null? pos)
	    (if (null? neg)
		sum
		(if (null? (cdr neg))
		    `(- ,sum ,(car neg))
		    `(- ,sum (+ ,@neg))))
	    (if (null? neg)
		`(+ ,sum ,@pos)
		(if (null? (cdr neg))
		    `(- (+ ,sum ,@pos) ,(car neg))
		    `(- (+ ,sum ,@pos) (+ ,@neg)))))))
  (let plp ((p pos) (sum :zero) (respos '()))
    (cond ((null? p)
	   (let nlp ((n neg) (sum sum) (resneg '()))
	     (cond ((null? n)
		    (make-answer sum
				 (reverse respos)
				 (reverse resneg)))
		   ((number? (car n))
		    (nlp (cdr n) (- sum (car n)) resneg))
		   (else
		    (nlp (cdr n) sum (cons (car n) resneg))))))
	  ((number? (car p))
	   (plp (cdr p) (+ sum (car p)) respos))
	  (else
	   (plp (cdr p) sum (cons (car p) respos))))))

(define (symb1:* a1 a2)
  (cond ((product? a1)
	 (cond ((product? a2)
		(mulup-args (append (operands a1) (operands a2)) '()))
	       ((quotient? a2)
		(if (null? (cdr (operands a2)))
		    (mulup-args (operands a1) (operands a2))
		    (mulup-args (append (operands a1)
                                        (list (car (operands a2))))
				(cdr (operands a2)))))
	       (else (mulup-args (append (operands a1) (list a2)) '()))))
	((quotient? a1)
	 (if (null? (cdr (operands a1)))
	     (cond ((product? a2) (mulup-args (operands a2) (operands a1)))
		   ((quotient? a2)
		    (if (null? (cdr (operands a2)))
			(mulup-args '() (append (operands a1) (operands a2)))
			(mulup-args (list (car (operands a2)))
				    (append (operands a1) (cdr (operands a2))))))
		   (else (mulup-args (list a2) (operands a1))))
	     (cond ((product? a2)
		    (mulup-args (append (list (car (operands a1))) (operands a2))
				(cdr (operands a1))))
		   ((quotient? a2)
		    (if (null? (cdr (operands a2)))
			(mulup-args (list (car (operands a1)))
				    (append (cdr (operands a1))
                                            (operands a2)))
			(mulup-args (list (car (operands a1))
                                          (car (operands a2)))
				    (append (cdr (operands a1))
                                            (cdr (operands a2))))))
		   (else (mulup-args (list (car (operands a1)) a2)
				     (cdr (operands a1)))))))
	(else
	 (cond ((product? a2)
		(mulup-args (append (list a1) (operands a2)) '()))
	       ((quotient? a2)
		(if (null? (cdr (operands a2)))
		    (mulup-args (list a1) (operands a2))
		    (mulup-args (append (list a1) (list (car (operands a2))))
				(cdr (operands a2)))))
	       (else (mulup-args (list a1 a2) '()))))))

(define (mulup-args pos neg)
  (define (make-product numfact factors)
    (if (null? factors)
        numfact
        (if (one? numfact)
            (if (null? (cdr factors))
                (car factors)
                `(* ,@factors))
            (if (null? (cdr factors))
                `(* ,numfact ,(car factors))
                `(* ,numfact ,@factors)))))
  (define (make-answer pfactor pos nfactor neg)
    (let ((num (make-product pfactor pos))
          (den (make-product nfactor neg)))
      (cond ((and (number? den) (zero? den))
             (error "zero divide in mulup-args"))
            ((and (number? num) (zero? num))
             :zero)
            ((and (number? den) (one? den))
             num)
            ((and (number? num) (number? den))
             (/ num den))
            ((equal? num den)
             :one)
            (else
             `(/ ,num ,den)))))
  (let plp ((p pos) (pfactor :one) (respos '()))
    (cond ((null? p)
	   (let nlp ((n neg) (nfactor :one) (resneg '()))
	     (cond ((null? n)
		    (make-answer pfactor
				 (reverse respos)
                                 nfactor
				 (reverse resneg)))
		   ((number? (car n))
		    (nlp (cdr n) (* nfactor (car n)) resneg))
		   (else
		    (nlp (cdr n) nfactor (cons (car n) resneg))))))
	  ((number? (car p))
	   (plp (cdr p) (* pfactor (car p)) respos))
	  (else
	   (plp (cdr p) pfactor (cons (car p) respos))))))
	   
(define (symb1:- a1 a2)
  (cond ((sum? a1)
	 (cond ((sum? a2)
		(addup-args (operands a1) (operands a2)))
	       ((difference? a2)
		(if (null? (cdr (operands a2)))
		    (addup-args (append (operands a1) (operands a2)) '())
		    (addup-args (append (operands a1) (cdr (operands a2)))
				(list (car (operands a2))))))
	       (else (addup-args (operands a1) (list a2)))))
	((difference? a1)
	 (if (null? (cdr (operands a1)))
	     (cond ((sum? a2)
                    (addup-args '() (append (operands a1) (operands a2))))
		   ((difference? a2)
		    (if (null? (cdr (operands a2)))
			(addup-args (operands a2) (operands a1))
			(addup-args (cdr (operands a2))
				    (append (operands a1)
					    (list (car (operands a2)))))))
		   (else (addup-args '() (append (operands a1) (list a2)))))
	     (cond ((sum? a2)
		    (addup-args (list (car (operands a1)))
				(append (cdr (operands a1)) (operands a2))))
		   ((difference? a2)
		    (if (null? (cdr (operands a2)))
			(addup-args (append (list (car (operands a1)))
                                            (operands a2))
				    (cdr (operands a1)))
			(addup-args (cons (car (operands a1))
                                          (cdr (operands a2)))
				    (append (cdr (operands a1))
					    (list (car (operands a2)))))))
		   (else (addup-args (list (car (operands a1)))
				     (append (cdr (operands a1)) (list a2)))))))
	(else
	 (cond ((sum? a2)
		(addup-args (list a1) (operands a2)))
	       ((difference? a2)
		(if (null? (cdr (operands a2)))
		    (addup-args (append (list a1) (operands a2)) '())
		    (addup-args (append (list a1) (cdr (operands a2)))
				(list (car (operands a2))))))
	       (else (addup-args (list a1) (list a2)))))))

(define (symb1:/ a1 a2)
  (cond ((product? a1)
	 (cond ((product? a2)
		(mulup-args (operands a1) (operands a2)))
	       ((quotient? a2)
		(if (null? (cdr (operands a2)))
		    (mulup-args (append (operands a1) (operands a2)) '())
		    (mulup-args (append (operands a1) (cdr (operands a2)))
				(list (car (operands a2))))))
	       (else (mulup-args (operands a1) (list a2)))))
	((quotient? a1)
	 (if (null? (cdr (operands a1)))
	     (cond ((product? a2)
                    (mulup-args '() (append (operands a1) (operands a2))))
		   ((quotient? a2)
		    (if (null? (cdr (operands a2)))
			(mulup-args (operands a2) (operands a1))
			(mulup-args (cdr (operands a2))
				    (append (operands a1)
					    (list (car (operands a2)))))))
		   (else (mulup-args '() (append (operands a1) (list a2)))))
	     (cond ((product? a2)
		    (mulup-args (list (car (operands a1)))
				(append (cdr (operands a1)) (operands a2))))
		   ((quotient? a2)
		    (if (null? (cdr (operands a2)))
			(mulup-args (append (list (car (operands a1)))
                                            (operands a2))
				    (cdr (operands a1)))
			(mulup-args (cons (car (operands a1))
                                          (cdr (operands a2)))
				    (append (cdr (operands a1))
					    (list (car (operands a2)))))))
		   (else (mulup-args (list (car (operands a1)))
				     (append (cdr (operands a1)) (list a2)))))))
	(else
	 (cond ((product? a2)
		(mulup-args (list a1) (operands a2)))
	       ((quotient? a2)
		(if (null? (cdr (operands a2)))
		    (mulup-args (append (list a1) (operands a2)) '())
		    (mulup-args (append (list a1) (cdr (operands a2)))
				(list (car (operands a2))))))
	       (else (mulup-args (list a1) (list a2)))))))

turn compose, arg-shift and arg-scale into generic functions

We have enough implementations now that this makes sense!

power-series
rational function
polynomial
function

all have sensible arg-shift and arg-scale, that keep their types around. Rational functions, series and polynomials can compose without needing to convert to functions too.

Add Quad (DoubleDouble) implementation, to boost precision without fully committing to Algebraic

Here's a nice library that does this now: https://github.com/Kerbaya/ieee754lib

Related Spire issue: typelevel/spire#236

In scmutils, this lives in numerics/quad/quad-flo.scm:

;;;     Fundamental quad manipulation
;;;        given double arguments.

;;; Knuth's theorem B provides a method of adding 2 floats
;;;  to produce a sum and an error estimate.
;;;   It takes 6 floating operations and no comparisions.

(define (+B u v k)			;k=(lambda (sum error) ...)
  (let ((s (flo:+ u v)))
    (let ((up (flo:- s v)))
      (let ((vpp (flo:- s up)))
        (k s (flo:+ (flo:- u up) (flo:- v vpp)))))))


(define (-B u v k)			;k=(lambda (sum error) ...)
  (let ((s (flo:- u v)))
    (let ((up (flo:+ s v)))
      (let ((vpp (flo:- s up)))
        (k s
           (flo:- (flo:- u up)
                  (flo:+ v vpp)))))))


;;; Multiplication of two floating operations, transcribed from
;;;  Quinn & Tremaine 1990.

(define IEEE-P 53)			;see standard.
(define tremaine-j
  (ceiling (/ IEEE-p 2)))
(define tremaine-c
  (exact->inexact (+ (expt 2 tremaine-j) 1)))

(define (*T u v k)			;k=(lambda (prod error) ...)
  (let* ((up (flo:* u tremaine-c))
         (u1 (flo:- u up))
         (u1 (flo:+ u1 up))
         (u2 (flo:- u u1))
         (vp (flo:* v tremaine-c))
         (v1 (flo:- v vp))
         (v1 (flo:+ v1 vp))
         (v2 (flo:- v v1))
         (prod (flo:* u v))
         (e1 (flo:* u1 v1))
         (e2 (flo:* u2 v2))
         (e12 (flo:* u1 v2))
         (e21 (flo:* u2 v1))
         (err (flo:- e1 prod))
         (err (flo:+ err e12))
         (err (flo:+ err e21))
         (perr (flo:+ err e2)))
    (k prod perr)))


(define (/T u v k)			;k=(lambda (prod error) ...)
  (let ((quot (flo:/ u v)))
    (*T quot v
        (lambda (rat err)
          (let* ((qerr (flo:- u rat))
                 (qerr (flo:/ (flo:- qerr err) v)))
            (k quot qerr))))))

;;;  Quad arithmetic -- a quad is a cons of the high and low part.

(define (quad:+ q1 q2)
  (let ((q1h (car q1))
        (q1l (cdr q1))
        (q2h (car q2))
        (q2l (cdr q2)))
    (+B q1h q2h
        (lambda (s e)
          (+B s (flo:+ (flo:+ q1l q2l) e)
              cons)))))

(define (quad:- q1 q2)
  (let ((q1h (car q1))
        (q1l (cdr q1))
        (q2h (car q2))
        (q2l (cdr q2)))
    (-B q1h q2h
        (lambda (s e)
          (+B s (flo:+ (flo:- q1l q2l) e)
              cons)))))

(define (quad:* q1 q2)
  (let ((q1h (car q1))
        (q1l (cdr q1))
        (q2h (car q2))
        (q2l (cdr q2)))
    (*T q1h q2h
        (lambda (p e)
          (+B p
              (flo:+ (flo:+ (flo:* q1l q2h) (flo:* q1h q2l)) e)
              cons)))))

(define (quad:/ q1 q2)
  (let ((q1h (car q1))
        (q1l (cdr q1))
        (q2h (car q2))
        (q2l (cdr q2)))
    (/T q1h q2h
        (lambda (q e)
          (+B q
              (flo:+ (flo:+ (flo:/ q1l q2h)
                            (flo:* -1.0 (flo:* q (/ q2l q2h))))
                     e)
              cons)))))

(define (->quad r)
  (cons (exact->inexact r) 0.0))

(define quad:0 (->quad 0.0))
(define quad:1 (->quad 1.0))
(define quad:2 (->quad 2.0))
(define quad:3 (->quad 3.0))
(define quad:4 (->quad 4.0))


(define *quad-epsilon*
  (flo:* *machine-epsilon*
         *machine-epsilon*))

(define (quad:square x)
  (quad:* x x))

(define (quad:sqrt x)
  (let ((xh (car x))
        (xl (cdr x)))
    (let lp ((guess (->quad (sqrt xh))))
      (let ((new (quad:/ (quad:+ guess (quad:/ x guess)) quad:2)))
        (if (= (car guess) (car new))
            new
            (lp new))))))


(define quad:sin-bugger-factor
  (expt 2 (* -1 IEEE-p)))

(define (quad:sin x)
  (let ((xh (quad:high x)))
    (if (< (flo:abs xh) quad:sin-bugger-factor)
        x
        (let ((y (quad:sin (quad:/ x quad:3))))
          (quad:- (quad:* quad:3 y)
                  (quad:* quad:4
                          (quad:* y
                                  (quad:* y y))))))))

(define (alternate-sine x)
  (let ((xsq (quad:* x x)))
    (let lp ((n 1) (xp x) (f quad:1) (sum quad:0))
      (let ((term (quad:/ xp f)))
        (if (< (abs (quad:high term)) *quad-epsilon*)
            sum
            (lp (+ n 2)
                (quad:* (->quad -1.0) (quad:* xsq xp))
                (quad:* f (->quad (* (+ n 1) (+ n 2))))
                (quad:+ sum term)))))))


(define quad:high car)
(define quad:low cdr)

#| probably wrong

(define rat:pi 314159265358979323846264338327950288419716939937510/100000000000000000000000000000000000000000000000000)

(define (quad:floor x)
(floor (quad:high x)))

(define (quad:truncate x)
(->quad (truncate (quad:high x))))  ;;; what if the error has an integer part?

(define continued-wallp true)

(define (best-fraction-from-quad x maxint)
(let* ((a0 (inexact->exact (floor (quad:high x))))
(r0 (quad:- x (quad:truncate x))))
(if (<= (quad:high r0) (* (abs (quad:high x)) *quad-epsilon*))
a0
(let* ((d (quad:/ quad:1 r0))
(a1 (inexact->exact (truncate (quad:high d)))))
(let lp ((s d) (last-p a0) (last-q 1) (p (+ (* a0 a1) 1)) (q a1))
(if continued-wallp
(write-line (list s last-p last-q p q)))
(if (> q maxint)
(/ last-p last-q)
(let ((r (quad:- s (quad:truncate s))))
(if (<= (quad:high r) (* (quad:high s) *quad-epsilon*))
(/ p q)
(let* ((new-d (quad:/ quad:1 r))
(a (inexact->exact (truncate (quad:high new-d)))))
(lp new-d
p
q
(+ last-p (* a p))
(+ last-q (* a q))))))))))))

|#

(define (cquad:+ z1 z2)
  (map quad:+ z1 z2))

(define (cquad:- z1 z2)
  (map quad:- z1 z2))

(define (cquad:* z1 z2)
  (list (quad:- (quad:* (car z1) (car z2))
                (quad:* (cadr z1) (cadr z2)))
        (quad:+ (quad:* (car z1) (cadr z2))
                (quad:* (cadr z1) (car z2)))))

(define (cquad:/ z1 z2)
  (let ((m (quad:+ (quad:* (car z2) (car z2))
                   (quad:* (cadr z2) (cadr z2)))))
    (list (quad:/ (quad:+ (quad:* (car z1) (car z2))
                          (quad:* (cadr z1) (cadr z2)))
                  m)
          (quad:/ (quad:- (quad:* (cadr z1) (car z2))
                          (quad:* (car z1) (cadr z2)))
                  m))))

(define (q->cq q)
  (list q quad:0))

(define (c->cq c)
  (list (->quad (real-part c))
        (->quad (imag-part c))))


(define (quad:magnitude z)
  (quad:sqrt (quad:+ (quad:* (car z) (car z)) (quad:* (cadr z) (cadr z)))))

;;; some other trig functions

(define (compute-quad:pi)
  (let lp ((guess (->quad pi)))
    (write-line guess)
    (let ((new-guess (quad:+ guess (quad:sin guess))))
      (if (< (abs (quad:high (quad:- new-guess guess))) *quad-epsilon*)
          new-guess
          (lp new-guess)))))

(define quad:pi '(3.141592653589793 . 1.2246467991473535e-16))
(define quad:pi/2 '(1.5707963267948966 . 6.123233995736767e-17))

(define (quad:cos x)
  (quad:sin (quad:+ x quad:pi/2)))

(define (quad:tan x) (quad:/ (quad:sin x) (quad:cos x)))

(define (quad:atan x)
  (let ((fp (quad:+ quad:1 (quad:* x x)))
        (guess (->quad (atan (quad:high x)))))
    (let lp ((guess guess))
      (let ((new-guess (quad:- guess (quad:/ (quad:- (quad:tan guess) x) fp))))
        (if (< (abs (quad:high (quad:- new-guess guess))) *quad-epsilon*)
            new-guess
            (lp new-guess))))))

generate a standalone REPL binary using GraalVM

Imagine an instantly-firing-up terminal calculator. We can do this GraalVM. Some resources for thinking this through:

https://github.com/lread/clj-graal-docs
For something even more ambitious, we might make a Terminal User Interface with Clojurescript. See these templates for beautiful examples.

port remaining ODE initial value solvers from scmutils

scmutils supports the following ODE solvers in its numerics package:

Bulirsch-Stoer
"quality controlled" runge kutta 4: https://people.cs.clemson.edu/~dhouse/courses/817/papers/adaptive-h-c16-2.pdf
quality-controlled Trapezoidal
quality-controlled "euler" (qcceuler?)
Explicit Gear

We currently only have support for Bulirsch-Stoer (the best, default!), via Colin Smith's "ODEX" library in Javascript and the Apache-Commons implementation on the JVM.

Fleshing this out with the care we used for the quadrature package would be brilliant.

here are the files:

add a "version" variable that the API can query from the REPL

It would be lovely if Nextjournal and friends had some sicmutils-version variable inside env that they could query to figure out what is actually running! I'm not sure of the best pattern here. Does anyone have any pointers?

Clojure

Here is how to do it on the JVM, modeled after *clojure.core/clojure-version*:

(let [version-string (slurp (.getResourceAsStream
                             (clojure.lang.RT/baseLoader)
                             "SICMUTILS_VERSION"))
      [_ major minor incremental qualifier snapshot]
      (re-matches
       #"(\d+)\.(\d+)\.(\d+)(?:-([a-zA-Z0-9_]+))?(?:-(SNAPSHOT))?"
       version-string)
      version {:major       (Integer/valueOf ^String major)
               :minor       (Integer/valueOf ^String minor)
               :incremental (Integer/valueOf ^String incremental)
               :qualifier   (if (= qualifier "SNAPSHOT") nil qualifier)}]
  (def ^:dynamic *sicmutils-version*
    (if (.contains version-string "SNAPSHOT")
      (assoc version :interim true)
      version)))

Use spec or similar with `defintegrator` to tag various methods as accepting open or closed intervals

The Quadrature interface contributed in sicmutils/sicmutils#129 lets you specify in the options map whether your interval is open, closed or half-open.

Certain integrators can only handle certain cases - the Midpoint and Milne rules, or the open version of Romberg, handle open cases, while all of the trapezoid-based methods only handle closed intervals.

"Adaptive" can take an open and closed integrator and make a new integrator that can handle any of these combinations.

Proposal

We might modify the defintegrator interface in sicmutils.numerical.quadrature.common to take a Clojure Spec declaring what kinds of intervals the defined integrator can handle.

We could also make a simpler version of adaptive that only splits a single time, to allow integrators like Romberg to handle all of the kinds of intervals (open-closed would split at an internal midpoint and attack the left and right with fully open and fully closed, respectively.)

Add generic <, <=, >, >= to match =

We need these to get proper comparisons working between numbers and literal-number and Differential instances. We have sicmutils.value/compare, so these should be simple to add and document.

I think they should replace the clojure.core versions in sicmutils.env.

Also: This takes us into territory well covered by clojure.algo.generic https://github.com/clojure/algo.generic/blob/master/src/main/clojure/clojure/algo/generic/comparison.clj. We should in fact add the full suite of operations, following their lead.

Document ideal interactive editor experience, setup to achieve it

The task here is to get a writeup of what we want interactive math and physics explorations, on their own and in the context of a textbook, to feel like with a Clojure library like sicmutils. We can start this in a discussion here, but the goal is to populate:

This wiki page https://github.com/sicmutils/sicmutils/wiki/Interactive-Editor-Setup
the official docs once we get something stable

Experimentally determine `p` in Richardson extrapolation

This MathOverflow comment mentions a nice way to do this:

Also, Sussman in "Abstraction in Numerical Methods" (bottom of page 6) notes:

By the way, if the sequence $p_1, p_2, p_3, ...$ is not known in advance... Alternatively, the appropriate exponent $p$ for a given column can be inferred numerically from the early terms of that column; we have done that in our library, but do not pursue it here.

I can't find the actual library he's mentioning, so we'll have to figure this out ourselves!

Port, optimize AD benchmarks from Siskind + Pearlmutter's work

@barak and @qobi presented on a series of AD benchmarks that they performed as part of their work on Stalingrad:

"Using Programming Language Theory to MakeAutomatic Differentiation Sound and Efficient" gives a nice overview of the benchmarks
this page stores the actual code
this presentation from NIPS 2017 presents the benchmarks

scmutils didn't do terribly well for a few reasons, one of which is no compilation of the derivatives... but no matter!

The task here is to get the SICMUtils code for these benchmarks tidied up, measured and going faster.

I took a pass at getting the code mechanically ported from the two straightforward benchmarks that worked in scmutils, particle-FF and saddle-FF. These are described in the first paper linked above.

NOTE: probabilistic-lambda-calculus, probabilistic-prolog and backprop all have forward-mode examples too that would work with some porting. Also once I get sicmutils/sicmutils#226 up and running we can get the rest of the benchmarks running here too. I know they'd be faster with a more intelligent compile implementation that could handle multiple arguments more gracefully. This would need to generate multiple nested functions in forward mode.

;; Run these from sicmutils.env... I have some defensive g/ prefixes that we probably don't need.

(defn multivariate-argmin [f x]
  (let [g (D f)]
    (loop [x   x
           fx  (f x)
           gx  (g/transpose (g x))
           eta 1e-5
           i   0]
      (cond (<= (compare (g/abs gx) 1e-5) 0)
            x

            (= i 10)
            (recur x fx gx (g/* 2.0 eta) 0)

            :else
		        (let [x-prime (g/- x (g/* eta gx))]
		          (if (<= (compare (g/abs (- x x-prime)) 1e-5) 0)
			          x
			          (let [fx-prime (f x-prime)]
			            (if (< (compare fx-prime fx) 0)
			              (recur x-prime fx-prime (g/transpose (g x-prime)) eta (+ i 1))
			              (recur x fx gx (/ eta 2.0) 0)))))))))

(defn naive-euler [w]
  (let [charges      [[10.0 (- 10.0 w)] [10.0 0.0]]
	      x-initial    [0.0 8.0]
	      xdot-initial [0.75 0.0]
	      delta-t      1e-1
	      p (fn [x]
	          (transduce (map (fn [c] (/ 1.0 (g/abs (- x c)))))
                       g/+
                       0.0
                       charges))]
    (loop [x    x-initial
           xdot xdot-initial]
      (let [[x' y' :as x-new] (g/+ x (g/* delta-t xdot))]
		    (if (g/negative? y')
          (let [delta-t-f (g// (g/- 0.0 (second x))
					                     (second xdot))
			          x-t-f     (g/+ x (g/* delta-t-f xdot))]
			      (g/square (first x-t-f)))
		      (let [xddot (g/* -1.0 (g/transpose ((D p) x)))]
            (recur x-new (g/+ xdot (g/* delta-t xddot)))))))))

(defn multivariate-argmax [f x]
  (multivariate-argmin (fn [x] (g/- (f x))) x))

(defn multivariate-max [f x]
  (f (multivariate-argmax f x)))

(defn particle-FF []
  (multivariate-argmin naive-euler 0.0))

(defn saddle-FF []
  (let [start [1.0 1.0]
	      f     (fn [[x1 y1] [x2 y2]]
	              (- (+ (g/square x1) (g/square y1))
                   (+ (g/square x2) (g/square y2))))
	      x1* (multivariate-argmin
	           (fn [x1]
	             (multivariate-max
                (fn [x2] (f x1 x2))
	              start))
	           start)
	      x2* (multivariate-argmax
	           (fn [x2] (f x1* x2)) start)]
    [x1* x2*]))

And the results (with nothing to compare it to, since this is a new M1 mac and I haven't run any other benchmarks here):

sicmutils.env> (time (saddle-FF))
"Elapsed time: 9141.064291 msecs"
[(up 8.246324826140356E-6 8.246324826140356E-6) (up 8.246324826140356E-6 8.246324826140356E-6)]

sicmutils.env> (time (particle-FF))
"Elapsed time: 28148.853958 msecs"
0.20719187464861197

generic support underneath Complex type

Right now, the Complex datatype is based on

Double in Clojure,
js/Number in Clojurescript due to the Complex.js dependency

josdejong/mathjs#694 has some discussion about the problems this has created in math.js. The big one is that you can't extend complex arithmetic to arbitrary precision.

We COULD probably do this if we rebuilt Complex in Clojurescript, on top of the generic arithmetic. (I think you might be able to get complex matrices for free??)

https://github.com/infusion/Complex.js/blob/master/complex.js is not that big of a library. porting it over to cljc, and backing it with generic operations, would be a really nice project.

Implement compose for rational functions, polynomials

rcf.scm has a compose function that converts to, roughly, this:

;; The following function still reads slightly mystical; this comes from rcf.scm
;; in the scmutils library. I'm sure with some staring and understanding we
;; could build a version that allows composition with any variable in the
;; left [[RationalFunction]] instance.

(defn compose
  "only plugs r2 in for the principal indeterminate.

  TODO we COULD do a version that composes with a different one??"
  [r1 r2]
  {:pre [(rational-function? r1)]}
  (if-not (rational-function? r2)
    (g/div (p/evaluate (bare-u r1) [r2])
           (p/evaluate (bare-v r1) [r2]))
    (let [nr1 (bare-u r1)
          nr2 (bare-u r2)
          dr1 (bare-v r1)
          dr2 (bare-v r2)
          dn  (p/degree nr1)
          dd  (p/degree dr1)
          narity (+ (p/arity dr1) 1)
          nnr1 (p/extend 1 (p/reciprocal nr1))
          ndr1 (p/extend 1 (p/reciprocal dr1))
          scales [(second (p/new-variables narity)) 1]
          pn (p/evaluate (p/reciprocal
                          (p/arg-scale nnr1 scales))
                         [nr2 dr2])
          pd (p/evaluate (p/reciprocal
                          (p/arg-scale ndr1 scales))
                         [nr2 dr2])]
      (cond (> dn dd) (g/div pn (p/mul (p/expt dr2 (- dn dd)) pd))
            (< dn dd) (g/div (p/mul (p/expt dr2 (- dd dn)) pn) pd)
            :else (g/div pn pd)))))

But it's a bit mystical and I don't see how it works yet. The task is to get this working, and get a similar version working for polynomials.

Add native Trig functions for js/BigInt, goog.math.{Long, Integer}

add trig functions for js/BigInt, goog.math.{Long, Integer}.

These work as of sicmutils/sicmutils#168 , but they work by converting these types to js/Number, which won't work if they're too large.

This is fine if we can first confine the numbers back down to a reasonable range.

See Roylance's "Expressing Mathematical Routines Constructively" for a scheme implementation that does this.

Continued fraction data type

@alexgian sent me a lovely paper by the legendary Bill Gosper paper on the continued fraction representation of real numbers.

This would be an AWESOME data type to add to SICMUtils, and a great target for a literate-style essay, like the numerical methods work or Power Series.

From the site:

"Multiprecision arithmetic algorithms usually represent real numbers as decimals, or perhaps as their base-2n analogues. But this representation has some puzzling properties. For example, there is no exact representation of even as simple a number as one-third. Continued fractions are a practical but little-known alternative.

"Continued fractions are a representation of the real numbers that are in many ways more mathematically natural than the usual decimal or binary representations. All rational numbers have simple representations, and so do many irrational numbers, such as sqrt(2) and e1. One reason that continued fractions are not often used, however, is that it's not clear how to involve them in basic operations like addition and multiplication. This was an unsolved problem until 1972, when Bill Gosper found practical algorithms for continued fraction arithmetic."

Flesh out literal fn syntax for structured args, values to match scmutils

We do support such literal functions, but with an undocumented syntax for domain and range. This should be harmonized with the book.

Add faster, sparse polynomial GCD implementation

The code for this lives in the simplify package, in

sparse-gcd.scm
sparse-interpolate.scm

Candidates:

"A Fast Parallel Sparse Polynomial GCD Algorithm" http://www.cecm.sfu.ca/CAG/papers/jscgcd2018.pdf
Zippel, "Probabilistic Algorithms for Sparse Polynomials" https://rjlipton.files.wordpress.com/2009/11/zippel.pdf

UPDATE: I started working on this a bit as part of the monster sicmutils/sicmutils#341 effort, then backed off. Here's the code I had started sketching to make this work:

:else clause in GCD dispatch...

(def ^:dynamic *euclid-breakpoint-arity* 3)
(def ^:dynamic *gcd-cut-losses* nil)

(or (with-limited-time [1.0 :seconds]
		  (fn [] (sparse-gcd u v)))

	  (with-limited-time [1.0 :seconds]
		  (fn [] (classical-gcd u v)))

	  (with-limited-time [100.0 :seconds]
		  (fn [] (sparse-gcd u v)))

	  (and (< arity *euclid-breakpoint-arity*)
		     (with-limited-time [100.0 :seconds]
		       (fn [] (classical-gcd u v))))
	  (sparse-gcd u v)

	  (if *gcd-cut-losses*
		  (or (with-limited-time *gcd-cut-losses*
			      (fn [] (classical-gcd u v)))
		      1)
		  (classical-gcd u v)))

Modularize sicmutils?

@littleredcomputer , curious for your thoughts here.

As I work through the cljs port sicmutils/sicmutils#40 , clear layers are emerging:

the generic arithmetic system built on types available in the OS, and the generics themselves
extensions to the arithmetic system, like Complex and friends, that rely on dependencies.
compound types like Matrix, Series, etc
the analyzer and simplifier
calculus modules
book examples

For some reason it's much more natural in other build systems like SBT to offer modules: https://github.com/twitter/algebird so maybe this is not a good idea here... but it does feel like it might be a nice move to explicitly break the library into these modules, so that someone could depend on JUST the numerical stuff, or JUST the automatic differentiation, without pulling in the rest of the dependencies for the book.

If we start adding better tools for visualization into the library, like a Quil dependency, I can see this becoming even more true.

EDIT: Here's an example, from the Clojurians slack, of how this might work: https://github.com/juxt/crux/blob/master/project.clj

Extend protocols to goog.math.Matrix, Vec2, Vec3

These come packaged with Google Closure; it would be great to have these available as another matrix library in SICMUtils. If we don't depend on the namespace by default, then anyone compiling Clojurescript won't see it show up in their output.

@littleredcomputer , wdyt... is this the sort of thing that would be better served in a new, small project? sicmutils-closure, to pick a very confusing name choice?

EDIT I think emmy.goog.math is a good namespace here.

Change to 'simplify automatically', add easy ways to prevent this from happening

This should be user-controllable, at the expression level and the system level.

You really just want simplification the FIRST time. Then if you call "factor" or something, you want that to not run through the simplifier again and UN-factor. So explicit "simplify" style calls should not trigger further simplification. How to make that happen?

Pros:

some simplifications already happen in the constructor. If you are going to turn them off full-stop, that's not possible now
simplified forms are the default for scmutils

Cons:

When you print something and see that it's very simple, then pass this thing as an argument, you might not realize that the UNSIMPLIFIED version is getting passed in, which might have bad effects.

Implement lcm, quotient and remainder for rational functions

These three are obviously missing generics; the implementations are not necessarily obvious!

Implement "Lazy Multivariate Higher-Order Forward-Mode AD" for cheap higher-order derivatives

Thanks @qobi for the reference.

This is the solution to the problem of an exponential number of derivatives being calculated for (((expt D n) f) 'x) == the size of the power set of a set of size n (see the comment below).

The task is to implement the method described in this paper: https://engineering.purdue.edu/~qobi/papers/popl2007a.pdf

A method is presented for computing all higher-order partial derivatives of a multivariate function R^n → R. This method works by evaluating the function under a nonstandard interpretation, lifting reals to multivariate power series. Multivariate power series, with potentially an infinite number of terms with nonzero coefficients, are represented using a lazy data structure constructed out of linear terms. A complete implementation of this method in SCHEME is presented, along with a straightforward exposition, based on Taylor expansions, of the method’s correctness.

And THEN to figure out how to expose this behavior inside the D operator as it exists currently.

References

Here is a presentation from Mu Wang on how to take advantage of symmetries in reverse mode: "Reverse Hessian Algorithm With Preaccumulation"
Another Mu Wang paper, "Capitalizing on live variables: new algorithms forefﬁcient Hessian computation via automatic differentiation"
Mu Wang’s ReverseAD software based on his thesis: https://github.com/wangmu0701/ReverseAD
Poster from the Matthew Johnson paper on taylor-mode in JAX
Associated paper from Matthew Johnson, with a blistering review (and great references!) from Siskind: "Taylor-Mode Automatic Differentiation for Higher-Order Derivatives in JAX"
here's a recent Master's thesis that implements the higher-order derivatives using symbolics.jl, we could copy this approach

Implement unit systems and SI Units, fold in constants from scmutils

Hey!

disclaimer: There are many moving pieces of sicmutils I don't understand; I find this to be quite difficult.

NOTE: this paper, which discusses "dimensional algebra": https://fermatslibrary.com/s/can-one-take-the-logarithm-or-the-sine-of-a-dimensioned-quantity-or-a-unit-dimensional-analysis-involving-transcendental-functions#email-newslettter

Also note the scmutils files that implement this live in scmutils/units:

I'd love to see first class unit support. I asked in the Clojurians Slack #sicmutils channel, and got some pointers from @sritchie. Messages pasted for context:

@teodorlu

Have you given any thought to putting SI units on numbers; allowing 1 km / 30 minutes = 0.5 km/h and similar calculations? You talked a bit about the extensibility of the number systems (multimethods for +, *, ...) on the scicloj meetup.

Or, one could define plain SI units symbolically, and derived SI units in terms of base SI units. Perhaps adding defmethods and a new type is the wrong way to go.

@sritchie

it’s a brilliant idea! In fact the original scmutils has a “units” package that I would love to mentor / help someone port over into sicmutils

https://github.com/Tipoca/scmutils/tree/master/src/units

it takes a while to learn how scmutils works so it would be a tough go on your own, but if you’re interested this would be an AWESOME project

https://github.com/Tipoca/scmutils/blob/master/src/units/with-units.scm#L243

that would let us put all these constants in! https://github.com/Tipoca/scmutils/blob/master/src/units/constants.scm#L108

Another note: the https://github.com/Tipoca/scmutils/tree/master/src/units[`units`] folder hosts a bunch of physical constants, plus arithmetic on units and the ability to work with numbers and expressions that have units attached.

Allow for equation representation (ie symbolic `=`)

@hcarvalhoaves and CTW in Clojurians Slack noted that this would be great.

Some tasks:

make a constructor for = and inequality expressions
give them a new type, so they don't dispatch to the same functions as ::v/numeric
get these rendering appropriately in ->infix etc

We might want inequalities represented with a new record type, rather than an expression with an embedded <=. But they should certainly freeze to that.

`g/partial-derivative` for structures shouldn't use fn impl for NON-structural inputs

The current implementation of g/partial-derivative for ::s/structure instances is the same implementation used by ::function.

The reason for this, I think (and the reason it works!) Is that all of the machinery that derivative.cljc develops to allow for structured inputs (Jacobian support) gets to work a SINGLE time and call the structure as a function for each perturbed input in the incoming structure.

The alternative would be to implement partial-derivative on structures like this:

(defmethod g/partial-derivative [::s/structure v/seqtype]
  [structure selectors]
  (s/mapr #(g/partial-derivative % selectors) structure))

This would fix two problems I've noticed:

if a structure contains a PowerSeries, or really any structure whose derivative we can take cheaply WITHOUT calling derivative and using forward-mode AD, we currently always bypass that implementation. a PowerSeries in particular can take its own derivative easily, and would stay a PowerSeries. The current implementation actually calls it, dropping it back down to Series (and not appropriately shifting all of the entries!!)
If you use an explicit IPersistentVector that you haven't turned into a structure, the whole current approach will fail because a vector implements IFn by looking up an entry, no applying its args to contained functions.

Compare these two cases for an example:

sicmutils.calculus.derivative> ((D [g/sin g/cos]) 5)
Execution error (IllegalArgumentException) at sicmutils.calculus.derivative/derivative$fn (derivative.cljc:191).
Key must be integer

sicmutils.calculus.derivative> ((D (s/up g/sin g/cos)) 5)
(up 0.28366218546322625 0.9589242746631385)

[type] Multivariate power series

We have univariate power series. Here's some prior art on multivariate power series: https://github.com/michibo/pwrsrs

Add Real implementation based on Boehm's implementation

Here's a nice paper we can add: https://blog.acolyer.org/2020/10/02/toward-an-api-for-the-real-numbers/

I THINK this is implemented in Spire: https://github.com/typelevel/spire/blob/master/core/src/main/scala/spire/math/Real.scala

Benchmark and (possibly) speed up term-rewriting

@kovasb has some excellent experiments in https://github.com/kovasb/combinator on how to do very fast term rewriting in clj/cljs, 2-3x faster than Mathematica is able to support. @littleredcomputer , you may have studied this already (please close the issue if so!) but I'd love your thoughts on the zipper-based implementation here, and whether it might give us some advantages over the current implementation.

Some other Clojure term-rewriting systems:

https://github.com/mschuene/stratege, based on Stratego/XT
https://github.com/noprompt/meander, with rewrite rules

add reader literal for Literal instances

Right now these print as the thing that's wrapped, leading to confusing error messages like

0 does not equal 0

in tests. If we make these print as #sicm/expression, it will all be much clearer. Here's what it will take:

;;print-method
(let [rep (str "#" (.-sym (.-type s))
               " " (.toString s))]
  (.write w ^String rep))

;; clojurescript version
(let [rep (str "#" (.-sym (type s))
               " " (pr-str expression))]
  (-write writer rep))

Plus an entry in data_readers.clj.

NOTE that we also want print-dup methods for some of the types!

Make sicmutils.expression.compile backend pluggable

Right now the compile functions use a dictionary called compiled-fn-whitelist to associate symbols with functions. After simplification and common subexpression elimination, the compiler swaps these functions back in for their associated symbols and calls eval on the result.

There are a few alternatives a user might want to the current compiled-fn-whitelist:

Swap the GENERIC functions back in, so you can perform a simplification but still stay in flexible, generic land
super fast math on the JVM, via Geex or Fastmath

Geex looks like it might be the easiest to integrate.

Implement Brent's Method for root finding

I'm making a ticket about this because during my deep dive on Brent's Method for minimization, I discovered that the algorithm described on the wiki page is in fact wrong! https://en.wikipedia.org/wiki/Brent%27s_method

So the task here is to:

ask me what's wrong with it (if I take this ticket this one's solved)
study the Brent's method for minimization code
implement the root finding version, with appropriate shared code!

Add "Algebraic" data type - rational + roots as exact, valid operations

Here's a Spire implementation: https://github.com/typelevel/spire/blob/master/core/src/main/scala/spire/math/Algebraic.scala
a nice post about how to PRINT these types: https://fredrikj.net/blog/2021/03/printing-algebraic-numbers/

Nice notes describing this and related papers:

/**
 * Algebraic provides an exact number type for algebraic numbers. Algebraic
 * numbers are roots of polynomials with rational coefficients. With it, we can
 * represent expressions involving addition, multiplication, division, n-roots
 * (eg. `sqrt` or `cbrt`), and roots of rational polynomials. So, it is similar
 * [[Rational]], but adds roots as a valid, exact operation. The cost is that
 * this will not be as fast as [[Rational]] for many operations.
 *
 * In general, you can assume all operations on this number type are exact,
 * except for those that explicitly construct approximations to an Algebraic
 * number, such as `toBigDecimal`.
 *
 * For an overview of the ideas, algorithms, and proofs of this number type,
 * you can read the following papers:
 *
 *  - "On Guaranteed Accuracy Computation." C. K. Yap.
 *  - "Recent Progress in Exact Geometric Computation." C. Li, S. Pion, and C. K. Yap.
 *  - "A New Constructive Root Bound for Algebraic Expressions" by C. Li & C. K. Yap.
 *  - "A Separation Bound for Real Algebraic Expressions." C. Burnikel, et al.
 */

add symb:elementary-access? from Kernel

This function checks if the access-chain terminates in another structure or a final symbolic thing. Stare at this more to understand what it is doing... the args list is a bit confusing. But check the call site to see.

This particular version assumes that we are dealing with a frozen structure.

This is used ONLY in rules.cljc. Look there in the partial derivative simplifications for where check occurs.

(define (symb:elementary-access? access-chain args)
  (define (sea chain thing)
    (cond ((and (not (null? chain))
		(pair? thing)
		(or (eq? (car thing) 'up)
		    (eq? (car thing) 'down)))
	   (sea (cdr chain)
		(list-ref (cdr thing) (car chain))))
	  ((null? chain)
	   (or (not (pair? thing))
	       (not (or (eq? (car thing) 'up)
			(eq? (car thing) 'down)))))
	  (else #f)))
  (if (= (length args) 1)
      (sea access-chain (car args))
      (sea (cdr access-chain)
	   (list-ref args (car access-chain)))))

Add generic `fmap` operation

If we do this, and take a tip from clojure.algo.generic, it will be easier to write generic functions that do stuff like replace-tag from sicmutils.differential.

https://github.com/clojure/algo.generic/blob/master/src/main/clojure/clojure/algo/generic/functor.clj

Many types implement something like this, thanks to @littleredcomputer . Getting them together under one roof is a nice move.

Make cljs arithmetic auto-promoting

to match the CLJ behavior. Here's an example of how clj does this:

TODO gotta get promoting arithmetic in cljs! here is how https://github.com/clojure/clojure/blob/master/src/jvm/clojure/lang/Numbers.java#L488 file a ticket.

Operators should return named functions on application

This is how to do this:

sicmutils.env> (freeze (with-meta (D (literal-function 'f)) {:name "cake"}))
"cake"
sicmutils.env> (with-meta (D (literal-function 'f)) {:name "cake"})
#function[clojure.lang.AFunction/1]

Considering D, `(partial :key)` to treat map-shaped args as structural inputs

This would let us do ((partial :key) f) {:key 10 :val 100}. I am thinking more about supporting maps on the input...

Of course this ticket also implies that D should return a Jacobian shaped like a map. That breaks the contract that the output is compatible for contraction with an increment in the input.

Do we restrict this behavior to Grad, which keeps orientations preserved?

Add suppress-args, convert tests

This is a nice system to take some derivation and return an unapplied functional version of the result. Here is the code, and some examples from Scheme.

Once we have this we'll need to

go find all the spots that use an anonymous present function in the tests
convert them to instead use suppress-arguments
probably add a with-suppressed-args version that does the same but only inside of a let binding. Or maybe just call it suppress...

(define *suppressed-argument-list* '())

(define *suppressed-argument-list-counter* 0)

(define (suppress-arguments arguments)
  (let ((n (+ (length *suppressed-argument-list*) 1)))
    (set! *suppressed-argument-list-counter*
	        (+ *suppressed-argument-list-counter* 1))
    (set! *suppressed-argument-list*
	        (cons (cons arguments
		                  (symbol "args."
			                        *suppressed-argument-list-counter*))
		            *suppressed-argument-list*))
    n))

(define (show-suppressed-arguments)
  (pp (map (lambda (al)
	     `(,(cdr al) = ,@(car al)))
	   *suppressed-argument-list*)))

(define (clear-arguments)
  (set! *suppressed-argument-list* '())
  (set! *suppressed-argument-list-counter* 0)
  0)

(define (arg-suppressor expression)
  (if (pair? expression)
      (let ((v (assoc (cdr expression) *suppressed-argument-list*)))
	      (if v
	          (list (arg-suppressor (car expression)) (cdr v))
	          (cons (arg-suppressor (car expression))
		              (arg-suppressor (cdr expression)))))
      expression))


;;; For example

(let ((t 't) (xy (up 'x 'y)) (uv (up 'r 's)))
  (* (((partial 2) Hp) (up t uv (- (((partial 2) F1) t xy uv))))
     (((partial 2) ((partial 1) F1)) 't xy uv)))
#|
(down
 (+
  (*
   (((partial 1 0) ((partial 2 1) F1)) t (up x y) (up r s))
   (((partial 2 1) Hp)
    (up t
	(up r s)
	(down (* -1 (((partial 2 0) F1) t (up x y) (up r s)))
	      (* -1 (((partial 2 1) F1) t (up x y) (up r s)))))))
  ...mess...)
 ...mess...)
|#

;;; We choose arguments to suppress:

(suppress-arguments '((up t
			  (up r s)
			  (down (* -1 (((partial 2 0) F1) t (up x y) (up r s)))
				(* -1 (((partial 2 1) F1) t (up x y) (up r s)))))))
#| 1 |#

(suppress-arguments '(t (up x y) (up r s)))
#| 2 |#


;;; Now look at the pretty result:

(let ((t 't) (xy (up 'x 'y)) (uv (up 'r 's)))
  (* (((partial 2) Hp) (up t uv (- (((partial 2) F1) t xy uv))))
     (((partial 2) ((partial 1) F1)) 't xy uv)))
#|
(down
 (+ (* (((partial 2 0) Hp) args.1) (((partial 1 0) ((partial 2 0) F1)) args.2))
    (* (((partial 2 1) Hp) args.1) (((partial 1 0) ((partial 2 1) F1)) args.2)))
 (+ (* (((partial 2 0) Hp) args.1) (((partial 1 1) ((partial 2 0) F1)) args.2))
    (* (((partial 2 1) Hp) args.1) (((partial 1 1) ((partial 2 1) F1)) args.2))))
|#

(show-suppressed-arguments)
((args.2 = t (up x y) (up r s))
 (args.1 = (up t
	             (up r s)
	             (down (* -1 (((partial 2 0) F1) t (up x y) (up r s)))
		                 (* -1 (((partial 2 1) F1) t (up x y) (up r s)))))))

Add `known-real?`, other assumptions to Literal

This is an idea from scmutils that is not actually used there. The idea is to add metadata to Literal instances that states more specific claims we know about those expressions.

Actually, while typing this I see that we can potentially add assumptions like positive? etc so forms as they're simplified. So this is a bigger idea.

But for now, consider metadata that would make a literal respond to real?, integer? etc, perhaps with a different v/kind.

Here are notes from a previous think on this:

Known Real?

known-real? is a predicate in numbers.scm - this is a system for registering symbols as "real". The symbolic environment (from numsymb.scm) uses this information in the call to g/conjugate for symbols and just mirrors the term back out. We could of course do this in the other complex calls.

We could ALSO take this unused idea from express.scm:

(define (make-real-literal expression)
  (let ((e (make-numerical-literal expression)))
    (add-property! e 'real #t)
    e))

and use the metadata field that's now in Literal. BUT thinking now, we can actually just attach metadata to symbols themselves and have them respond directly also to known-real?.

UPDATE: symbolic expressions and symbols now both respond to with-meta. All we need here is:

implement the known-real? function
add functions that will allow us to add arbitrary tags to expressions, and decide on a convention here for :real? true and metadata like this.

From numbers.scm:

(define *known-reals* '())

(define (known-real? z)
  (cond ((structure? z)
         (s:forall known-real? z))
        ((matrix? z)
         (let ((m (m:num-rows matrix))
               (n (m:num-cols matrix))
               (mat (matrix->array z)))
           (let rowlp ((i 0))
             (if (fix:= i m)
                 #t
                 (let collp ((j 0))
                   (if (fix:= j n)
                       (rowlp (fix:+ i 1))
                       (if (known-real? (array-ref mat i j))
                           (collp (fix:+ j 1))
                           #f)))))))
        ((differential? z)
         (for-all? (differential->terms z)
                   (lambda (term)
                     (known-real? (differential-coefficient term)))))
        (else
         (there-exists? *known-reals*
                        (lambda (w)
                          (or (equal? w z)
                              (let ((diff
                                     (ignore-errors
                                      (lambda ()
                                        (simplify (g:- w z))))))
                                (and (not (condition? diff))
                                     (exact-zero? diff)))))))))

;;; Permanent declaration

(define (declare-known-reals . stuff)
  (set! *known-reals* (list-union stuff *known-reals*)))

(define (declare-unknown-reals . stuff)
  (set! *known-reals* (list-difference *known-reals* stuff)))


;;; Temporary declaration

(define (with-known-reals stuff thunk)
  (fluid-let ((*known-reals*
               (list-union stuff *known-reals*)))
    (thunk)))

Add Gauss quadrature rules, Gauss-Kronrod implementation

I have my eyes on https://github.com/JuliaApproximation/FastGaussQuadrature.jl as a nice target to clone here. This would be an awesome addition to our evenly-spaced-points methods we have now!

sicmutils.function/memoize should work on operators, other HOF

Currently this code just works on functions. If you memoize something else, it will get wrapped in a function.

Make memoize generic, and install it for:

operators
structures, matrices, anything that implements using fmap
in Clojurescript, install separately for MetaFn and Fn etc

Then do some perf testing to see if this slows things down! Maybe we have a function-specific memoize and then a dispatching one too, so if you KNOW you are using a function you can specialize.

Replace numsymb implementations with pattern matching versions

For example, derivative becomes:

(r/choice
 (r/rule (D ?f) r/=> ((expt D 2) ?f))
 (r/rule ((expt D ?n) ?f) r/=> ((expt D (? #(inc (% '?n)))) ?f))
 (r/rule ?f r/=> (D ?f)))

This is faster and more readable!

Add multiple-argument gcd, lcm, patterned on `g/*` and friends

This is how the operation is defined in Scheme:

(define (g:gcd:n args)
  (cond ((null? args) :zero)
        ((null? (cdr args)) (car args))
        (else
         (let lp
             ((as (cddr args))
              (ans (g:gcd:bin (car args) (cadr args))))
           (cond ((null? as) ans)
                 ((g:one? ans) ans)
                 (else
                  (lp (cdr as) (g:gcd:bin ans (car as)))))))))

The task here would be to

rename the generics now to something like g/-gcd and g/-lcm
expose a multi-arg g/gcd, g/lcm that wraps them and does this trick with both functions.

interval arithmetic / ball arithmetic implementation

I saw a reference to "Arb: Efficient Arbitrary-Precision Midpoint-Radius Interval Arithmetic" on hackernews; it would be awesome to get proper interval arithmetic in to sicmutils.

More references:

blog from the author
https://www.texmacs.org/joris/ball/ball.html is a monstrously good resource with tons of background.
existing Clojure impl of interval arithmetic https://github.com/arrdem/imprecise

Extend generics to more fast / generic matrix protocols like core.matrix, gl-matrix

It would be really nice if it were possible to build our matrix operations on top of Clojure's core.matrix protocols.

This may be too much! But a big advantage here is that we could plug in to very fast matrix math for our state evolution.

The other possible matrix target to support is https://github.com/uncomplicate/neanderthal. https://github.com/whilo/denisovan exists, if someone wants to use neanderthal as a core.matrix implementation... and of course they can use neanderthal directly here. What's important is that the matrix representation be extensible.

This symbolic manipulation library claims core.matrix support: https://github.com/clojure-numerics/expresso and might be worth looking at.

Clojurescript, JS

I mentioned goog.math.Matrix in sicmutils/placeholder#82 , but I'm not sure this is terribly fast. gl-matrix is very quick, on the other hand. It would be nice to let users do the FAST stuff they want to, and then interact with other sicmutils data types via the generics. If only as a template for how folks can fuse into the ecosystem.

implement sd/IPerturbed for RationalFunction instances

This is hard because we'll need to detect differentials in the input and respond there, if we want to do this. It all works fine IF you restrict everything to the numerator...

sd/IPerturbed
  (perturbed? [_]
    (sd/perturbed? u))

  (replace-tag [this old new]
    (RationalFunction. arity
                       (sd/replace-tag u old new)
                       v
                       m))

  (extract-tangent [this tag]
    (RationalFunction. arity
                       (sd/extract-tangent u tag)
                       v
                       m))

But this is not very genuine.

Here's a test that shows how this can get wacky:

(testing "IPerturbed"
    (letfn [(f [x]
              (rf/make
               (p/make [1 2 (g/square x) 3])
               (p/make [5 (g/square x)])))
            (g [x y]
              (let [n (g/+ 1 (g/* 2 y)
                           (g/* (g/square x) (g/square y))
                           (g/* 3 (g/cube y)))
                    d (g/+ 5 (g/* (g/square x) y))]
                (g// n d)))]
      (((D f) 'x) 'y)
      (((cd/partial 0) g) 'x 'y)))