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:

Status badges

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)

toy_variety_fix.patch (3.0 KB) - added by malb 12 years ago.
mpoly_factor_proof.patch (1.1 KB) - added by malb 12 years ago.
fix typo
variety_CC.patch (3.4 KB) - added by john_perry 12 years ago.
variety_CC2.patch (881 bytes) - added by john_perry 12 years ago.
apply over variety_CC.patch

Download all attachments as: .zip

Change History (32)

comment:1 Changed 12 years ago by malb

The patch fixes the exception, however:

sage: sage: R.<x,y> = CC[]
sage: sage: I = R.ideal(x^2+y^2-1,x*y-1)
sage: sage: I.variety()
verbose 0 (1735: multi_polynomial_ideal.py, variety) Warning: falling back to very slow toy implementation.
[{y: -1.00000000000000}, {y: 0}, {y: 1.00000000000000}]

which is certainly the wrong answer (x is missing), thus there seems to be another bug.

comment:2 Changed 12 years ago by malb

It seems toy_variety does not switch to 'lex' when it should?

Changed 12 years ago by malb

comment:3 Changed 12 years ago by malb

  • 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 AlexGhitza

  • 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 mabshoff

  • Milestone changed from sage-3.4.2 to sage-4.0

Wrong milestone - better luck in Sage 4.0.

Cheers,

Michael

Changed 12 years ago by malb

fix typo

comment:6 Changed 12 years ago by malb

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 john_perry

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 john_perry

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 john_perry

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 john_perry

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:11 Changed 12 years ago by john_perry

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 malb

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 john_perry

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 malb

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 john_perry

comment:15 Changed 12 years ago by john_perry

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 malb

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 malb

Ah, the doctest needs a :: prepended otherwise, it won't be run.

comment:18 Changed 12 years ago by john_perry

which doctest? are you talking about line 1848 of multi_polynomial_ideal.py?

comment:19 Changed 12 years ago by malb

Yes, the doctest starting at 1851.

Changed 12 years ago by john_perry

apply over variety_CC.patch

comment:20 Changed 12 years ago by john_perry

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 malb

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 john_perry

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 malb

No, if your Linux is 32-bit then your Sage is 32-bit :)

comment:24 Changed 12 years ago by john_perry

still no positive review then? :-( what's still needed?

comment:25 Changed 12 years ago by malb

  • 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 mvngu

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 malb

Indeed, only variety_CC.patch and variety_CC2.patch should be applied.

comment:28 Changed 12 years ago by mvngu

  • Authors set to John Perry
  • 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:

  1. variety_CC.patch
  2. variety_CC2.patch
Note: See TracTickets for help on using tickets.