#453 algsys fails in simple case

[Stavros Macrakis]

eqs:
[ v^2 - 2*v + u^2 + 10*u + 1,
v^2 - 10*v + u^2 + 2*u + 1 ]

solve(eqs,[u,v]) => []

But in fact eqs is solvable:

eqs1: [ eqs[1], eqs[1]-eqs[2] ];
solve(eqs1, [u,v]) =>
[[u = -(sqrt(34) + 6) / 2, v = (sqrt(2)*sqrt(17) + 6) /
2],
[u = (sqrt(34) - 6) / 2, v = -(sqrt(2)*sqrt(17) - 6) / 2]]

That solution is correct, as you can verify with subst
followed by radcan.

For that matter, eliminate also works.

Maxima 5.9.0/W2k
tested with both gcd:subres and gcd:spmod

Found starting with a problem report submitted by Kirk
Lancaster (27 Nov 2003).

[Barton Willis]

Logged In: YES
user_id=895922

Maxima can't even solve a triangular form of these equations:

(%o5) [2*v^2-12*v+1,v^2-2*v+u^2+10*u+1]
(%i6) algsys(%,[v,u]);
(%o6) []

Eliminating v^2 in the second equation allows Maxima to finish:

(%i7) [%o5[1],ratsubst((12*v-1)/2,v^2,%o5[2])];
(%o7) [2*v^2-12*v+1,(8*v+2*u^2+20*u+1)/2]
(%i8) algsys(%,[v,u]);
(%o8) [[v=-(sqrt(2)*sqrt(17)-6)/2,u=-(sqrt(34)+14)/2],...]

[David Billinghurst]

