#4252 algsys only returns 2 (of 3) solutions of a cubic

[David Billinghurst]

algsys only finds two solutions of the cubic x^3+x^2-x+%pi-1. Solve returns three, as expected.

Found while investigating bug [#4247]

Maxima 5.47.0 https://maxima.sourceforge.io
using Lisp CLISP 2.49+ (2010-07-17)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) display2d:false$
(%i2) eq1:x^3+x^2-x+%pi-1$
(%i3) algsys([eq1],[x]);
(%o3) [[x = -((2^(2/3)*(3^(9/2)*sqrt(27*%pi^2-32*%pi)+729*%pi-432)^(1/3)
          -3*2^(2/3)*(3^(3/2)*sqrt(27*%pi^2-32*%pi)-27*%pi+16)^(1/3)+6)
          /18)],
       [x = -(((2^(2/3)*sqrt(3)*%i-2^(2/3))
          *(3^(9/2)*sqrt(27*%pi^2-32*%pi)+729*%pi-432)^(1/3)
          +(2^(2/3)*3^(3/2)*%i+3*2^(2/3))
           *(3^(3/2)*sqrt(27*%pi^2-32*%pi)-27*%pi+16)^(1/3)+12)
          /36)]]
(%i4) solve(eq1,x);
(%o4) [x = (-1/2-(sqrt(3)*%i)/2)*((sqrt(%pi)*sqrt(27*%pi-32))/(2*3^(3/2))
                                 +(-(3*(%pi-1))-1)/6+(-1)/27)
                                 ^(1/3)
         +(4*((sqrt(3)*%i)/2+(-1)/2))/(9
                                      *((sqrt(%pi)*sqrt(27*%pi-32))
                                       /(2*3^(3/2))
                                       +(-(3*(%pi-1))-1)/6+(-1)/27)
                                       ^(1/3))+(-1)/3,
       x = ((sqrt(3)*%i)/2+(-1)/2)*((sqrt(%pi)*sqrt(27*%pi-32))/(2*3^(3/2))
                                   +(-(3*(%pi-1))-1)/6+(-1)/27)
                                   ^(1/3)
         +(4*(-1/2-(sqrt(3)*%i)/2))/(9*((sqrt(%pi)*sqrt(27*%pi-32))
                                       /(2*3^(3/2))
                                       +(-(3*(%pi-1))-1)/6+(-1)/27)
                                       ^(1/3))+(-1)/3,
       x = ((sqrt(%pi)*sqrt(27*%pi-32))/(2*3^(3/2))
         +(-(3*(%pi-1))-1)/6+(-1)/27)
         ^(1/3)
         +4/(9*((sqrt(%pi)*sqrt(27*%pi-32))/(2*3^(3/2))
               +(-(3*(%pi-1))-1)/6+(-1)/27)
               ^(1/3))+(-1)/3]

[David Billinghurst]

Deep in the bowels of algsys, function callsolve returns 3 solutions. One of them gets lost on the way back.

