Opened 14 years ago
Closed 14 years ago
#4965 closed defect (fixed)
[with patch, positive review] Z/nZ[x] via FLINT's zmod_poly
Reported by: | malb | Owned by: | somebody |
---|---|---|---|
Priority: | major | Milestone: | sage-3.3 |
Component: | basic arithmetic | Keywords: | |
Cc: | burcin, was | Merged in: | |
Authors: | Reviewers: | ||
Report Upstream: | Work issues: | ||
Branch: | Commit: | ||
Dependencies: | Stopgaps: |
Description
The attached patches wrap FLINT's Z/nZ[x] for word sized n. This provides a considerable speed-up (~ 20x) for univariate polynomial arithmetic over these rings and also cleans up the code using the polynomial_template.pxi
mechanism. For that polynomial_template.pxi
was also cleaned-up which should make it more suitable for other backend implementations.
Attachments (4)
Change History (44)
Changed 14 years ago by
comment:1 Changed 14 years ago by
Speed-Up on sage.math
Old
sage: P.<x> = PolynomialRing(GF(7)) sage: type(x) <type 'sage.rings.polynomial.polynomial_modn_dense_ntl.Polynomial_dense_mod_p'> sage: f = P.random_element(100) sage: g = P.random_element(100) sage: %timeit f*g 1000 loops, best of 3: 445 µs per loop
New
sage: P.<x> = PolynomialRing(GF(7)) sage: type(x) <type 'sage.rings.polynomial.polynomial_zmod_flint.Polynomial_zmod_flint'> sage: f = P.random_element(100) sage: g = P.random_element(100) sage: %timeit f*g 100000 loops, best of 3: 7.92 µs per loop
comment:2 Changed 14 years ago by
At least these two doctests fail.
sage -t "devel/sage/sage/schemes/hyperelliptic_curves/hyperelliptic_padic_field.py" sage -t "devel/sage/sage/schemes/elliptic_curves/ell_number_field.py"
I would greatly appreciate input/pointers where to look for the bug. Once I figured out which of the functions I've touched returns wrong/unexpected answers, I can fix it. However, as of now, I don't know which one it is.
comment:3 Changed 14 years ago by
More details:
sage -t ".../sage/schemes/hyperelliptic_curves/hyperelliptic_padic_field.py" ********************************************************************** File ".../sage/schemes/hyperelliptic_curves/hyperelliptic_padic_field.py", line 198: sage: w.coleman_integral(P, Q) Expected: O(11^6) Got: 6 + 7*11 + 8*11^2 + 7*11^3 + 3*11^4 + 10*11^5 + O(11^6) ********************************************************************** File ".../sage/schemes/hyperelliptic_curves/hyperelliptic_padic_field.py", line 200: sage: C.coleman_integral(x*w, P, Q) Expected: O(11^6) Got: 3 + 11 + 6*11^2 + 2*11^3 + 8*11^5 + O(11^6) ********************************************************************** File ".../sage/schemes/hyperelliptic_curves/hyperelliptic_padic_field.py", line 202: sage: C.coleman_integral(x^2*w, P, Q) Expected: 3*11 + 2*11^2 + 7*11^3 + 2*11^4 + 10*11^5 + O(11^6) Got: 5 + 11 + 8*11^2 + 11^3 + 7*11^4 + 6*11^5 + O(11^6) ********************************************************************** File ".../sage/schemes/hyperelliptic_curves/hyperelliptic_padic_field.py", line 213: sage: w.integrate(P, Q), (x*w).integrate(P, Q) Expected: (O(71^4), O(71^4)) Got: (40 + 4*71 + 21*71^2 + 39*71^3 + O(71^4), 2 + 58*71 + 27*71^2 + 49*71^3 + O(71^4)) ********************************************************************** File ".../sage/schemes/hyperelliptic_curves/hyperelliptic_padic_field.py", line 216: sage: w.integrate(P, R) Expected: 42*71 + 63*71^2 + 55*71^3 + O(71^4) Got: 14 + 27*71 + 46*71^2 + 10*71^3 + O(71^4)
comment:4 Changed 14 years ago by
File ".../sage/schemes/elliptic_curves/ell_number_field.py", line 728: sage: EllipticCurve([1 + a , -1 + a , 1 + a , -11 + a , 5 -9*a ]).conductor() Exception raised: Traceback (most recent call last): File ".../sage/schemes/elliptic_curves/ell_number_field.py", line 747, in conductor for d in self.local_data()],\ File ".../sage/schemes/elliptic_curves/ell_number_field.py", line 396, in local_data return [self._get_local_data(pr, proof) for pr in primes] File ".../sage/schemes/elliptic_curves/ell_number_field.py", line 418, in _get_local_data self._local_data[P] = EllipticCurveLocalData(self, P, proof) File ".../schemes/elliptic_curves/ell_local_data.py", line 140, in __init__ self._Emin, ch, self._val_disc, self._fp, self._KS, self._cp, self._split = self._tate(proof) File ".../sage/schemes/elliptic_curves/ell_local_data.py", line 520, in _tate r = proot(-b6, 3) File ".../sage/schemes/elliptic_curves/ell_local_data.py", line 451, in <lambda> proot = lambda x,e: F.lift(F(x).nth_root(e, extend = False, all = True)[0]) File "integer_mod.pyx", line 855, in sage.rings.integer_mod.IntegerMod_abstract.nth_root (sage/rings/integer_mod.c:7712) File "polynomial_element.pyx", line 3715, in sage.rings.polynomial.polynomial_element.Polynomial.roots (sage/rings/polynomial/polynomial_element.c:24924) File "polynomial_element.pyx", line 2096, in sage.rings.polynomial.polynomial_element.Polynomial.factor (sage/rings/polynomial/polynomial_element.c:15667) ValueError: factorization of 0 not defined
comment:5 follow-up: ↓ 6 Changed 14 years ago by
I'll try to look at the second of those since it seems to crop up in code which I wrote. But I em extremely busy at the moment so cannot promise anything very quickly.
comment:6 in reply to: ↑ 5 Changed 14 years ago by
Replying to cremona:
I'll try to look at the second of those since it seems to crop up in code which I wrote. But I em extremely busy at the moment so cannot promise anything very quickly.
Thanks!
comment:7 Changed 14 years ago by
- Cc burcin added
comment:8 Changed 14 years ago by
The updated patch
- addresses a couple of remarks by Bill Hart off-list
- changes the celement definition (following Burcin's alternative implementation)
Burcin also wrapped zmod_poly_t
two months ago and our patches match a lot. Thus, copyright should be shared. The only reason I'm using my patches as a basis is because I know them better. Btw. Burcin's patches rearrange a fair amount of stuff in polynomial_ring.pyx
which I don't. He also splits the polynomial_template.pxi
in polynomial_template.pxi
and polynomial_template_field.pxi
.
comment:9 Changed 14 years ago by
"Long" Doctest Failures
this one is easy:
sage -t --long "devel/sage/sage/schemes/elliptic_curves/padics.py" ... sage.rings.polynomial.polynomial_zmod_flint.Polynomial_zmod_flint._derivative ... sage.libs.ntl.ntl_lzz_pContext.ntl_zz_pContext ValueError: Modulus (=2384185791015625) is too big
File ".../sage/schemes/elliptic_curves/ell_number_field.py", line 728: sage: EllipticCurve([1 + a , -1 + a , 1 + a , -11 + a , 5 -9*a ]).conductor() ... ValueError: factorization of 0 not defined
sage -t --long "devel/sage/sage/schemes/elliptic_curves/sha_tate.py" ... sage: EllipticCurve('858k2').sha().an_padic(7) #rank 0, non trivial sha (long ... sage.libs.ntl.ntl_lzz_pContext.ntl_zz_pContext ValueError: Modulus (=558545864083284007) is too big
sage -t --long "devel/sage/sage/schemes/hyperelliptic_curves/hyperelliptic_padic_field.py" ********************************************************************** sage: w.coleman_integral(P, Q) Expected: O(11^6) Got: 6 + 7*11 + 8*11^2 + 7*11^3 + 3*11^4 + 10*11^5 + O(11^6) ********************************************************************** sage: C.coleman_integral(x*w, P, Q) Expected: O(11^6) Got: 3 + 11 + 6*11^2 + 2*11^3 + 8*11^5 + O(11^6) ********************************************************************** sage: C.coleman_integral(x^2*w, P, Q) Expected: 3*11 + 2*11^2 + 7*11^3 + 2*11^4 + 10*11^5 + O(11^6) Got: 5 + 11 + 8*11^2 + 11^3 + 7*11^4 + 6*11^5 + O(11^6) ********************************************************************** sage: w.integrate(P, Q), (x*w).integrate(P, Q) Expected: (O(71^4), O(71^4)) Got: (40 + 4*71 + 21*71^2 + 39*71^3 + O(71^4), 2 + 58*71 + 27*71^2 + 49*71^3 + O(71^4)) ********************************************************************** sage: w.integrate(P, R) Expected: 42*71 + 63*71^2 + 55*71^3 + O(71^4) Got: 14 + 27*71 + 46*71^2 + 10*71^3 + O(71^4) **********************************************************************
comment:10 Changed 14 years ago by
- Milestone changed from sage-3.4.1 to sage-3.4
The updated zmod_poly.patch
fixes the two
ValueError: Modulus (...) is too big
errors. So we're back to square one and these two fail:
sage -t "devel/sage/sage/schemes/hyperelliptic_curves/hyperelliptic_padic_field.py" sage -t "devel/sage/sage/schemes/elliptic_curves/ell_number_field.py"
comment:11 Changed 14 years ago by
The patch update adds a simple zn_poly interface to Polynomial_zmod_flint
.
Performance
def f(p,n): P = PolynomialRing(GF(next_prime(p)),'x') f = P.random_element(n) g = P.random_element(n) t0 = cputime() r0 = f*g t0 = cputime(t0) t1 = cputime() r1 = f._mul_zn_poly(g) t1 = cputime(t1) assert(r0 == r1) return p,n,t0,t1 for i in range(21): f(2**47,2**i)
returns
# (140737488355328, 1, 0.0, 0.0) # (140737488355328, 2, 0.0, 0.0) # (140737488355328, 4, 0.00099999999999766942, 0.0) # (140737488355328, 8, 0.0, 0.0) # (140737488355328, 16, 0.0, 0.0) # (140737488355328, 32, 0.0059990000000027521, 0.0) # (140737488355328, 64, 0.0, 0.0) # (140737488355328, 128, 0.0, 0.0) # (140737488355328, 256, 0.0, 0.0) # (140737488355328, 512, 0.0, 0.00099999999999766942) # (140737488355328, 1024, 0.00099999999999766942, 0.0) # (140737488355328, 2048, 0.0020000000000024443, 0.0019989999999978636) # (140737488355328, 4096, 0.0049989999999979773, 0.005000000000002558) # (140737488355328, 8192, 0.010998000000000729, 0.011997999999998399) # (140737488355328, 16384, 0.023995999999996798, 0.023997000000001378) # (140737488355328, 32768, 0.050992000000000814, 0.052991999999996153) # (140737488355328, 65536, 0.1149820000000048, 0.10598499999999689) # (140737488355328, 131072, 0.29195599999999189, 0.21996599999999944) # (140737488355328, 262144, 0.6119070000000022, 0.45393199999999467) # (140737488355328, 524288, 1.5217689999999919, 1.0278430000000043) # (140737488355328, 1048576, 3.1365230000000111, 2.0966819999999871)
comment:12 Changed 14 years ago by
Here's what is causing the problem in the ell_number_field.py test:
If F is a field of 3 elements, b=F(0), then b.nth_root(0) gives an error:
sage: K.<a>=NumberField(x^2-x+3) sage: F=K.residue_field(a); F; F.order() Residue field of Fractional ideal (a) 3 sage: b=F(0); b 0 sage: b.nth_root(3) ValueError: factorization of 0 not defined
The nth_root() function being called here is at fault, but that's as far as I got.
comment:13 Changed 14 years ago by
PS To the above: The nth_root() function should definitely catch x=0 as a special trivial case. But here that would just hide the problem, that calling roots() on the polynomial x^3
over the field F causes a crash. I have not time to dig deeper: I would expect that first one would take gcd(x^3,x^3-x)
...
comment:14 Changed 14 years ago by
Thanks a lot John! I now have a function to debug which should be sufficient to track down the bug.
comment:15 Changed 14 years ago by
Thanks to John's comment we're now down to one doctest failure:
sage -t --long ...sage/schemes/hyperelliptic_curves/hyperelliptic_padic_field.py
********************************************************************** File ".../sage/schemes/hyperelliptic_curves/hyperelliptic_padic_field.py", line 198: sage: w.coleman_integral(P, Q) Expected: O(11^6) Got: 6 + 7*11 + 8*11^2 + 7*11^3 + 3*11^4 + 10*11^5 + O(11^6) ********************************************************************** File ".../sage/schemes/hyperelliptic_curves/hyperelliptic_padic_field.py", line 200: sage: C.coleman_integral(x*w, P, Q) Expected: O(11^6) Got: 3 + 11 + 6*11^2 + 2*11^3 + 8*11^5 + O(11^6) ********************************************************************** File ".../sage/schemes/hyperelliptic_curves/hyperelliptic_padic_field.py", line 202: sage: C.coleman_integral(x^2*w, P, Q) Expected: 3*11 + 2*11^2 + 7*11^3 + 2*11^4 + 10*11^5 + O(11^6) Got: 5 + 11 + 8*11^2 + 11^3 + 7*11^4 + 6*11^5 + O(11^6) ********************************************************************** File ".../sage/schemes/hyperelliptic_curves/hyperelliptic_padic_field.py", line 213: sage: w.integrate(P, Q), (x*w).integrate(P, Q) Expected: (O(71^4), O(71^4)) Got: (40 + 4*71 + 21*71^2 + 39*71^3 + O(71^4), 2 + 58*71 + 27*71^2 + 49*71^3 + O(71^4)) ********************************************************************** File ".../sage/schemes/hyperelliptic_curves/hyperelliptic_padic_field.py", line 216: sage: w.integrate(P, R) Expected: 42*71 + 63*71^2 + 55*71^3 + O(71^4) Got: 14 + 27*71 + 46*71^2 + 10*71^3 + O(71^4)
comment:16 Changed 14 years ago by
I'm glad that this bug made me look at that code, since what I found was that the only situations where I was calling that nth_root() function were for square roots in residue fields of characteristic 2 and cube roots in char. 3. In both cases there is a special p'th root function implemented by jhpalmieri in #4553 which I reviewed and which could be used instead, which should be quicker. So this bug will have had a useful side-effect or two...
comment:17 follow-ups: ↓ 21 ↓ 22 Changed 14 years ago by
Some comments:
1) Is there any reason len was removed from GF(2)[x]? If we're going to take it out, it should at lest be deprecated, but I think it's useful to have. 2) It is preferable to use .pxd files rather than .pxi files.
I'm still looking into why the doctest failure above, and nothing obvious is coming to mind. Falling back to ntl polynomials gives the same wrong answer, so it must be something different in the way polynomials are interfaced/used.
Changed 14 years ago by
comment:18 follow-up: ↓ 19 Changed 14 years ago by
See attached.
Perhaps __rshift__
and __lshift__
should call shift, it seems redundant to have code in all three. Also, I bet ntl and flint can shift faster than a multiply and/or divide.
comment:19 in reply to: ↑ 18 Changed 14 years ago by
Perhaps
__rshift__
and__lshift__
should call shift, it seems redundant to have code in all three. Also, I bet ntl and flint can shift faster than a multiply and/or divide.
I'd expect them to detect shifts, but I might be wrong. Should there be a celement_rshift
/celement_lshift
?
comment:20 Changed 14 years ago by
Malb, my wrapper provides support for the shifts, you can just copy it from there.
I think keeping a generic implementation in the Polynomial_template class, and overriding it if there is support for faster methods is good enough. Though maybe separating the user interface and functionality into 2 different functions (e.g. shift_left(), _shift_left_c()) might be better, so that the user interface and docstrings are consistent for subclasses.
comment:21 in reply to: ↑ 17 Changed 14 years ago by
Replying to robertwb:
Some comments:
1) Is there any reason len was removed from GF(2)[x]? If we're going to take it out, it should at lest be deprecated, but I think it's useful to have.
len shouldn't have been there in the first place because it is inconsistent with the rest of the univariate polynomials.
2) It is preferable to use .pxd files rather than .pxi files.
Okay, I'll change that then.
comment:22 in reply to: ↑ 17 Changed 14 years ago by
Replying to robertwb:
I'm still looking into why the doctest failure above, and nothing obvious is coming to mind. Falling back to ntl polynomials gives the same wrong answer, so it must be something different in the way polynomials are interfaced/used.
I can't reproduce the behaviour you're describing. Maybe you missed that there are two places where Polynomial_zmod_poly
is constructed (modn and modp). It fails if I leave the modn stuff in.
comment:23 Changed 14 years ago by
it seems like Robert's patch fixes the issue (I just realized now). I'm running -t --long
now to see if everything is alright. Then I'll switch stuff around to pxd and do the shift business.
comment:24 Changed 14 years ago by
malb@sage:~/sage-3.2.2/devel/sage$ sage -tp 10 --long sage/ ... All tests passed!
comment:25 Changed 14 years ago by
malb: Yes, I was missing the fact that I had to change modn and modp, but in any case I found the bug in the end.
Regarding len, I think it would be useful to have for all univariate polynomials, but it's not a big deal.
Regarding shifts, I think we already have too many functions, and shouldn't be introducing even more. I think there should be a celement_shift (which takes positive and negative values), and the template __lshift__, __rshift__
, and shift
all call this. It's probably the exception to not have specialized shift methods, and in this case one would manually implement them in celement_shift.
comment:26 Changed 14 years ago by
I've made the shifting issue a new ticket: #4982
comment:27 Changed 14 years ago by
So I need to switch from pxi to pxd for
polynomial_template.pxi
ntl_GF2X_linkage.pxi
flint_zmod_poly_linkage.pxi
and the ticket is ready for review. Or am I missing something?
comment:28 Changed 14 years ago by
Sorry for not being clear. The files you listed:
- polynomial_template.pxi
- flint_zmod_poly_linkage.pxi
are the ones that should be .pxi (as they're actually included). The files that should be .pxd are
- flint.pxi
- ntl_interface.pxi
- zmod_poly.pxi
as they are declaration files. An example is worth a thousand words, I'll attach a patch (which you can fold into yours if you want).
comment:29 Changed 14 years ago by
OK, I attached a patch. This should probably be folded into zmod_poly.patch
comment:30 Changed 14 years ago by
I had a lock and the direct use of ling in the FLINT headers is a bad thing on LLP platforms, i.e. Solaris, for at least efficiency reasons. I will complain about this to Bill and see what happens, but I guess this should not prevent this patch from going in for now.
Cheers,
Michael
comment:31 Changed 14 years ago by
- Summary changed from [with patch, needs work] Z/nZ[x] via FLINT's zmod_poly to [with patch, needs review] Z/nZ[x] via FLINT's zmod_poly
to apply:
polynomial_template.patch
zmod_poly.patch
comment:32 Changed 14 years ago by
- Milestone changed from sage-3.4 to sage-3.3
For the record: Doctests on sage.math pass.
Robert: If humanly possible can you review this in the next couple hours? Then this patch will make it into 3.3.alpha0.
Cheers,
Michael
comment:33 Changed 14 years ago by
The code looks good, and seems to work for me.
One failure in the doctests:
sage -t "sage/rings/polynomial//polynomial_zmod_flint.pyx" ********************************************************************** File "/Users/robert/sage/sage-3.1.3/devel/sage-trac4965/sage/rings/polynomial/polynomial_zmod_flint.pyx", line 236: sage: f*g == f._mul_zn_poly(g) Exception raised: Traceback (most recent call last): File "/Users/robert/sage/current/local/bin/ncadoctest.py", line 1231, in run_one_test self.run_one_example(test, example, filename, compileflags) File "/Users/robert/sage/current/local/bin/sagedoctest.py", line 38, in run_one_example OrigDocTestRunner.run_one_example(self, test, example, filename, compileflags) File "/Users/robert/sage/current/local/bin/ncadoctest.py", line 1172, in run_one_example compileflags, 1) in test.globs File "<doctest __main__.example_8[5]>", line 1, in <module> f*g == f._mul_zn_poly(g)###line 236: sage: f*g == f._mul_zn_poly(g) AttributeError: 'sage.rings.polynomial.polynomial_modn_dense_ntl.Po' object has no attribute '_mul_zn_poly' ********************************************************************** 1 items had failures: 1 of 10 in __main__.example_8 ***Test Failed*** 1 failures. For whitespace errors, see the file /Users/robert/sage/current/tmp/.doctest_polynomial_zmod_flint.py
should be an easy fix.
comment:34 Changed 14 years ago by
I cannot reproduce your failure neither could mabshoff (from the above statement). Maybe there was a problem applying the patches?
comment:35 Changed 14 years ago by
I bet I have a different cutoff between flint and ntl for a 32-bit machine.
comment:36 Changed 14 years ago by
- Cc was added
You are of course right. I changed the doctest to use next_prime(GF(2^31))
which should work on 32-bit systems. However, this raises the question whether we want to be consistent across all platforms for the cutoff: w.r.t. Matrix_modn_dense
mabshoff and wstein indicated their desire to do so in another ticket.
comment:37 follow-up: ↓ 38 Changed 14 years ago by
Looks like you have to go down to 2^{30. }
sage: P.<x> = PolynomialRing(GF(next_prime(2^31))) sage: type(x) <type 'sage.rings.polynomial.polynomial_modn_dense_ntl.Polynomial_dense_mod_p'> sage: P.<x> = PolynomialRing(GF(next_prime(2^30))) sage: type(x) <type 'sage.rings.polynomial.polynomial_zmod_flint.Polynomial_zmod_flint'>
Personally, I think it's a bad idea to impose a 32-b it cutoff for 64-bit platforms.
comment:38 in reply to: ↑ 37 Changed 14 years ago by
comment:39 Changed 14 years ago by
- Summary changed from [with patch, needs review] Z/nZ[x] via FLINT's zmod_poly to [with patch, positive review] Z/nZ[x] via FLINT's zmod_poly
That resolved all of my concerns, positive review.
comment:40 Changed 14 years ago by
- Resolution set to fixed
- Status changed from new to closed
Merged polynomial_template.patch and zmod_poly.patch in Sage 3.3.alpha1
apply this first