250

# Lisp as Renaissance Workshop: A Lispy Tour through Mathematical Physics

These slides are from my keynote at the 2022 European Lisp Symposium in Porto, Portugal: https://european-lisp-symposium.org/2022/index.html

The source code for these slides and the demos in the talk live at https://github.com/sritchie/programming-2022.

Abstract:
Lisp is an exquisite medium for the communication of computational ideas. From our most accurate observations of physical reality up through chemistry, biology, and cognition, the universe seems to be computing itself; modeling and simulating these systems in machines has led to incredible technological wealth.

Deep principles and beautiful abstractions seem to drive these systems, but they have always been hard to discover; and we are floundering at the computational frontiers of intelligence, synthetic biology and control systems for our climate. The only way to push forward is to build powerful tools that can communicate and teach.

This talk will take a tour through SICMUtils, a Lisp system designed as a workshop for conducting serious work in mathematical physics and sharing those explorations in a deeply interactive, multiplayer way. The library’s growth parallels our human scientific history; hopefully tools like this will help us write the next chapter. March 23, 2022

## Transcript

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

2. SICMUTILS
4

3. [(+ (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

4. (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

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

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

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

8. SCMUTILS BY GJS
11

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

12

10. SICM AND FDG

14

16

12. DISAPPOINTMENT @ X!
17

13. 19TH CENTURY SCIENTIFIC
COMMUNICATION
19

14. HISTORY OF VECTOR ANALYSIS
20

15. 21

16. CODE AS COMMUNICATION
23

17. 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

18. BRENT'S BOOK, 1973
25

19. FORTRAN: NUMERICAL RECIPES, 1986
26

20. 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::declval3
{
4
BOOST_MATH_STD_USING
5
bits = (std::min)(policies::digits >
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

21. 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::declval3
{
4
BOOST_MATH_STD_USING
5
bits = (std::min)(policies::digits >
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

22. 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::declval3
{
4
BOOST_MATH_STD_USING
5
bits = (std::min)(policies::digits >
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::declval3
{
4
BOOST_MATH_STD_USING
5
bits = (std::min)(policies::digits >
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

23. 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::declval3
{
4
BOOST_MATH_STD_USING
5
bits = (std::min)(policies::digits >
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::declval3
{
4
BOOST_MATH_STD_USING
5
bits = (std::min)(policies::digits >
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::declval3
{
4
BOOST_MATH_STD_USING
5
bits = (std::min)(policies::digits >
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

24. 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::declval3
{
4
BOOST_MATH_STD_USING
5
bits = (std::min)(policies::digits >
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::declval3
{
4
BOOST_MATH_STD_USING
5
bits = (std::min)(policies::digits >
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::declval3
{
4
BOOST_MATH_STD_USING
5
bits = (std::min)(policies::digits >
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::declval3
{
4
BOOST_MATH_STD_USING
5
bits = (std::min)(policies::digits >
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

25. 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

26. 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

27. ACTUAL CORE IDEA
30

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

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

30. 34

31. LOGO'S MICROWORLDS
35

32. 36

33. NOTEBOOKS, MATHEMATICA

37

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

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

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

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

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

39. WHY CARE?
41

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

41. WHAT TO DO?
44

42. SICM AND FDG AS CLUES

45

43. 46

44. EXPLANATION AS SIDE EFFECT
47

45. 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

46. EULER-LAGRANGE EQUATIONS
d
dt
∂L
∂ ˙
q

∂L
∂q
= 0
50

47. EULER-LAGRANGE EQUATIONS
"What could this expression possibly mean?"
d
dt
∂L
∂ ˙
q

∂L
∂q
= 0
50

48. EXPAND:
d
dt
∂L
∂ ˙
q

∂L
∂q
= 0
51

49. 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

50. 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

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

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

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

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

55. 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

56. (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

57. 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

58. 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

59. 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

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

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

62. HOW TO GET INVOLVED? WHAT'S NEXT?
SICM and FDG: Executable Textbooks
Full library as essays
Collaborative editing, simulation