Opened 5 years ago
Closed 2 years ago
#16516 closed enhancement (fixed)
Faster roots computation for sparse polynomials over ZZ
Reported by:  bruno  Owned by:  

Priority:  major  Milestone:  sage8.0 
Component:  commutative algebra  Keywords:  roots, sparse polynomial 
Cc:  Merged in:  
Authors:  Bruno Grenet  Reviewers:  Vincent Delecroix, Travis Scrimshaw, Jeroen Demeyer 
Report Upstream:  N/A  Work issues:  
Branch:  686d2b3 (Commits)  Commit:  686d2b366a2726e78de271e379f0e49209fb7c8e 
Dependencies:  Stopgaps: 
Description (last modified by )
Algorithms exist for the computation of the roots of sparse polynomials, which are much faster than the "generic" algorithms which work for dense polynomials.
In this ticket, I implement one of these algorithms, for integer roots of sparse polynomials over ZZ.
Attachments (5)
Change History (74)
comment:1 Changed 5 years ago by
 Description modified (diff)
 Summary changed from Faster roots computation for sparse polynomials over ZZ to Faster roots computation for sparse polynomials over ZZ and QQ
comment:2 Changed 5 years ago by
 Branch set to u/bruno/faster_roots_computation_for_sparse_polynomials_over_zz
comment:3 Changed 5 years ago by
 Commit set to 920c16e8952cb1f885c1fd1f3bcdfd5d572c17fe
 Description modified (diff)
 Status changed from new to needs_review
 Summary changed from Faster roots computation for sparse polynomials over ZZ and QQ to Faster roots computation for sparse polynomials over ZZ
comment:4 Changed 5 years ago by
 Commit changed from 920c16e8952cb1f885c1fd1f3bcdfd5d572c17fe to e725544b76363ea517a262419a092f9414c1b4a4
comment:5 Changed 5 years ago by
Hi Bruno,
It seems reasonable to have it in integer_ring.pyx
. And I think that it would be good to allow it for dense polynomial as well. You certainly know what is the good measure of sparseness... you could implement a method .density()
of .sparness_XYZ()
on polynomials that would help you to find the good threshold for the different algorithms.
For better code, please:
 remove the trailing whitespaces
 there is no need for
;
at the end of each line  for list update it is better to use
.append(x)
for adding one element and.extend(iterator)
for extending with several elements not ring is None
>ring is not None
(theis not
is a Python operator)
And you must provide timings to see whether it is better with your code for small input. You can use the command timeit
as in
sage: R.<T> = ZZ[] sage: p = R.random_element(degree=10) sage: timeit("p.roots()") 625 loops, best of 3: 223 µs per loop
It would be nice to make the following works:
sage: R.<T> = PolynomialRing(ZZ, sparse=True) sage: p = R.random_element(degree=10) sage: p.roots(QQbar) Traceback (most recent call last): ... AttributeError: 'Polynomial_generic_sparse' object has no attribute 'gcd'
It will be a nice thing to have in Sage!
More remarks to come, Vincent
comment:6 Changed 5 years ago by
 Commit changed from e725544b76363ea517a262419a092f9414c1b4a4 to f344cdc5b4a9e316a52f5372a30a3db1ac5043ab
Branch pushed to git repo; I updated commit sha1. New commits:
f344cdc  Cosmestic changes + add conditions on the constant coefficient

comment:7 Changed 5 years ago by
Hi Vincent,
Thanks for the comments. I pushed a new version taking them into account, especially for the "better code" part.
For roots over QQbar, I think it is not really related to this ticket. It is related to #15790 (implementation of GCD for sparse polynomials). I'd like to work in this other ticket too.
Now for the most complicate part (to my mind), timings. The issue is the following: Whether the new algorithm is better than the previous one is not a matter of density. For instance, if you have p = 3  x + x^7 * q
where q
is a very highdegree dense polynomial, integer roots of p
are common roots of 3x
and q
. And the same thing may happen in the middle of a polynomial. In a sense, the only way to detect such things is to use the loop in the algorithm. Of course, for lowdegree random polynomials, this will make the algorithm slower.
For the discussion, I refer to the previous algorithm
self._roots_from_factorization(...)
as the dense algorithm, and the new algorithm I implemented as the sparse algorithm.
I did tests:
 In the first one, the roots are computed using the dense algorithm and then the sparse algorithm, and I compute the ratio of the time it takes. The ration is
