janverschelde / PHCpack

The primary source code repository for PHCpack, a software package to solve polynomial systems with homotopy continuation methods.
http://www.phcpack.org
GNU General Public License v3.0
59 stars 21 forks source link

sets.embed gives wrong embedding #59

Closed wilrop closed 2 years ago

wilrop commented 2 years ago

I've found a bug in the sets.embed function where it will give a wrong embedding. Consider the short piece of code underneath:

from phcpy.sets import embed
from phcpy.solver import solve

start_polys = ['(5 * p1_0*p2_0 + 9 * p1_0*p2_2 + 5 * p1_2*p2_0 + 1 * p1_2*p2_2) - (5 * p1_0*p2_0 + 1 * p1_0*p2_2 + 2 * p1_2*p2_0 + 9 * p1_2*p2_2);', 'p0_1 + p0_2 - 1;', '(6 * p0_1*p2_0 + 7 * p0_1*p2_2 + 4 * p0_2*p2_0 + 3 * p0_2*p2_2) - (9 * p0_1*p2_0 + 4 * p0_1*p2_2 + 8 * p0_2*p2_0 + 4 * p0_2*p2_2);', 'p1_0 + p1_2 - 1;', '(5 * p0_1*p1_0 + 7 * p0_1*p1_2 + 2 * p0_2*p1_0 + 6 * p0_2*p1_2) - (4 * p0_1*p1_0 + 4 * p0_1*p1_2 + 9 * p0_2*p1_0 + 3 * p0_2*p1_2);', 'p2_0 + p2_2 - 1;']
embedded_polys = embed(6, 0, start_polys)
correct_embedded = ['8*p1_0*p2_2 + 3*p2_0*p1_2 - 8*p2_2*p1_2;', 'p0_1 + p0_2 - 1;', '-3*p2_0*p0_1 - 4*p2_0*p0_2 + 3*p2_2*p0_1 - p2_2*p0_2;', 'p1_0 + p1_2 - 1;', 'p1_0*p0_1 - 7*p1_0*p0_2 + 3*p1_2*p0_1 + 3*p1_2*p0_2;', 'p2_0 + p2_2 - 1;']

print(start_polys)
print(embedded_polys)

solutions = solve(start_polys, verbose=False)
for solution in solutions:
    print(solution)

print("next")
emb_sol = solve(embedded_polys, verbose=False)
for solution in emb_sol:
    print(solution)

print("next") 
correct_emb_sol = solve(correct_embedded, verbose=False)
for solution in correct_emb_sol:
    print(solution)

This gives as an output:

['(5 * p1_0*p2_0 + 9 * p1_0*p2_2 + 5 * p1_2*p2_0 + 1 * p1_2*p2_2) - (5 * p1_0*p2_0 + 1 * p1_0*p2_2 + 2 * p1_2*p2_0 + 9 * p1_2*p2_2);', 'p0_1 + p0_2 - 1;', '(6 * p0_1*p2_0 + 7 * p0_1*p2_2 + 4 * p0_2*p2_0 + 3 * p0_2*p2_2) - (9 * p0_1*p2_0 + 4 * p0_1*p2_2 + 8 * p0_2*p2_0 + 4 * p0_2*p2_2);', 'p1_0 + p1_2 - 1;', '(5 * p0_1*p1_0 + 7 * p0_1*p1_2 + 2 * p0_2*p1_0 + 6 * p0_2*p1_2) - (4 * p0_1*p1_0 + 4 * p0_1*p1_2 + 9 * p0_2*p1_0 + 3 * p0_2*p1_2);', 'p2_0 + p2_2 - 1;']
['8*p1_0*p2_2 + 3*p2_0*p1_2 - 8*p2_2*p1_2;', 'p0_1 + p0_2 - 1;', ' - -3*p2_0*p0_1 - 4*p2_0*p0_2 + 3*p2_2*p0_1 - p2_2*p0_2;', 'p1_0 + p1_2 - 1;', 'p1_0*p0_1 - 7*p1_0*p0_2 + 3*p1_2*p0_1 + 3*p1_2*p0_2;', 'p2_0 + p2_2 - 1;']
t :  1.00000000000000E+00   0.00000000000000E+00
m : 1
the solution for t :
 p1_0 :  4.76373451010818E-01   9.58807317440962E-94
 p2_0 :  1.93967840523004E-01  -5.75284390464577E-93
 p2_2 :  8.06032159476996E-01   4.11036594881613E-93
 p1_2 :  5.23626548989182E-01  -1.91761463488192E-93
 p0_1 :  4.62802478382692E-01  -1.91761463488192E-93
 p0_2 :  5.37197521617308E-01   1.91761463488192E-93
