Ticket #11782 (closed defect: fixed)

Opened 22 months ago

Last modified 19 months ago

Bug in discriminant of polynomials over Z/nZ with n composite

Reported by: johanbosman Owned by: AlexGhitza
Priority: major Milestone: sage-4.8
Component: algebra Keywords: polynomial, discriminant, integer mod ring
Cc: Work issues:
Report Upstream: N/A Reviewers: Julian Rueth
Authors: Johan Bosman Merged in: sage-4.8.alpha4
Dependencies: Stopgaps:

Description (last modified by saraedum) (diff)

The following behaviour is inconsistent:

sage: f = ZZ[x](2*x^3 + x^2 + x)
sage: f.discriminant()
-7
sage: GF(3)[x](f).discriminant()
2
sage: ZZ.quo(9)[x](f).discriminant()
0

And the following raises an error:

sage: ZZ.quo(9)[x](3*x^2 + 3*x + 3).discriminant()
...
ZeroDivisionError: Inverse does not exist.

Apply

  1. 11782_use_sylvester.patch Download
  2. 11782_review.patch Download

to the Sage library.

Attachments

11782_use_sylvester.patch Download (2.5 KB) - added by johanbosman 22 months ago.
11782_review.patch Download (2.5 KB) - added by saraedum 19 months ago.
review patch

Change History

comment:1 Changed 22 months ago by johanbosman

The first problem is caused by the fact that FLINT doesn't compute resultants correctly. This is also mentioned in FLINT's documentation: see p. 58 of  http://www.flintlib.org/flint-1.6.pdf: n must be prime to compute resultants correctly.

The second problem is of course caused by the fact that the leading coefficient is a zero divisor.

Changed 22 months ago by johanbosman

comment:2 Changed 22 months ago by johanbosman

  • Status changed from new to needs_review

comment:3 Changed 19 months ago by saraedum

  • Status changed from needs_review to needs_info

The docstring of discriminant() says

           Uses the identity `R_n(f) := (-1)^(n (n-1)/2) R(f,
           f') a_n^(n-k-2)`, where `n` is the degree of self,
           `a_n` is the leading coefficient of self, `f'`
           is the derivative of `f`, and `k` is the degree
           of `f'`.

So is it correct that if a_n is not invertible this should read

a_n^(2+k-n) R_n(f) = (-1)^(n (n-1)/2) R(f,f')

If I understand your patch correctly, you're assuming 2+k-n=1 in the lines

            mat[0, 0] = self.base_ring()(1)
            mat[n - 1, 0] = self.base_ring()(n)

since you're only dividing by a_n. But k could be different from n-1.

I don't know much about discriminants/resultants, especially if the polynomials are defined over rings with zerodivisors. So maybe I'm just misunderstanding the definition/implementation here.

comment:4 Changed 19 months ago by johanbosman

  • Status changed from needs_info to needs_review
  • Authors set to Johan Bosman

That code is reached whenever 

an = self[n]**(n - k - 2)

raises a ZeroDivisionError?, which means n-k-2 < 0 thus n-k <= 1 (and thus n-1 <= k). Now, since k is the degree of the derivative of f and n is the degree of f, we have k <= n-1. This two inequalities together yield k = n-1. ;).

comment:5 Changed 19 months ago by johanbosman

"This" => "These" in the last sentence (is it possible to edit your own comments?)

comment:6 Changed 19 months ago by saraedum

I see. I had assumed that the catch block was a full alternative implementation of the try block. You're of course right, it doesn't have to be.

comment:7 Changed 19 months ago by saraedum

  • Reviewers set to Julian Rueth
  • Description modified (diff)

I'm a bit unhappy with catching a ZeroDivisionError. The error could be raised in the self.resultant(d) call and would then be silently ignored — and the result could be wrong in that case. So I made the check more explicit in my review patch.

Also I added some comments and an explicit test to the flint polynomials.

If you don't mind these changes, feel free to set it to positive review.

comment:8 Changed 19 months ago by johanbosman

  • Status changed from needs_review to needs_work

If a_n is not a unit, but also not a zero divisor, it still makes sense to invert it, so your "if" condition is too strict.

What do you think about the following solution?

try:         
    an = self[n]**(n - k - 2) 
except ZeroDivisionError: 
    mat = self.sylvester_matrix(d)            
    mat[0, 0] = self.base_ring()(1)
    mat[n - 1, 0] = self.base_ring()(n)
    return u * mat.determinant()
else:
    return self.base_ring()(u * self.resultant(d) * an)

Changed 19 months ago by saraedum

review patch

comment:9 Changed 19 months ago by saraedum

  • Status changed from needs_work to needs_review

You're right. I tried to follow your suggestion in the new review patch.

comment:10 Changed 19 months ago by johanbosman

  • Status changed from needs_review to positive_review

Okäse!

comment:11 Changed 19 months ago by jdemeyer

  • Status changed from positive_review to closed
  • Resolution set to fixed
  • Merged in set to sage-4.8.alpha4
Note: See TracTickets for help on using tickets.