>1
if the sparse algorithm is worst than the dense one. The left column is the degree, and then 10 tests are performed with random polynomial of the given degree. As you can see, the sparse algorithm is often worst, by a factor between 1.5 and 3, and is sometimes better by a factor of 10 to 25.  In the second one, I added a condition of the sparsity (i.e. number of nonzero terms) of the polynomial:
if p.degree()  p.sparsity()  p.valuation() < 10
, I use the dense algorithm. The timings are then almost equal for both algorithm for small degrees, but sometimes (not as often as for the first test) the new algorithm is much better. For large degrees, we get the same kind of results as in the first test.
As a consequence of these tests, I do not know how the algorithm should be chosen. Do you have any idea (or any idea of which other tests I should perform to be able to chose)?
Finally, I also use the fact that a nonzero root of an integer polynomial should divide its lowestdegree coefficient (constant coefficient if the valuation is zero) to accelerate the computations. On the one hand, if the absolute value cst_coeff
of the constant coefficient is small, the dense algorithm can compute the product q
of the (xi)*(x+i)
for 1 <= i <= cst_coeff
and compute the roots of the gcd
of self
and q
instead of q
. My tests show that it really improves the computation time. (Note yet that because of #15790, this cannot work for the moment with sparse polynomials.) On the other hand, in the sparse algorithm, it is useless to compute the split of self
into a list polys
when the cst_coeff
is 1
.
Tests
First test:
2 [1.80, 2.15, 2.27, 2.14, 1.80, 1.37, 1.37, 2.10, 2.19, 2.19] 4 [2.17, 1.32, 1.41, 2.12, 2.22, 1.89, 2.01, 1.43, 1.85, 2.11] 6 [2.14, 1.38, 0.968, 1.13, 0.466, 2.12, 1.88, 2.15, 1.87, 1.96] 8 [2.07, 1.93, 2.02, 1.83, 2.08, 1.83, 2.57, 1.76, 1.82, 1.95] 10 [1.77, 2.06, 2.10, 2.13, 2.16, 1.78, 0.709, 2.26, 2.16, 2.02] 12 [1.76, 1.79, 0.671, 1.14, 1.26, 2.12, 2.00, 1.97, 1.78, 2.68] 14 [1.66, 0.651, 1.66, 2.07, 2.31, 0.590, 1.13, 0.930, 2.04, 1.66] 16 [2.13, 1.63, 2.24, 0.586, 0.802, 0.772, 1.71, 1.67, 0.592, 2.06] 18 [0.540, 0.667, 2.11, 1.70, 2.19, 1.99, 1.99, 1.98, 0.758, 1.07] 20 [0.568, 1.64, 1.01, 2.17, 1.64, 0.460, 2.07, 2.81, 1.83, 0.501] 22 [0.392, 2.18, 1.63, 1.50, 1.66, 1.95, 1.77, 0.394, 1.53, 2.08] 24 [0.545, 2.86, 1.64, 1.72, 1.94, 2.19, 1.66, 1.56, 0.351, 0.427] 26 [0.381, 1.51, 0.318, 0.394, 1.99, 0.467, 1.52, 1.89, 1.61, 1.53] 28 [1.55, 1.59, 1.95, 0.320, 2.08, 1.48, 1.99, 1.85, 1.47, 1.51] 30 [2.01, 2.02, 2.04, 1.47, 2.02, 0.267, 1.47, 0.339, 2.04, 2.02] 32 [2.01, 1.99, 1.99, 1.99, 2.01, 2.01, 2.01, 2.01, 0.259, 0.288] 34 [2.01, 2.00, 2.01, 0.287, 2.03, 2.01, 0.398, 0.245, 2.01, 1.98] 36 [0.237, 2.00, 2.02, 1.99, 2.02, 2.01, 2.01, 2.01, 1.94, 0.262] 38 [2.00, 0.256, 0.283, 2.01, 2.00, 0.242, 2.01, 2.01, 2.01, 0.235] 40 [1.99, 2.01, 1.98, 2.00, 2.00, 2.01, 2.00, 2.01, 2.01, 2.00] 42 [2.01, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 0.218, 1.99, 2.01] 44 [2.00, 2.00, 0.224, 2.00, 2.00, 2.00, 0.597, 2.00, 2.00, 0.203] 46 [0.174, 2.17, 0.170, 2.22, 0.555, 1.84, 0.645, 1.96, 0.188, 2.01] 48 [1.95, 2.03, 1.98, 2.03, 1.92, 1.94, 2.09, 2.00, 2.00, 2.02] 50 [1.97, 1.97, 2.03, 1.99, 2.00, 0.671, 2.00, 2.00, 0.155, 0.149] 52 [1.99, 0.172, 2.00, 2.01, 1.99, 0.173, 0.141, 2.01, 2.00, 1.99] 54 [1.99, 0.119, 2.00, 2.00, 1.99, 0.125, 0.154, 1.99, 1.99, 2.07] 56 [1.99, 2.00, 0.424, 1.98, 2.01, 1.99, 2.04, 2.01, 2.00, 2.00] 58 [2.00, 1.99, 0.590, 2.02, 0.123, 0.122, 2.01, 1.98, 2.00, 1.99] 60 [2.00, 1.99, 0.100, 2.02, 2.05, 2.02, 1.98, 1.95, 0.117, 1.96] 62 [0.0989, 1.97, 1.98, 0.126, 2.00, 1.98, 0.0992, 0.108, 0.120, 1.99] 64 [1.98, 1.98, 0.0956, 0.142, 0.0896, 2.02, 0.110, 0.106, 2.02, 2.00] 66 [2.01, 2.00, 0.210, 2.00, 2.00, 2.00, 1.98, 1.97, 1.99, 0.115] 68 [2.01, 0.0913, 2.00, 2.00, 2.01, 1.98, 1.99, 0.0813, 2.02, 1.99] 70 [0.0866, 0.0826, 2.00, 2.00, 2.00, 0.0900, 2.00, 1.99, 2.00, 0.105] 72 [2.00, 1.99, 0.268, 2.00, 2.00, 2.00, 2.00, 1.99, 2.00, 2.00] 74 [2.00, 0.430, 2.00, 2.00, 2.00, 0.0803, 2.00, 0.0924, 2.03, 2.02] 76 [1.97, 1.99, 1.98, 0.111, 2.00, 2.00, 2.00, 2.01, 0.0905, 0.0919] 78 [2.00, 2.02, 2.00, 1.99, 0.0716, 2.01, 2.00, 2.00, 2.00, 2.00] 80 [1.99, 0.0748, 2.00, 0.437, 1.99, 2.01, 2.00, 0.0731, 2.00, 2.01] 82 [2.00, 1.99, 0.0549, 0.0575, 2.01, 1.99, 0.0642, 0.0519, 1.99, 2.00] 84 [2.00, 0.0547, 1.99, 2.02, 2.01, 2.01, 1.99, 0.0712, 2.00, 2.01] 86 [1.97, 2.02, 2.00, 2.00, 2.01, 2.00, 2.00, 2.00, 0.0562, 2.01] 88 [2.00, 0.0534, 0.0605, 2.00, 2.00, 0.0567, 0.0709, 1.99, 2.00, 2.00] 90 [0.0553, 0.0496, 2.01, 1.99, 2.00, 2.00, 1.98, 0.0476, 2.01, 0.0541] 92 [1.98, 1.99, 2.00, 0.0545, 2.00, 2.00, 2.00, 1.99, 2.01, 2.01] 94 [1.99, 0.0425, 0.0579, 1.96, 2.01, 0.250, 0.106, 0.0561, 1.97, 2.00] 96 [0.0520, 2.00, 1.96, 2.00, 2.00, 0.0508, 2.00, 0.0480, 1.99, 2.01] 98 [0.0458, 0.0515, 0.0415, 2.01, 1.99, 2.01, 0.0510, 0.0526, 0.0455, 2.00]
Second test
2 [1.02, 1.02, 1.02, 1.02, 1.01, 1.01, 1.02, 1.02, 1.01, 1.02] 4 [1.02, 1.02, 1.02, 1.02, 1.02, 1.02, 1.02, 1.02, 1.02, 1.02] 6 [1.01, 1.03, 1.02, 1.03, 1.02, 1.02, 1.02, 1.02, 1.02, 1.02] 8 [1.01, 1.01, 1.02, 1.01, 1.01, 1.01, 1.02, 1.02, 1.01, 1.02] 10 [1.01, 1.02, 1.01, 1.01, 1.01, 1.01, 1.02, 1.02, 1.01, 1.02] 12 [1.01, 1.01, 1.02, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01] 14 [1.01, 1.01, 1.01, 1.02, 1.01, 1.01, 1.01, 1.02, 1.01, 1.01] 16 [1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01] 18 [1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01] 20 [1.01, 1.01, 1.01, 1.01, 1.00, 1.01, 1.01, 1.01, 1.01, 1.01] 22 [1.00, 1.00, 1.00, 1.01, 1.01, 1.01, 1.01, 1.00, 1.01, 1.00] 24 [1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.00, 1.00] 26 [0.996, 0.966, 1.00, 1.01, 1.01, 1.01, 1.00, 1.00, 1.00, 1.01] 28 [1.01, 1.00, 1.00, 1.01, 1.00, 1.01, 1.01, 1.00, 1.00, 1.00] 30 [1.00, 1.01, 1.00, 1.01, 1.01, 1.01, 1.00, 1.01, 1.01, 1.00] 32 [1.00, 1.06, 1.01, 1.00, 1.00, 1.00, 1.01, 1.00, 1.00, 1.00] 34 [1.00, 1.00, 1.01, 1.00, 1.00, 1.01, 1.00, 1.01, 1.00, 1.01] 36 [1.00, 1.01, 1.01, 1.01, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00] 38 [1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.358] 40 [1.01, 1.00, 1.00, 2.01, 1.00, 1.00, 1.00, 1.01, 1.00, 1.00] 42 [1.00, 1.00, 1.01, 1.00, 2.01, 1.00, 1.00, 1.00, 1.00, 1.00] 44 [1.00, 1.00, 1.00, 1.01, 1.00, 1.00, 1.00, 1.00, 1.01, 1.00] 46 [0.150, 1.00, 1.00, 1.00, 1.00, 1.00, 0.196, 1.00, 0.188, 2.00] 48 [0.153, 2.01, 1.01, 2.00, 1.01, 2.00, 1.00, 1.00, 1.00, 1.00] 50 [1.00, 1.00, 1.00, 1.00, 2.01, 1.00, 1.00, 0.163, 1.00, 1.99] 52 [1.00, 1.00, 1.00, 1.00, 0.166, 1.02, 1.00, 1.00, 1.01, 1.00] 54 [1.00, 0.166, 1.00, 1.00, 1.00, 2.01, 2.00, 0.153, 2.00, 0.992] 56 [1.00, 2.00, 1.00, 1.00, 2.01, 0.134, 2.00, 2.00, 0.999, 1.00] 58 [0.604, 0.423, 0.132, 2.00, 1.99, 2.00, 1.00, 1.00, 2.01, 0.119] 60 [1.00, 1.00, 0.117, 2.00, 1.00, 2.00, 2.01, 2.00, 0.119, 1.00] 62 [2.00, 0.151, 2.00, 2.00, 1.00, 0.139, 1.00, 2.00, 2.01, 0.106] 64 [1.99, 2.00, 0.103, 1.99, 2.00, 2.00, 0.120, 1.00, 1.00, 1.00] 66 [2.00, 2.00, 1.99, 2.00, 2.00, 2.00, 2.00, 2.01, 2.00, 1.00] 68 [1.99, 1.99, 1.00, 2.00, 2.00, 1.00, 2.00, 2.00, 0.950, 1.99] 70 [1.98, 2.00, 2.01, 1.00, 1.00, 0.0791, 2.00, 1.00, 0.118, 0.0964] 72 [0.429, 2.00, 2.00, 1.99, 1.98, 2.00, 1.99, 2.00, 0.128, 2.01] 74 [0.0890, 0.0795, 1.99, 2.00, 2.01, 2.00, 1.97, 0.0945, 2.00, 1.00] 76 [1.99, 1.99, 2.01, 2.01, 2.02, 2.01, 0.0878, 2.00, 0.0740, 2.01] 78 [1.00, 2.00, 2.03, 2.00, 1.00, 2.00, 0.0617, 0.0812, 2.00, 1.00] 80 [0.0651, 0.0687, 1.98, 2.00, 1.98, 2.00, 2.00, 0.0696, 0.0791, 0.0770] 82 [2.00, 0.0788, 2.00, 1.99, 0.0828, 2.00, 2.02, 2.00, 0.0669, 2.01] 84 [1.02, 1.99, 2.00, 2.01, 0.318, 2.00, 2.00, 1.99, 2.01, 2.01] 86 [2.00, 0.0669, 2.00, 2.00, 2.00, 1.98, 0.0551, 2.00, 2.00, 0.135] 88 [0.0540, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 1.99, 2.00, 0.0682] 90 [2.00, 2.02, 1.91, 2.01, 2.01, 2.00, 2.00, 2.00, 2.00, 2.00] 92 [2.00, 1.98, 2.00, 1.99, 1.99, 2.00, 2.00, 2.00, 2.00, 2.00] 94 [1.00, 0.0574, 2.00, 1.00, 2.00, 2.00, 0.0594, 2.00, 2.00, 2.00] 96 [2.00, 2.00, 2.00, 2.00, 0.0420, 2.00, 2.00, 2.00, 2.01, 2.00] 98 [0.0425, 0.0565, 2.01, 2.00, 2.00, 1.99, 1.99, 2.00, 0.0462, 2.00]
comment:8 Changed 5 years ago by
comment:9 Changed 5 years ago by
 Status changed from needs_review to needs_work
Hi,
Global remarks:
1) You seem to have a global problem with the difference between ==
and is
so have a look at pythonisoperator on stackoverflow.
2) Comment your code when it is complicated like
if cst_coeff is not ZZ(1): i_min=0 polys=[] for i in xrange(1,k): if e[i]e[i1] > c_max.nbits(): polys.append(R(p[ e[i_min]:e[i] ].shift(e[i_min]))) i_min=i c_max=c[i].abs() else: c_max=max(c[i].abs(),c_max) polys.append(R(p[ e[i_min]:1+e[k1] ].shift(e[i_min])))
open softwares should also be readable softwares
3) You have syntax error in the documentation (which will create error when you try to build the documentation with "make doc"):
``algorithm''  the algorithm to use
should be
``algorithm``  the algorithm to use
ie open and close with back quotes.
Specific ones:
4) In the method sparsity that you implemented in "polynomial_element.pyx", the variable c is not initialized. So
sage: K.<x>=QQ[] sage: (x^7 + x^3 + 1).sparsity() 32665
Moreover, the following test is not safe at all
if l.pop() is not zero
you can not believe that the zero always occupy the same memory. For instance
sage: K = QQ['x'] sage: K.zero().constant_coefficient() is QQ.zero() False
Moreover, using try except is totally useless in that case... You might be inspired by the implementation of coefficients
. Hopefully, I am not the one who teach you programming at school ;) Please test this method with other base rings (at least QQ, QQbar, ZZ/nZZ GF(p), GF(p^{n)).
}
5) Are you sure that the term sparsity is standard? I would rather go for something more explicit like "num_nonzero_coefficients" or something similar. I hoped to find a similar method in matrices but did not find it.
6) In your function _roots_univariate_polynomial
, there is no gain in using xrange
instead of range
. But there will be a big one if you define k
as an int!
6) Using the gcd
from sage.rings.arith
(in the line p=gcd(p,q)
) is slow compared to p = p.gcd(q)
. (and do not forget to remove the import)
Vincent Vincent
comment:10 Changed 5 years ago by
7) the method .extend()
of lists do not return anything as you can see
sage: [0].extend([1])
so there is a problem in the case algorithm "dense_with_gcd" in your function _roots_univariate_polynomial
.
8) How can you divide by the content??
sage: R.<x> = PolynomialRing(ZZ, sparse=True) sage: p = (x^7 + x^3 + 1) sage: p.content() Principal ideal (1) of Integer Ring sage: (x^7 + x^3 + 1) / p.content()
EDIT: my mistake. This is a bug in the code of polynomials: the method content must either return an ideal or a number but not one or the other depending on the implementation...
Nevertheless, you should test all cases of your code within the documentation!!
Vincent
comment:11 Changed 5 years ago by
For the content issue, I opened a thread on sagedevel.
comment:12 Changed 5 years ago by
 Commit changed from f344cdc5b4a9e316a52f5372a30a3db1ac5043ab to be6deca10dd8d25ace0e9d821c6cc5551ca8d47f