== err :  3.753E-16 = rco :  1.601E-01 = res :  1.665E-16 =
t :  1.00000000000000E+00   0.00000000000000E+00
m : 1
the solution for t :
 p1_0 : -4.91998451010818E-01   1.37158486953873E-99
 p2_0 :  7.80027431344608E-01   2.28597478256455E-100
 p2_2 :  2.19972568655392E-01  -2.28597478256455E-100
 p1_2 :  1.49199845101082E+00  -1.82877982605164E-99
 p0_1 :  2.01219752161731E+00   0.00000000000000E+00
 p0_2 : -1.01219752161731E+00  -4.57194956512910E-100
== err :  1.244E-16 = rco :  6.570E-02 = res :  8.327E-17 =
next
t :  1.00000000000000E+00   0.00000000000000E+00
m : 1
the solution for t :
 p1_0 : -2.94117647058824E-01   0.00000000000000E+00
 p2_2 :  2.34042553191489E-01  -2.41291762006314E-86
 p2_0 :  7.65957446808511E-01   1.60861174670876E-86
 p1_2 :  1.29411764705882E+00   0.00000000000000E+00
 p0_1 :  2.52500000000000E+00   1.28688939736701E-85
 p0_2 : -1.52500000000000E+00  -1.28688939736701E-85
== err :  8.816E-17 = rco :  5.665E-02 = res :  5.551E-17 =
t :  1.00000000000000E+00   0.00000000000000E+00
m : 1
the solution for t :
 p1_0 :  5.00000000000000E-01   0.00000000000000E+00
 p2_2 :  1.00000000000000E+00   6.07716335728627E-64
 p2_0 :  2.79168267377749E-18  -6.48230758110536E-64
 p1_2 :  5.00000000000000E-01   0.00000000000000E+00
 p0_1 :  5.00000000000000E-01  -1.21543267145725E-63
 p0_2 :  5.00000000000000E-01   1.21543267145725E-63
== err :  1.642E-47 = rco :  1.675E-01 = res :  2.410E-63 =
next
t :  1.00000000000000E+00   0.00000000000000E+00
m : 1
the solution for t :
 p1_0 :  4.76373451010818E-01   1.02951151789361E-84
 p2_2 :  8.06032159476996E-01   0.00000000000000E+00
 p2_0 :  1.93967840523004E-01   0.00000000000000E+00
 p1_2 :  5.23626548989182E-01  -7.72133638420204E-85
 p0_1 :  4.62802478382692E-01  -2.05902303578721E-84
 p0_2 :  5.37197521617308E-01   2.05902303578721E-84
== err :  3.753E-16 = rco :  1.601E-01 = res :  1.110E-16 =
t :  1.00000000000000E+00   0.00000000000000E+00
m : 1
the solution for t :
 p1_0 : -4.91998451010818E-01   5.14755758946803E-85
 p2_2 :  2.19972568655392E-01   0.00000000000000E+00
 p2_0 :  7.80027431344608E-01   0.00000000000000E+00
 p1_2 :  1.49199845101082E+00  -5.14755758946803E-85
 p0_1 :  2.01219752161731E+00   5.14755758946803E-85
 p0_2 : -1.01219752161731E+00  -5.14755758946803E-85
== err :  2.909E-16 = rco :  6.570E-02 = res :  3.608E-16 =

Observe that because the original polynomial system is already square (six variables and six polynomials), the embedding should not change it. I would expect this to be the intended behaviour when calling embed with the topdim set to 0. However, as you can see, the solutions for the first system are different from the second. The problem is in the third entry of the embedded polynomials: ' - -3*p2_0*p0_1 - 4*p2_0*p0_2 + 3*p2_2*p0_1 - p2_2*p0_2;' This reduction is not correct. The minus at the beginning should not be there. I have removed this minus sign in the correct_embedded system which then gives the same solutions as the original system.

For now, I'm working around this simply by doing a check on the length of the input system and not embedding it when it is already square. However, I'm unsure if the bug also manifests itself in cases where it is actually necessary to embed the system.

janverschelde commented 2 years ago

Thanks for reporting this bug. The superfluous minus for negative integers in the output of the embed function no longer happens.

wilrop commented 2 years ago

I have verified this and can confirm that it also works on my system. Thanks for solving it! You can close this issue if you like.

janverschelde commented 2 years ago

Thanks for checking. I will close this issue.