Diagonal solver precision bug

[wilrop]

There is a bug in the diagonal solver that makes the problem fail when going to precisions dd and qd. Below is a minimal working example of the bug in practice.

from phcpy.solver import solve
from phcpy.sets import embed
from phcpy.diagonal import diagonal_solver as diagsolve

polys1 = ['+ 3 * a - 2 * b;', '+ c + d - 1;']
polys2 = ['+ 12 * c - 4 * d;', '+ a + b - 1;']

precision = 'dd'

comb = polys1 + polys2
sols = solve(comb, verbose=False, precision=precision)
for sol in sols:
    print(sol)

print('-------------------')

nvar = 4
dim = nvar - len(polys1)
emb1 = embed(nvar, dim, polys1, precision=precision)
emb2 = embed(nvar, dim, polys2, precision=precision)

emb1[0] = '+ (a + b + c + d) - (a + b + c + d) + ' + emb1[0]
emb2[0] = '+ (a + b + c + d) - (a + b + c + d) + ' + emb2[0]

sols1 = solve(emb1, verbose=False, precision=precision)
sols2 = solve(emb2, verbose=False, precision=precision)

out_polys, out_sols = diagsolve(nvar, dim, emb1, sols1, dim, emb2, sols2, verbose=False, prc=precision)

for sol in out_sols:
    print(sol)

The system computes the correct solution when using the black box solver for all precisions. However, when splitting it up into two distinct systems and computing the intersection of witness sets, the solution is not found under the higher precisions. The solution can still be found when using precision d.

I'm not sure if it is relevant but in both higher precisions, the err and res values are essentially zero while rco is 1. When you set the precision to d this is not the case.