comment:13 Changed 5 years ago by
Below are my answers to Vincent's comments:
1) ==
vs is
: Should be OK now
2) Comment the code: Done
3) Syntax error in the doc: Done. "make doc" works now.
4) Method sparsity
in "polynomial_element.pyx": I totally changed the method since it was obviously badly written. I wanted to avoid constructing a list of nonzero monomials to get a better complexity but this was not a good idea! I hope the new (much simpler!) version is fine. I add quite a lot of tests in this function (as well as in others).
5) Name of the method sparsity
: I would say that this name is fairly standard. It may be used more often to speak about "the fact of being sparse or not" rather than "the number of nonzero monomials", but still you can find this term in several papers (cf for instance http://www.math.purdue.edu/~sbasu/ccc04.ps or http://eccc.hpiweb.de/report/2014/056/download/). Further, I find the term nicer than using a periphrasis but I let you tell me what you think is best. I don't mind changing the term if it appears not to be appropriate. I looked at similar terms in matrix_sparse.pyx and related files. I appears that num_nonzero
is used a bit, and sparsity
is only used in comments to speak about "the fact of being sparse or not".
6) range
vs xrange
, and def k
as int
: I defined k
(as well as other variables i
and j
) as int
, but I let xrange
rather than range
. I did not notice the promised big gain though.
7) .extend()
does not return anything: I note. Yet I removed the case "dense_with_gcd" since I have not been able to build a case where this is faster than the algorithm "dense".
8) content
issue: I let a distinction between sparse and dense polynomials until the issue is solved.
comment:14 Changed 5 years ago by
 Status changed from needs_work to needs_review
