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: sage-8.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 bruno)

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)

dense_pol_test.txt (4.1 KB) - added by bruno 5 years ago.
dense_pol_test.2.txt (4.1 KB) - added by bruno 5 years ago.
dense_pol_no_test.txt (4.3 KB) - added by bruno 5 years ago.
sparse_pol_no_test.txt (4.3 KB) - added by bruno 5 years ago.
sparse_pol_test.txt (4.1 KB) - added by bruno 5 years ago.

Download all attachments as: .zip

Change History (74)

comment:1 Changed 5 years ago by bruno

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

  • Branch set to u/bruno/faster_roots_computation_for_sparse_polynomials_over_zz

comment:3 Changed 5 years ago by bruno

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

  • Commit changed from 920c16e8952cb1f885c1fd1f3bcdfd5d572c17fe to e725544b76363ea517a262419a092f9414c1b4a4

Branch pushed to git repo; I updated commit sha1. New commits:

e725544New function for computing roots of sparse polynomials over ZZ

comment:5 Changed 5 years ago by vdelecroix

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 (the is 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 git

  • Commit changed from e725544b76363ea517a262419a092f9414c1b4a4 to f344cdc5b4a9e316a52f5372a30a3db1ac5043ab

Branch pushed to git repo; I updated commit sha1. New commits:

f344cdcCosmestic changes + add conditions on the constant coefficient

comment:7 Changed 5 years ago by bruno

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 high-degree dense polynomial, integer roots of p are common roots of 3-x 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 low-degree 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 lowest-degree 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 (x-i)*(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 bruno

  • Authors set to Bruno Grenet

comment:9 Changed 5 years ago by vdelecroix

  • 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 python-is-operator 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[i-1] > 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[k-1] ].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(pn)).

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 vdelecroix

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

Last edited 5 years ago by vdelecroix (previous) (diff)

comment:11 Changed 5 years ago by vdelecroix

For the content issue, I opened a thread on sage-devel.

comment:12 Changed 5 years ago by git

  • Commit changed from f344cdc5b4a9e316a52f5372a30a3db1ac5043ab to be6deca10dd8d25ace0e9d821c6cc5551ca8d47f

Branch pushed to git repo; I updated commit sha1. New commits:

7fec476Slighlty improved the algorithm + add a algorithm="sparse" option
e89350eChange sparsity function for dense polys + add doctests
5c4c70fDeclaration of int variables
be6decaclean up _roots_univariate_polynomial

comment:13 Changed 5 years ago by bruno

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.hpi-web.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 bruno

  • Status changed from needs_work to needs_review

comment:15 Changed 5 years ago by vbraun_spam

  • Milestone changed from sage-6.3 to sage-6.4

comment:16 Changed 5 years ago by vdelecroix

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 git

  • Commit changed from be6deca10dd8d25ace0e9d821c6cc5551ca8d47f to 0efc4ea81dc2c47bb9b0246f37424e296bb8ee42

Branch pushed to git repo; I updated commit sha1. New commits:

e4c386emerge 6.3
0efc4eatrac #16516: cosmetic

comment:18 Changed 5 years ago by bruno

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) (files sparse_pol_...) and PolynomialRing(ZZ,sparse=False) (files dense_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):

  1. the degree of the polynomials (here 2);
  2. list of 10 ratios corresponding to 10 independent tests;
  3. minimum of the list;
  4. average of the list;
  5. 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 bruno

Changed 5 years ago by bruno

Changed 5 years ago by bruno

Changed 5 years ago by bruno

Changed 5 years ago by bruno

comment:19 Changed 5 years ago by bruno

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 follow-up: Changed 5 years ago by vdelecroix

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 non-zero. 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 bruno

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 low-degree monomial has a zero coefficient: In such case, the sparse algorithm splits the input polynomial into a low-degree 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 speed-up 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
Last edited 5 years ago by bruno (previous) (diff)

comment:22 follow-up: Changed 5 years ago by jdemeyer

Wouldn't "density" be a better word than "sparsity"?

comment:23 Changed 5 years ago by jdemeyer

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".

Last edited 5 years ago by jdemeyer (previous) (diff)

comment:24 Changed 5 years ago by jdemeyer

  • Status changed from needs_review to needs_work

comment:25 Changed 5 years ago by git

  • Commit changed from 0efc4ea81dc2c47bb9b0246f37424e296bb8ee42 to 711e73c530ce0e56beb99fc2cddd95e3901d7385

Branch pushed to git repo; I updated commit sha1. New commits:

711e73cAutomatic choice of algorithm + replace "return None" by an exception + doc

comment:26 in reply to: ↑ 22 ; follow-up: Changed 5 years ago by bruno

  • 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 follow-up: Changed 5 years ago by jdemeyer

  • Status changed from needs_review to needs_work
  1. Remove commented-out code: #or p.degree() - p.sparsity() -p.valuation() > 10:
  2. Add doctests for the various algorithm= parameters, including algorithm="foobar"
  3. Remove if hasattr(K, '_roots_univariate_polynomial'): and change except NotImplementedError to except (AttributeError, NotImplementedError):

