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