comment:15 Changed 5 years ago by
 Milestone changed from sage6.3 to sage6.4
comment:16 Changed 5 years ago by
Hi Bruno,
I started modifying the code by myself. For that I opened a branch at u/vdelecroix/16516
. You can fetch from this one and adopt my commit if you like it (using git pull
inside your branch). I describe the changes below
0) I update above sage.6.3
1) Be careful, you use too often the "is" where "==" is needed: first example strings are not unique representation
sage: "dense" is "dense" True sage: other_dense = "denseX"[:1] sage: "dense" == other_dense True sage: "dense" is other_dense False
second example: you wrote
if g.is_zero(): # could also use the slightly faster # "g is R.zero_element()" since it is cached
But the zero of R is not unique
sage: (ZZ(1)  ZZ(1)) is ZZ(0) False
(though in that case, I only had to remove the comment)
2) pep8 recommands to have spaces x = y
instead of x=y
3) It is not needed to copy/paste doctests (that would imply that the code is tested twice and hence multiply the test timings by two). I modified the documentation of sparsity
appropriately.
I guess this is soon finished. Could you redo the exact timings you did?
Vincent
comment:17 Changed 5 years ago by
 Commit changed from be6deca10dd8d25ace0e9d821c6cc5551ca8d47f to 0efc4ea81dc2c47bb9b0246f37424e296bb8ee42
comment:18 Changed 5 years ago by
Hi Vincent,
I accept your changes, thank you for the work!
For the timings, I did the following tests, with the results attached as .txt
files. I tested p.roots(algorithm="...")
for random polynomials p
of degree 2 to 98, where the algorithm is either "sparse"
or "dense"
. I computed the ratio between the two timings (smaller than 1 says "sparse"
is better). There are four different files:
 I ran the tests either with a check
p.degree()  p.sparsity()  p.valuation() < 10
or without this test: The files are named..._test.txt
in the first case and..._no_test.txt
in the second one.  I ran the tests with polynomials in
PolynomialRing(ZZ,sparse=True)
(filessparse_pol_...
) andPolynomialRing(ZZ,sparse=False)
(filesdense_pol_...
).
To redo the same tests, here is the code I used:
nb = 3 for degre in xrange(2,100,2): liste = [] for _ in xrange(10): p = R.random_element(degre) t1 = timeit("p.roots(algorithm=\"sparse\")",seconds=True,number=nb) t2 = timeit("p.roots(algorithm=\"dense\")",seconds=True,number=nb) liste.append(QQ(t1/t2).n(digits=3)) print degre, liste, min(liste), sum(liste)/len(liste), max(liste)
To read the results, the lines looks like:
2 [0.939, 0.392, 0.523, 1.51, 1.17, 1.25, 1.17, 0.517, 0.503, 0.855] 0.392 0.884 1.51
where the 5 elements are (in order):
 the degree of the polynomials (here 2);
 list of 10 ratios corresponding to 10 independent tests;
 minimum of the list;
 average of the list;
 maximum of the list.
My conclusion is that we shall always use the sparse algorithm (at least for sparsely represented polynomials) without making any test such as the one I did in the timings. Do you agree?
Note: I did not check this really but the better timings may quite often come from the test on the constant coefficient ("is it ±1?") that I do but is not done in the other roots
function. In random polynomials, the constant coefficient is 1 with probability 2/5.
Changed 5 years ago by
Changed 5 years ago by
Changed 5 years ago by
Changed 5 years ago by
Changed 5 years ago by
comment:19 Changed 5 years ago by
A comment to build on my final note in the previous comment:
First, what is important is not only the constant coefficient but the lowest degree nonzero coefficient, which equals ±1 with probability 1/2 for R.random_polynomial(degree)
for R
some polynomial ring. I made once again the same tests, but with the 5th line replaced by
p = R.random_element(degre,1000,1000)
such that the lowest degree nonzero coefficient in ±1 with probability 1/999 (if I am not mistaken). Note also that the polynomials are very dense since the coefficient 0 appears only with probability 1/1999. The results are of course different (I put just a very few of them to get an idea):
Without test:
2 [1.25, 1.20, 0.939, 1.13, 1.14, 1.16, 1.15, 1.16, 1.14, 1.13] 0.939 1.14 1.25 4 [1.14, 1.04, 1.16, 0.964, 1.22, 1.14, 1.14, 1.13, 1.14, 1.13] 0.964 1.12 1.22 20 [1.04, 1.09, 1.05, 1.04, 1.07, 1.05, 1.03, 1.04, 1.04, 1.06] 1.03 1.05 1.09 30 [1.03, 1.03, 1.02, 1.08, 1.03, 1.03, 1.04, 1.08, 1.04, 1.03] 1.02 1.04 1.08 40 [1.01, 1.03, 1.01, 1.02, 1.09, 1.04, 1.04, 1.02, 1.02, 1.02] 1.01 1.03 1.09
With test:
2 [1.34, 1.17, 1.18, 1.04, 1.04, 1.11, 1.07, 1.03, 1.02, 1.09] 1.02 1.11 1.34 4 [1.04, 1.02, 1.03, 0.832, 1.02, 1.08, 1.04, 1.08, 1.04, 0.995] 0.832 1.02 1.08 20 [0.985, 1.14, 1.01, 1.01, 1.02, 0.994, 0.982, 0.934, 1.15, 1.19] 0.934 1.04 1.19 30 [0.989, 1.02, 1.01, 1.01, 0.991, 1.00, 1.01, 1.00, 1.02, 0.973] 0.973 1.00 1.02 40 [1.01, 1.02, 1.01, 1.00, 1.01, 1.00, 1.01, 0.999, 1.00, 0.998] 0.998 1.01 1.02
comment:20 followup: ↓ 21 Changed 5 years ago by
Hello Bruno,
In dense_pol_test, the entries looks like being either 1 or 0.001 in magnitude. Do you know what makes the sparse algorithm better in the second case?
How did you generate your random polynomials? The output of random_polynomial
is by definition a dense polynomial, i.e. all coefficients are nonzero. In order to test your algorithm you must look at sparse polynomial, i.e. when only 10% of the coefficients are nonzero. The datastructure (i.e. PolynomialRing(sparse=False)
vs PolynomialRing(sparse=True)
) should not play an important role in that game.
As you can see
sage: %runfile random_polynomial.py # home made random_integer_polynomial sage: p = random_integer_polynomial(degree=100, density=0.1) sage: p.sparsity() 7 sage: timeit("p.roots(algorithm='dense')") 125 loops, best of 3: 5.59 ms per loop sage: timeit("p.roots(algorithm='sparse')") 625 loops, best of 3: 242 µs per loop
(be careful, I had to change a line in the code in order to be able to use dense/sparse alternative here).
Vincent
comment:21 in reply to: ↑ 20 Changed 5 years ago by
Salut Vincent,
In dense_pol_test, the entries looks like being either 1 or 0.001 in magnitude. Do you know what makes the sparse algorithm better in the second case?
Actually, there are two possibilities which both occur for the sparse algorithm to beat the dense one:
 Either some lowdegree monomial has a zero coefficient: In such case, the sparse algorithm splits the input polynomial into a lowdegree polynomial and another one (of larger degree), then take their gcd and compute the roots. This is faster than the dense algorithm and here this really exploit the sparsity of the polynomial to speedup the computation. For instance, if the constant coefficient is nonzero and the linear coefficient is zero, the sparse algorithm boils down to computing the roots of a monomial which is of course very fast.
 Or, and this does not exploit sparsity at all, the constant coefficient is ±1: In this case, there is a kind of trivial optimization in my code to only test for the possible roots 1, 0 and 1. This is sensible in the tests I made since the constant coefficient (or more accurately the lowest degree monomial) is 1 with probability 1/2.
