For the new year, my number one wish is a high-quality polynomial
system solver for Maxima.
First, I'll describe a bug I found in algsys and give a possible
fix; second I'll list the steps I've made at getting a new Maxima
polynomial solver and ask for your advice.
Here is the bug:
------------------------------------------------------------------
GCL (GNU Common Lisp) Version(2.4.0) Sun Oct 14 17:17:45 CDT 2001
Licensed under GNU Library General Public License
Contains Enhancements by W. Schelter
Maxima 5.6 Sun Oct 14 17:17:07 CDT 2001 (with enhancements by W. Schelter).
Licensed under the GNU Public License (see file COPYING)
/* Let m = matrix([a,b],[c,d]). Then m^2=id implies the following equations. */
(C1) eqs : [b*c+a^2-1, b*d + a*b,c*d + a*c, d^2 + b*c - 1];
(D1) [b*c+a^2-1,b*d+a*b,c*d+a*c,d^2+b*c-1]
(C2) algsys(eqs,[a,b,c,d]);
Error: Caught fatal error [memory may be damaged]
Fast links are on: do (si::use-fast-links nil) for debugging
Error signalled by CATCH.
Broken at MACSYMA-TOP-LEVEL. Type :H for Help.
--------------------------------------------------------------------
The culprit may be the function punivarp (defined in algsys.lisp).
I'm guessing that punivarp(poly) is supposed to evaluate to true if
poly is a univariate polynomial in the (global) variable *ivar*. Further,
the argument poly is supposed to be in CRE form. (Or is it similar to
CRE form, but specialized to polynomials?) For the function call
in (c2), punivarp gets called with zero as its argument. It's clear
that punivarp will misbehave when called with zero as its argument:
(defun punivarp (poly)
(do ((l (cdr poly) (cddr l)))
((null l) t)
(or (numberp (cadr l))
(and (eq (caadr l) *ivar*)
(punivarp (cadr l)))
(return nil))))
Changing punivarp to
(defun punivarp (poly)
(or ($constantp poly)
(do ((l (cdr poly) (cddr l)))
((null l) t)
(or (numberp (cadr l))
(and (eq (caadr l) *ivar*)
(punivarp (cadr l)))
(return nil)))))
eliminates the bug:
(C5) load("/home/barton/mac2/maxima-5.6/src/algsys.lisp")$
Loading /home/barton/mac2/maxima-5.6/src/algsys.lisp
Finished loading /home/barton/mac2/maxima-5.6/src/algsys.lisp
(C6) algsys(eqs,[a,b,c,d]);
(D6) [[a = 1,b = %R1,c = 0,d = -1],[a = -1,b = %R2,c = 0,d = 1],
[a = -1,b = 0,c = 0,d = -1],[a = 1,b = 0,c = 0,d = 1]]
The bug is gone, but the solution set isn't complete; the equations
have a two parameter family of solutions. With algsys_uses_triangsys
set to false, Macsyma 422 returns the solution set identical to
(d6).
In an effort to get a better polynomial solver, I have:
1. Tried to use the Maxima grobner code; it wouldn't load.
2. I asked Gail Zacharias if she had grobner code that she could
donate to Maxima. She doesn't have a usable copy and doubts that
she has the legal right to donate it.
3. I asked Dongming Wang if he had Macsyma code for his elimination
method that he could donate. He told me that some time ago
he gave Macsyma Inc some advice on a Macsyma implementation but that
he didn't have any Macsyma code.
4. Dan Stanger converted a Maple to MuPad translator to a Maple to
Maxima translator. From the Maple share library, he translated
a polynomial solver into 7,000 lines of Maxima. I wrote a bit
of supporting Lisp code for polynomial degree, etc. (My code
is, I suspect, far from optimal; I don't know enough about CRE
expressions to do these things efficiently.)
Advice?
Barton