Skip to content

Instantly share code, notes, and snippets.

@cswiercz
Created January 11, 2016 17:13
Show Gist options
  • Save cswiercz/c136eda374d04a9d2172 to your computer and use it in GitHub Desktop.
Save cswiercz/c136eda374d04a9d2172 to your computer and use it in GitHub Desktop.
Pure-Sage Integral Basis timings
f = -y^8 + x^6 + 2*x^5
Fri Jan 8 11:53:49 2016 test.profile
10111400 function calls (10092038 primitive calls) in 18.657 seconds
Ordered by: internal time
List reduced from 1618 to 20 due to restriction <20>
ncalls tottime percall cumtime percall filename:lineno(function)
23382 1.540 0.000 4.408 0.000 qqbar.py:7433(_interval_fast)
161795/161708 1.430 0.000 1.631 0.000 complex_interval_field.py:364(__call__)
24650 1.230 0.000 2.266 0.000 {method '_rational_' of 'sage.rings.number_field.number_field_element.NumberFieldElement' objects}
101293 0.973 0.000 2.212 0.000 polynomial_ring.py:318(_element_constructor_)
616/280 0.787 0.001 1.033 0.004 qqbar.py:2444(super_poly)
2749578 0.674 0.000 0.674 0.000 {isinstance}
57756 0.607 0.000 1.239 0.000 multi_polynomial_ring.py:184(__call__)
58491 0.546 0.000 0.564 0.000 {method '_coerce_' of 'sage.structure.parent_old.Parent' objects}
1648 0.492 0.000 5.643 0.003 puiseux_series_ring_element.py:291(_add_)
158911 0.392 0.000 5.264 0.000 qqbar.py:3151(__init__)
50759 0.381 0.000 1.318 0.000 rings.py:708(__getitem__)
23662 0.377 0.000 1.246 0.000 {method 'polynomial' of 'sage.rings.number_field.number_field_element.NumberFieldElement' objects}
1504 0.359 0.000 3.669 0.002 puiseux_series_ring_element.py:318(_mul_)
51054 0.344 0.000 0.802 0.000 polynomial_ring_constructor.py:47(PolynomialRing)
16163 0.332 0.000 1.322 0.000 laurent_series_ring.py:273(_element_constructor_)
18622 0.291 0.000 2.664 0.000 qqbar.py:2778(an_muldiv_element)
8364 0.267 0.000 2.980 0.000 {method 'laurent_polynomial' of 'sage.rings.laurent_series_ring_element.LaurentSeries' objects}
36347 0.242 0.000 0.297 0.000 number_field.py:6660(_coerce_non_number_field_element_in)
115705/112301 0.240 0.000 1.066 0.000 qqbar.py:4158(__eq__)
24235 0.239 0.000 0.928 0.000 power_series_ring.py:647(_element_constructor_)
Fri Jan 8 11:53:49 2016 test.profile
10111400 function calls (10092038 primitive calls) in 18.657 seconds
Ordered by: cumulative time
List reduced from 1618 to 20 due to restriction <20>
ncalls tottime percall cumtime percall filename:lineno(function)
1 0.001 0.001 18.664 18.664 integralbasis.py:171(integral_basis)
1 0.001 0.001 18.656 18.656 integralbasis.py:217(_integral_basis_monic)
7 0.033 0.005 15.365 2.195 integralbasis.py:265(compute_bd)
1072/688 0.012 0.000 9.877 0.014 {sum}
736 0.071 0.000 7.186 0.010 integralbasis.py:302(<genexpr>)
49657 0.126 0.000 6.819 0.000 qqbar.py:3291(_mul_)
1648 0.492 0.000 5.643 0.003 puiseux_series_ring_element.py:291(_add_)
37432 0.204 0.000 5.448 0.000 multi_polynomial_element.py:220(_mul_)
150897 0.128 0.000 5.312 0.000 qqbar.py:4018(__init__)
158911 0.392 0.000 5.264 0.000 qqbar.py:3151(__init__)
23382 1.540 0.000 4.408 0.000 qqbar.py:7433(_interval_fast)
1504 0.359 0.000 3.669 0.002 puiseux_series_ring_element.py:318(_mul_)
8360 0.042 0.000 3.326 0.000 puiseux_series_ring_element.py:143(__init__)
2 0.000 0.000 3.276 1.638 integralbasis.py:129(compute_series_truncations)
8364 0.267 0.000 2.980 0.000 {method 'laurent_polynomial' of 'sage.rings.laurent_series_ring_element.LaurentSeries' objects}
27631 0.115 0.000 2.718 0.000 qqbar.py:7136(__new__)
18622 0.291 0.000 2.664 0.000 qqbar.py:2778(an_muldiv_element)
24650 1.230 0.000 2.266 0.000 {method '_rational_' of 'sage.rings.number_field.number_field_element.NumberFieldElement' objects}
101293 0.973 0.000 2.212 0.000 polynomial_ring.py:318(_element_constructor_)
18 0.023 0.001 1.882 0.105 integralbasis.py:320(solve_coefficient_system)
b =
[1, y, y^2/x, y^3/x, y^4/x^2, y^5/x^3, y^6/x^3, y^7/x^4]
@cswiercz
Copy link
Author

