@@ -350,9 +350,9 @@ class RiemannSurface():
350
350
by the monomials of degree up to `d-3`.
351
351
352
352
- ``integration_method`` -- (default: ``'rigorous'``). String specifying the
353
- integration method to use when calculating the integrals of differentials.
353
+ integration method to use when calculating the integrals of differentials.
354
354
The options are ``'heuristic'`` and ``'rigorous'``, the latter of
355
- which is often the most efficient.
355
+ which is often the most efficient.
356
356
357
357
EXAMPLES::
358
358
@@ -386,19 +386,19 @@ class RiemannSurface():
386
386
sage: all( len(T.minpoly().roots(K)) > 0 for T in A)
387
387
True
388
388
389
- The ``'heuristic'`` integration method uses the method ``integrate_vector``
390
- defined in ``sage.numerical.gauss_legendre`` to compute integrals of differentials.
391
- As mentioned there, this works by iteratively doubling the number of nodes
389
+ The ``'heuristic'`` integration method uses the method ``integrate_vector``
390
+ defined in ``sage.numerical.gauss_legendre`` to compute integrals of differentials.
391
+ As mentioned there, this works by iteratively doubling the number of nodes
392
392
used in the quadrature, and uses a heuristic based on the rate at which the
393
393
result is seemingly converging to estimate the error. The ``'rigorous'``
394
- method uses results from [Neu2018]_, and bounds the algebraic integrands on
394
+ method uses results from [Neu2018]_, and bounds the algebraic integrands on
395
395
circular domains using Cauchy's form of the remainder in Taylor approximation
396
396
coupled to Fujiwara's bound on polynomial roots (see Bruin-DisneyHogg-Gao,
397
- in preparation). Note this method of bounding on circular domains is also
398
- implemented in :meth:`_compute_delta`. The net result of this bounding is
397
+ in preparation). Note this method of bounding on circular domains is also
398
+ implemented in :meth:`_compute_delta`. The net result of this bounding is
399
399
that one can know (an upper bound on) the number of nodes required to achieve
400
- a certain error. This means that for any given integral, assuming that the
401
- same number of nodes is required by both methods in order to achieve the
400
+ a certain error. This means that for any given integral, assuming that the
401
+ same number of nodes is required by both methods in order to achieve the
402
402
desired error (not necessarily true in practice), approximately half
403
403
the number of integrand evaluations are required. When the required number
404
404
of nodes is high, e.g. when the precision required is high, this can make
@@ -408,7 +408,7 @@ class RiemannSurface():
408
408
if the computation is 'fast', the heuristic method may outperform the
409
409
rigorous method, but for slower computations the rigorous method can be much
410
410
faster::
411
-
411
+
412
412
sage: f = z*w^3+z^3+w
413
413
sage: p = 53
414
414
sage: Sh = RiemannSurface(f, prec=p, integration_method='heuristic')
@@ -417,18 +417,18 @@ class RiemannSurface():
417
417
sage: import time
418
418
sage: nodes.cache.clear()
419
419
sage: ct = time.time()
420
- sage: Rh = Sh.riemann_matrix()
420
+ sage: Rh = Sh.riemann_matrix()
421
421
sage: ct1 = time.time()-ct
422
422
sage: nodes.cache.clear()
423
423
sage: ct = time.time()
424
- sage: Rr = Sr.riemann_matrix()
424
+ sage: Rr = Sr.riemann_matrix()
425
425
sage: ct2 = time.time()-ct
426
426
sage: ct2/ct1 # random
427
427
1.2429363969691192
428
428
429
429
This disparity in timings can get increasingly worse, and testing has shown
430
430
that even for random quadrics the heuristic method can be as bad as 30 times
431
- slower.
431
+ slower.
432
432
433
433
TESTS:
434
434
@@ -1549,27 +1549,27 @@ def _bounding_data(self, differentials):
1549
1549
INPUT:
1550
1550
1551
1551
- ``differentials`` -- list. A list of polynomials in ``self._R`` giving
1552
- the numerators of the differentials, as per the output of
1553
- :meth:`cohomology_basis`.
1552
+ the numerators of the differentials, as per the output of
1553
+ :meth:`cohomology_basis`.
1554
1554
1555
1555
OUTPUT:
1556
1556
1557
- A tuple ``(CCzg, [(g, dgdz, F, a0_info), ...])`` where each element of
1557
+ A tuple ``(CCzg, [(g, dgdz, F, a0_info), ...])`` where each element of
1558
1558
the list corresponds to an element of ``differentials``. Introducing the
1559
- notation ``RBzg = PolynomialRing(self._R, ['z','g'])`` and
1559
+ notation ``RBzg = PolynomialRing(self._R, ['z','g'])`` and
1560
1560
``CCzg = PolynomialRing(self._CC, ['z','g'])``, we have that:
1561
1561
1562
- - ``g`` is the full rational function in ``self._R.fraction_field()``
1562
+ - ``g`` is the full rational function in ``self._R.fraction_field()``
1563
1563
giving the differential,
1564
1564
- ``dgdz`` is the derivative of ``g`` with respect to ``self._R.gen(0)``
1565
- , written in terms of ``self._R.gen(0)`` and ``g``, hence laying in
1565
+ , written in terms of ``self._R.gen(0)`` and ``g``, hence laying in
1566
1566
``RBzg``,
1567
- - ``F`` is the minimal polynomial of ``g`` over ``self._R.gen(0)``,
1567
+ - ``F`` is the minimal polynomial of ``g`` over ``self._R.gen(0)``,
1568
1568
laying in the polynomial ring ``CCzg``,
1569
- - ``a0_info`` is a tuple ``(lc, roots)`` where ``lc`` and ``roots`` are
1569
+ - ``a0_info`` is a tuple ``(lc, roots)`` where ``lc`` and ``roots`` are
1570
1570
the leading coefficient and roots of the polynomial in ``CCzg.gen(0)``
1571
- that is the coefficient of the term of ``F`` of highest degree in
1572
- ``CCzg.gen(1)``.
1571
+ that is the coefficient of the term of ``F`` of highest degree in
1572
+ ``CCzg.gen(1)``.
1573
1573
1574
1574
EXAMPLES::
1575
1575
@@ -1592,8 +1592,8 @@ def _bounding_data(self, differentials):
1592
1592
-0.500000000000000 + 0.866025403784439*I]))])
1593
1593
1594
1594
"""
1595
- # This copies previous work by NB, outputting the zipped list required
1596
- # for a certified line integral.
1595
+ # This copies previous work by NB, outputting the zipped list required
1596
+ # for a certified line integral.
1597
1597
RB = self ._R .base_ring ()
1598
1598
P = PolynomialRing (RB , 'Z' )
1599
1599
k = P .fraction_field ()
@@ -1602,30 +1602,30 @@ def _bounding_data(self, differentials):
1602
1602
L = k .extension (fZW , 'Wb' )
1603
1603
dfdw_L = self ._dfdw (P .gen (0 ), L .gen (0 ))
1604
1604
integrand_list = [h / self ._dfdw for h in differentials ]
1605
- # minpoly_univ gives the minimal polynomial for h, in variable x, with
1605
+ # minpoly_univ gives the minimal polynomial for h, in variable x, with
1606
1606
# coefficients given by polynomials in P (i.e. rational polynomials in Z).
1607
1607
minpoly_univ = [(h (P .gen (0 ), L .gen (0 ))/ dfdw_L ).minpoly ().numerator ()
1608
1608
for h in differentials ]
1609
1609
RBzg = PolynomialRing (RB , ['z' , 'g' ])
1610
- # The following line changes the variables in these minimal polynomials
1610
+ # The following line changes the variables in these minimal polynomials
1611
1611
# as Z -> z, x -> G, then evaluates at G = QQzg.gens(1) ( = g )
1612
1612
RBzgG = PolynomialRing (RBzg , 'G' )
1613
1613
minpoly_list = [RBzgG ([c (RBzg .gen (0 )) for c in list (h )])(RBzg .gen (1 ))
1614
1614
for h in minpoly_univ ]
1615
1615
# h(z,g)=0 --> dg/dz = - dhdz/dhdg
1616
1616
dgdz_list = [- h .derivative (RBzg .gen (0 ))/ h .derivative (RBzg .gen (1 ))
1617
1617
for h in minpoly_list ]
1618
-
1618
+
1619
1619
CCzg = PolynomialRing (self ._CC , ['z' ,'g' ])
1620
1620
CCminpoly_list = [CCzg (h ) for h in minpoly_list ]
1621
-
1621
+
1622
1622
a0_list = [P (h .leading_coefficient ()) for h in minpoly_univ ]
1623
1623
# Note that because the field over which the Riemann surface is defined
1624
- # is embedded into CC, it has characteristic 0, and so we know the
1624
+ # is embedded into CC, it has characteristic 0, and so we know the
1625
1625
# irreducible factors are all separable, i.e. the roots have multiplicity
1626
- # one.
1626
+ # one.
1627
1627
a0_info = [(self ._CC (a0 .leading_coefficient ()),
1628
- flatten ([self ._CCz (F ).roots (multiplicities = False )* m
1628
+ flatten ([self ._CCz (F ).roots (multiplicities = False )* m
1629
1629
for F , m in a0 .factor ()]))
1630
1630
for a0 in a0_list ]
1631
1631
return CCzg , list (zip (integrand_list , dgdz_list , CCminpoly_list , a0_info ))
@@ -1635,10 +1635,10 @@ def rigorous_line_integral(self, upstairs_edge, differentials, bounding_data):
1635
1635
Perform vectorized integration along a straight path.
1636
1636
1637
1637
Using the error bounds for Gauss-Legendre integration found in [Neu2018]_
1638
- and a method for bounding an algebraic integrand on a circular domains
1639
- using Cauchy's form of the remainder in Taylor approximation coupled to
1640
- Fujiwara's bound on polynomial roots (see Bruin-DisneyHogg-Gao, in
1641
- preparation), this method calculates (semi-)rigorously the integral of a
1638
+ and a method for bounding an algebraic integrand on a circular domains
1639
+ using Cauchy's form of the remainder in Taylor approximation coupled to
1640
+ Fujiwara's bound on polynomial roots (see Bruin-DisneyHogg-Gao, in
1641
+ preparation), this method calculates (semi-)rigorously the integral of a
1642
1642
list of differentials along an edge of the upstairs graph.
1643
1643
1644
1644
INPUT:
@@ -1651,7 +1651,7 @@ def rigorous_line_integral(self, upstairs_edge, differentials, bounding_data):
1651
1651
the equation defining the Riemann surface.
1652
1652
1653
1653
- ``bounding_data`` -- tuple containing the data required for bounding
1654
- the integrands. This should be in the form of the output from
1654
+ the integrands. This should be in the form of the output from
1655
1655
:meth:`_bounding_data`.
1656
1656
1657
1657
OUTPUT:
@@ -1677,81 +1677,81 @@ def rigorous_line_integral(self, upstairs_edge, differentials, bounding_data):
1677
1677
1678
1678
.. NOTE::
1679
1679
1680
- Uses data that ``homology_basis`` initializes.
1680
+ Uses data that ``homology_basis`` initializes.
1681
1681
1682
1682
Note also that the data of the differentials is contained within
1683
- ``bounding_data``. It is, however, still advantageous to have this
1683
+ ``bounding_data``. It is, however, still advantageous to have this
1684
1684
be a separate argument, as it lets the user supply a fast-callable
1685
- version of the differentials, to significantly speed up execution
1686
- of the integrand calls, and not have to re-calculate these
1685
+ version of the differentials, to significantly speed up execution
1686
+ of the integrand calls, and not have to re-calculate these
1687
1687
fast-callables for every run of the function. This is also the benefit
1688
- of representing the differentials as a polynomial over a known
1689
- common denominator.
1688
+ of representing the differentials as a polynomial over a known
1689
+ common denominator.
1690
1690
1691
1691
.. TODO::
1692
1692
1693
1693
Note that bounding_data contains the information of the integrands,
1694
1694
so one may want to check for consistency between ``bounding_data``
1695
- and ``differentials``. If so one would not want to do so at the
1696
- expense of speed.
1695
+ and ``differentials``. If so one would not want to do so at the
1696
+ expense of speed.
1697
1697
1698
- Moreover, the current implementation bounds along a line by
1698
+ Moreover, the current implementation bounds along a line by
1699
1699
splitting it up into segments, each of which can be covered entirely
1700
- by a single circle, and then placing inside that the ellipse
1700
+ by a single circle, and then placing inside that the ellipse
1701
1701
required to bound as per [Neu2018]_. This is reliably more efficient
1702
- than the heuristic method, especially in poorly-conditioned cases
1702
+ than the heuristic method, especially in poorly-conditioned cases
1703
1703
where discriminant points are close together around the edges, but
1704
1704
in the case where the branch locus is well separated, it can require
1705
1705
slightly more nodes than necessary. One may want to include a method
1706
- here to transition in this regime to an algorithm that covers the
1707
- entire line with one ellipse, then bounds along that ellipse with
1708
- multiple circles.
1706
+ here to transition in this regime to an algorithm that covers the
1707
+ entire line with one ellipse, then bounds along that ellipse with
1708
+ multiple circles.
1709
1709
"""
1710
- # Note that this, in its current formalism, makes no check that bounding
1711
- # data at all corresponds to the differentials given. The onus is then
1710
+ # Note that this, in its current formalism, makes no check that bounding
1711
+ # data at all corresponds to the differentials given. The onus is then
1712
1712
# on the design of other functions which use it.
1713
-
1714
- # CCzg is required to be known as we need to know the ring which the minpolys lie in.
1713
+
1714
+ # CCzg is required to be known as we need to know the ring which the minpolys lie in.
1715
1715
CCzg , bounding_data_list = bounding_data
1716
-
1716
+
1717
1717
i0 , _ = upstairs_edge [0 ]
1718
1718
i1 , _ = upstairs_edge [1 ]
1719
1719
z0 = self ._vertices [i0 ]
1720
1720
z1 = self ._vertices [i1 ]
1721
1721
zwt , z1_minus_z0 = self .make_zw_interpolator (upstairs_edge )
1722
-
1722
+
1723
1723
# list of (centre, radius) pairs that still need to be processed
1724
1724
ball_stack = [(self ._RR (1 / 2 ), self ._RR (1 / 2 ), 0 )]
1725
- alpha = self ._RR (912 / 1000 )
1726
- # alpha set manually for scaling purposes. Basic benchmarking shows
1727
- # that ~0.9 is a sensible value.
1725
+ alpha = self ._RR (912 / 1000 )
1726
+ # alpha set manually for scaling purposes. Basic benchmarking shows
1727
+ # that ~0.9 is a sensible value.
1728
1728
E_global = self ._RR (2 )** (- self ._prec + 3 )
1729
1729
1730
- # Output will iteratively store the output of the integral.
1730
+ # Output will iteratively store the output of the integral.
1731
1731
V = VectorSpace (self ._CC , len (differentials ))
1732
1732
output = V (0 )
1733
1733
1734
- # The purpose of this loop is as follows: We know we will be using
1734
+ # The purpose of this loop is as follows: We know we will be using
1735
1735
# Gauss-Legendre quadrature to do the integral, and results from [Neu2018]_
1736
- # tell us an upper bound on the number of nodes required to achieve a
1737
- # given error bound for this quadrature, provided we have a bound for
1738
- # the integrand on a certain ellipse in the complex plane. The method
1736
+ # tell us an upper bound on the number of nodes required to achieve a
1737
+ # given error bound for this quadrature, provided we have a bound for
1738
+ # the integrand on a certain ellipse in the complex plane. The method
1739
1739
# developed by Bruin and Gao that uses Cauchy and Fujiwara can bound an
1740
1740
# algebraic integrand on a circular region. Hence we need a way to change
1741
- # from bounding with an ellipse to bounding with a circle. The size of
1742
- # these circles will be constrained by the distance to the nearest point
1743
- # where the integrand blows up, i.e. the nearest branchpoint. Basic
1744
- # benchmarking showed that it was in general a faster method to split
1745
- # the original line segment into multiple smaller line segments, and
1741
+ # from bounding with an ellipse to bounding with a circle. The size of
1742
+ # these circles will be constrained by the distance to the nearest point
1743
+ # where the integrand blows up, i.e. the nearest branchpoint. Basic
1744
+ # benchmarking showed that it was in general a faster method to split
1745
+ # the original line segment into multiple smaller line segments, and
1746
1746
# compute the contribution from each of the line segments bounding with
1747
1747
# a single circle, the benefits mainly coming when the curve is poorly
1748
- # conditioned s.t. the branch points are close together. The following
1749
- # loop does exactly this, repeatedly bisecting a segment if it is not
1748
+ # conditioned s.t. the branch points are close together. The following
1749
+ # loop does exactly this, repeatedly bisecting a segment if it is not
1750
1750
# possible to cover it entirely in a ball which encompasses an appropriate
1751
- # ellipse.
1751
+ # ellipse.
1752
1752
def local_N (ct , rt ):
1753
1753
cz = (1 - ct )* z0 + ct * z1 # This is the central z-value of our ball.
1754
- distances = [(cz - b ).abs () for b in self .branch_locus ]
1754
+ distances = [(cz - b ).abs () for b in self .branch_locus ]
1755
1755
rho_z = min (distances )
1756
1756
rho_t = rho_z / (z1_minus_z0 ).abs ()
1757
1757
rho_t = alpha * rho_t + (1 - alpha )* rt # sqrt(rho_t*rt) could also work
@@ -1765,11 +1765,11 @@ def local_N(ct, rt):
1765
1765
n = minpoly .degree (CCzg .gen (1 ))
1766
1766
ai_new = [(minpoly .coefficient ({CCzg .gen (1 ):i }))(z = cz + self ._CCz .gen (0 ))
1767
1767
for i in range (n )]
1768
- ai_pos = [self ._RRz ([c .abs () for c in h .list ()])
1768
+ ai_pos = [self ._RRz ([c .abs () for c in h .list ()])
1769
1769
for h in ai_new ]
1770
1770
m = [a (rho_z )/ z_1 for a in ai_pos ]
1771
1771
l = len (m )
1772
- M_tilde = 2 * max ((m [i ].abs ())** (1 / self ._RR (l - i ))
1772
+ M_tilde = 2 * max ((m [i ].abs ())** (1 / self ._RR (l - i ))
1773
1773
for i in range (l ))
1774
1774
cg = g (cz ,cw )
1775
1775
cdgdz = dgdz (cz ,cg )
@@ -1786,15 +1786,15 @@ def local_N(ct, rt):
1786
1786
1787
1787
if not lN :
1788
1788
cz = (1 - ct )* z0 + ct * z1
1789
- distances = [(cz - b ).abs () for b in self .branch_locus ]
1789
+ distances = [(cz - b ).abs () for b in self .branch_locus ]
1790
1790
rho_z = min (distances )
1791
1791
rho_t = rho_z / (z1_minus_z0 ).abs ()
1792
1792
1793
1793
if rho_t <= rt :
1794
1794
ball_stack .append ((ncts [0 ], nrt , 0 ))
1795
1795
ball_stack .append ((ncts [1 ], nrt , 0 ))
1796
1796
continue
1797
-
1797
+
1798
1798
lN = local_N (ct , rt )
1799
1799
1800
1800
nNs = [local_N (nct , nrt ) for nct in ncts ]
@@ -1834,7 +1834,7 @@ def matrix_of_integral_values(self, differentials, integration_method="heuristic
1834
1834
- ``differentials`` -- a list of polynomials.
1835
1835
1836
1836
- ``integration_method`` -- (default: ``'heuristic'``). String specifying
1837
- the integration method to use. The options are ``'heuristic'`` and
1837
+ the integration method to use. The options are ``'heuristic'`` and
1838
1838
``'rigorous'``.
1839
1839
1840
1840
OUTPUT:
@@ -1885,7 +1885,7 @@ def normalize_pairs(L):
1885
1885
raise ValueError ("Invalid integration method" )
1886
1886
1887
1887
integral_dict = {edge : line_int (edge ) for edge in occurring_edges }
1888
-
1888
+
1889
1889
rows = []
1890
1890
for cycle in cycles :
1891
1891
V = VectorSpace (self ._CC , len (differentials )).zero ()
@@ -1928,7 +1928,7 @@ def period_matrix(self):
1928
1928
3
1929
1929
1930
1930
One can check that the two methods give similar answers::
1931
-
1931
+
1932
1932
sage: from sage.schemes.riemann_surfaces.riemann_surface import RiemannSurface
1933
1933
sage: R.<x,y> = QQ[]
1934
1934
sage: f = y^2 - x^3 + 1
0 commit comments