Opened 12 years ago
Closed 12 years ago
#5958 closed defect (fixed)
[with patch, positive review] MPolynomial_polydict.factor() should accept proof parameter
Reported by: | malb | Owned by: | malb |
---|---|---|---|
Priority: | major | Milestone: | sage-4.1.2 |
Component: | commutative algebra | Keywords: | singular, factor |
Cc: | john_perry | Merged in: | Sage 4.1.2.alpha0 |
Authors: | John Perry | Reviewers: | Martin Albrecht |
Report Upstream: | Work issues: | ||
Branch: | Commit: | ||
Dependencies: | Stopgaps: |
Description
The parameter should be ignored, but for compatibility it is necessary.
E.g. this should work:
sage: R.<x,y> = CC[] sage: I = R.ideal(x^2+y^2-1,x*y-1) sage: I.variety()
but it raises an except in 3.4.1.
Attachments (4)
Change History (32)
comment:1 Changed 12 years ago by
comment:2 Changed 12 years ago by
It seems toy_variety
does not switch to 'lex' when it should?
Changed 12 years ago by
comment:3 Changed 12 years ago by
- Summary changed from [with patch, needs work] MPolynomial_polydict.factor() should accept proof parameter to [with patch, needs review] MPolynomial_polydict.factor() should accept proof parameter
The attached patch switches to lex before calling toy_variety.triangular_factorization
comment:4 Changed 12 years ago by
- Summary changed from [with patch, needs review] MPolynomial_polydict.factor() should accept proof parameter to [with patch, needs work] MPolynomial_polydict.factor() should accept proof parameter
Two issues:
- in the first patch, "proofably" should be "provably"
- I'm getting some slight numerical noise in the doctests from the second patch:
sage -t "devel/sage-main/sage/rings/polynomial/multi_polynomial_ideal.py" ********************************************************************** File "/opt/sage-devel/devel/sage-main/sage/rings/polynomial/multi_polynomial_ideal.py", line 1908: sage: for v in I.variety(): print v Expected: verbose 0 (...: multi_polynomial_ideal.py, variety) Warning: falling back to very slow toy implementation. {y: -0.866025403784439 - 0.500000000000000*I, x: -0.866025403784438 + 0.500000000000000*I} {y: -0.866025403784438 + 0.500000000000000*I, x: -0.866025403784438 - 0.499999999999999*I} {y: 0.866025403784438 + 0.500000000000001*I, x: 0.866025403784440 - 0.499999999999999*I} {y: 0.866025403784439 - 0.500000000000000*I, x: 0.866025403784439 + 0.500000000000000*I} Got: verbose 0 (1735: multi_polynomial_ideal.py, variety) Warning: falling back to very slow toy implementation. {y: -0.866025403784439 - 0.500000000000000*I, x: -0.866025403784439 + 0.500000000000001*I} {y: -0.866025403784439 + 0.500000000000000*I, x: -0.866025403784439 - 0.500000000000000*I} {y: 0.866025403784438 + 0.500000000000000*I, x: 0.866025403784438 - 0.499999999999999*I} {y: 0.866025403784439 - 0.500000000000000*I, x: 0.866025403784438 + 0.500000000000001*I} ********************************************************************** 1 items had failures: 1 of 36 in __main__.example_34 ***Test Failed*** 1 failures. For whitespace errors, see the file /opt/sage-devel/tmp/.doctest_multi_polynomial_ideal.py [12.8 s] exit code: 1024 ---------------------------------------------------------------------- The following tests failed: sage -t "devel/sage-main/sage/rings/polynomial/multi_polynomial_ideal.py" Total time for all tests: 12.8 seconds
comment:5 Changed 12 years ago by
- Milestone changed from sage-3.4.2 to sage-4.0
Wrong milestone - better luck in Sage 4.0.
Cheers,
Michael
comment:6 Changed 12 years ago by
I fixed the "proofable" vs. "provable" typo. However, I am a bit clueless how to deal with the numerical noise. The ...
works for most of it but not for 0.49999999999
vs. {{{0.500000000}}. Ideas?
comment:7 Changed 12 years ago by
Martin,
Sorry for the late reply, but according to Lazard's paper, the algorithm I used ("Triangular") does NOT need a lexicographic term ordering. On p. 124, beginning of Section 6:
These algorithms do not depend explicitly on the ordering; however, they are mainly designed for degree orderings for which the base is more easily obtained; they work also for lexicographic ordering, but, for them, we dispose of a structure theorem which permit us to provide a better algorithm (section 8).
The algorithm I implemented was Triangular (line 6). So the switch the lexicographic is unnecessary.
I don't know how to fix the problem with complex coefficients, though.
comment:8 Changed 12 years ago by
Incidentally, I'll look at the problem with the missing x as soon as I finish downloading & compiling the latest Sage. (4.1.1?)
comment:9 Changed 12 years ago by
I have found two bugs. One of them is not solvable from my end, and possibly not at all (others who know more should comment).
The first one: elim_pol is not always computing the correct polynomial. After removing the switch to lexicographic order, change line 358 (?) of toy_variety.py from
for each in xrange(len(coeffs)-1):
to
for each in xrange(len(coeffs)):
I have no idea why I had it count until the next-to-last element of coeffs; all of them are necessary. I can submit a patch if you like.
However: even after this fix, a problem remains: line 292 (?) triangular_factorization tries to compute a Groebner basis of an ideal whose generators should have a common solution. This is the source of the 1.0 appearing in the triangular variety. Unfortunately, the computed basis is 1.0, suggesting that the ideal has no common solution! I have an idea why this is happening, but I can't yet say for sure.
comment:10 Changed 12 years ago by
To follow up, the problem occurs when computing the Groebner basis over CC.
I'll use the example given:
The Groebner basis computed is [y^{3 + x - y, x}2 + y^{2 - 1.00000000000000, x*y - 1.00000000000000]. }
The result p from elim_pol is y^{4 - y}2 + 1.0. (This reflects the bugfix I identified above; it used to return y^{4-y}2.)
The first factor q of p is y - 0.866025403784439 - 0.500000000000000*I.
The reduction of B modulo q gives us
[x - 0.866025403784438 + 0.500000000000001*I,
x^{2 - 0.499999999999999 + 0.866025403784439*I, (0.866025403784439 + 0.500000000000000*I)*x - 1.00000000000000] }
This should be a consistent system: the first polynomial is a factor of the second, and the solution to the third is nearly the same as the solution to the first: 0.866025403784438 - 0.500000000000001*I vs. 0.866025403784438 - 0.500000000000000*I. This appears to be a roundoff/floating point error.
The system above should produce a Groebner basis with one polynomial, but it returns [1.00000000000000] instead. This is why nothing is coming back for x. Anyone know how to fix it?
comment:11 Changed 12 years ago by
Sorry for the repost, but I had some superscript typos in the previou.
To follow up, the problem occurs when computing the Groebner basis over CC.
I'll use the example given:
The Groebner basis computed is [y3 + x - y, x2 + y2 - 1.00000000000000, x*y - 1.00000000000000].
The result p from elim_pol is y4 - y2 + 1.0. (This reflects the bugfix I identified above; it used to return y4-y2.)
The first factor q of p is y - 0.866025403784439 - 0.500000000000000*I.
The reduction of B modulo q gives us
[x - 0.866025403784438 + 0.500000000000001*I,
x2 - 0.499999999999999 + 0.866025403784439*I, (0.866025403784439 + 0.500000000000000*I)*x - 1.00000000000000]
This should be a consistent system: the first polynomial is a factor of the second, and the solution to the third is nearly the same as the solution to the first: 0.866025403784438 - 0.500000000000001*I vs. 0.866025403784438 - 0.500000000000000*I. This appears to be a roundoff/floating point error.
The system above should produce a Groebner basis with one polynomial, but it returns [1.00000000000000] instead. This is why nothing is coming back for x. Anyone know how to fix it?
comment:12 Changed 12 years ago by
Well, we should probably just disallow computing a GB over CC anyway, it doesn't really make sense because CC is not realy a field.
comment:13 Changed 12 years ago by
You mean CC as implemented, not CC in theory?
I could insert a test for whether the field is CC, and if so raise an exception. That said, the roots we're looking for are algebraic, and so can be described symbolically; i.e., without floating point. So if the polynomial starts from QQ, one could in theory construct the roots. Should the exception just reject the user's input, advise the user to try an extension field, or try to construct one itself?
Is this something I should ask about on sage-devel?
comment:14 Changed 12 years ago by
Yes, I mean fixed precision floating point numbers. We could just print a warning? Even if the solutions are algebraic the rounding errors can cause a zero to look like a one.
Changed 12 years ago by
comment:15 Changed 12 years ago by
The uploaded patch should incorporate Martin's changes and should also provide a warning if the field is CC. I'd still like to find a way to make it work in CC, though.
comment:16 Changed 12 years ago by
The patch looks good, applies cleanly against 4.1.1 and doctests pass on 64-bit Linux (sage.math). There might be a numerical noise issue (see above) which needs to be addressed.
comment:17 Changed 12 years ago by
Ah, the doctest needs a ::
prepended otherwise, it won't be run.
comment:18 Changed 12 years ago by
which doctest? are you talking about line 1848 of multi_polynomial_ideal.py?
comment:19 Changed 12 years ago by
Yes, the doctest starting at 1851.
comment:20 Changed 12 years ago by
Okay, double colons on that file. They aren't prepended to the example, though; they appear at the end of line 1849 in the same way as at the end of 1831.
PLEASE tell me this is okay now. Fixing this drove me crazy with hg...
comment:21 Changed 12 years ago by
Yes, thats correct this way. However, I expect the doctest will now fail on 32-bit systems (it wasn't run before).
comment:22 Changed 12 years ago by
I believe that my system is running as a 32-bit system. I know that Linux is running as 32-bit; does sage run differently?
comment:23 Changed 12 years ago by
No, if your Linux is 32-bit then your Sage is 32-bit :)
comment:24 Changed 12 years ago by
still no positive review then? :-( what's still needed?
comment:25 Changed 12 years ago by
- Summary changed from [with patch, needs work] MPolynomial_polydict.factor() should accept proof parameter to [with patch, positive review] MPolynomial_polydict.factor() should accept proof parameter
Sorry, here we go.
comment:26 Changed 12 years ago by
Which patches should be applied and in what order? The patch variety_CC.patch
contains all changes made in mpoly_factor_proof.patch
. Does that mean mpoly_factor_proof.patch
can be left out?
comment:27 Changed 12 years ago by
Indeed, only variety_CC.patch
and variety_CC2.patch
should be applied.
comment:28 Changed 12 years ago by
- Merged in set to Sage 4.1.2.alpha0
- Resolution set to fixed
- Reviewers set to Martin Albrecht
- Status changed from new to closed
John -- Please remember to put a sensible commit message in your patch. It's recommended that you reference the ticket number in the commit message. Merged patches in this order:
variety_CC.patch
variety_CC2.patch
The patch fixes the exception, however:
which is certainly the wrong answer (
x
is missing), thus there seems to be another bug.