The above is integral basis profile information in a pure-Sage version on the example f = -y**8 + x**6 + 2*x**5. Some observations:

  • a lot of time is spent doing Puiseux series arithmetic
  • much of that time is done in computing sum(a[j]*b[j](xx,rxi) for j in range(d). Note that in the master branch this is done using a separate, cached function.

Cythonizing Puiseux series seems to be a natural next step.

@cswiercz
Copy link
Author

After caching the results of sum(a[j]*b[j](xx,rxi) for j in range(d) into an evaluate_A function (similar to the Sympy-version of the code) we get the following timings. It's about a 40% performance improvemnt.

Note that with caching enabled the bottleneck is now mostly within QQbar arithmetic. My guess is that Cythonization of Puiseux series will still help tremendously. Perhaps Cythonization of integral basis calculations would help, as well, but I doubt it.

f = -y^8 + x^6 + 2*x^5

Fri Jan  8 13:00:54 2016    test.profile

         5963719 function calls (5941535 primitive calls) in 11.089 seconds

   Ordered by: internal time
   List reduced from 1257 to 20 due to restriction <20>

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
148855/148769    1.067    0.000    1.220    0.000 complex_interval_field.py:364(__call__)
    12467    0.721    0.000    2.234    0.000 qqbar.py:7433(_interval_fast)
    10805    0.528    0.000    0.971    0.000 {method '_rational_' of 'sage.rings.number_field.number_field_element.NumberFieldElement' objects}
 1017/329    0.521    0.001    0.714    0.002 qqbar.py:2444(super_poly)
    51314    0.501    0.000    0.946    0.000 polynomial_ring.py:318(_element_constructor_)
  1604883    0.403    0.000    0.403    0.000 {isinstance}
    31590    0.240    0.000    0.845    0.000 rings.py:708(__getitem__)
    99748    0.230    0.000    2.732    0.000 qqbar.py:3151(__init__)
    32005    0.224    0.000    0.519    0.000 polynomial_ring_constructor.py:47(PolynomialRing)
     1392    0.221    0.000    0.889    0.001 qqbar.py:6740(_complex_refine_interval)
    18946    0.215    0.000    0.451    0.000 multi_polynomial_ring.py:184(__call__)
       18    0.197    0.011    0.197    0.011 {method 'nffactor' of 'sage.libs.pari.gen.gen_auto' objects}
    12796    0.194    0.000    0.668    0.000 {method 'polynomial' of 'sage.rings.number_field.number_field_element.NumberFieldElement' objects}
    19680    0.185    0.000    0.191    0.000 {method '_coerce_' of 'sage.structure.parent_old.Parent' objects}
       14    0.173    0.012    0.173    0.012 {method 'polredbest' of 'sage.libs.pari.gen.gen_auto' objects}
    21885    0.173    0.000    0.206    0.000 number_field.py:6660(_coerce_non_number_field_element_in)
     3861    0.146    0.000    0.316    0.000 {method '_rational_' of 'sage.rings.number_field.number_field_element_quadratic.NumberFieldElement_quadratic' objects}
     1362    0.145    0.000    0.199    0.000 {sage.rings.polynomial.refine_root.refine_root}
      432    0.140    0.000    1.738    0.004 puiseux_series_ring_element.py:291(_add_)
      288    0.118    0.000    0.156    0.001 {method 'minpoly' of 'sage.rings.number_field.number_field_element.NumberFieldElement_absolute' objects}


Fri Jan  8 13:00:54 2016    test.profile

         5963719 function calls (5941535 primitive calls) in 11.089 seconds

   Ordered by: cumulative time
   List reduced from 1257 to 20 due to restriction <20>

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.001    0.001   11.098   11.098 integralbasis.py:171(integral_basis)
        1    0.001    0.001   11.096   11.096 integralbasis.py:217(_integral_basis_monic)
        7    0.048    0.007    8.197    1.171 integralbasis.py:265(compute_bd)
    32592    0.076    0.000    2.915    0.000 qqbar.py:3291(_mul_)
        2    0.000    0.000    2.893    1.446 integralbasis.py:129(compute_series_truncations)
    99748    0.230    0.000    2.732    0.000 qqbar.py:3151(__init__)
    69848    0.062    0.000    2.605    0.000 qqbar.py:4018(__init__)
     5270    0.019    0.000    2.279    0.000 qqbar.py:4042(__cmp__)
    12467    0.721    0.000    2.234    0.000 qqbar.py:7433(_interval_fast)
45553/42050    0.012    0.000    2.147    0.000 {cmp}
       18    0.024    0.001    2.016    0.112 integralbasis.py:331(solve_coefficient_system)
      144    0.050    0.000    1.953    0.014 puiseux_series_ring_element.py:359(_cmp_)
    11325    0.065    0.000    1.813    0.000 multi_polynomial_element.py:220(_mul_)
      432    0.140    0.000    1.738    0.004 puiseux_series_ring_element.py:291(_add_)
 2491/889    0.004    0.000    1.638    0.002 qqbar.py:3678(exactify)
    12140    0.040    0.000    1.612    0.000 multi_polynomial_element.py:159(__cmp__)
    12140    0.063    0.000    1.554    0.000 {method 'compare' of 'sage.rings.polynomial.polydict.PolyDict' objects}
     4004    0.021    0.000    1.545    0.000 qqbar.py:3379(__hash__)
    14666    0.095    0.000    1.518    0.000 qqbar.py:7136(__new__)
  652/326    0.048    0.000    1.512    0.005 {method 'roots' of 'sage.rings.polynomial.polynomial_element.Polynomial' objects}


b =
[1, y, y^2/x, y^3/x, y^4/x^2, y^5/x^3, y^6/x^3, y^7/x^4]

@cswiercz
Copy link
Author

Python Puiseux series timings:

sage: P.<x> = PuiseuxSeriesRing(QQbar, 'x')
sage: L.<t> = LaurentSeriesRing(QQbar, 't')
sage: p = 2*x**(-2) + 1 + x + 2*x**2 + 5*x**5
sage: l = 2*t**(-2) + 1 + t + 2*t**2 + 5*t**5
sage: timeit('p+p')
125 loops, best of 3: 1.52 ms per loop
sage: timeit('l+l')
625 loops, best of 3: 52.6 µs per loop
sage: timeit('p*p')
125 loops, best of 3: 2.29 ms per loop
sage: timeit('l*l')
625 loops, best of 3: 649 µs per loop
sage: timeit('p**3')
125 loops, best of 3: 1.7 ms per loop
sage: timeit('l**3')
125 loops, best of 3: 1.7 ms per loop

@cswiercz
Copy link
Author

Timings post Cythonization (first round):

f = -y^8 + x^6 + 2*x^5

Fri Jan 15 09:28:15 2016    test.profile

         3825380 function calls (3808255 primitive calls) in 7.606 seconds

   Ordered by: internal time
   List reduced from 1239 to 20 due to restriction <20>

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
107410/107324    0.746    0.000    0.852    0.000 complex_interval_field.py:364(__call__)
  881/261    0.577    0.001    0.789    0.003 qqbar.py:2444(super_poly)
     6539    0.322    0.000    1.094    0.000 qqbar.py:7433(_interval_fast)
  1008583    0.273    0.000    0.273    0.000 {isinstance}
    25858    0.264    0.000    0.428    0.000 polynomial_ring.py:318(_element_constructor_)
     1392    0.247    0.000    0.995    0.001 qqbar.py:6740(_complex_refine_interval)
       18    0.210    0.012    0.210    0.012 {method 'nffactor' of 'sage.libs.pari.gen.gen_auto' objects}
       14    0.190    0.014    0.190    0.014 {method 'polredbest' of 'sage.libs.pari.gen.gen_auto' objects}
     3341    0.182    0.000    0.336    0.000 {method '_rational_' of 'sage.rings.number_field.number_field_element.NumberFieldElement' objects}
    67904    0.168    0.000    1.543    0.000 qqbar.py:3151(__init__)
     3861    0.162    0.000    0.351    0.000 {method '_rational_' of 'sage.rings.number_field.number_field_element_quadratic.NumberFieldElement_quadratic' objects}
    19518    0.161    0.000    0.580    0.000 rings.py:708(__getitem__)
     1362    0.156    0.000    0.204    0.000 {sage.rings.polynomial.refine_root.refine_root}
    15045    0.153    0.000    0.181    0.000 number_field.py:6660(_coerce_non_number_field_element_in)
    19933    0.149    0.000    0.365    0.000 polynomial_ring_constructor.py:47(PolynomialRing)
        7    0.143    0.020    5.342    0.763 integralbasis.py:264(compute_bd)
11095/8133    0.111    0.000    0.119    0.000 {hash}
     7828    0.102    0.000    0.251    0.000 multi_polynomial_ring.py:184(__call__)
     6800    0.097    0.000    0.355    0.000 {method 'polynomial' of 'sage.rings.number_field.number_field_element.NumberFieldElement' objects}
     7202    0.081    0.000    0.913    0.000 qqbar.py:7136(__new__)


Fri Jan 15 09:28:15 2016    test.profile

         3825380 function calls (3808255 primitive calls) in 7.606 seconds

   Ordered by: cumulative time
   List reduced from 1239 to 20 due to restriction <20>

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.001    0.001    7.616    7.616 integralbasis.py:170(integral_basis)
        1    0.001    0.001    7.614    7.614 integralbasis.py:216(_integral_basis_monic)
        7    0.143    0.020    5.342    0.763 integralbasis.py:264(compute_bd)
       18    0.027    0.001    2.277    0.127 integralbasis.py:330(solve_coefficient_system)
        2    0.000    0.000    2.266    1.133 integralbasis.py:128(compute_series_truncations)
     2134    0.015    0.000    2.238    0.001 qqbar.py:4042(__cmp__)
 1759/565    0.003    0.000    1.714    0.003 qqbar.py:3678(exactify)
  652/326    0.055    0.000    1.705    0.005 {method 'roots' of 'sage.rings.polynomial.polynomial_element.Polynomial' objects}
     4004    0.023    0.000    1.705    0.000 qqbar.py:3379(__hash__)
      326    0.021    0.000    1.689    0.005 complex_roots.py:147(complex_roots)
  418/225    0.010    0.000    1.596    0.007 qqbar.py:8072(exactify)
    67904    0.168    0.000    1.543    0.000 qqbar.py:3151(__init__)
     2788    0.011    0.000    1.539    0.001 multi_polynomial_element.py:159(__cmp__)
     2788    0.019    0.000    1.523    0.001 {method 'compare' of 'sage.rings.polynomial.polydict.PolyDict' objects}
        4    0.028    0.007    1.436    0.359 puiseux.py:746(xseries)
    38276    0.034    0.000    1.368    0.000 qqbar.py:4018(__init__)
     8457    0.027    0.000    1.148    0.000 qqbar.py:3327(_add_)
     6539    0.322    0.000    1.094    0.000 qqbar.py:7433(_interval_fast)
  342/332    0.001    0.000    1.033    0.003 {map}
       32    0.000    0.000    1.025    0.032 puiseux.py:792(<lambda>)


b =
[1, y, y^2/x, y^3/x, y^4/x^2, y^5/x^3, y^6/x^3, y^7/x^4]

@cswiercz
Copy link
Author

Here are some timings from abelfunctions.differentials_numerators. I thought that most of the performance issues were with the use of Groeber bases but it seems that half of the time spent is still in integral_basis. Darn...

f = 2*x^3*y^5 + x^2*y^6 - 1

Wed Jan 20 14:26:52 2016    test.profile

         11629538 function calls (11531472 primitive calls) in 19.429 seconds

   Ordered by: internal time
   List reduced from 1813 to 20 due to restriction <20>

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
 1441/513    1.263    0.001    1.681    0.003 qqbar.py:2444(super_poly)
157171/156892    1.204    0.000    1.364    0.000 complex_interval_field.py:364(__call__)
    20696    0.943    0.000    3.142    0.000 qqbar.py:7433(_interval_fast)
  2908264    0.715    0.000    0.715    0.000 {isinstance}
    77239    0.704    0.000    1.270    0.000 polynomial_ring.py:318(_element_constructor_)
45997/44379    0.496    0.000    1.518    0.000 multi_polynomial_ring.py:184(__call__)
    12749    0.488    0.000    1.058    0.000 {method '_rational_' of 'sage.rings.number_field.number_field_element_quadratic.NumberFieldElement_quadratic' objects}
     9122    0.483    0.000    0.894    0.000 {method '_rational_' of 'sage.rings.number_field.number_field_element.NumberFieldElement' objects}
   216072    0.481    0.000    4.364    0.000 qqbar.py:3151(__init__)
     2769    0.426    0.000    4.501    0.002 multi_polynomial_element.py:105(__call__)
    56625    0.421    0.000    1.600    0.000 rings.py:708(__getitem__)
    59423    0.401    0.000    1.977    0.000 polynomial_ring_constructor.py:47(PolynomialRing)
       37    0.350    0.009    0.350    0.009 {method 'nffactor' of 'sage.libs.pari.gen.gen_auto' objects}
35766/26445    0.335    0.000    0.363    0.000 {hash}
     3912    0.323    0.000    0.327    0.000 dynamic_class.py:328(dynamic_class_internal)
45746/45716    0.317    0.000    0.456    0.000 {method '_coerce_' of 'sage.structure.parent_old.Parent' objects}
       31    0.312    0.010    0.312    0.010 {method 'polredbest' of 'sage.libs.pari.gen.gen_auto' objects}
    33282    0.309    0.000    0.366    0.000 number_field.py:6660(_coerce_non_number_field_element_in)
    21209    0.276    0.000    1.007    0.000 {method 'polynomial' of 'sage.rings.number_field.number_field_element.NumberFieldElement' objects}
115997/115267    0.269    0.000    0.606    0.000 multi_polynomial_element.py:313(__init__)


Wed Jan 20 14:26:52 2016    test.profile

         11629538 function calls (11531472 primitive calls) in 19.429 seconds

   Ordered by: cumulative time
   List reduced from 1813 to 20 due to restriction <20>

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.006    0.006   19.447   19.447 differentials.py:149(differentials_numerators)
        2    0.000    0.000   16.673    8.336 integralbasis.py:170(integral_basis)
        2    0.001    0.000   16.663    8.331 integralbasis.py:216(_integral_basis_monic)
       14    0.118    0.008   12.071    0.862 integralbasis.py:264(compute_bd)
    13029    0.069    0.000    4.900    0.000 qqbar.py:3379(__hash__)
     2769    0.426    0.000    4.501    0.002 multi_polynomial_element.py:105(__call__)
   216072    0.481    0.000    4.364    0.000 qqbar.py:3151(__init__)
   184851    0.136    0.000    4.168    0.000 qqbar.py:4018(__init__)
        4    0.001    0.000    3.866    0.966 integralbasis.py:128(compute_series_truncations)
2881/2341    0.084    0.000    3.863    0.002 {sum}
      224    0.001    0.000    3.823    0.017 integralbasis.py:325(evaluate_A)
       33    0.044    0.001    3.727    0.113 integralbasis.py:330(solve_coefficient_system)
     1120    0.163    0.000    3.406    0.003 integralbasis.py:327(<genexpr>)
    23866    0.075    0.000    3.406    0.000 qqbar.py:3327(_add_)
4212/1294    0.007    0.000    3.299    0.003 qqbar.py:3678(exactify)
    20696    0.943    0.000    3.142    0.000 qqbar.py:7433(_interval_fast)
  844/455    0.018    0.000    3.040    0.007 qqbar.py:8072(exactify)
    21871    0.241    0.000    2.632    0.000 qqbar.py:7136(__new__)
       16    0.055    0.003    2.584    0.162 puiseux.py:746(xseries)
    47206    0.094    0.000    2.164    0.000 qqbar.py:3291(_mul_)


d =
[y, x*y^3, x*y^4]

@cswiercz
Copy link
Author

All of this is moot, for the moment, because I found out that Singular implements an integral basis algorithm over Q[x,y] that is (a) fast and (b) works.

Based on van Hoeij's paper it seems like his algorithm allows for arbitrary (at least, algebraic) extensions of Q as base rings. So there is potentially good reason to keep this version around.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment