bug(polys): fix multivariate square free factorisation #26514
Conversation
|
✅ Hi, I am the SymPy bot. I'm here to help you write a release notes entry. Please read the guide on how to write release notes. Your release notes are in good order. Here is what the release notes will look like:
This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.13. Click here to see the pull request description that was parsed.
Update The release notes on the wiki have been updated. |
sympy/polys/sqfreetools.py
Outdated
| g, p, q = dmp_inner_gcd(p, h, u, K) | ||
|
|
||
| if all or dmp_degree(g, u) > 0: | ||
| result.append((g, i)) |
There was a problem hiding this comment.
Should the i here be the dmp_degree(g,u)? The abstractions are losing me a bit, but I know that dmp_sqf_list returns the powers of the factors as the 2nd element of the tuple whereas here the indexed variable i is being put in that position.
There was a problem hiding this comment.
The idea here is that we have a polynomial p = p1 * p2**2 * p3**3 * ... and we want to find p1, p2, etc which are coprime. The way we do that is by differentiating the polynomial and computing the gcd like g = gcd(p, p'). This works because if p2**2 divides p then p2 must divide both p and p' (you can check by the product rule). Therefore p2 must divide g. However p2 will not divide p''. We can find p2 then by observing that it is a factor of p and p' but not of p''. What matters for the multiplicity of the factor is not the degree of g or the degree of p2 but rather how many times we have differentiated i.e. i.
The actual (Yun's) algorithm is a little more complicated than this but it also described on Wikipedia:
https://en.wikipedia.org/wiki/Square-free_polynomial#Yun's_algorithm
What is not mentioned there but is mentioned in Yun's paper is that for multivariate polynomials like Q[x,y,z] we need to think of them as being like p(z) in Q[x,y][z] where the coefficients of p(z) (elements of Q[x,y]) should be coprime. In other words p(z) should be primitive in the sense of a univariate polynomial being primitive. This is the difference between dmp_ground_primitive (make primitive in Q[x,y,z]) and dmp_primitive (make primitive in Q[x,y][z]). With dmp_ground_primitive the extracted content is an element of Q but with dmp_primitive it is an element of Q[x,y]. The change I made here is to extract that content in Q[x,y] and then recursively compute its square free factorisation.
What is potentially worth changing is that this version of dmp_sqf_list does not produce the factors in order whereas the previous code does:
In [1]: R, x, y = ring(['x', 'y'], QQ)
In [2]: p = (x+y)**2*(x-1)*(y-1)*(x-2)**2*(y-2)**2
In [3]: R.dmp_sqf_list(p)
Out[3]: (MPQ(1,1), [(x - 1, 1), (x**2 + x*y - 2*x - 2*y, 2), (y - 1, 1), (y - 2, 2)])With master that would be:
In [5]: R.dmp_sqf_list(p)
Out[5]: (MPQ(1,1), [(x - 1, 1), (x**2 + x*y - 2*x - 2*y, 2)])The master result is incorrect but it always has a well define format with ascending powers because results are appended as i increases.
The higher level sqf function somehow accommodates for this but Poly.sqf does not:
In [1]: p = (x+y)**2*(x-1)*(y-1)*(x-2)**2*(y-2)**2
In [2]: sqf(expand(p))
Out[2]:
2
⎛ 2 2 2 2 ⎞
(x⋅y - x - y + 1)⋅⎝x ⋅y - 2⋅x + x⋅y - 4⋅x⋅y + 4⋅x - 2⋅y + 4⋅y⎠
In [4]: Poly(p).sqf_list()
Out[4]: (1, [(Poly(x - 1, x, y, domain='ZZ'), 1), (Poly(x**2 + x*y - 2*x - 2*y, x, y, domain='ZZ'), 2), (Poly(y - 1, x, y, domain='ZZ'), 1), (Poly(y - 2, x, y, domain='ZZ'), 2)])I'll patch this up in dmp_sqf_list.
|
Benchmark results from GitHub Actions Lower numbers are good, higher numbers are bad. A ratio less than 1 Significantly changed benchmark results (PR vs master) Significantly changed benchmark results (master vs previous release) | Change | Before [2487dbb5] | After [b9bafe06] | Ratio | Benchmark (Parameter) |
|----------|----------------------|---------------------|---------|----------------------------------------------------------------------|
| - | 67.3±1ms | 43.6±0.3ms | 0.65 | integrate.TimeIntegrationRisch02.time_doit(10) |
| - | 66.2±0.3ms | 43.6±0.4ms | 0.66 | integrate.TimeIntegrationRisch02.time_doit_risch(10) |
| + | 18.5±0.4μs | 30.9±0.3μs | 1.67 | integrate.TimeIntegrationRisch03.time_doit(1) |
| - | 5.26±0.07ms | 2.83±0.01ms | 0.54 | logic.LogicSuite.time_load_file |
| - | 73.0±0.4ms | 28.8±0.1ms | 0.39 | polys.TimeGCD_GaussInt.time_op(1, 'dense') |
| - | 25.9±0.1ms | 16.9±0.1ms | 0.66 | polys.TimeGCD_GaussInt.time_op(1, 'expr') |
| - | 73.6±0.7ms | 29.0±0.1ms | 0.39 | polys.TimeGCD_GaussInt.time_op(1, 'sparse') |
| - | 256±1ms | 126±0.9ms | 0.49 | polys.TimeGCD_GaussInt.time_op(2, 'dense') |
| - | 258±1ms | 125±0.3ms | 0.48 | polys.TimeGCD_GaussInt.time_op(2, 'sparse') |
| - | 661±6ms | 373±0.5ms | 0.56 | polys.TimeGCD_GaussInt.time_op(3, 'dense') |
| - | 661±4ms | 375±2ms | 0.57 | polys.TimeGCD_GaussInt.time_op(3, 'sparse') |
| - | 499±5μs | 283±1μs | 0.57 | polys.TimeGCD_LinearDenseQuadraticGCD.time_op(1, 'dense') |
| - | 1.79±0.02ms | 1.05±0.01ms | 0.59 | polys.TimeGCD_LinearDenseQuadraticGCD.time_op(2, 'dense') |
| - | 5.82±0.1ms | 3.04±0.02ms | 0.52 | polys.TimeGCD_LinearDenseQuadraticGCD.time_op(3, 'dense') |
| - | 460±7μs | 232±2μs | 0.5 | polys.TimeGCD_QuadraticNonMonicGCD.time_op(1, 'dense') |
| - | 1.51±0.01ms | 674±3μs | 0.45 | polys.TimeGCD_QuadraticNonMonicGCD.time_op(2, 'dense') |
| - | 4.96±0.2ms | 1.67±0ms | 0.34 | polys.TimeGCD_QuadraticNonMonicGCD.time_op(3, 'dense') |
| - | 376±3μs | 208±2μs | 0.55 | polys.TimeGCD_SparseGCDHighDegree.time_op(1, 'dense') |
| - | 2.43±0.01ms | 1.24±0.01ms | 0.51 | polys.TimeGCD_SparseGCDHighDegree.time_op(3, 'dense') |
| - | 10.2±0.06ms | 4.36±0.03ms | 0.43 | polys.TimeGCD_SparseGCDHighDegree.time_op(5, 'dense') |
| - | 368±4μs | 171±3μs | 0.47 | polys.TimeGCD_SparseNonMonicQuadratic.time_op(1, 'dense') |
| - | 2.52±0.01ms | 907±2μs | 0.36 | polys.TimeGCD_SparseNonMonicQuadratic.time_op(3, 'dense') |
| - | 9.70±0.1ms | 2.69±0.02ms | 0.28 | polys.TimeGCD_SparseNonMonicQuadratic.time_op(5, 'dense') |
| - | 1.04±0.01ms | 422±5μs | 0.41 | polys.TimePREM_LinearDenseQuadraticGCD.time_op(3, 'dense') |
| - | 1.71±0.02ms | 500±2μs | 0.29 | polys.TimePREM_LinearDenseQuadraticGCD.time_op(3, 'sparse') |
| - | 5.85±0.03ms | 1.81±0.02ms | 0.31 | polys.TimePREM_LinearDenseQuadraticGCD.time_op(5, 'dense') |
| - | 8.60±0.04ms | 1.50±0.01ms | 0.17 | polys.TimePREM_LinearDenseQuadraticGCD.time_op(5, 'sparse') |
| - | 282±0.9μs | 65.3±0.5μs | 0.23 | polys.TimePREM_QuadraticNonMonicGCD.time_op(1, 'sparse') |
| - | 3.44±0.03ms | 402±1μs | 0.12 | polys.TimePREM_QuadraticNonMonicGCD.time_op(3, 'dense') |
| - | 3.94±0.02ms | 278±1μs | 0.07 | polys.TimePREM_QuadraticNonMonicGCD.time_op(3, 'sparse') |
| - | 7.15±0.06ms | 1.28±0.01ms | 0.18 | polys.TimePREM_QuadraticNonMonicGCD.time_op(5, 'dense') |
| - | 8.83±0.06ms | 833±2μs | 0.09 | polys.TimePREM_QuadraticNonMonicGCD.time_op(5, 'sparse') |
| - | 5.00±0.02ms | 2.97±0.01ms | 0.59 | polys.TimeSUBRESULTANTS_LinearDenseQuadraticGCD.time_op(2, 'sparse') |
| - | 12.2±0.2ms | 6.56±0.07ms | 0.54 | polys.TimeSUBRESULTANTS_LinearDenseQuadraticGCD.time_op(3, 'dense') |
| - | 22.1±0.1ms | 8.98±0.04ms | 0.41 | polys.TimeSUBRESULTANTS_LinearDenseQuadraticGCD.time_op(3, 'sparse') |
| - | 5.22±0.01ms | 864±5μs | 0.17 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(1, 'sparse') |
| - | 12.5±0.05ms | 7.04±0.03ms | 0.56 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(2, 'sparse') |
| - | 102±0.6ms | 25.9±0.09ms | 0.25 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(3, 'dense') |
| - | 166±0.8ms | 53.4±0.1ms | 0.32 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(3, 'sparse') |
| - | 177±1μs | 112±0.6μs | 0.63 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(1, 'dense') |
| - | 363±2μs | 215±0.8μs | 0.59 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(1, 'sparse') |
| - | 4.33±0.05ms | 838±7μs | 0.19 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(3, 'dense') |
| - | 5.30±0.03ms | 381±2μs | 0.07 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(3, 'sparse') |
| - | 20.7±0.1ms | 2.84±0.01ms | 0.14 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(5, 'dense') |
| - | 22.9±0.1ms | 628±3μs | 0.03 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(5, 'sparse') |
| - | 482±3μs | 136±1μs | 0.28 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(1, 'sparse') |
| - | 4.74±0.06ms | 618±2μs | 0.13 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(3, 'dense') |
| - | 5.31±0.06ms | 140±1μs | 0.03 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(3, 'sparse') |
| - | 13.2±0.06ms | 1.30±0.01ms | 0.1 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(5, 'dense') |
| - | 14.1±0.07ms | 145±1μs | 0.01 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(5, 'sparse') |
| - | 131±0.3μs | 74.4±2μs | 0.57 | solve.TimeMatrixOperations.time_rref(3, 0) |
| - | 252±1μs | 87.7±0.4μs | 0.35 | solve.TimeMatrixOperations.time_rref(4, 0) |
| - | 24.2±0.1ms | 10.1±0.05ms | 0.42 | solve.TimeSolveLinSys189x49.time_solve_lin_sys |
| - | 28.0±0.08ms | 15.3±0.07ms | 0.55 | solve.TimeSparseSystem.time_linsolve_Aaug(20) |
| - | 54.4±0.1ms | 25.1±0.3ms | 0.46 | solve.TimeSparseSystem.time_linsolve_Aaug(30) |
| - | 28.0±0.2ms | 15.2±0.04ms | 0.54 | solve.TimeSparseSystem.time_linsolve_Ab(20) |
| - | 53.9±0.3ms | 24.7±0.09ms | 0.46 | solve.TimeSparseSystem.time_linsolve_Ab(30) |
Full benchmark results can be found as artifacts in GitHub Actions |
|
Thanks |
References to other Issues or PRs
Partial fix for gh-26497
Brief description of what is fixed or changed
Other comments
Release Notes
sqfwas fixed. Previously incorrect results were returned in some cases when computing the square free factorisation of multivariate polynomials.