algsys can now solve the triangular form ot the equations after fixing [bugs:#3208]

(%i1) eqs:[2*v^2-12*v+1,v^2-2*v+u^2+10*u+1]$

(%i2) algsys(eqs,[v,u]);

(%o2) [[v = -(sqrt(2)*sqrt(17)-6)/2,u = -(sqrt(34)+14)/2],
       [v = (sqrt(2)*sqrt(17)+6)/2,u = (sqrt(34)-14)/2],
       [v = (sqrt(2)*sqrt(17)+6)/2,u = -(sqrt(34)+6)/2],
       [v = -(sqrt(2)*sqrt(17)-6)/2,u = (sqrt(34)-6)/2]]

[David Billinghurst]

I have traced this with git HEAD (after recent commits [2986fc] and [640ca7] ) and reproduced the process manually.

algsys rejects the first solution as 4sqrt(34)+(sqrt(2)sqrt(17)-6)^2/2+2(sqrt(2)sqrt(17)-6)-23 isn't simplified to zero by new function simplify-after-solve, which use sqrtdenest and ratsimp() with algebraic=true.

Even sqrt(34)-sqrt(2)*sqrt(17) defies ratsimp().

rootscontract() will simplify this this. What else does it do?

(%i2) eq1:y^2-2*y+x^2+10*x+1
(%i3) eq2:y^2-10*y+x^2+2*x+1

/* Eliminate y */

(%i4) r:resultant(eq1,eq2,y)
                                    2
(%o4)                        64 (2 x  + 12 x + 1)
(%i5) rf:factor(r)
                                    2
(%o5)                        64 (2 x  + 12 x + 1)
(%i6) rf:inpart(rf,2)
                                   2
(%o6)                           2 x  + 12 x + 1

/* Solve resultant for x */

(%i7) X:solve(rf,x)
                           sqrt(34) + 6      sqrt(34) - 6
(%o7)               [x = - ------------, x = ------------]
                                2                 2
(%i8) X:reverse(X)
                         sqrt(34) - 6        sqrt(34) + 6
(%o8)               [x = ------------, x = - ------------]
                              2                   2

/* Substitute solution X[1] into eq1 and eq2 */

(%i9) eq11:subst(X[1],eq1)
                                                           2
                2                            (sqrt(34) - 6)
(%o9)          y  - 2 y + 5 (sqrt(34) - 6) + --------------- + 1
                                                    4
(%i10) ev(eq11:ratsimp(2*eq11),algebraic = true)
                            2
(%o10)                   2 y  - 4 y + 4 sqrt(34) - 23
(%i11) eq21:subst(X[1],eq2)
                                                       2
                   2                     (sqrt(34) - 6)
(%o11)            y  - 10 y + sqrt(34) + --------------- - 5
                                                4
(%i12) ev(eq21:ratsimp(2*eq21),algebraic = true)
                            2
(%o12)                   2 y  - 20 y - 4 sqrt(34) + 25

/* Check resultant - should = 0 */

(%i13) r2:resultant(eq11,eq21,y)
(%o13)                                 0

/* Subst X[1] into eq2 and solve for y */
(%i14) Y1:solve(eq21,y)
              sqrt(2) sqrt(4 sqrt(34) + 25) - 10
(%o14) [y = - ----------------------------------, 
                              2
                                            sqrt(2) sqrt(4 sqrt(34) + 25) + 10
                                        y = ----------------------------------]
                                                            2
(%i15) Y1:sqrtdenest(Y1)
                                   3/2
              sqrt(2) (sqrt(17) + 2   ) - 10
(%o15) [y = - ------------------------------, 
                            2
                                                                     3/2
                                                sqrt(2) (sqrt(17) + 2   ) + 10
                                            y = ------------------------------]
                                                              2
(%i16) ev(Y1:ratsimp(Y1),algebraic = true)
                   sqrt(2) sqrt(17) - 6      sqrt(2) sqrt(17) + 14
(%o16)      [y = - --------------------, y = ---------------------]
                            2                          2

/* Substitute soln Y[1] into eq1 */

(%i17) eq111:subst(Y1[1],eq11)
                                          2
                    (sqrt(2) sqrt(17) - 6)
(%o17) 4 sqrt(34) + ----------------------- + 2 (sqrt(2) sqrt(17) - 6) - 23
                               2
(%i18) ev(eq111:ratsimp(eq111),algebraic = true)
                                        5/2
(%o18)                    4 sqrt(34) - 2    sqrt(17)

/* PROBLEM: eq111 == 0, but algsys dosn't know */

(%i19) float(eq111)
(%o19)                      - 7.105427357601002e-15

/* Found a solution but reject it */

(%i20) s:[X[1],Y1[1]]
                     sqrt(34) - 6        sqrt(2) sqrt(17) - 6
(%o20)          [x = ------------, y = - --------------------]
                          2                       2
(%i21) ev([eq1,eq2],s)
                                         2                         2
                           (sqrt(34) - 6)    (sqrt(2) sqrt(17) - 6)
(%o21) [5 (sqrt(34) - 6) + --------------- + -----------------------
                                  4                     4
                                                  2                         2
                                    (sqrt(34) - 6)    (sqrt(2) sqrt(17) - 6)
 + sqrt(2) sqrt(17) - 5, sqrt(34) + --------------- + -----------------------
                                           4                     4
 + 5 (sqrt(2) sqrt(17) - 6) - 5]

(%i22) ev(ratsimp(%),algebraic = true)
                          3/2            3/2
(%o22)     [2 sqrt(34) - 2    sqrt(17), 2    sqrt(17) - 2 sqrt(34)]
(%i23) rootscontract(%)
(%o23)                              [0, 0]

[David Billinghurst]

Fixed in commit [cb3dc1]. Following tests added to testsuite.

/* SF bug #453 algsys fails in simple case */
algsys([v^2-2*v+u^2+10*u+1,v^2-10*v+u^2+2*u+1],[v,u]);
[[v = -((sqrt(34)-6)/2), u =   (sqrt(2)*sqrt(17)-6)/2],
 [v =   (sqrt(34)+6)/2,  u = -((sqrt(2)*sqrt(17)+6)/2)]];

algsys([2*v^2-12*v+1,v^2-2*v+u^2+10*u+1],[v,u]);

[[v = -((sqrt(2)*sqrt(17)-6)/2), u = -((sqrt(34)+14)/2)],
 [v =   (sqrt(2)*sqrt(17)+6)/2,  u =   (sqrt(34)-14)/2],
 [v =   (sqrt(2)*sqrt(17)+6)/2,  u = -((sqrt(34)+6)/2)],
 [v = -((sqrt(2)*sqrt(17)-6)/2), u =   (sqrt(34)-6)/2]];