## #453 algsys fails in simple case

None
closed
5
2016-10-06
2003-11-28
No

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

solve(eqs,[u,v]) =&gt; []

But in fact eqs is solvable:

eqs1: [ eqs[1], eqs[1]-eqs[2] ];
solve(eqs1, [u,v]) =&gt;
[[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

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).

## Discussion

• Robert Dodier - 2006-07-13
• labels: --> 460522

• Robert Dodier - 2006-09-09
• labels: 460522 --> Lisp Core - Solving equations

• Barton Willis - 2006-09-15

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],...]

• 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]]
```

#### Related

Last edit: David Billinghurst 2016-10-06
• assigned_to: David Billinghurst
• Group: --> None

• 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]
```

#### Related

• status: open --> closed

• 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]];
```