Maxima 5.47.0 https://maxima.sourceforge.io
using Lisp CLISP 2.49+ (2010-07-17)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) display2d:false$
(%i2) eq1:x^3+x^2-x+%pi-1$
(%i3) ?trace(?callsolve)$
;; Tracing function CALLSOLVE.
(%i4) algsys([eq1],[x]);
1. Trace: (CALLSOLVE '((#:X16628 3 1 2 1 1 -1 0 (#:%PI16629 1 1 0 -1)) . #:X16628))
1. Trace: CALLSOLVE ==>
((((MEQUAL) $X
   ((MTIMES SIMP) ((RAT SIMP) -1 36)
    ((MPLUS SIMP) 12
     ((MTIMES SIMP)
      ((MPLUS SIMP) ((MTIMES SIMP) 3 ((MEXPT SIMP) 2 ((RAT SIMP) 2 3)))
       ((MTIMES SIMP) ((MEXPT SIMP) 2 ((RAT SIMP) 2 3))
        ((MEXPT SIMP) 3 ((RAT SIMP) 3 2)) $%I))
      ((MEXPT SIMP)
       ((MPLUS SIMP) 16 ((MTIMES SIMP RATSIMP) -27 $%PI)
        ((MTIMES SIMP) ((MEXPT SIMP) 3 ((RAT SIMP) 3 2))
         ((MEXPT SIMP)
          ((MPLUS SIMP) ((MTIMES SIMP RATSIMP) -32 $%PI)
           ((MTIMES SIMP) 27 ((MEXPT SIMP RATSIMP) $%PI 2)))
          ((RAT SIMP) 1 2))))
       ((RAT SIMP) 1 3)))
     ((MTIMES SIMP)
      ((MPLUS SIMP) ((MTIMES SIMP) -1 ((MEXPT SIMP) 2 ((RAT SIMP) 2 3)))
       ((MTIMES SIMP) ((MEXPT SIMP) 2 ((RAT SIMP) 2 3))
        ((MEXPT SIMP) 3 ((RAT SIMP) 1 2)) $%I))
      ((MEXPT SIMP)
       ((MPLUS SIMP) -432 ((MTIMES SIMP RATSIMP) 729 $%PI)
        ((MTIMES SIMP) ((MEXPT SIMP) 3 ((RAT SIMP) 9 2))
         ((MEXPT SIMP)
          ((MPLUS SIMP) ((MTIMES SIMP RATSIMP) -32 $%PI)
           ((MTIMES SIMP) 27 ((MEXPT SIMP RATSIMP) $%PI 2)))
          ((RAT SIMP) 1 2))))
       ((RAT SIMP) 1 3)))))))
 (((MEQUAL) $X
   ((MTIMES SIMP) ((RAT SIMP) -1 18)
    ((MPLUS SIMP) 6
     ((MTIMES SIMP) -3 ((MEXPT SIMP) 2 ((RAT SIMP) 2 3))
      ((MEXPT SIMP)
       ((MPLUS SIMP) 16 ((MTIMES SIMP RATSIMP) -27 $%PI)
        ((MTIMES SIMP) ((MEXPT SIMP) 3 ((RAT SIMP) 3 2))
         ((MEXPT SIMP)
          ((MPLUS SIMP) ((MTIMES SIMP RATSIMP) -32 $%PI)
           ((MTIMES SIMP) 27 ((MEXPT SIMP RATSIMP) $%PI 2)))
          ((RAT SIMP) 1 2))))
       ((RAT SIMP) 1 3)))
     ((MTIMES SIMP) ((MEXPT SIMP) 2 ((RAT SIMP) 2 3))
      ((MEXPT SIMP)
       ((MPLUS SIMP) -432 ((MTIMES SIMP RATSIMP) 729 $%PI)
        ((MTIMES SIMP) ((MEXPT SIMP) 3 ((RAT SIMP) 9 2))
         ((MEXPT SIMP)
          ((MPLUS SIMP) ((MTIMES SIMP RATSIMP) -32 $%PI)
           ((MTIMES SIMP) 27 ((MEXPT SIMP RATSIMP) $%PI 2)))
          ((RAT SIMP) 1 2))))
       ((RAT SIMP) 1 3)))))))
 (((MEQUAL) $X
   ((MTIMES SIMP) ((RAT SIMP) -1 18)
    ((MPLUS SIMP) 6
     ((MTIMES SIMP) -3 ((MEXPT SIMP) 2 ((RAT SIMP) 2 3))
      ((MEXPT SIMP)
       ((MPLUS SIMP) 16 ((MTIMES SIMP RATSIMP) -27 $%PI)
        ((MTIMES SIMP) ((MEXPT SIMP) 3 ((RAT SIMP) 3 2))
         ((MEXPT SIMP)
          ((MPLUS SIMP) ((MTIMES SIMP RATSIMP) -32 $%PI)
           ((MTIMES SIMP) 27 ((MEXPT SIMP RATSIMP) $%PI 2)))
          ((RAT SIMP) 1 2))))
       ((RAT SIMP) 1 3)))
     ((MTIMES SIMP) ((MEXPT SIMP) 2 ((RAT SIMP) 2 3))
      ((MEXPT SIMP)
       ((MPLUS SIMP) -432 ((MTIMES SIMP RATSIMP) 729 $%PI)
        ((MTIMES SIMP) ((MEXPT SIMP) 3 ((RAT SIMP) 9 2))
         ((MEXPT SIMP)
          ((MPLUS SIMP) ((MTIMES SIMP RATSIMP) -32 $%PI)
           ((MTIMES SIMP) 27 ((MEXPT SIMP RATSIMP) $%PI 2)))
          ((RAT SIMP) 1 2))))
       ((RAT SIMP) 1 3))))))))
(%o4) [[x = -((2^(2/3)*(3^(9/2)*sqrt(27*%pi^2-32*%pi)+729*%pi-432)^(1/3)
          -3*2^(2/3)*(3^(3/2)*sqrt(27*%pi^2-32*%pi)-27*%pi+16)^(1/3)+6)
          /18)],
       [x = -(((2^(2/3)*sqrt(3)*%i-2^(2/3))
          *(3^(9/2)*sqrt(27*%pi^2-32*%pi)+729*%pi-432)^(1/3)
          +(2^(2/3)*3^(3/2)*%i+3*2^(2/3))
           *(3^(3/2)*sqrt(27*%pi^2-32*%pi)-27*%pi+16)^(1/3)+12)
          /36)]]

[David Billinghurst]

OK. The issue is with the second solution found by solve.

(%i1) display2d:false$
(%i2) expr:((sqrt(3)*%i)/2+(-1)/2)*((sqrt(%pi)*sqrt(27*%pi-32))/(2*3^(3/2))
                        +(-(3*(%pi-1))-1)/6+(-1)/27)
                        ^(1/3)
 +(4*(-1/2-(sqrt(3)*%i)/2))/(9*((sqrt(%pi)*sqrt(27*%pi-32))/(2*3^(3/2))
                               +(-(3*(%pi-1))-1)/6+(-1)/27)
                               ^(1/3))+(-1)/3$
(%i3) rectform(float(expr));
(%o3) 0.8945133161601935*%i+0.5099711601368959

The lisp simplifying function within algsys simplify-after-subst converts this to the third solution found by solve. One of these two solutions is then removed as a duplicate solution by function CONDENSESOLNL, called from ALGSYS0.

(%i4) r2:?simplify\-after\-subst(expr);
(%o4) -((2^(2/3)*(3^(9/2)*sqrt(27*%pi^2-32*%pi)+729*%pi-432)^(1/3)
 -3*2^(2/3)*(3^(3/2)*sqrt(27*%pi^2-32*%pi)-27*%pi+16)^(1/3)+6)
 /18)
(%i5)  rectform(float(r2));
(%o5) -2.0199423202737896

The function simplify-after-subst was introduced to avoid exactly this sort of nonsense with complex roots of polynomials.

[David Billinghurst]

Picking apart the solution finds two subexpressions with the same problem.

Expressions p1 and p2 below have different floating point values after processing with lisp function simplify-after-subst

(%i2) display2d:false$
(%i3) s2:((sqrt(3)*%i)/2+(-1)/2)*((sqrt(%pi)*sqrt(27*%pi-32))/(2*3^(3/2))
                                 +(-(3*(%pi-1))-1)/6+(-1)/27)
                                 ^(1/3)
         +(4*(-1/2-(sqrt(3)*%i)/2))/(9*((sqrt(%pi)*sqrt(27*%pi-32))
                                       /(2*3^(3/2))
                                       +(-(3*(%pi-1))-1)/6+(-1)/27)
                                       ^(1/3))+(-1)/3$

/* Here is the problem. simplify\-after\-subst changes value */
(%i4) rectform(float(s2));
(%o4) 0.8945133161601935*%i+0.5099711601368959;
(%i5) t2:simplify\-after\-subst(s2);
(%o5) -((2^(2/3)*(3^(9/2)*sqrt(27*%pi^2-32*%pi)+729*%pi-432)^(1/3)
 -3*2^(2/3)*(3^(3/2)*sqrt(27*%pi^2-32*%pi)-27*%pi+16)^(1/3)+6)
 /18)
(%i6) rectform(float(t2));
(%o6) -2.0199423202737896

/* Pick s2 apart as sum p1 + p2 + p3 */
(%i7) part(s2,0);
(%o7) "+"
(%i8) p1:part(s2,1);
(%o8) ((sqrt(3)*%i)/2-1/2)*((sqrt(%pi)*sqrt(27*%pi-32))/(2*3^(3/2))
                           +(-(3*(%pi-1))-1)/6-1/27)
                           ^(1/3)
(%i9) p2:part(s2,2);
(%o9) (4*(-((sqrt(3)*%i)/2)-1/2))/(9*((sqrt(%pi)*sqrt(27*%pi-32))/(2*3^(3/2))
                                     +(-(3*(%pi-1))-1)/6-1/27)
                                     ^(1/3))
(%i10) p3:part(s2,3);
(%o10) -(1/3)
(%i11) s2-(p1+p2+p3);
(%o11) 0

/* p1 and p2 have the same issue */
(%i12) rectform(float(p1));
(%o12) 0.16342849479571833-0.28306645639069*%i
(%i13) simplify\-after\-subst(p1);
(%o13) (3^(3/2)*sqrt(27*%pi^2-32*%pi)-27*%pi+16)^(1/3)/(3*2^(1/3))
(%i14) rectform(float(%));
(%o14) -0.32685698959143794

(%i15) rectform(float(p2));
(%o15) 1.1775797725508834*%i+0.6798759986745109
(%i16) simplify\-after\-subst(p2);
(%o16) -((3^(9/2)*sqrt(27*%pi^2-32*%pi)+729*%pi-432)^(1/3)/(9*2^(1/3)))
(%i17) rectform(float(%));
(%o17) -1.3597519973490189

[David Billinghurst]

simplify-after-subst tries two approaches
1) ratsimp
2) if $constantp(e) sqrtdenest + rectform + rootscontract + ratsimp
then takes smallest expression