Tests
I redid some tests with polynomials of smaller density. My code for producing random polynomial of given density and degree is as follows:
def random_element_with_density(R,deg,dens): K = R.base_ring() p = {} c = K.random_element() while c == 0: c = K.random_element() p[deg]=c while len(p) < dens * deg: e = ZZ.random_element(deg) if p.has_key(e): continue c = K.random_element() if c == 0: continue p[e] = c return R(p)
First, I did the same tests as before but with density 1/10. (The results are inverse of before, that is >1
means the sparse algorithm is faster. Moreover, all the tests are made without any condition of the form if p.degree()  p.sparsity()  p.valuation() < 10
.) I display only some results for concision:
5 [3.04, 3.42, 3.59, 3.77, 3.36, 3.54, 3.44, 4.01, 3.47, 3.88] 3.04 3.55 4.01 25 [10.2, 1.44, 2.16, 4.05, 4.42, 5.22, 3.71, 3.65, 1.53, 1.48] 1.44 3.79 10.2 45 [11.1, 1.89, 10.6, 5.16, 5.64, 16.3, 7.58, 12.4, 2.35, 12.0] 1.89 8.50 16.3 65 [10.8, 20.9, 21.3, 22.6, 12.4, 22.3, 10.9, 7.10, 25.1, 5.49] 5.49 15.9 25.1 85 [40.5, 12.0, 40.3, 41.9, 23.7, 30.2, 21.5, 5.25, 14.4, 4.98] 4.98 23.5 41.9
To come up with the second source of improvement (which has nothing to do with sparsity), I did the same test after removing the condition if not cst_coeff.abs().is_one()
:
5 [3.15, 3.37, 3.79, 3.05, 3.48, 3.43, 3.42, 3.43, 3.38, 3.45] 3.05 3.39 3.79 25 [3.16, 1.84, 2.62, 1.48, 1.53, 1.37, 2.63, 1.73, 1.82, 3.64] 1.37 2.18 3.64 45 [4.15, 12.5, 5.01, 12.4, 5.31, 4.47, 5.57, 8.59, 5.49, 4.65] 4.15 6.82 12.5 65 [24.1, 13.8, 15.2, 11.0, 6.08, 13.7, 23.8, 13.0, 1.62, 14.2] 1.62 13.7 24.1 85 [29.2, 29.5, 18.0, 19.4, 30.0, 10.1, 25.8, 21.3, 7.81, 30.6] 7.81 22.2 30.6
As you can see, there is still an improvement though a lesser one.
Finally, I did the same comparisons of algorithms with different densities. Each line corresponds to a degree and each column to a density. Each entry in the following tables is an average of 30 tests for the corresponding degree and density. The first table is with the "trivial optimization" and the second one without. You can see that without the optimization, the sparse algorithm becomes worse when the density is close to 1.
0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1  

5  3.37  3.40  2.33  2.44  1.96  1.79  2.35  2.00  1.73  1.95 
10  3.34  2.92  2.40  2.56  1.97  2.53  2.11  2.47  1.98  1.96 
15  3.25  3.13  2.67  2.79  2.63  2.63  2.69  2.17  2.34  1.96 
20  4.44  3.60  3.32  3.53  2.67  2.68  2.90  3.46  2.59  3.32 
25  3.57  2.95  4.40  2.93  3.29  3.22  3.29  3.20  2.05  3.00 
30  3.67  5.22  6.28  5.85  5.34  6.13  4.70  4.53  4.88  4.51 
35  5.25  6.17  7.34  8.25  6.99  6.53  6.45  7.68  4.66  4.99 
40  5.38  7.95  8.17  9.27  7.89  7.81  7.55  6.98  5.27  7.17 
45  7.45  10.1  9.40  9.36  9.17  8.86  9.73  7.97  6.85  5.92 
50  8.01  11.4  12.4  10.6  10.8  10.1  8.77  10.9  9.08  9.51 
55  10.9  13.5  12.9  14.0  13.4  12.8  12.0  10.0  10.4  7.01 
60  10.4  14.4  15.5  16.5  15.1  17.8  11.1  12.5  10.7  9.76 
65  14.2  15.2  16.2  18.6  16.3  14.5  13.9  16.1  11.9  12.2 
70  16.3  20.4  21.7  18.6  19.1  19.1  17.2  16.1  14.4  13.7 
75  18.6  22.2  25.5  23.5  24.2  17.3  24.1  18.9  16.8  19.8 
80  24.5  23.2  25.5  25.9  25.8  26.5  15.5  23.3  14.5  20.3 
85  21.3  30.0  28.5  32.4  29.4  24.8  21.5  22.4  25.6  28.9 
90  30.7  32.5  30.7  31.7  33.1  31.7  22.2  23.3  17.8  18.7 
95  30.6  34.6  33.2  36.1  37.5  31.6  30.3  28.2  28.1  18.0 
0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1  

