Partial solution for bug 1452341: solve(x^(5/2)+1,x)
Subject: Partial solution for bug 1452341: solve(x^(5/2)+1,x)
From: Richard Fateman
Date: Thu, 25 May 2006 07:39:10 -0700
I'm pretty sure that you will not be able to confirm all the roots
that are produced by solve, by using ratsimp.
(Sorry).
I think that Mathematica has a message that says essentially that
"there may be roots that are not in this list" and (perhaps optionally,
if asked to check the solutions) "not all the solutions have been
confirmed".
There are heuristics like numerical evaluation, giving random values
to any extra parameters. But these are iffy. Doing exact substitution
in finite fields, when possible, removes floating-point roundoff error;
there are papers on hash-coding of symbolic expressions...
RJF
----- Original Message -----
From: "Raymond Toy" <raymond.toy at ericsson.com>
To: <maxima at math.utexas.edu>
Sent: Thursday, May 25, 2006 7:22 AM
Subject: Partial solution for bug 1452341: solve(x^(5/2)+1,x)
>
> In this bug, solve(x^(5/2)+1,x) produces some roots that aren't
> actually roots. It does this by first solving y^5+1 = 0, and solving
> x^(1/2) = y for each of the roots found previously and returns each of
> these solutions. But some of these are not actually roots.
>
> The code that implements this is in solvespec and solvespec1 in
> solve.lisp. I've modified this code to verify each of the roots is
> actually a root of the original equation, and discards any that are
> not.
>
> However, I'm not sure if the substitution of the root will always
> produce 0. I just run $ratsimp to try to simplify the expression if
> necessary. If the result is not zero, that root is discarded, even if
> it's because ratsimp can't simplify the reesult to zero.
>
> What do people think about this? Some roots might be missed, and you
> won't even know about it.
>
> Ray
>
> P.S.
>
> Here is the replacement function if anyone wants to test this out:
>
> (defun solvespec (exp $%emode)
> (prog (a b c)
> (setq a (pdis (caddr exp)))
> (setq c (pdis (list (car exp) 1 1)))
> (cond ((null (cdddr exp))
> (return (solve c *var (times (cadr exp) mult)))))
> (setq b (pdis (pminus (cadddr (cdr exp)))))
> ;; First, save *roots, in case there's something there already.
> ;; Then rebind it to nil, so we know exactly the roots we've
> ;; found for this equation. (Is it really necessary to do this?
> ;; Anyway, this should be very safe.)
> (let ((save-roots *roots))
> (let ((*roots nil))
> ;; Solve a*f(x)^n+b.
> (solvespec1 c
> (simpnrt (div* b a) (cadr exp))
> (make-rat 1 (cadr exp))
> (cadr exp))
> ;; Now look at all the roots we found and see if they really
> ;; are roots.
> (do ((r *roots (cddr r)))
> ((null r))
> ;; Substitute the root that we found into the equation and
> ;; see if it evaluates to zero. If not, it's not a root.
> ;; Is calling ratsimp good enough? We assume it is here.
> (let ((val ($ratsimp
> (sub (mul a (power (subst (third (car r)) *var c)
> (cadr exp)))
> b))))
> (when (alike1 val 0)
> ;; Yay, it's a root. Push it onto save-roots. FIXME:
> ;; But what if it's not a root? Should we push that
> ;; onto *fail? What if it really is a root and ratsimp
> ;; can't simplify it to zero? Then what?
> (push (second r) save-roots)
> (push (first r) save-roots)))))
> (setf *roots save-roots))))
>
> _______________________________________________
> Maxima mailing list
> Maxima at math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima
>