src/psolve.lisp

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