rectform is causing problems. p1 and p2 are expressions from above

(%i2) display2d:false$
(%i3) p1:((sqrt(3)*%i)/2-1/2)*((sqrt(%pi)*sqrt(27*%pi-32))/(2*3^(3/2))
                              +(-(3*(%pi-1))-1)/6+(-1)/27)
                              ^(1/3)$
(%i4) rectform(float(p1));
(%o4) 0.16342849479571833-0.28306645639069*%i
(%i5) r1:rectform(p1);
(%o5) -(-((sqrt(%pi)*sqrt(27*%pi-32))/(2*3^(3/2)))-(-(3*(%pi-1))-1)/6+1/27)^(1
                                                                            /3)
(%i6) rectform(float(%));
(%o6) -0.32685698959143666

(%i7) p2:(4*(-((sqrt(3)*%i)/2)-1/2))/(9*((sqrt(%pi)*sqrt(27*%pi-32))
                                        /(2*3^(3/2))
                                        +(-(3*(%pi-1))-1)/6+(-1)/27)
                                        ^(1/3))$
(%i8) rectform(float(p2));
(%o8) 1.1775797725508834*%i+0.6798759986745109
(%i9) r2:rectform(p2);
(%o9) -(4/(9*(-((sqrt(%pi)*sqrt(27*%pi-32))/(2*3^(3/2)))
             -(-(3*(%pi-1))-1)/6+1/27)
             ^(1/3)))