5  3.47  3.32  1.88  1.79  1.39  1.65  1.60  1.56  1.16  1.00 
10  3.32  2.62  2.00  1.80  2.30  1.88  1.78  1.30  1.24  1.17 
15  3.09  2.28  2.18  2.28  2.11  1.71  1.43  1.43  1.25  1.14 
20  3.27  3.02  2.82  2.48  2.43  2.20  1.76  1.74  1.45  0.948 
25  3.99  3.50  2.77  3.40  3.05  3.00  2.09  2.13  1.30  1.14 
30  3.39  5.01  5.22  4.48  4.49  4.79  3.78  1.74  1.24  1.18 
35  5.26  5.09  6.10  5.89  6.02  5.07  4.38  2.51  1.76  0.972 
40  4.98  6.62  6.31  7.26  7.11  4.87  3.70  4.71  1.37  1.65 
45  5.97  7.89  9.24  8.37  7.81  6.92  4.84  2.76  1.41  1.44 
50  7.72  8.55  9.28  10.2  8.48  8.50  5.52  5.66  2.61  0.983 
55  10.3  11.6  11.5  11.5  11.2  10.1  6.66  2.55  1.58  0.986 
60  11.5  13.3  13.5  13.3  13.1  11.8  5.60  5.63  2.91  0.986 
65  14.3  14.3  16.2  17.6  14.6  13.0  7.44  4.32  6.36  1.83 
70  13.7  16.9  17.0  19.1  13.6  13.3  9.85  6.61  3.10  0.991 
75  15.8  19.8  18.2  20.9  18.7  15.2  13.4  5.46  2.69  1.01 
80  19.2  22.6  25.7  20.2  21.9  18.5  9.91  9.94  5.85  0.992 
85  24.1  22.4  27.2  24.5  23.7  15.5  17.1  7.90  7.02  0.988 
90  24.5  29.8  28.7  29.4  27.6  17.2  11.9  8.36  6.91  2.20 
95  26.2  29.5  32.2  28.5  21.3  23.4  15.5  9.56  3.82  0.991 
comment:22 followup: ↓ 26 Changed 5 years ago by
Wouldn't "density" be a better word than "sparsity"?
comment:23 Changed 5 years ago by
You really should document this:
if ring is not self and ring is not None: return None
and related
 return K._roots_univariate_polynomial(self, ring=ring, multiplicities=multiplicities, algorithm=algorithm) + ret = K._roots_univariate_polynomial(self, ring=ring, multiplicities=multiplicities, algorithm=algorithm) + if ret is not None: return ret
You could also seriously consider changing return None
to the more Pythonic raise NotImplementedError
.
When a method has an algorithm
argument, you should document the possible values for algorithm
and say what the default is. Also, you need to doctest those various options for algorithm
and reject bad values such as algorithm="foobar"
.
comment:24 Changed 5 years ago by
 Status changed from needs_review to needs_work
comment:25 Changed 5 years ago by
 Commit changed from 0efc4ea81dc2c47bb9b0246f37424e296bb8ee42 to 711e73c530ce0e56beb99fc2cddd95e3901d7385
Branch pushed to git repo; I updated commit sha1. New commits:
711e73c  Automatic choice of algorithm + replace "return None" by an exception + doc

comment:26 in reply to: ↑ 22 ; followup: ↓ 28 Changed 5 years ago by
 Status changed from needs_work to needs_review
Replying to jdemeyer:
Wouldn't "density" be a better word than "sparsity"?
I have the feeling that "density" is for a ratio. In my case, what I'd like to denote is the number of terms, not the number of terms divided by some quantity such as the "maximal number of possible terms" (if such a thing can be defined).
I've tried to take your other comments into account. Please let me know!
comment:27 followup: ↓ 29 Changed 5 years ago by
 Status changed from needs_review to needs_work
 Remove commentedout code:
#or p.degree()  p.sparsity() p.valuation() > 10:
 Add doctests for the various
algorithm=
parameters, includingalgorithm="foobar"
 Remove
if hasattr(K, '_roots_univariate_polynomial'):
and changeexcept NotImplementedError
toexcept (AttributeError, NotImplementedError):
comment:28 in reply to: ↑ 26 Changed 5 years ago by
Replying to bruno:
Replying to jdemeyer:
Wouldn't "density" be a better word than "sparsity"?
I have the feeling that "density" is for a ratio.
And I have the feeling that "sparsity" should be a larger value if the polynomial is more sparse. In my intuition, the "sparsity" of a polynomial like x^10 + x^9 + x^8 + ... + x + 1
would be 0 since it's absolutely not sparse.
How about number_of_terms()
?
comment:29 in reply to: ↑ 27 Changed 5 years ago by
Replying to jdemeyer:
And I have the feeling that "sparsity" should be a larger value if the polynomial is more sparse. In my intuition, the "sparsity" of a polynomial like
x^10 + x^9 + x^8 + ... + x + 1
would be 0 since it's absolutely not sparse.
Right, this is actually misleading.
How about
number_of_terms()
?
Sounds good. Vincent proposed num_nonzero_coefficients()
that I found a bit long. I think the nonzero
part can safely be removed. Maybe I would opt for number_of_coefficients()
since one gets the nonzero coefficients of a polynomial using coefficients()
.
Replying to jdemeyer:
 Remove commentedout code:
#or p.degree()  p.sparsity() p.valuation() > 10:
Actually, I'd like to decomment this code (maybe after some modification) based on the tests I ran. Even if a polynomial is densely represented, it may be actually quite sparse and then the sparse algorithm can be much better. The difficulty here is to find the right condition. We could also always use the sparse algorithm (which falls back to the dense algorithm when the input is not sparse enough) but this adds some overhead. This is especially sensible for lowdegree polynomials where the computation is very fast.
 Add doctests for the various
algorithm=
parameters, includingalgorithm="foobar"
 Remove
if hasattr(K, '_roots_univariate_polynomial'):
and changeexcept NotImplementedError
toexcept (AttributeError, NotImplementedError):
I'll do that.
comment:30 Changed 5 years ago by
 Commit changed from 711e73c530ce0e56beb99fc2cddd95e3901d7385 to 4b15fc26b7ea40fcc9a897c1875fde28726a38ef
comment:31 Changed 5 years ago by
Why overcomplicate your doctest (see last commit) with a try
/except
? You can simply doctest the exception...
comment:32 Changed 5 years ago by
I've made the changes suggested.
The delicate point, from the beginning, is to automatically choose which algorithm to use. The choice I made is to always go through the sparse algorithm when no algorithm is specified. As you can see on the file extensive_tests.sobj I attached,¹ the proportion of polynomials for which the sparse algorithm is slower than the dense algorithm is very small. Actually, I made the timings a second time with the exact same polynomials, and I finally found only 9 polynomials (out of 24900) for which the sparse algorithm is slower with a ratio >1.2.
¹ The object is a list of triples (p,t1,t2)
where p
is a polynomial, t1
the timing (in seconds) with the dense algorithm and t2
the timing with the sparse algorithm.
[edit: answer to the newest comment of jdemeyer] What do you mean by this?
comment:33 Changed 5 years ago by
 Commit changed from 4b15fc26b7ea40fcc9a897c1875fde28726a38ef to 26a3fbb893bc728964e20e78e39e18156720149d