comment:28 in reply to: ↑ 26 Changed 5 years ago by jdemeyer

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 bruno

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:

  1. Remove commented-out code: #or p.degree() - p.sparsity() -p.valuation() > 10:

Actually, I'd like to de-comment 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 low-degree polynomials where the computation is very fast.

  1. Add doctests for the various algorithm= parameters, including algorithm="foobar"
  2. Remove if hasattr(K, '_roots_univariate_polynomial'): and change except NotImplementedError to except (AttributeError, NotImplementedError):

I'll do that.

comment:30 Changed 5 years ago by git

  • Commit changed from 711e73c530ce0e56beb99fc2cddd95e3901d7385 to 4b15fc26b7ea40fcc9a897c1875fde28726a38ef

Branch pushed to git repo; I updated commit sha1. New commits:

89a867bReplace sparsity by number_of_terms
fccdd25remove has_attr and replaced by an except AttributeError
2f8b2a0By default, sparse algorithm
4b15fc2doctest possible values for algorithm=...

comment:31 Changed 5 years ago by jdemeyer

Why overcomplicate your doctest (see last commit) with a try/except? You can simply doctest the exception...

comment:32 Changed 5 years ago by bruno

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?

Last edited 5 years ago by bruno (previous) (diff)

comment:33 Changed 5 years ago by git

  • Commit changed from 4b15fc26b7ea40fcc9a897c1875fde28726a38ef to 26a3fbb893bc728964e20e78e39e18156720149d

Branch pushed to git repo; I updated commit sha1. New commits:

26a3fbbbetter doctest

comment:34 Changed 5 years ago by bruno

  • Status changed from needs_work to needs_review

Ok, I found how one is supposed to doctest exceptions, sorry not to have checked this before!

comment:35 Changed 5 years ago by jdemeyer

The sentence "Default is based on whether p.parent() is sparse or not." should be fixed then.

comment:36 Changed 5 years ago by git

  • Commit changed from 26a3fbb893bc728964e20e78e39e18156720149d to 52fdc7876281d9086057d1cf985f7ae2c21df095

Branch pushed to git repo; I updated commit sha1. New commits:

52fdc78Update doc

comment:37 Changed 5 years ago by jdemeyer

The code cannot be interrupted. Please have a look at http://www.sagemath.org/doc/developer/coding_in_cython.html#using-sig-check (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 follow-up: Changed 5 years ago by jdemeyer

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

Note: in Python, you should use != instead of <>

comment:40 in reply to: ↑ 38 Changed 5 years ago by bruno

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 git

  • Commit changed from 52fdc7876281d9086057d1cf985f7ae2c21df095 to d72a8f1554235a389d8c42e278d3779895d568b5

Branch pushed to git repo; I updated commit sha1. New commits:

eb64fb5Add sig_on() & sig_off()'s
a056a34Replace <> by !=
0ddbcf2Remove test on cst coeff
8b029c5Add shortcut for totally dense polynomials
d72a8f1Automatic choice of alg: sparse for degree >100, dense otherwise

comment:42 Changed 4 years ago by bruno

  • 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 (x-1)^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 slow-down that exists in some examples using the "sparse" algorithm rather than the "dense" algorithm becomes negligible from around degree 100, with never more that 10% slow down, and almost always <1% (and of course quite often there is a speed-up!). Note that I also add a shortcut for completely dense polynomials (with density 1) 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 git

  • Commit changed from d72a8f1554235a389d8c42e278d3779895d568b5 to f97f22bd0b549ff3f0e9406ab10e4ce989e0f6cd

Branch pushed to git repo; I updated commit sha1. New commits:

f97f22bMerge 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 chapoton

  • Status changed from needs_review to needs_work

4 failings doctests, see patchbot report

comment:45 Changed 4 years ago by git

  • Commit changed from f97f22bd0b549ff3f0e9406ab10e4ce989e0f6cd to e9b274f3c413f21738ad8d231323911ac5087f4b

Branch pushed to git repo; I updated commit sha1. New commits:

5516ab1Merge branch 'u/bruno/faster_roots_computation_for_sparse_polynomials_over_zz' of git://trac.sagemath.org/sage into develop
e9b274fReplace zero_element by zero

comment:46 Changed 4 years ago by bruno

  • Status changed from needs_work to needs_review

Should be OK now.

comment:47 follow-up: Changed 4 years ago by chapoton

  • Status changed from needs_review to needs_work

one failing doctest

comment:48 Changed 4 years ago by git

  • Commit changed from e9b274f3c413f21738ad8d231323911ac5087f4b to 80be713ae2673ca1d769cc5f0ceb7b4a2153dd9c

Branch pushed to git repo; I updated commit sha1. New commits:

