The specificity of Lisp is, maybe, best visible on problems that require processing of data naturally represented as S-expressions. An example of such problem is "Probability that random propositional formula is tautology." In this post, the variation of classical example - symbolic differentiation - is presented: newLISP (and it is very similar in other Lisps) program that computes the tangent of the graph of the function f of single variable, defined as composition of elementary functions (+, -, sin, cos ...) in given point x0. The tangent of the function f in x0 is defined as linear function

y(x) = a·(x - x0) + b,

where a = f '(x0) and b = f(x0). The program consists of:

  • The function tangent that accepts two arguments: the function f in form of lambda list and the floating point number x0. The function tangent analyzes the code of the function f, calculates the values a and b according to mathematical definition of tangent (see above) and constructs the code of the function

    (lambda(x)(+. (*. a (-. x x0)) b))

    where symbols a, b and x0 are replaced with their respective values; that function is then returned as result of the function tangent. The implementation is very short, but it delegates main part of the work to the function d.
  • The function d accepts two arguments: formula and variable. It returns the expression obtained by symbolic differentiation of the formula for given variable. It is main function in this program, typical recursive "code walking" through formula. Simple "domain specific language" is used for descriptions of the rules  for symbolic differentiation that should be applied, for example

    (+. (*. df g) (*. f dg))

    for multiplication. The application of  rules for symbolic differentiation frequently give result that can be simplified. For instance, derivation of (*. 2 x) is

    (+. (*. 0 x) (*. 2 1))

    so function simplify is called when appropriate.
  • The function simplify accepts only one argument: symbolic expression formula. It analyzes and simplify formula, again, using typical recursive "code walking". For instance, S-expression (- x x) is simplified to 0. There are infinitely many possible rules for simplification; the function performs only few of these. The function uses eval at one place: if some expression contains only operators and constants, but not variables then the easisest way to simplify it is to eval it. (One might speculate that simplify is generalized eval.) Simplification is, actually, not necessary for computation of the tangent. However, its use in the context of symbolic differentiation is reasonable.
Few operators from my library are used; function expand// i.e. parallel expand, fexpr println= for convenient printing and floating point arithmetic operators +., -., *. and /., synonymous for built-in add, sub, mul and div. These functions are not essential.

The graph of the complicated function and tangent on the curve suggests that program, generally, works.

(setf f0 (lambda(x)
           (+. (sin (*. 12 x))
               (cos (*. 32 x))
               (tan (*. x 1.4))
               (asin x)
               (acos  x)
               (atan x)
               (*. x (cos (/. 7 x)))
               (sqrt (/. 9 x))
               (pow x x)
               (*. x (sinh x))
               (*. x (cosh x))
               (asinh x)
               (sin (acosh (+. x 1)))
               (atanh x))))
(setf x0 0.4)
(println= (tangent f0 x0))
; (tangent f0 x0)=(lambda (x) (+. (*. -21.491 (-. x 0.4)) 10.252))


Finally, the code:
(setf [print.supressed] true [println.supressed] true)
(load "http://www.instprog.com/Instprog.default-library.lsp")
(setf [print.supressed] nil [println.supressed] nil)

