LISP AS RENAISSANCE WORKSHOP A Lispy Tour through Mathematical Physics , Mentat Collective Sam Ritchie 2

SICMUTILS 4

[(+ (sin 'x) 'x) (+ (sin 12) 2)] [(+ (sin x) x) 1.4634270819995652] (tex$$ (up (square (cos (* 't 'phi))) (floor (* 4 'zeta)))) ⎛ ⎝ cos 2 (ϕ t) ⌊4 ζ⌋ ⎞ ⎠ 5

(let [f (literal-function 'f (-> (UP Real Real) Real))] (tex$$ ((D f) (up 'alpha_1 'alpha_2)))) ⎡ ⎣ ∂0 f ∂1 f ⎛ ⎝ ⎛ ⎝ α 1 α 2 ⎞ ⎠ ⎞ ⎠ ⎛ ⎝ ⎛ ⎝ α 1 α 2 ⎞ ⎠ ⎞ ⎠ ⎤ ⎦ 6

❤️ OPEN SOURCE ❤️ https://github.com/sicmutils/sicmutils

https://github.com/sritchie/programming-2022 7

AGENDA Why build another CAS? Code as Communication? Literate Systems (and their problems) Why it matters What to do? + SICMUtils 9

SCMUTILS BY GJS 11

"programs must be written for people to read, and only incidentally for machines to execute." ~ Hal Abelson, Structure and Interpretation of Computer Programs 12

SICM AND FDG 14

MOTIVATION, GOOGLE X 16

DISAPPOINTMENT @ X! 17

19TH CENTURY SCIENTIFIC COMMUNICATION 19

HISTORY OF VECTOR ANALYSIS 20

21

CODE AS COMMUNICATION 23

NUMERICAL CODE AS COMMUNICATION? "We personally like Brent's algorithm for univariate minimization, as found on pages 79-80 of his book 'Algorithms for Minimization Without Derivatives'. It is pretty reliable and pretty fast, but we cannot explain how it works." ~ scmutils, refman.txt 24

BRENT'S BOOK, 1973 25

FORTRAN: NUMERICAL RECIPES, 1986 26

C++: BOOST, 2006 template 1 std::pair brent_find_minima(F f, T min, T max, int bits 2 noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval > 6 T tolerance = static_cast(ldexp(1.0, 1-bits)); 7 T x; // minima so far 8 T w; // second best point 9 T v; // previous value of w 10 T u; // most recent evaluation point 11 T delta; // The distance moved in the last step 12 T delta2; // The distance moved in the step before last 13 T fu, fv, fw, fx; // function evaluations at u, v, w, x 14 T mid; // midpoint of min and max 15 T fract1 fract2; // minimal relative movement in x 16 27

C++: BOOST, 2006 template 1 std::pair brent_find_minima(F f, T min, T max, int bits 2 noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval > 6 T tolerance = static_cast(ldexp(1.0, 1-bits)); 7 T x; // minima so far 8 T w; // second best point 9 T v; // previous value of w 10 T u; // most recent evaluation point 11 T delta; // The distance moved in the last step 12 T delta2; // The distance moved in the step before last 13 T fu, fv, fw, fx; // function evaluations at u, v, w, x 14 T mid; // midpoint of min and max 15 T fract1 fract2; // minimal relative movement in x 16 } 62 else 63 { 64 // golden section: 65 delta2 = (x >= mid) ? min - x : max - x; 66 delta = golden * delta2; 67 } 68 // update current position: 69 u = (fabs(delta) >= fract1) ? T(x + delta) : (delta > 0 70 fu = f(u); 71 if(fu <= fx) 72 { 73 // good new point is an improvement! 74 // update brackets: 75 if(u >= x) 76 i 77 27

C++: BOOST, 2006 template 1 std::pair brent_find_minima(F f, T min, T max, int bits 2 noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval > 6 T tolerance = static_cast(ldexp(1.0, 1-bits)); 7 T x; // minima so far 8 T w; // second best point 9 T v; // previous value of w 10 T u; // most recent evaluation point 11 T delta; // The distance moved in the last step 12 T delta2; // The distance moved in the step before last 13 T fu, fv, fw, fx; // function evaluations at u, v, w, x 14 T mid; // midpoint of min and max 15 T fract1 fract2; // minimal relative movement in x 16 template 1 std::pair brent_find_minima(F f, T min, T max, int bits 2 noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval > 6 T tolerance = static_cast(ldexp(1.0, 1-bits)); 7 T x; // minima so far 8 T w; // second best point 9 T v; // previous value of w 10 T u; // most recent evaluation point 11 T delta; // The distance moved in the last step 12 T delta2; // The distance moved in the step before last 13 T fu, fv, fw, fx; // function evaluations at u, v, w, x 14 T mid; // midpoint of min and max 15 T fract1 fract2; // minimal relative movement in x 16 ; x = u; 83 fv = fw; 84 fw = fx; 85 fx = fu; 86 } 87 else 88 { 89 // Oh dear, point u is worse than what we have alrea 90 // even so it *must* be better than one of our endpo 91 if(u < x) 92 min = u; 93 else 94 max = u; 95 if((fu <= fw) || (w == x)) 96 { 97 // however it is at least second best: 98 27

C++: BOOST, 2006 template 1 std::pair brent_find_minima(F f, T min, T max, int bits 2 noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval > 6 T tolerance = static_cast(ldexp(1.0, 1-bits)); 7 T x; // minima so far 8 T w; // second best point 9 T v; // previous value of w 10 T u; // most recent evaluation point 11 T delta; // The distance moved in the last step 12 T delta2; // The distance moved in the step before last 13 T fu, fv, fw, fx; // function evaluations at u, v, w, x 14 T mid; // midpoint of min and max 15 T fract1 fract2; // minimal relative movement in x 16 template 1 std::pair brent_find_minima(F f, T min, T max, int bits 2 noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval > 6 T tolerance = static_cast(ldexp(1.0, 1-bits)); 7 T x; // minima so far 8 T w; // second best point 9 T v; // previous value of w 10 T u; // most recent evaluation point 11 T delta; // The distance moved in the last step 12 T delta2; // The distance moved in the step before last 13 T fu, fv, fw, fx; // function evaluations at u, v, w, x 14 T mid; // midpoint of min and max 15 T fract1 fract2; // minimal relative movement in x 16 template 1 std::pair brent_find_minima(F f, T min, T max, int bits 2 noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval > 6 T tolerance = static_cast(ldexp(1.0, 1-bits)); 7 T x; // minima so far 8 T w; // second best point 9 T v; // previous value of w 10 T u; // most recent evaluation point 11 T delta; // The distance moved in the last step 12 T delta2; // The distance moved in the step before last 13 T fu, fv, fw, fx; // function evaluations at u, v, w, x 14 T mid; // midpoint of min and max 15 T fract1 fract2; // minimal relative movement in x 16 w = u; 100 fv = fw; 101 fw = fu; 102 } 103 else if((fu <= fv) || (v == x) || (v == w)) 104 { 105 // third best: 106 v = u; 107 fv = fu; 108 } 109 } 110 111 }while(--count); 112 113 max iter -= count; 114 27

C++: BOOST, 2006 template 1 std::pair brent_find_minima(F f, T min, T max, int bits 2 noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval > 6 T tolerance = static_cast(ldexp(1.0, 1-bits)); 7 T x; // minima so far 8 T w; // second best point 9 T v; // previous value of w 10 T u; // most recent evaluation point 11 T delta; // The distance moved in the last step 12 T delta2; // The distance moved in the step before last 13 T fu, fv, fw, fx; // function evaluations at u, v, w, x 14 T mid; // midpoint of min and max 15 T fract1 fract2; // minimal relative movement in x 16 template 1 std::pair brent_find_minima(F f, T min, T max, int bits 2 noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval > 6 T tolerance = static_cast(ldexp(1.0, 1-bits)); 7 T x; // minima so far 8 T w; // second best point 9 T v; // previous value of w 10 T u; // most recent evaluation point 11 T delta; // The distance moved in the last step 12 T delta2; // The distance moved in the step before last 13 T fu, fv, fw, fx; // function evaluations at u, v, w, x 14 T mid; // midpoint of min and max 15 T fract1 fract2; // minimal relative movement in x 16 template 1 std::pair brent_find_minima(F f, T min, T max, int bits 2 noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval > 6 T tolerance = static_cast(ldexp(1.0, 1-bits)); 7 T x; // minima so far 8 T w; // second best point 9 T v; // previous value of w 10 T u; // most recent evaluation point 11 T delta; // The distance moved in the last step 12 T delta2; // The distance moved in the step before last 13 T fu, fv, fw, fx; // function evaluations at u, v, w, x 14 T mid; // midpoint of min and max 15 T fract1 fract2; // minimal relative movement in x 16 template 1 std::pair brent_find_minima(F f, T min, T max, int bits 2 noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval > 6 T tolerance = static_cast(ldexp(1.0, 1-bits)); 7 T x; // minima so far 8 T w; // second best point 9 T v; // previous value of w 10 T u; // most recent evaluation point 11 T delta; // The distance moved in the last step 12 T delta2; // The distance moved in the step before last 13 T fu, fv, fw, fx; // function evaluations at u, v, w, x 14 T mid; // midpoint of min and max 15 T fract1 fract2; // minimal relative movement in x 16 fw = fu; 102 } 103 else if((fu <= fv) || (v == x) || (v == w)) 104 { 105 // third best: 106 v = u; 107 fv = fu; 108 } 109 } 110 111 }while(--count); 112 113 max_iter -= count; 114 115 return std::make_pair(x, fx); 116 } 117 27

PYTHON: SCIPY, 2001 def optimize(self): 1 # set up for optimization 2 func = self.func 3 xa, xb, xc, fa, fb, fc, funcalls = self.get_bracket_info() 4 _mintol = self._mintol 5 _cg = self._cg 6 ################################# 7 #BEGIN CORE ALGORITHM 8 ################################# 9 x = w = v = xb 10 fw = fv = fx = func(*((x,) + self.args)) 11 if (xa < xc): 12 a = xa 13 b = xc 14 else: 15 a = xc 16 28

SCHEME: 1987 ;;; Brent's algorithm for univariate minimization -- transcrib 1 ;;; pages 79-80 of his book "Algorithms for Minimization Witho 2 3 (define (brent-min f a b eps) 4 (let ((a (min a b)) (b (max a b)) 5 (maxcount 100) 6 (small-bugger-factor *sqrt-machine-epsilon*) 7 (g (/ (- 3 (sqrt 5)) 2)) 8 (d 0) (e 0) (old-e 0) (p 0) (q 0) (u 0) (fu 0)) 9 (let* ((x (+ a (* g (- b a)))) 10 (fx (f x)) 11 (w x) (fw fx) (v x) (fv fx)) 12 (let loop ((count 0)) 13 (if (> count maxcount) 14 (list 'maxcount x fx count) ;failed to converge 15 (let* ((tol (+ (* eps (abs x)) small-bugger-factor)) 16 29

ACTUAL CORE IDEA 30

SO WHAT? What's wrong with this? How can we do better? 31

EXISTING IDEAS Literate Programming LOGO's Microworlds Notebooks, Mathematica 33

34

LOGO'S MICROWORLDS 35

36

NOTEBOOKS, MATHEMATICA 37

WHY AREN'T THESE WORKING? Literate Programming == one-way Science is Multiplayer Dynamic is Too Seductive Real Work doesn't happen here 39

WHY AREN'T THESE WORKING? Literate Programming == one-way Science is Multiplayer Dynamic is Too Seductive Real Work doesn't happen here 39

WHY AREN'T THESE WORKING? Literate Programming == one-way Science is Multiplayer Dynamic is Too Seductive Real Work doesn't happen here 39

WHY AREN'T THESE WORKING? Literate Programming == one-way Science is Multiplayer Dynamic is Too Seductive Real Work doesn't happen here 39

WHY AREN'T THESE WORKING? Literate Programming == one-way Science is Multiplayer Dynamic is Too Seductive Real Work doesn't happen here 39

WHY CARE? 41

"FUTURE OF EDUCATION" Mythic Understanding Romantic Understanding Philosophic Understanding Ironic Understanding 42

WHAT TO DO? 44

SICM AND FDG AS CLUES 45

46

EXPLANATION AS SIDE EFFECT 47

EULER-LAGRANGE EQUATIONS S[q](t a , t b ) = D ∫ t b t a L(t, q(t), Dq(t))dt D ∫ t b t a L(t, q(t), Dq(t))dt = 0 49

EULER-LAGRANGE EQUATIONS d dt ∂L ∂ ˙ q − ∂L ∂q = 0 50

EULER-LAGRANGE EQUATIONS "What could this expression possibly mean?" d dt ∂L ∂ ˙ q − ∂L ∂q = 0 50

EXPAND: d dt ∂L ∂ ˙ q − ∂L ∂q = 0 51

EXPAND: d dt ∂L ∂ ˙ q − ∂L ∂q = 0 d dt ∂L(t, q, ˙ q) ∂ ˙ q − ∂L(t, q, ˙ q) ∂q = 0 ⎛ ⎝ q=w(t) ˙ q= dw(t) dt ⎞ ⎠ q=w(t) ˙ q= dw(t) dt 51

OKAY, FINE d dt ∂L ∂ ˙ q − ∂L ∂q = 0 d dt ((∂ 2 L) (t, w(t), d dt w(t))) −(∂ 1 L) (t, w(t), d dt w(t)) = 0 52

SUBSTITUTIONS (Df)(t) = d dx f(x) x=t Γ[w](t) = (t, w(t), d dt w(t)) 53

d dt ((∂ 2 L) (Γ[w](t))) − (∂ 1 L) (Γ[w](t)) = 0 54

d dt ((∂ 2 L) (Γ[w](t))) − (∂ 1 L) (Γ[w](t)) = 0 D ((∂ 2 L) ∘ (Γ[w])) − (∂ 1 L) ∘ (Γ[w]) = 0 54

d dt ((∂ 2 L) (Γ[w](t))) − (∂ 1 L) (Γ[w](t)) = 0 D ((∂ 2 L) ∘ (Γ[w])) − (∂ 1 L) ∘ (Γ[w]) = 0 (defn Lagrange-equations [L] 1 (fn [w] 2 (- (D (comp ((partial 2) L) (Gamma w))) 3 (comp ((partial 1) L) (Gamma w))))) 4 54

d dt ((∂ 2 L) (Γ[w](t))) − (∂ 1 L) (Γ[w](t)) = 0 D ((∂ 2 L) ∘ (Γ[w])) − (∂ 1 L) ∘ (Γ[w]) = 0 (defn Lagrange-equations [L] 1 (fn [w] 2 (- (D (comp ((partial 2) L) (Gamma w))) 3 (comp ((partial 1) L) (Gamma w))))) 4 (fn [w] (- (D (comp ((partial 2) L) (Gamma w))) (comp ((partial 1) L) (Gamma w))))) (defn Lagrange-equations [L] 1 2 3 4 54

(let [L (- (comp (literal-function 'T) velocity) (comp (literal-function 'V) coordinate)) w (literal-function 'w)] (tex$$ (((Lagrange-equations L) w) 't))) D 2 w (t) D 2 T (Dw (t)) + DV (w (t)) 55

SIMPLE HARMONIC MOTION (- (* 1/2 m (square v)) (defn L-harmonic 1 "Returns a Lagrangian of a simple harmonic oscillator (mass-s 2 3 m is the mass and k is the spring constant used in Hooke's la 4 [m k] 5 (fn [[_ q v]] 6 7 (* 1/2 k (square q))))) 8 (defn proposed-solution [t] (* 'a (cos (+ (* 'omega t) 'phi)))) (let [L (L-harmonic 'm 'k) w proposed-solution] (tex$$ (((Lagrange-equations L) w) 't))) −a m ω 2 cos (ω t + ϕ) + a k cos (ω t + ϕ) 56

SIMPLE HARMONIC MOTION (- (* 1/2 m (square v)) (defn L-harmonic 1 "Returns a Lagrangian of a simple harmonic oscillator (mass-s 2 3 m is the mass and k is the spring constant used in Hooke's la 4 [m k] 5 (fn [[_ q v]] 6 7 (* 1/2 k (square q))))) 8 (* 1/2 k (square q))))) (defn L-harmonic 1 "Returns a Lagrangian of a simple harmonic oscillator (mass-s 2 3 m is the mass and k is the spring constant used in Hooke's la 4 [m k] 5 (fn [[_ q v]] 6 (- (* 1/2 m (square v)) 7 8 (defn proposed-solution [t] (* 'a (cos (+ (* 'omega t) 'phi)))) (let [L (L-harmonic 'm 'k) w proposed-solution] (tex$$ (((Lagrange-equations L) w) 't))) −a m ω 2 cos (ω t + ϕ) + a k cos (ω t + ϕ) 56

HOOKE'S LAW (let [L (L-harmonic 'm 'k) w (literal-function 'w)] (tex$$ (((Lagrange-equations L) w) 't))) k w (t) + m D 2 w (t) 57

CLERK DEMO Let's convince ourselves that that is true, by doing the spring demo in the browser. 59

COMMUNITY Sussman et. al. Colin Smith Clerk: Martin Kavalar, Jack Rusher, Nextjournal team SCI: Michiel Borkent (@borkdude) Mathbox: Chris Chudzicki, Steven Wittens 61

HOW TO GET INVOLVED? WHAT'S NEXT? SICM and FDG: Executable Textbooks Full library as essays Collaborative editing, simulation Please Steal! 63

THANKS! Sam Ritchie, Mentat Collective Slides, Demos live at @sritchie https://github.com/sritchie/programming-2022 65

No content