1b327dbMerge branch 'u/bruno/faster_roots_computation_for_sparse_polynomials_over_zz' of git://trac.sagemath.org/sage into t16516_roots_sparse_polys
80be713Add missing sig_off()

comment:49 in reply to: ↑ 47 Changed 4 years ago by bruno

  • Status changed from needs_work to needs_review

Replying to chapoton:

one failing doctest

One sig_off() was missing, fixed.

comment:50 follow-up: Changed 3 years ago by chapoton

  • Status changed from needs_review to needs_work

1 failing doctest, due to deprecation of polynomial slicing

comment:51 Changed 3 years ago by git

  • Commit changed from 80be713ae2673ca1d769cc5f0ceb7b4a2153dd9c to d47a525a9c316d1b84c3323fb11c46a88803cafb

Branch pushed to git repo; I updated commit sha1. This was a forced push. Last 10 new commits:

763a898trac #16516: cosmetic
f7b499816516: Automatic choice of algorithm + replace "return None" by an exception + doc
d06e3e516516: remove has_attr and replaced by an except AttributeError
971f70e16516: doctest possible values for algorithm=...
865929916516: Add sig_on() & sig_off()'s
8554f3e16516: Replace <> by !=
be39af016516: 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 bruno

  • 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:53 Changed 3 years ago by chapoton

  • Status changed from needs_review to needs_work

no longer applies

comment:54 Changed 2 years ago by git

  • Commit changed from d47a525a9c316d1b84c3323fb11c46a88803cafb to 48d4fc9e73c4e227711e6a656782791489eeef37

Branch pushed to git repo; I updated commit sha1. New commits:

48d4fc9Merge branch 'develop' into 16516_faster_roots_ZZ

comment:55 Changed 2 years ago by bruno

  • Milestone changed from sage-6.4 to sage-8.0
  • Status changed from needs_work to needs_review

Applies again, needs review again!

comment:56 follow-ups: Changed 2 years ago by tscrim

Minor comments:

  • Mark the docstring with r""" and
    -Returns 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 vdelecroix

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 git

  • Commit changed from 48d4fc9e73c4e227711e6a656782791489eeef37 to df27d498e2f52fc0c3d902f5cca73e235188f923

Branch pushed to git repo; I updated commit sha1. New commits:

1604cc916516: Move reference to global reference file
080a61716516: Better docstring + cosmetic change
6bdb88f16516: Better use of sig_on/sig_off
2a96eb216516: xrange -> range
df27d4916516: Better primitive part

comment:59 in reply to: ↑ 56 Changed 2 years ago by bruno

Replying to tscrim:

Minor comments:

  • Mark the docstring with r""" and
    -Returns 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 Ctrl-C. There are still a fairly large number of them.

comment:60 Changed 2 years ago by vdelecroix

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 vdelecroix

  • Status changed from needs_review to needs_work

comment:62 Changed 2 years ago by tscrim

  • Branch changed from u/bruno/faster_roots_computation_for_sparse_polynomials_over_zz to public/polynomials/faster_roots_sparse_zz-16516
  • 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 in Polynomial instead of self.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:

f9f4ceaMerge branch 'u/bruno/faster_roots_computation_for_sparse_polynomials_over_zz' of git://trac.sagemath.org/sage into public/polynomials/faster_roots_sparse_zz-16516
16f7c09It is faster to check hasattr than to catch the exception.
e7c0ad3Remove sig_on()/sig_off() because calling Python functions.
9b9bd39Use parent(p) instead of p.parent().
6b7116bSpeed improvements for generic sparse polynomials.
aa6ffb9Using self._parent instead of self.parent().
2290039Doing some reorganization and getting some more speed out of the factorizing.

comment:63 Changed 2 years ago by bruno

-        # 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 514-bit RSA integer).

By the way, thanks for cleaning my code!

comment:64 Changed 2 years ago by git

  • Commit changed from 22900398ba60e1cb1a2fec3131518c191dc6e3ac to 74dd4b2a2eeb879f74e8f170a92d90aacb496896

Branch pushed to git repo; I updated commit sha1. New commits:

74dd4b216516: Remove content before factorization

comment:65 Changed 2 years ago by bruno

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.

Last edited 2 years ago by bruno (previous) (diff)

comment:66 Changed 2 years ago by tscrim

Then I am ready to set this to a positive review if there are no other objections.

comment:67 Changed 2 years ago by git

  • Commit changed from 74dd4b2a2eeb879f74e8f170a92d90aacb496896 to 686d2b366a2726e78de271e379f0e49209fb7c8e

Branch pushed to git repo; I updated commit sha1. New commits:

53edcc0Merge branch 'public/polynomials/faster_roots_sparse_zz-16516' of git://trac.sagemath.org/sage into public/polynomials/faster_roots_sparse_zz-16516
686d2b3Update comment to reflect the code change.

comment:68 Changed 2 years ago by tscrim

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

  • Branch changed from public/polynomials/faster_roots_sparse_zz-16516 to 686d2b366a2726e78de271e379f0e49209fb7c8e
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.