Fix typo in display-html-help
[maxima.git] / src / psolve.lisp
blob3098af447b7706feac9f7b6fcda79c4b0715ce62
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
11 (in-package :maxima)
13 (macsyma-module psolve)
15 (declare-top (special mult *roots *failures))
16 (declare-top (special expsumsplit *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*))
27 (defmvar flag4 nil)
29 (defun solvecubic (x)
30 (prog (s1 a0 a1 a2 discr lcoef adiv3 omega^2 pdiv3 qdiv-2
31 omega y1 u y2)
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))
43 ;; coefficients now a0,a1,a2, and leading coef 1.
45 (setf a2 (rdis a2) a1 (rdis a1) a0 (rdis a0))
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 ))))
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) )))
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)))
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))))))))))
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
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))
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)))
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))
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
159 (return (mapc #'(lambda (j) (solve3 j mult)) z1))))
161 ;;; SOLVES RESOLVENT CUBIC EQUATION
162 ;;; GENERATED FROM QUARTIC
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)
175 (when (member 0 *roots :test #'equal) (return 0))
176 (return (caddar (cdr (reverse *roots))))))
178 (declare-top (unspecial mult))