Branch pushed to git repo; I updated commit sha1. New commits:
26a3fbb  better doctest

comment:34 Changed 5 years ago by
 Status changed from needs_work to needs_review
Ok, I found how one is supposed to doctest exception
s, sorry not to have checked this before!
comment:35 Changed 5 years ago by
The sentence "Default is based on whether p.parent()
is sparse or not." should be fixed then.
comment:36 Changed 5 years ago by
 Commit changed from 26a3fbb893bc728964e20e78e39e18156720149d to 52fdc7876281d9086057d1cf985f7ae2c21df095
Branch pushed to git repo; I updated commit sha1. New commits:
52fdc78  Update doc

comment:37 Changed 5 years ago by
The code cannot be interrupted. Please have a look at http://www.sagemath.org/doc/developer/coding_in_cython.html#usingsigcheck (or better yet: $SAGE_ROOT/src/doc/output/html/en/developer/coding_in_cython.html
built from git develop branch) on how to fix this.
comment:38 followup: ↓ 40 Changed 5 years ago by
 Status changed from needs_review to needs_work
This is much slower than it should be, why is it?
sage: z=(polygen(ZZ)1)^2000; timeit('z.roots(algorithm="sparse")') 5 loops, best of 3: 2.83 s per loop
comment:39 Changed 5 years ago by
Note: in Python, you should use !=
instead of <>
comment:40 in reply to: ↑ 38 Changed 5 years ago by
Replying to jdemeyer:
This is much slower than it should be, why is it?
sage: z=(polygen(ZZ)1)^2000; timeit('z.roots(algorithm="sparse")') 5 loops, best of 3: 2.83 s per loop
It's a severe flaw in my implementation, though I do not know yet where it comes from. I'll try to find it out...
comment:41 Changed 4 years ago by
 Commit changed from 52fdc7876281d9086057d1cf985f7ae2c21df095 to d72a8f1554235a389d8c42e278d3779895d568b5
comment:42 Changed 4 years ago by
 Status changed from needs_work to needs_review
I changed the following aspects of my code:
 Replace each
<>
by!=
 Add
sig_on
/sig_off
's  Remove the tests on the constant coefficient: This was a hack to speed up some computations (when the constant coefficient is
±1
), but this leads to the slow down noticed by jdemeyer for the roots of(x1)^2000
. It had nothing to do with the"sparse"
algorithm. (More on this in the note below.)  I changed the automatic choice of algorithm: Based on extensive tests (not exhaustive though, by definition), it appears that the slowdown that exists in some examples using the
"sparse"
algorithm rather than the"dense"
algorithm becomes negligible from around degree100
, with never more that 10% slow down, and almost always <1% (and of course quite often there is a speedup!). Note that I also add a shortcut for completely dense polynomials (with density1
) when we are sure the"dense"
algorithm will be used as a fallback to avoid some useless (though not expensive) computations.
Almost unrelated note on the "constant coefficient trick": As written above, I have removed some tests that were made on the constant coefficient. The trick I used is that the constant coefficient is the product of the roots, thus if a polynomial has a "small" constant coefficient, it can be interesting to use some kind of brute force algorithm to discover these roots, or to first take a gcd before computing the roots. I removed the trick since it has some drawbacks. If it is more correctly written, it may become useful but in any case, it should be the object of another ticket. Maybe some day...
comment:43 Changed 4 years ago by
 Commit changed from d72a8f1554235a389d8c42e278d3779895d568b5 to f97f22bd0b549ff3f0e9406ab10e4ce989e0f6cd
Branch pushed to git repo; I updated commit sha1. New commits:
f97f22b  Merge branch 'u/bruno/faster_roots_computation_for_sparse_polynomials_over_zz' of git://trac.sagemath.org/sage into ticket/16516/roots_sparse_polys

comment:44 Changed 4 years ago by
 Status changed from needs_review to needs_work
4 failings doctests, see patchbot report
comment:45 Changed 4 years ago by
 Commit changed from f97f22bd0b549ff3f0e9406ab10e4ce989e0f6cd to e9b274f3c413f21738ad8d231323911ac5087f4b
comment:47 followup: ↓ 49 Changed 4 years ago by
 Status changed from needs_review to needs_work
one failing doctest
comment:48 Changed 4 years ago by
 Commit changed from e9b274f3c413f21738ad8d231323911ac5087f4b to 80be713ae2673ca1d769cc5f0ceb7b4a2153dd9c
comment:49 in reply to: ↑ 47 Changed 4 years ago by
 Status changed from needs_work to needs_review
comment:50 followup: ↓ 52 Changed 3 years ago by
 Status changed from needs_review to needs_work
1 failing doctest, due to deprecation of polynomial slicing
comment:51 Changed 3 years ago by
 Commit changed from 80be713ae2673ca1d769cc5f0ceb7b4a2153dd9c to d47a525a9c316d1b84c3323fb11c46a88803cafb
Branch pushed to git repo; I updated commit sha1. This was a forced push. Last 10 new commits:
763a898  trac #16516: cosmetic

f7b4998  16516: Automatic choice of algorithm + replace "return None" by an exception + doc

d06e3e5  16516: remove has_attr and replaced by an except AttributeError

971f70e  16516: doctest possible values for algorithm=...

8659299  16516: Add sig_on() & sig_off()'s

8554f3e  16516: Replace <> by !=

be39af0  16516: Add shortcut for totally dense polynomials

abbfd18  #16516: replace zero_element by zero

96f40ff  #16516: remove trailing whitespace

d47a525  #16516: handle deprecation of polynomial slicing + cosmetic

comment:52 in reply to: ↑ 50 Changed 3 years ago by
 Status changed from needs_work to needs_review
Replying to chapoton:
1 failing doctest, due to deprecation of polynomial slicing
Corrected.
In the same time, I rebased on 7.2beta4 because this branch defined a method number_of_terms
which has actually be introduced in #18617 (merged into 7.1beta6 or something like that). This implied some changes (including removing the now duplicate method number_of_terms
). I also took the opportunity to perform some cosmetic changes.
comment:54 Changed 2 years ago by
 Commit changed from d47a525a9c316d1b84c3323fb11c46a88803cafb to 48d4fc9e73c4e227711e6a656782791489eeef37
Branch pushed to git repo; I updated commit sha1. New commits:
48d4fc9  Merge branch 'develop' into 16516_faster_roots_ZZ

comment:55 Changed 2 years ago by
 Milestone changed from sage6.4 to sage8.0
 Status changed from needs_work to needs_review
Applies again, needs review again!
comment:56 followups: ↓ 57 ↓ 59 Changed 2 years ago by
Minor comments:
 Mark the docstring with