(%i10) rectform(float(%));
(%o10) -1.3597519973490217

This behavior is unexpected. The manual says

 -- Function: rectform (<expr>)

     Returns an expression ‘a + b %i’ equivalent to <expr>, such that
     <a> and <b> are purely real.

     Example:

          (%i1) rectform(sqrt(2)*%e^(%i*%pi/4));
          (%o1)                        %i + 1
          (%i2) rectform(sqrt(b^2+a^2)*%e^(%i*atan2(b, a)));
          (%o2)                       %i b + a
          (%i3) rectform(sqrt(5)*%e^(%i*atan(2)));
          (%o3)                       2 %i + 1

[David Billinghurst]

After removing the call to $rectform from SIMPLIFY-AFTER-SUBST algsys returns three solutions.

(%i9) eq1:x^3+x^2-x+%pi-1$
(%i10) sa:algsys([eq1],[x]);
(%o10) [[x = -((2^(2/3)*(3^(9/2)*sqrt(27*%pi^2-32*%pi)+729*%pi-432)^(1/3)
           -3*2^(2/3)*(3^(3/2)*sqrt(27*%pi^2-32*%pi)-27*%pi+16)^(1/3)+6)
           /18)],
        [x = ((2^(2/3)*sqrt(3)*%i+2^(2/3))
           *(3^(9/2)*sqrt(27*%pi^2-32*%pi)+729*%pi-432)^(1/3)
           +(2^(2/3)*3^(3/2)*%i-3*2^(2/3))
            *(3^(3/2)*sqrt(27*%pi^2-32*%pi)-27*%pi+16)^(1/3)-12)
           /36],
        [x = -(((2^(2/3)*sqrt(3)*%i-2^(2/3))
           *(3^(9/2)*sqrt(27*%pi^2-32*%pi)+729*%pi-432)^(1/3)
           +(2^(2/3)*3^(3/2)*%i+3*2^(2/3))
            *(3^(3/2)*sqrt(27*%pi^2-32*%pi)-27*%pi+16)^(1/3)+12)
           /36)]]

