Rename lapack.texi to lapack.texi.m4

[maxima.git] / src / csimp.lisp

;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;     The data in this file contains enhancements.                   ;;;;;
;;;                                                                    ;;;;;
;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
;;;     All rights reserved                                            ;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;     (c) Copyright 1980 Massachusetts Institute of Technology         ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(in-package :maxima)

(macsyma-module csimp)

(declare-top (special $exponentialize
                      var genvar $ratprint
                      errorsw))

(load-macsyma-macros rzmac)

(loop for (a b) on
       '(%sin %asin %cos %acos %tan %atan
         %cot %acot %sec %asec %csc %acsc
         %sinh %asinh %cosh %acosh %tanh %atanh
         %coth %acoth %sech %asech %csch %acsch)
       by #'cddr
       do  (putprop a b '$inverse) (putprop b a '$inverse))

(defmfun $demoivre (exp)
  (let ($exponentialize nexp)
    (cond ((atom exp) exp)
          ((and (eq (caar exp) 'mexpt) (eq (cadr exp) '$%e)
                (setq nexp (demoivre (caddr exp))))
           nexp)
          (t (recur-apply #'$demoivre exp)))))

(defun demoivre (l)
  (cond ($exponentialize
         (merror (intl:gettext "demoivre: 'demoivre' and 'exponentialize' cannot both be true.")))
        (t (setq l (islinear l '$%i))
           (and l (not (equal (car l) 0))
                (m* (m^ '$%e (cdr l))
                    (m+ (list '(%cos) (car l))
                        (m* '$%i (list '(%sin) (car l)))))))))

;; If expr is of the form a*var1+b where a is freeof var1
;; then (a . b) is returned else nil.
(defun islinear (expr var1)
  (let ((a (let ((*islinp* t))
             (sdiff expr var1))))
    (if (freeof var1 a)
        (cons a (maxima-substitute 0 var1 expr)))))

(defmfun $partition (e var1)
  (let ((e (mratcheck e))
        (var1 (getopr var1)))
    (cond
      (($listp e)
       (do ((l (cdr e) (cdr l)) (l1) (l2) (x))
           ((null l) (list '(mlist simp)
                           (cons '(mlist simp) (nreverse l1))
                           (cons '(mlist simp) (nreverse l2))))
         (setq x (mratcheck (car l)))
         (cond ((free x var1) (setq l1 (cons x l1)))
               (t (setq l2 (cons x l2))))))

      ((mplusp e)
       (destructuring-bind (constant . variable) (partition e var1 0)
         (list '(mlist simp) constant variable)))

      ((mtimesp e)
       (destructuring-bind (constant . variable) (partition e var1 1)
         (list '(mlist simp) constant variable)))

      (t
       (merror
        (intl:gettext
         "partition: first argument must be a list or '+' or '*' expression; found ~M") e)))))

;; Apply PREDICATE to each element in the sequence SEQ and return
;;
;;   (VALUES YES NO),
;;
;; where YES and NO are lists consisting of elements for which PREDICATE is true
;; or false, respectively.
(defun partition-by (predicate seq)
  (let ((yes) (no))
    (map nil
         (lambda (x)
           (if (funcall predicate x)
               (push x yes)
               (push x no)))
         seq)
    (values yes no)))

;; Partition an expression, EXP, into terms that contain VAR1 and terms that
;; don't contain VAR1. EXP is either considered as a product or a sum (for which
;; you should pass K = 1 or 0, respectively).

(defun partition (exp var1 k)
  (let ((op (if (= k 0) 'mplus 'mtimes))
        (exp-op (unless (atom exp) (caar exp))))
    (cond
      ;; Exit immediately if exp = var1
      ((alike1 exp var1) (cons k exp))

      ;; If we have the wrong operation then the expression is either entirely
      ;; constant or entirely variable.
      ((not (eq exp-op op))
       (if (freeof var1 exp)
           (cons exp k)
           (cons k exp)))

      ;; Otherwise, we want to partition the arguments into constant and
      ;; variable.
      (t
       (multiple-value-bind (constant variable)
           (partition-by (lambda (x) (freeof var1 x)) (cdr exp))
         (cons (cond ((null constant) k)
                     ((null (cdr constant)) (car constant))
                     (t (simplifya
                         (cons (list op) (nreverse constant)) t)))
               (cond ((null variable) k)
                     ((null (cdr variable)) (car variable))
                     (t (simplifya
                         (cons (list op) (nreverse variable)) t)))))))))

;;To use this INTEGERINFO and *ASK* need to be special.
;;(defun integerpw (x)
;; ((lambda (*ask*)
;;    (integerp10 (ssimplifya (sublis '((z** . 0) (*z* . 0)) x))))
;;  t))

;;(defun integerp10 (x)
;; ((lambda (d)
;;   (cond ((or (null x) (not (free x '$%i))) nil)
;;       ((mnump x) (integerp x))
;;       ((setq d (assolike x integerinfo)) (eq d 'yes))
;;       (*ask* (setq d (cond ((integerp x) 'yes) (t (needinfo x))))
;;              (setq integerinfo (cons (list x d) integerinfo))
;;              (eq d 'yes))))
;; nil))

(setq var (make-symbol "foo"))

(defun numden (e)
  (prog (varlist)
     (setq varlist (list var))
     (newvar (setq e (fmt e)))
     (setq e (cdr (ratrep* e)))
     (setq dn*
           (simplifya (pdis (ratdenominator e))
                      nil))
     (setq nn*
           (simplifya (pdis (ratnumerator e))
                      nil)))
  (values nn* dn*))

;; Like NUMDEN but dependency on VAR is explicit.  Use this instead of
;; NUMDEN if possible.  Also, we do not set the special variables DN*
;; and NN*.  The callers are responsible for doing that, now.
(defun numden-var (e var1)
  (let ((varlist (list var1))
        dn nn)
     (newvar (setq e (fmt e)))
     (setq e (cdr (ratrep* e)))
     (setq dn
           (simplifya (pdis (ratdenominator e))
                      nil))
     (setq nn
           (simplifya (pdis (ratnumerator e))
                      nil))
     (values nn dn)))

(defun fmt (exp)
  (let (nn*)
    (cond ((atom exp) exp)
          ((mnump exp) exp)
          ((eq (caar exp) 'mexpt)
           (cond ((and (mnump (caddr exp))
                       (eq ($sign (caddr exp)) '$neg))
                  (list '(mquotient)
                        1
                        (cond ((equal (caddr exp) -1)
                               (fmt (cadr exp)))
                              (t (list (list (caar exp))
                                       (fmt (cadr exp))
                                       (timesk -1 (caddr exp)))))))
                 ((atom (caddr exp))
                  (list (list (caar exp))
                        (fmt (cadr exp))
                        (caddr exp)))
                 ((and (mtimesp (setq nn* (sratsimp (caddr exp))))
                       (mnump (cadr nn*))
                       (equal ($sign (cadr nn*)) '$neg))
                  (list '(mquotient)
                        1
                        (list (list (caar exp))
                              (fmt (cadr exp))
                              (cond ((equal (cadr nn*) -1)
                                     (cons '(mtimes)
                                           (cddr nn*)))
                                    (t (neg nn*))))))
                 ((eq (caar nn*) 'mplus)
                  (fmt (spexp (cdr nn*) (cadr exp))))
                 (t (cons (ncons (caar exp))
                          (mapcar #'fmt (cdr exp))))))
          (t (cons (delsimp (car exp)) (mapcar #'fmt (cdr exp)))))))

(defun spexp (expl dn*)
  (cons '(mtimes) (mapcar #'(lambda (e) (list '(mexpt) dn* e)) expl)))

(defun subin (y x)
  (cond ((not (among var x)) x)
        (t (maxima-substitute y var x))))

;; Like SUBIN but dependency on VAR is explicit.  Use this instead
;; when possible.
(defun subin-var (y x ivar)
  (cond ((not (among ivar x)) x)
        (t (maxima-substitute y ivar x))))

;; Right-hand side (rhs) and left-hand side (lhs) of binary infix expressions.
;; These are unambiguous for relational operators, some other built-in infix operators,
;; and user-defined infix operators (declared by the infix function).

;; a - b and a / b are somewhat problematic, since subtraction and division are not
;; ordinarily represented as such (rather a - b = a + (-1)*b and a / b = a * b^(-1)).
;; Also, - can be unary. So let's not worry about - and / .

;; Other problematic cases: The symbols $< $<= $= $# $>= $> have a LED property,
;; but these symbols never appear in expressions returned by the Maxima parser;
;; MLESSP, MLEQP, MEQUAL etc are substituted. So ignore those symbols here.
(let
  ((relational-ops
     ;; < <= = # equal notequal >= >
     '(mlessp mleqp mequal mnotequal $equal $notequal mgeqp mgreaterp
       %mlessp %mleqp %mequal %mnotequal %equal %notequal %mgeqp %mgreaterp))

   (other-infix-ops
     ;; := ::= : :: ->
     '(mdefine mdefmacro msetq mset marrow
       %mdefine %mdefmacro %msetq %mset %marrow)))

  (defmfun $rhs (rel)
     (if (atom rel)
       0
       (if (or (member (caar rel) (append relational-ops other-infix-ops) :test #'eq)
               ;; This test catches user-defined infix operators.
               (eq (get (caar rel) 'led) 'parse-infix))
         (caddr rel)
         0)))

  (defmfun $lhs (rel)
     (if (atom rel)
       rel
       (if (or (member (caar rel) (append relational-ops other-infix-ops) :test #'eq)
               ;; This test catches user-defined infix operators.
               (eq (get (caar rel) 'led) 'parse-infix))
         (cadr rel)
         rel))))

(defun ratgreaterp (x y)
  (cond ((and (ratnump x) (ratnump y))
         (great x y))
        ((equal ($asksign (m- x y)) '$pos))))

;; Simplify the exponential function of the type exp(p/q*%i*%pi+x) using the
;; periodicity of the exponential function and special values for rational
;; numbers with a denominator q = 2, 3, 4, or 6. e is the argument of the 
;; exponential function. For float and bigfloat numbers in the argument e only
;; simplify for an integer representation or a half integral value.
;; The result is an exponential function with a simplified argument.
(defun %especial (e)
  (prog (varlist y k kk j ans $%emode $ratprint genvar)
     (let (($keepfloat nil) ($float nil))
       (unless (setq y (pip ($ratcoef e '$%i))) (return nil))
       ;; Subtract the term y*%i*%pi from the expression e.
       (setq k ($expand (add e (mul -1 '$%pi '$%i y)) 1))
       ;; This is a workaround to get the type (integer, float, or bigfloat)
       ;; of the expression. kk must evaluate to 1, 1.0, or 1.0b0.
       ;; Furthermore, if e is nonlinear, kk does not simplify to a number ONE.
       ;; Because of this we do not simplify something like exp((2+x)^2*%i*%pi)
       (setq kk (div (sub ($expand e) k) (mul '$%i '$%pi y)))
       ;; Return if kk is not an integer or kk is ONE, but y not an integer
       ;; or a half integral value.
       (if (not (or (integerp kk)
                    (and (onep1 kk)
                         (integerp (add y y)))))
           (return nil))
       (setq j (trigred y))
       (setq ans (spang1 j t)))
     (cond ((among '%sin ans)
            (cond ((equal y j) (return nil))
                  ((zerop1 k)
                   ;; To preverse the type we add k into the result.
                   (return (power '$%e (mul '$%pi '$%i (add k j)))))
                  (t 
                    ;; To preserve the type we multiply kk into the result.
                    (return 
                      (power '$%e (add (mul kk k) (mul kk '$%pi '$%i j))))))))
     (setq y (spang1 j nil))
     ;; To preserve the type we multiply kk into the result.
     (return (mul (power '$%e (mul kk k)) (add y (mul '$%i ans))))))

(defun trigred (r)
  (prog (m n eo flag)
     (cond ((numberp r) (return (cond ((even r) 0) (t 1)))))
     (setq m (cadr r))
     (cond ((minusp m) (setq m (- m)) (setq flag t)))
     (setq n (caddr r))
     loop (cond ((> m n)
                 (setq m (- m n))
                 (setq eo (not eo))
                 (go loop)))
     (setq m (list '(rat)
                   (cond (flag (- m)) (t m))
                   n))
     (return (cond (eo (addk m (cond (flag 1) (t -1))))
                   (t m)))))

(defun polyinx (exp x ind)
  (prog (genvar varlist var $ratfac)
     (setq var x)
     (cond ((numberp exp)(return t))
           ((polyp exp)
            (cond (ind (go on))
                  (t (return t))))
           (t (return nil)))
     on (setq genvar nil)
     (setq varlist (list x))
     (newvar exp)
     (setq exp (cdr (ratrep* exp)))
     (cond
       ((or (numberp (cdr exp))
            (not (eq (car (last genvar)) (cadr exp))))
        (setq x (pdis (cdr exp)))
        (return (cond ((eq ind 'leadcoef)
                       (div* (pdis (caddr (car exp))) x))
                      (t (setq exp (car exp))
                         (div* (cond ((atom exp) exp)
                                     (t
                                      (pdis (list (car exp)
                                                  (cadr exp)
                                                  (caddr exp)))))
                               x))
                      ))))))

(defun polyp (a)
  (cond ((atom a) t)
        ((member (caar a) '(mplus mtimes) :test #'eq)
         (every #'polyp (cdr a)))
        ((eq (caar a) 'mexpt)
         (cond ((free (cadr a) var)
                (free (caddr a) var))
               (t (and (integerp (caddr a))
                       (> (caddr a) 0)
                       (polyp (cadr a))))))
        (t (andmapcar #'(lambda (subexp)
                          (free subexp var))
                      (cdr a)))))

;; Like polyp but takes an extra arg for the variable.
(defun polyp-var (a var2)
  (cond ((atom a) t)
        ((member (caar a) '(mplus mtimes) :test #'eq)
         (every #'(lambda (p)
                    (polyp-var p var2))
                (cdr a)))
        ((eq (caar a) 'mexpt)
         (cond ((free (cadr a) var2)
                (free (caddr a) var2))
               (t (and (integerp (caddr a))
                       (> (caddr a) 0)
                       (polyp-var (cadr a) var2)))))
        (t (andmapcar #'(lambda (subexp)
                          (free subexp var2))
                      (cdr a)))))

(defun pip (e)
  (prog (varlist d c)
     (newvar e)
     (cond ((not (member '$%pi varlist :test #'eq)) (return nil)))
     (setq varlist '($%pi))
     (newvar e)
     (let (($ratfac nil))
       ;; non-nil $ratfac changes form of CRE
       (setq e (cdr (ratrep* e))))
     (setq d (cdr e))
     (cond ((not (atom d)) (return nil))
           ((equal e '(0 . 1))
            (setq c 0)
            (go loop)))
     (setq c (ptterm (cdar e) 1))
     loop (cond ((atom c)
                 (cond ((equal c 0) (return nil))
                       ((equal 1 d) (return c))
                       (t (return (list '(rat) c d))))))
     (setq c (ptterm (cdr c) 0))
     (go loop)))

(defun spang1 (j ind)
  (prog (ang ep $exponentialize $float $keepfloat)
     (cond ((floatp j) (setq j (maxima-rationalize j))
            (setq j (list '(rat simp) (car j) (cdr j)))))
     (setq ang j)
     (cond
       (ind nil)
       ((numberp j)
        (cond ((zerop j) (return 1)) (t (return -1))))
       (t (setq j
                (trigred (add2* '((rat simp) 1 2)
                                (list (car j)
                                      (- (cadr j))
                                      (caddr j)))))))
     (cond ((numberp j) (return 0))
           ((mnump j) (setq j (cdr j))))
     (return
       (cond ((equal j '(1 2)) 1)
             ((equal j '(-1 2)) -1)
             ((or (equal j '(1 3))
                  (equal j '(2 3)))
              (div ($sqrt 3) 2))
             ((or (equal j '(-1 3))
                  (equal j '(-2 3)))
             (div ($sqrt 3) -2))
             ((or (equal j '(1 6))
                  (equal j '(5 6)))
              '((rat simp) 1 2))
             ((or (equal j '(-1 6))
                  (equal j '(-5 6)))
              '((rat simp) -1 2))
             ((or (equal j '(1 4))
                  (equal j '(3 4)))
              (div 1 ($sqrt 2)))
             ((or (equal j '(-1 4))
                  (equal j '(-3 4)))
             (div -1 ($sqrt 2)))
             (t (cond ((mnegp ang)
                       (setq ang (timesk -1 ang) ep t)))
                (setq ang (list '(mtimes simp)
                                ang
                                '$%pi))
                (cond (ind (cond (ep (list '(mtimes simp)
                                           -1
                                           (list '(%sin simp) ang)))
                                 (t (list '(%sin simp) ang))))
                      (t (list '(%cos simp) ang))))))))

(defun archk (a b v)
  (simplify
   (cond ((and (equal a 1) (equal b 1)) v)
         ((and (equal b -1) (equal 1 a))
          (list '(mtimes) -1 v))
         ((equal 1 b)
          (list '(mplus) '$%pi (list '(mtimes) -1 v)))
         (t (list '(mplus) v (list '(mtimes) -1 '$%pi))))))

(defun genfind (h v)
;;; finds gensym coresponding to v h
  (do ((varl (caddr h) (cdr varl))
       (genl (cadddr h) (cdr genl)))
;;;is car of rat form
      ((eq (car varl) v) (car genl))))