This iPython notebook contains the example code presented in the paper Computing the Riemann Constant Vector by Bernard Deconinck, Matthew S. Patterson, and Chris Swierczewski.
First, we construct the Riemann surface, $X$, defined by the plane algebraic curve
$$ f(x,y) = x^2 y^3 - x^4 + 1. $$We computationally determine its genus as well as a basis $\{\tilde{\omega}_1, \ldots, \tilde{\omega}_g\}$ of Abelian differentials of the first kind, $\Omega_X^1$, defined on the surface $X$.
In [1]:
from sympy.abc import x,y
from abelfunctions import (RiemannSurface, RiemannTheta, Jacobian,
AbelMap, RiemannConstantVector, puiseux)
f = x**2*y**3 - x**4 + 1
X = RiemannSurface(f,x,y)
g = X.genus()
omega = X.holomorphic_differentials()
print 'differentials:'
for omegai in omega:
print omegai
print 'genus:', g
Next, we compute a canonical basis of cycles $\{a_1, \ldots, a_g, b_1, \ldots, b_g\}$ of the first homology group $H_1(X, \mathbb{Z})$. That is, every cycle on $X$ can be written as a linear combination of these cycles.
In [2]:
%matplotlib inline
a = X.a_cycles()
b = X.b_cycles()
# use 256 points to plot the x- and y-parts of the path
figx = a[0].plot_x(256, color='blue', linewidth=2, linestyle='dashed')
figy = a[0].plot_y(256, color='green', linewidth=2)
Using these data we can compute the period matrix of the surface:
$$ \tau = [A \; B] \in \mathbb{C}^{g \times 2g} $$where
$$ A_{ij} = \oint_{a_j} \tilde{\omega}_i, \quad \text{and} \quad B_{ij} = \oint_{a_j} \tilde{\omega}_i $$If using a normalized basis of differentials, $\{\omega_1, \ldots, \omega_g\}$, the period matrix is of the form
$$ \tau = [I \; \Omega] \in \mathbb{C}^{g \times 2g}. $$This is equivalent to setting $\Omega = A^{-1}B$ and, in fact, since Abelfunctions returns a non-normalized basis this is how the matrix $\Omega$ is computed. We computationally verify that $\Omega$ is a "Riemann matrix": a symmetric matrix with positive definite imaginary part.
In [3]:
# for brevity, we only print the first four significant digits
import numpy
numpy.set_printoptions(precision=4, suppress=True)
tau = X.period_matrix()
A = tau[:g,:g]
B = tau[:g,g:]
Omega = X.riemann_matrix() # returns A**(-1)*B
print A
print
print B
print
print Omega
In [4]:
symmetric_error = numpy.linalg.norm(Omega - Omega.T)
imag_part_evals = numpy.linalg.eigvals(Omega.imag)
print 'error:', symmetric_error
print 'evals:', imag_part_evals
Places $P \in X$ are represented using a Puiseux series. If $x$ and $y$ are the affine coordinates of the underlying curve $f$ then we write
$$ P = \begin{cases} x(t) = \alpha + \lambda t^e \\ y(t) = \sum_i \beta_i t^{n_i} \end{cases} $$In the case when $e = 1$ (an "unramified place") and $x \neq \infty$ a place is synonymous with a tuple $(\alpha, \beta) \in \mathbb{C}^2$ such that $f(\alpha, \beta) = 0$.
A divisor $\mathcal{D}$ is a formal sum of places, $\mathcal{D} = \sum_i n_i P_i$, where each of the $n_i$ are the multiplicities of the places $P_i$ in the divisor. Below we compute local representation of several places on the Riemann surface and construct some divisors from them.
In [5]:
places_above_zero = X(0)
print places_above_zero
In [6]:
print 'Places:'
places_above_two = X(2)
for P in places_above_two:
print P
print 'Puiseux:'
series_at_two = puiseux(f,x,y,2)
for p in series_at_two:
print p
In [7]:
P = places_above_zero[0]
Q = places_above_two[0]
D = 3*P + Q
print 'Divisor:', D
print 'Degree:', D.degree
In [8]:
D2 = sum(places_above_two)
print D2
Given a fixed base place $P_0 \in X$, the Abel Map $A : X \to J(X)$ is defined
$$ A(P) = \left( \int_{P_0}^P \omega_1, \ldots, \int_{P_0}^P \omega_g \right) $$where the path of integration from $P_0$ to $P$ is the same for each Abelian differential of the first kind. The Abel map is linear over divisors. That is, if $\mathcal{D} = \sum_i n_iP_i$ then
$$ A(\mathcal{D}) = \sum_i n_i A(P_i). $$When we want to change or make the back point explicit we write $A(P_0,\mathcal{D})$. Below we evaluate the Abel map at the places and divisors constructed above.
In [9]:
J = Jacobian(X) # reduces vectors modulo lattice ZZ^g + Omega ZZ^g
z1 = AbelMap(P) # Abel map from P0 to P
z2 = AbelMap(Q) # Abel map from P0 to Q
z3 = AbelMap(P,Q) # Abel map from P to Q
print z1
print z2
print z3
# numerically verify that A(P,Q) = A(P0,Q) - A(P0,P)
print numpy.linalg.norm( J((z2-z1) - z3) )
In [10]:
w = AbelMap(D)
print w
# verify that w = 3*z1 + z2 mod period lattice
z = J(3*z1 + z2)
print numpy.linalg.norm(w-z)
Let $P = (x(t), y(t))$ be the Puiseux series representation of a place $P \in X$. Given a meromorphic differential $\nu$ we can write the "localization of $\nu$ at $P$" as a Laurent series in the local parameter $t$:
$$ \nu\big|_P = (ct^n + \cdots) dt. $$Then the "valuation" of $\nu$ at $P$ is
$$ \text{val}(\nu,P) = n. $$The valuation divisor $(\nu) \in \text{Div}(X)$ of $\nu$ is the divisor
$$ (\nu) = \sum_{P \in X} \text{val}(\nu,P) P. $$A divisor is canonical if and only if it's the valuation divisor of a holomorphic differential on $X$. All canonical divisors are of degree $2g-2$. Below, we compute a canonical divisor for each holmomrphic differential basis element $\tilde{\omega}_i$ and verify that their degree is $6$.
In [11]:
C0 = omega[0].valuation_divisor()
for place,multiplicity in C0:
print multiplicity, place
print 'Degree:', C0.degree
In [12]:
C1 = omega[1].valuation_divisor()
for place,multiplicity in C1:
print multiplicity, place
print 'Degree:', C1.degree
In [13]:
C2 = omega[2].valuation_divisor()
for place,multiplicity in C2:
print multiplicity, place
print 'Degree:', C2.degree
In [14]:
C3 = omega[3].valuation_divisor()
for place,multiplicity in C3:
print multiplicity, place
print 'Degree:', C3.degree
The Riemann constant vector satisfies the following two theorems:
Theorem 1: $\mathcal{C}$ is a canonical divisor if any only if $K(P_0) \equiv -A(P_0,\mathcal{C})$.
Theorem 2: $\theta(W,\Omega) = 0$ if and only if $W = A(P_0,\mathcal{D}) + K(P_0)$ where $\mathcal{D}$ is a degree $g-1$ divisor.
We compute $K$ below and verify that these two theorems are satisfied.
In [15]:
K = RiemannConstantVector # alias the RCV function for brevity
P0 = X.base_place()
print K(P0)
In [16]:
z = J(2*K(P0) + AbelMap(C3))
print z
In [17]:
W = K(P0)
v = RiemannTheta.oscillatory_part(W,Omega)
print abs(v)
In [18]:
D = sum(places_above_two)
W = J(AbelMap(D) + K(P0))
v = RiemannTheta.oscillatory_part(W, Omega)
print abs(v)
In [19]:
import matplotlib.pyplot as plt
from abelfunctions.riemann_constant_vector import initialize_half_lattice_vectors
In [20]:
h = initialize_half_lattice_vectors(X)
V = 0.5*AbelMap(C3)
oscillatory_parts = [
RiemannTheta.oscillatory_part(J(hi.T-V),Omega)
for hi in h
]
oscillatory_magnitudes = sorted(map(abs,oscillatory_parts))
In [21]:
%matplotlib inline
fig = plt.figure()
ax = fig.add_subplot(111)
ax.semilogy(range(256), oscillatory_magnitudes,
linestyle='', marker='d',markerfacecolor='grey',markersize=8)
ax.axis([-32,256+32,1e-9,10])
ax.set_xticks([-32] + [32*k for k in range(10)])
ax.set_xticklabels([''] + [str(32*k) for k in range(9)] + [''])
ax.grid(True)
fig.set_figwidth(9)
fig.set_figheight(4)
fig.show()
In [ ]: