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)) 180