Mathematical Symbolic Differentiation in Guile Scheme

;;; Differentiation in code ;;; ======================= ;;; ;;; The following code calculates what is taught in a first semester ;;; university math course on Differential Calculus along with a ;;; routine to simplify s-expression based mathematical expressions. ;;; ;;; This post is to show how little code it takes to compute the ;;; mathematics taught in the first semester of the typical university ;;; math curriculum. The interesting part is the differentiate ;;; function which is about 40 lines of code. ;;; ;;; To use this, you would load this file into your favorite scheme ;;; implementation and from the repl call derivative, which will ;;; calculate the derivative of the given function and simplify the ;;; result. ;;; ;;; Here is a sample session using guile: ;;; ;;; $ guile ;;; guile> (load "mathx.scm") ;;; guile> (simplify '(* (/ (* x x) (+ (+ 2 3) (* 4 3))) y)) ;;; (* (/ (^ x 2) 17) y) ;;; guile> (derivative (^ x 2) x) ;;; (* 2 x) ;;; guile> (differentiate '(^ x 2) 'x) ;;; (* (* 2 (^ x 1)) 1) ;;; guile> (simplify (differentiate '(+ (^ x 2) (^ y 2)) 'y)) ;;; (* 2 y) ;;; guile> (simplify (differentiate '(* (^ x 2) (^ y 2)) 'y)) ;;; (* (^ x 2) (* 2 y)) ;;; guile> (simplify (differentiate '(sin (+ (^ x 2) (^ y 2))) 'x)) ;;; (* (cos (+ (^ x 2) (^ y 2))) (* 2 x)) ;;; guile> (simplify '(/ x (+ (+ 2 3) (* 4 3)))) ;;; (/ x 17) ;;; guile> (simplify '(* (/ (* x x) (+ (+ 2 3) (* 4 3))) y)) ;;; (* (/ (^ x 2) 17) y) ;;; guile> (derivative (^ %e (cos (^ x 2))) x) ;;; (* (^ %e (cos (^ x 2))) (- (* (sin (^ x 2)) (* 2 x)))) ;;; ;;; No guarantees this is bug free or suitable for any purpose. ;;; ;;; Author: Burton Samograd <burtonjohnsamograd@protonmail.com> ;;; License: AGPL ;;; Copyleft 2012 (define *constants* '(%e %pi)) (define (constant? x) (member x *constants*)) (define (flatten x) (cond ((null? x) '()) ((not (pair? x)) (list x)) (else (append (flatten (car x)) (flatten (cdr x)))))) (define (error . rest) (display "error: ") (display rest) (newline)) (define (smaller l1 l2) ;; return the smaller list, with smaller being the one with the less ;; nodes when flattened, or l1 if there is a tie (define (smaller-helper fl1 fl2) (if (not (or (null? fl1) (null? fl2))) (smaller-helper (cdr fl1) (cdr fl2)) (if (null? fl1) l1 l2))) (smaller-helper (flatten l1) (flatten l2))) (define (simplify-1-arg-expr expr) (case (car expr) ((-) (if (number? (cadr expr)) (* -1 (cadr expr)) `(- ,(simplify (cadr expr))))) ((ln) expr) (else expr))) (define (simplify-commutative? expr l r) ;; TODO: I have a feeling that this logic operation could be cleared up ;; if expr matches (+ (+ ..) (+ ..)), (+ x (+ ..)) or (+ (+ ..) x) (or (and (list? l) (list? r) (eq? (car expr) (car l)) (eq? (car expr) (car r))) (and (list? r) (not (list? l)) (eq? (car r) (car expr))) (and (list? l) (not (list? r)) (eq? (car l) (car expr))))) (define (simplify-commutative expr l r) ;; try to see if a commutative reordering will result in a smaller epxression (let ((op (car expr))) (cond ((not (list? r)) (smaller expr (simplify `(,op (,op ,r ,(cadr l)) ,(caddr l))))) ((not (list? l)) (smaller expr (simplify `(,op (,op ,l ,(cadr r)) ,(caddr r))))) (else ;; both l and r are lists ;; this is slow, but it checks all possible reductions and gives the best ;;(+ (+ x 1) (+ x 2)) ;;(+ (+ ll lr) (+ rl rr)) ;;rl lr ll rr ;;ll rr rl lr ;;rr lr rl ll ;;ll rl lr rr (let* ((rl (cadr r)) (rr (caddr r)) (ll (cadr l)) (lr (caddr l))) (smaller expr (smaller (smaller `(,op ,(simplify `(,op ,rl ,lr)) ,(simplify `(,op ,ll ,rr))) `(,op ,(simplify `(,op ,ll ,rr)) ,(simplify `(,op ,rl ,lr)))) (smaller `(,op ,(simplify `(,op ,rr ,lr)) ,(simplify `(,op ,rl ,ll))) `(,op ,(simplify `(,op ,ll ,rl)) ,(simplify `(,op ,lr ,rr)))))))) ))) (define (simplify-* expr l r) (cond ((and (number? l) (number? r)) (* l r)) ((eq? l 1) r) ((eq? r 1) l) ((or (eq? l 0) (eq? r 0)) 0) ((equal? l r) `(^ x 2)) ((simplify-commutative? expr l r) (simplify-commutative expr l r)) (else expr))) (define (simplify-/ expr l r) (cond ((and (number? l) (number? r)) (/ l r)) ((eq? r 1) l) ((eq? l 0) 0) ((eq? r 0) (error "Divide by zero in expr: " expr)) ((equal? l r) 1) ;; sin x / cos x = tan x ((and (list? l) (list? r) (eq? (car l) 'sin) (eq? (car r) 'cos) (equal? (cadr l) (cadr r))) (simplify `(tan ,(cadr l)))) ;; n / sin x = n * csc x ((and (list? r) (eq? (car r) 'sin)) (simplify `(* ,l (csc ,(cadr r))))) ;; n / cos x = n * sec x ((and (list? r) (eq? (car r) 'cos)) (simplify `(* ,l (sec ,(cadr r))))) (else expr))) (define (simplify-+ expr l r) (cond ((and (number? l) (number? r)) (+ l r)) ((eq? l 0) r) ((eq? r 0) l) ((equal? l r) `(* 2 ,l)) ((simplify-commutative? expr l r) (simplify-commutative expr l r)) (else expr))) (define (simplify-- expr l r) (cond ((eq? l 0) `(- ,r)) ((eq? r 0) l) ((equal? l r) 0) ((and (number? l) (number? r)) (- l r)) (else expr))) (define (simplify-^ expr l r) (cond ((eq? l 0) 0) ((eq? l 1) 1) ((eq? r 0) 1) ((eq? r 1) l) ((and (number? l) (number? r)) (exp (* (log l) r))) ;; maybe bad (else expr))) (define (simplify-2-arg-expr expr) (let ((l (simplify (cadr expr))) (r (simplify (caddr expr)))) ;; (display l) (display " ") (display r) (newline) (set! expr (list (car expr) l r)) (case (car expr) ((*) (simplify-* expr l r)) ((/) (simplify-/ expr l r)) ((+) (simplify-+ expr l r)) ((-) (simplify-- expr l r)) ((^) (simplify-^ expr l r))))) (define (simplify expr) ; (display "reducing: ") (display expr) (newline) (if (list? expr) (let ((narg (length (cdr expr)))) (if (= narg 1) (simplify-1-arg-expr expr) (simplify-2-arg-expr expr))) expr)) (define (f-of? expr x) ;; is expr a function of x? look for presence of x in flattened expr (cond ((null? expr) #f) ((not (pair? expr)) (eq? expr x)) (else (or (f-of? (car expr) x) (f-of? (cdr expr) x))))) (define (csc x) (/ 1 (sin x))) (define (sec x) (/ 1 (cos x))) (define (cot x) (/ 1 (tan x))) (define (ln x) (log x)) (define (differentiate expr v) ;(display "-> ") (display expr) (newline) (if (and (list? expr) (f-of? expr v)) (let ((op (car expr)) (narg (- (length expr) 1))) (case narg ((1) ; single argument functions (let ((x (cadr expr))) (case op ((-) `(- ,(differentiate x v))) ((ln) `(/ ,(differentiate x v) ,x)) ((sin) `(* (cos ,x) ,(differentiate x v))) ((cos) `(- (* (sin ,x) ,(differentiate x v)))) ((tan) `(* (^ (sec ,x) 2) ,(differentiate x v))) ((cot) `(- (* (^ (csc ,x) 2) ,(differentiate x v)))) ((sec) `(* (* (sec ,x) (tan ,x)) ,(differentiate x v))) ((csc) `(- (* (* (csc ,x) (cot ,x)) ,(differentiate x v)))) ((arcsin asin) `(/ ,(differentiate x v) (^ (- 1 (^ ,x 2)) (/ 1 2) ))) ((arccos acos) `(- ,(differentiate `(asin ,x) v))) ((arctan atan) `(/ ,(differentiate x v)) (- 1 (^ ,x 2))) ((arccot acot) `(- ,(differentiate `(atan ,x) v))) ((arcsec asec) `(/ ,(differentiate x v) (* (abs ,x) (^ (- (^ ,x 2 ) 1) 1/2)))) ((arccsc acsc) `(- ,(differentiate `(asec ,x) v)))))) ((2) ; two argument functions (let ((l (simplify (cadr expr))) (r (simplify (caddr expr)))) (case op ((+ -) `(,(car expr) ,(differentiate l v) ,(differentiate r v))) ((*) `(+ (* ,l ,(differentiate r v)) (* ,r ,(differentiate l v)))) ((/) `(/ (- (* ,l ,(differentiate r v)) (* ,r ,(differentiate l v))) (* ,r ,r))) ((^) (case l ((%e) `(* ,expr ,(differentiate r v))) (else `(* (* ,r (^ ,l ,(simplify `(- ,r 1)))) ,(differentiate l v))))) ))))) (if (and (symbol? expr) (eq? expr v)) 1 0))) (define-macro (derivative f d) `(simplify (differentiate ',f ',d)))