r"""
andReturns the root of the univariate polynomial ``p''. +Return the root of the univariate polynomial ``p``.
`_roots_from_factorization`
should either be:meth:`_roots_from_factorization`
or``_roots_from_factorization``
as it is not latex code.ZZ
>`\ZZ`
. No semicolons needed and our error messages star with a lowercase to match python:
raise ValueError("Roots of 0 are not defined"); raise ValueError("roots of 0 are not defined")
 Jeroen can probably explain it better than me, but he has told me you should try to avoid wrapping any Python code in
sig_on
/sig_off
as it is not needed. Since almost all of this looks like Python or Cython files that should be interruptible, I don't think many of them are necessary. You might want to ask him about this though.  You should move the reference to the master reference file.
comment:57 in reply to: ↑ 56 Changed 2 years ago by
Replying to tscrim:
 Jeroen can probably explain it better than me, but he has told me you should try to avoid wrapping any Python code in
sig_on
/sig_off
as it is not needed. Since almost all of this looks like Python or Cython files that should be interruptible, I don't think many of them are necessary. You might want to ask him about this though.
Even better, RTFM https://readthedocs.org/projects/cysignals/
comment:58 Changed 2 years ago by
 Commit changed from 48d4fc9e73c4e227711e6a656782791489eeef37 to df27d498e2f52fc0c3d902f5cca73e235188f923
comment:59 in reply to: ↑ 56 Changed 2 years ago by
Replying to tscrim:
Minor comments:
 Mark the docstring with
r"""
andReturns the root of the univariate polynomial ``p''. +Return the root of the univariate polynomial ``p``.`_roots_from_factorization`
should either be:meth:`_roots_from_factorization`
or``_roots_from_factorization``
as it is not latex code.ZZ
>`\ZZ`
. No semicolons needed and our error messages star with a lowercase to match python:
raise ValueError("Roots of 0 are not defined"); raise ValueError("roots of 0 are not defined") You should move the reference to the master reference file.
All of these have been taken into account.
 Jeroen can probably explain it better than me, but he has told me you should try to avoid wrapping any Python code in
sig_on
/sig_off
as it is not needed. Since almost all of this looks like Python or Cython files that should be interruptible, I don't think many of them are necessary. You might want to ask him about this though.
I tried to put sig_on
/sig_off
only when needed, performing a number of tests of CtrlC
. There are still a fairly large number of them.
comment:60 Changed 2 years ago by
Bruno, your sig_on/sig_off
handling is terribly wrong. For example in def _roots_univariate_polynomial
you have
sig_on() roots = p._roots_from_factorization(p.factor(), multiplicities) sig_off()
There is a lot of Python code between this sig_on
and sig_off
(here I count 4). Basically, the code between sig_on
sand sig_off
should not involve any Python object (this is not formally true but a good advice to follow). Moreover, in this particular case interruption should work perfectly fine without sig_on/sig_off
.
The only cases where you should use it is
 in a portion of your code where you use only C objects and that can potentially takes time
cdef int n sig_on() while some_condition: n = n + 3 sig_off()
 a call to a C function that can be long
my_arg1 = my_py_object._my_C_attribute1 my_arg2 = my_py_object._my_C_attribute2 sig_on() my_c_function(my_arg1, my_arg2) sig_off()
comment:61 Changed 2 years ago by
 Status changed from needs_review to needs_work
comment:62 Changed 2 years ago by
 Branch changed from u/bruno/faster_roots_computation_for_sparse_polynomials_over_zz to public/polynomials/faster_roots_sparse_zz16516
 Commit changed from df27d498e2f52fc0c3d902f5cca73e235188f923 to 22900398ba60e1cb1a2fec3131518c191dc6e3ac
 Reviewers set to Vincent Delecroix, Travis Scrimshaw
 Status changed from needs_work to needs_review
Agg...trac logged me out before I posted my message. Now I don't want to rewrite my long message.
Short version:
 I removed
sig_on
/sig_off
everywhere because it was calling Python functions (which should by C/Cython functions).  I fixed a bug when it was fully dense but had a positive valuation (see added doctest).
 I reorganized the ordering of the code so special cases are more quickly handled.
 Got some speed out of the generic polynomial methods.
 Used
self._parent
inPolynomial
instead ofself.parent()
because the former is relatively much faster. We can't do this in the special Python methods because of how Cython handles them.
New commits:
f9f4cea  Merge branch 'u/bruno/faster_roots_computation_for_sparse_polynomials_over_zz' of git://trac.sagemath.org/sage into public/polynomials/faster_roots_sparse_zz16516

16f7c09  It is faster to check hasattr than to catch the exception.

e7c0ad3  Remove sig_on()/sig_off() because calling Python functions.

9b9bd39  Use parent(p) instead of p.parent().

6b7116b  Speed improvements for generic sparse polynomials.

aa6ffb9  Using self._parent instead of self.parent().

2290039  Doing some reorganization and getting some more speed out of the factorizing.

comment:63 Changed 2 years ago by
 # The dense algorithm is to compute the roots from the factorization + # The dense algorithm is to compute the roots from the factorization. + # It is faster to let the factorization take care of the content if algorithm == "dense": return p._roots_from_factorization(p.factor(), multiplicities)
That is not true! The problem is that the factorization factors the content, which may take a very long time of course. As a consequence, there is a failing test in src/sage/rings/polynomial/polynomial_modn_dense_ntl.pyx
: The function small_roots
uses root finding over ZZ
at some point and the test at line 475 hangs because of the factorization of the content (which is a 514bit RSA integer).
By the way, thanks for cleaning my code!
comment:64 Changed 2 years ago by
 Commit changed from 22900398ba60e1cb1a2fec3131518c191dc6e3ac to 74dd4b2a2eeb879f74e8f170a92d90aacb496896
Branch pushed to git repo; I updated commit sha1. New commits:
74dd4b2  16516: Remove content before factorization

comment:65 Changed 2 years ago by
I added a few lines to remove the content before calling _roots_from_factorization
in the dense algorithm. This adds some code duplication but keeps this content computation at the latest. I'll wait the patchbots verdicts but the tests pass AFAICT.
comment:66 Changed 2 years ago by
Then I am ready to set this to a positive review if there are no other objections.
comment:67 Changed 2 years ago by
 Commit changed from 74dd4b2a2eeb879f74e8f170a92d90aacb496896 to 686d2b366a2726e78de271e379f0e49209fb7c8e
comment:68 Changed 2 years ago by
 Reviewers changed from Vincent Delecroix, Travis Scrimshaw to Vincent Delecroix, Travis Scrimshaw, Jeroen Demeyer
 Status changed from needs_review to positive_review
My last commit is a merge and removing part of the comment that wasn't true. Since the patchbot was (essentially) green on the previous branch and my change was trivial, I'm setting this to a positive review since there were no objections.
comment:69 Changed 2 years ago by
 Branch changed from public/polynomials/faster_roots_sparse_zz16516 to 686d2b366a2726e78de271e379f0e49209fb7c8e
 Resolution set to fixed
 Status changed from positive_review to closed
Branch pushed to git repo; I updated commit sha1. New commits:
New function for computing roots of sparse polynomials over ZZ