Fix factor-of-2 in PolyOptimizer quadratic coupling encoding

[gcattan]

Both quadratic branches of _encode() wrote half the coupling weight they should. The comments ("factor-of-2 symmetry accounted for") show the reasoning: the author assumed a symmetric pair contributes 2*J[i,j] to the energy and pre-divided by two.

It does not. Solvers read J through the 1/2 * s.J.s convention -- see the energy readout in CaSuDaSolver.solve(), and note that field = m @ J is exactly the gradient of that form -- under which the (i,j) and (j,i) entries and the leading 1/2 cancel, so a symmetric pair contributes J[i,j] undivided.

h was already correct, so the encoded landscape was not a rescaling of the target polynomial but a distortion of it: the linear part kept its weight while the quadratic part lost half of its own. Minimising -0.75x - 0.75y + xy over {0,1}^2 flattens the true optima (1,0)/(0,1) at -0.75 into a three-way tie with (1,1) at -0.5. Across 40 random polynomials, 37 encodings failed to be an affine image of their polynomial; doubling J fixes all 40 exactly.

p-kit already had the right convention elsewhere: TransverseFieldIsing.from_qubo uses -Q/4 for the pair coefficient under the identical (1+s)/2 mapping. The two converters disagreed by exactly this factor of 2.

test_quadratic_J had pinned the buggy value, having been written from the implementation rather than the mathematics, which is how this survived. It now asserts the derived value. The new tests check that J/h reproduce the polynomial up to an additive constant -- comparing only argmin is unsafe, since a mis-scaled coupling produces ties that an argmin assertion passes on an arbitrary tie-break.

DocplexOptimizer inherits the fix.

[toncho11]

Yes, looks good.