src/psolve.lisp

   1 ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
   2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   3 ;;;     The data in this file contains enhancments.                    ;;;;;
   4 ;;;                                                                    ;;;;;
   5 ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
   6 ;;;     All rights reserved                                            ;;;;;
   7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   8 ;; ** (c) Copyright 1980 Massachusetts Institute of Technology **
   9 ;; whole file revised to avoid conflict with CRE forms. 4/27/2016 Richard Fateman
  10
  11 (in-package :maxima)
  12
  13 (macsyma-module psolve)
  14
  15 (declare-top (special mult *roots *failures $solvefactors))
  16 (declare-top (special expsumsplit $dispflag checkfactors *g
  17                       $algebraic equations ;List of E-labels
  18                       *power *varb *flg $derivsubst
  19                       $%emode genvar genpairs varlist broken-not-freeof
  20                       mult    ;Some crock which tracks multiplicities.
  21                       *roots ;alternating list of solutions and multiplicities
  22                       *failures ;alternating list of equations and multiplicities
  23                       *myvar $listconstvars
  24                       *has*var *var $dontfactor
  25                       $keepfloat $ratfac
  26                       xm* xn* mul*))
  27
  28 (defmvar flag4 nil)
  29
  30 (defun solvecubic (x)
  31   (prog (s1 a0 a1 a2 discr lcoef adiv3 omega^2 pdiv3 qdiv-2
  32          omega y1 u y2)
  33
  34      (setq x (cdr x))
  35      (setq lcoef (cadr x))
  36      (setq adiv3
  37          (mul
  38                  '((rat) -1 3)
  39                  (rdis (setq a2 (ratreduce (ptterm x 2)
  40                                            lcoef)))))
  41      (setq a1 (ratreduce (ptterm x 1) lcoef))
  42      (setq a0 (ratreduce (ptterm x 0) lcoef))
  43
  44     ;;      coefficients now a0,a1,a2,  and leading coef 1.
  45
  46    (setf a2 (rdis a2) a1 (rdis a1) a0 (rdis a0))
  47
  48
  49
  50    (setq s1 (mul' ((rat) 1 2) '$%i (power 3 '((rat) 1 2))))
  51    (setq omega (add '((rat) -1 2)  s1)
  52          omega^2 (add '((rat) -1 2) (mul -1 s1)))
  53      (setq pdiv3
  54          (add (mul a1 '( (rat) 1  3))
  55               (mul  (power      a2 2) '((rat) -1 9))))
  56      (and (not (equal pdiv3 0)) (go harder))
  57      (setq y1
  58          (mul
  59              '((rat) 1 3)
  60           (add
  61            (simpnrt  (setq y2 (add (power a2 3)
  62                                    (mul -27 a0)))
  63                      3)                 ; cube root
  64            (mul -1 a2 ))))
  65
  66      (and flag4 (return (solve3 y1 mult)))
  67    (setq y2 (simpnrt  (mul  y2 '((rat) 1  27)) 3))
  68      (return (mapc #'(lambda (j) (solve3 j mult))
  69                    (list y1
  70                         (add (mul omega y2) adiv3)
  71                         (add (mul omega^2 y2) adiv3))))
  72      harder
  73      (setq qdiv-2
  74          (add (mul (add (mul a1 a2) (mul -3  a0))
  75                    '((rat) 1 6))
  76               (mul (power a2 3) '((rat) -1  27) )))
  77
  78      (cond ((equal qdiv-2 0)
  79             (setq u (simpnrt pdiv3 2))
  80             (setq y1 adiv3))
  81          (t (setq discr (add
  82                          (power pdiv3 3)
  83                          (power qdiv-2 2)))
  84
  85               (cond ((equal discr 0)
  86                      (setq u (simpnrt qdiv-2 3)))
  87                     (t (setq discr (simpnrt discr 2))
  88                        (and (complicated discr)
  89                             (setq discr (adispline discr)))
  90                      (setq u (power
  91                               (add
  92                                                      qdiv-2
  93                                                      discr)
  94                               '((rat) 1 3)))
  95                        (and (complicated u)
  96                             (setq u (adispline u)))))))
  97      (if (equal u 0) (merror (intl:gettext "SOLVECUBIC: arithmetic overflow.")))
  98      (or y1
  99        (setq y1 (add adiv3  u  (mul  -1 pdiv3 (power u -1)))))
 100      (return
 101        (cond (flag4 (solve3 y1 mult))
 102              (t (mapc
 103                  #'(lambda (j) (solve3 j mult))
 104                  (list y1
 105                     (add  adiv3 (mul omega u)  (mul -1 pdiv3 omega^2 (power u -1)))
 106                     (add  adiv3 (mul omega^2 u)(mul -1 pdiv3 omega   (power u -1))))))))))
 107
 108 (defun solvequartic (x)
 109   (prog (a0 a1 a2 b1 b2 b3 b0 lcoef z1 r tr1 tr2 d d1 e sqb3)
 110      (setq x (cdr x) lcoef (cadr x))
 111    (setq b3 (rdis(ratreduce (ptterm x 3) lcoef)))
 112    (setq b2 (rdis(ratreduce (ptterm x 2) lcoef)))
 113    (setq b1 (rdis(ratreduce (ptterm x 1) lcoef)))
 114    (setq b0 (rdis(ratreduce (ptterm x 0) lcoef)))
 115    (setq a2 (mul -1 b2))
 116    (setq a1 (sub (mul b1 b3) (setq a0 (mul b0 4))))
 117      (setq a0
 118          (sub (sub (mul b2 a0) (mul (setq sqb3(power b3 2)) b0 )) (power b1 2)))
 119    (setq tr2   (mul'((rat) 1 4)
 120                 (sub (sub (mul b3  b2 4)
 121                           (mul 8 b1))
 122                      (mul sqb3 b3 )) ))
 123      (setq z1 (resolvent a2 a1 a0))
 124      (setq r
 125          (add
 126                           z1
 127           (sub (mul sqb3 '((rat) 1 4))
 128                                         b2)))
 129      (setq r (simpnrt r 2))
 130      (and (equal r 0) (go l0))
 131      (and (complicated r) (setq r (adispline r)))
 132      (and (complicated tr2) (setq tr2 (adispline tr2)))
 133      (setq tr1
 134          (sub
 135           (sub (mul sqb3 '((rat) 1  2))
 136                b2)
 137            z1))
 138      (and (complicated tr1) (setq tr1 (adispline tr1)))
 139      (setq tr2 (div* tr2 r))
 140      (go lb1)
 141    l0   (setq d1 (simpnrt (add (power z1 2) (mul -4 b0)) 2))
 142    (setq tr2 (mul 2 d1))
 143      (and (complicated tr2) (setq tr2 (adispline tr2)))
 144    (setq tr1 (sub (mul sqb3 '((rat) 3 4))  (mul b2 2)))
 145      (and (complicated tr1) (setq tr1 (adispline tr1)))
 146   lb1
 147      (setq d (div (power (add tr1 tr2) '((rat simp) 1 2)) 2))
 148      (setq e (div (power (sub tr1 tr2) '((rat simp) 1 2)) 2))
 149      (and (complicated d) (setq d (adispline d)))
 150      (and (complicated e) (setq e (adispline e)))
 151    (setq a2 (mul b3 '((rat) -1 4)))
 152      (setq a1 (div* r 2))
 153
 154      (setq z1
 155      (list (add a2 a1 d)                ;1
 156            (add a2 a1 (mul -1 d))       ;2
 157            (add a2 (mul -1 a1) e)       ;3
 158            (add a2 (mul -1 a1) (mul -1 e)) ;4
 159                ))
 160      (return (mapc #'(lambda (j) (solve3 j mult)) z1))))
 161
 162 ;;; SOLVES RESOLVENT CUBIC EQUATION
 163 ;;; GENERATED FROM QUARTIC
 164
 165 (defun resolvent (a2 a1 a0)
 166   (prog (*roots flag4 *failures $solvefactors yy) ;undoes binding in
 167      (setq flag4 t $solvefactors t)            ;algsys
 168     (setf yy (gensym))
 169      (solve (add
 170              (power yy 3)
 171              (mul a2   (power yy 2))
 172              (mul a1  yy)
 173              a0)
 174             yy
 175             1)
 176      (when (member 0 *roots :test #'equal) (return 0))
 177      (return (caddar (cdr (reverse *roots))))))
 178
 179 (declare-top (unspecial mult))