(define (tangent f x0)
   (letn((variable           (first (first f)))
         (expression         (last f))
         (derived-expression (d expression variable))
         (a                  (eval (expand derived-expression
                                           (list (list variable x0)))))
         (b                  (f x0)))
       (expand// '(lambda(x)(+. (*. a (-. x x0)) b))
                  'a 'b 'x0)))

(define (d formula variable)
      ((= formula variable) 1)
      ((atom? formula) 0)
      ((list? formula)
       (letn((operator (first formula))
             (operands (rest formula))
                 (letn((flatexpr (flat expr))
                       (f  (if (find 'f  flatexpr)(operands 0)))
                       (df (if (find 'df flatexpr)(d f variable)))
                       (g  (if (find 'g  flatexpr)(operands 1)))
                       (dg (if (find 'dg flatexpr)(d g variable))))
                   (expand// expr 'f 'df 'g 'dg)))))

         (case operator 

           (+. (cons '+. (map (lambda(op)(d op variable)) operands)))
           (-. (cons '-. (map (lambda(op)(d op variable)) operands)))

           (*. (case (length operands)
                 (1    (lexpand 'df))
                 (2    (lexpand '(+. (*. df g) (*. f dg))))
                 (true (d (list '*. (first operands)
                                    (cons '*. (rest operands)))

           (/. (case (length operands)
                 (1    (d (list '/. 1 (first operands)) variable)) 
                 (2    (lexpand '(/. (-. (*. df g) (*. f dg)) (*. g g))))
                 (true (d (list '/.  (first operands)
                                     (cons '*. (rest operands)))

           (pow (d (lexpand '(exp (*. g (log f)))) variable))
           (exp (lexpand '(*. f df)))
           (log (if (= (length operands) 1)
                    (lexpand '(/. df f))
                    (d (lexpand '(/. (log f) (log g))) variable)))

           (sqrt  (lexpand '(*. 0.5 df (/. 1 (sqrt f)))))
           (sin   (lexpand '(*. (cos f) df)))
           (cos   (lexpand '(*. (-. (sin f)) df)))
           (tan   (lexpand '(/. df (pow (cos f) 2))))
            (asin (lexpand '(/. df (sqrt (-. 1 (*. f f))))))
           (acos  (lexpand '(-. (/. df (sqrt (-. 1 (*. f f)))))))
           (atan  (lexpand '(/. df (+. 1 (*. f f)))))
           (sinh  (lexpand '(*. (cosh f) df)))
           (cosh  (lexpand '(*. (sinh f) df)))
           (tanh  (lexpand '(*. (-. 1 (pow (tanh f) 2)) df)))
           (asinh (lexpand '(/. df (sqrt (+.(*. f f) 1)))))
           (acosh (lexpand '(/. df (sqrt (-. (*. f f) 1)))))
           (atanh (lexpand '(/. df (-. 1 (*. f f)))))


(define (simplify formula)
    ((atom? formula) formula)
    ((list? formula)
     (letn((operator (first formula))
           (operands (map simplify (rest formula)))
           (formula (cons operator operands)))

         ; if all operands are constants, then 
         ; simplified formula is evaluated formula
         ((for-all number? operands)(eval formula))
         ; (*. x), (+. x) => x
         ((and (or (= operator '*.) (= operator '+.)) 
               (= (length operands) 1))
           (first operands))
         ; (*. ... 0 ...) => 0  
         ((and (= operator '*.) (find 0 operands)) 0)
         ; (*. ... 1 ...) => (*. ...)
         ((and (= operator '*.) (find 1 operands))
          (simplify (clean (curry = 1) formula)))
         ; (+. ... 0 ...) => 0
         ((and (= operator '+.) (find 0 operands))
          (simplify (clean zero? formula)))
         ; (-. (-. ...)) => ...
         ((match '(-. (-. ?)) formula)
           (last (last formula)))
         ; (-. minuend ...)
         ((and (= operator '-.) (> (length operands) 1))
          (letn((minuend (first operands))
                (subtrahends (rest operands))
                (subtrahend  (simplify (cons '+. subtrahends))))
            (cond ((zero? minuend)        (simplify (list '-. 
                  ((zero? subtrahend)     minuend)
                  ((= minuend subtrahend) 0)
                  (true                   (list '-. minuend 
         ; (/. (/. ...))
         ((match '(/. (/. ?)) formula) (last (last formula)))
         ; (/. dividend ...)
         ((and (= operator '/.) (> (length operands) 1))
          (letn((dividend (first operands))
                (divisors (rest operands))
                (divisor  (simplify (cons '*. divisors))))
            (cond ((zero? dividend)     0)
                  ((= divisor 1)        dividend)
                  ((= divisor -1)       (simplify (list '-. dividend)))
                  ((= dividend divisor) 1)
                  (true                 (list '/. dividend divisor)))))
         (true formula)))))) 

No comments:

Post a Comment