My notes in algsys.lisp say

;; $ratsimp with algebraic:true can transform
;;     sqrt(2)*sqrt(-1/(sqrt(3)*%i+1)) => (sqrt(3)*%i)/2+1/2
;; but $rectform is required for
;;     sqrt(sqrt(3)*%i-1)) => (sqrt(3)*%i)/sqrt(2)+1/sqrt(2)
;; and $rootscontract is required for
;;     sqrt(34)-sqrt(2)*sqrt(17) => 0
[...]
;; Rectform does more than is wanted.  A function that denests and
;; rationalizes nested complex radicals would be better.

It may be time to extend sqrtdenest or raddenest to expressions of the form sqrt(sqrt(3)*%i-1)).

[David Billinghurst]

I haven't found too many references to denesting complex expressions . There is

https://brownmath.com/alge/nestrad.htm

[David Billinghurst]

The sqrtdenest enhancement is now bug [#4255]. It is worthwhile enhancement independently of this

[David Billinghurst]

Three solutions are returned when domain:complex .

(%i1) display2d:false$
(%i2) eq1:x^3+x^2-x+%pi-1$
(%i3) domain:complex$
(%i4) algsys([eq1],[x]);
(%o4) [[x = (4*(54/(9*sqrt((27*%pi^2-32*%pi)/3)-27*%pi+16))^(1/3)
          +9*((9*sqrt((27*%pi^2-32*%pi)/3)-27*%pi+16)/54)^(1/3)-3)
          /9],
       [x = -((4*(-(54/(9*sqrt((27*%pi^2-32*%pi)/3)-27*%pi+16)))^(1/3)
          +9*(-((9*sqrt((27*%pi^2-32*%pi)/3)-27*%pi+16)/54))^(1/3)+3)
          /9)],
       [x = ((4*sqrt(3)*%i-4)*(54/(9*sqrt((27*%pi^2-32*%pi)/3)-27*%pi+16))
                              ^(1/3)
          +(-(3^(5/2)*%i)-9)*((9*sqrt((27*%pi^2-32*%pi)/3)-27*%pi+16)/54)
                             ^(1/3)-6)
          /18]]

This is understandable as the roots contain cube roots of negative real numbers, which have different values depending on setting of domain.

[David Billinghurst]

See the thread https://sourceforge.net/p/maxima/mailman/maxima-discuss/thread/CAMTHLKj1zk9G8H%3Dfpqi5Mx63jfXt4tH7FSHXQsvM7ubNCh1h1Q%40mail.gmail.com/#msg58730878

In particular https://sourceforge.net/p/maxima/mailman/message/58731017/ where Stavros says

domain:real/complex sounds like a really powerful global setting -- am I
working in the reals or in the complex numbers?

But in fact, it has a very narrow meaning -- it (together with m1pbranch and
radexpand) simply controls the simplification of (x^a)^b in certain cases.

It doesn't make sense for it to modify the behavior of
rectform/polarform/....

Note that domain:complex affects simplification, so it may be necessary to set domain:complex much higher up the algsys call tree. What does solve do? It seems to return the correct solution to algsys, which then breaks stuff.