Computing the Riemann Constant Vector
=====================================

We use the Python software package [`abelfunctions`](http://abelfunctions.cswiercz.info) to compute with Riemann surfaces obtained via the desingularization and compactification of complex plane algebraic curves. `abelfunctions` is written by Chris Swierczewski ([cswiercz@uw.edu](mailto:cswiercz@uw.edu)).

Throughout this notebook we use the genus 4 curve

$$
    f(x,y) = x^2y^3 - x^4 + 1
$$

as an example.

The primary object in `abelfunctions` is "`RiemannSurface`".


```
%matplotlib inline

from sympy.abc import x,y
from abelfunctions import (RiemannSurface, RiemannTheta, Jacobian,
                           RiemannConstantVector, puiseux, AbelMap)

f = x**2*y**3 - x**4 + 1
X = RiemannSurface(f,x,y)

print X
```

Demo - Homology Basis
=====================

A canonical basis of cycles
$$
    H_1(X,\mathbb{Z}) = \{a_1, \ldots, a_g, b_1, \ldots, b_g\}
$$
is computed using

    X.a_cycles()
    
and
    
    X.b_cycles().


```
a = X.a_cycles()
b = X.b_cycles()
```


```
# plot the x- and y-projections of the first a-cycle using
# 256 interpolating points
figx = a[0].plot_x(256, color='b', linestyle='dashed', linewidth=2)
figy = a[0].plot_y(256, color='g', linewidth=2)
```


```
# plot in 3d with the path parameter, t, on the z-axis
#
# note that there are no self-intersections in the path
elev = 50     # 50,-130
azim = -130

figx_3d = a[0].plot3d_x(256, color='blue', linewidth=2)
ax = figx_3d.gca()
ax.view_init(elev=elev, azim=azim)

figy_3d = a[0].plot3d_y(256, color='green', linewidth=2)
ax = figy_3d.gca()
ax.view_init(elev=elev, azim=azim)
```

Demo - Basis of Abelian Differentials of the First Kind
=======================================================

**Note**: `abelfunctions` does not necessarily compute a *normalized* basis of holomorphic one-forms.


```
omega = X.holomorphic_differentials()
for omegai in omega:
    print omegai
```

We can plot the value of the differentials along Riemann surface paths.
* blue: real part
* red: imaginary part

Here we plot 

$$
\omega_1 = \frac{1}{3x^2y^2} \; dx
$$

along the first $a$-cycle.


```
fig = omega[0].plot(a[0],256)
```

Demo - Riemann / Period Matrices
================================

Integrating each Abelian differential of the first kind basis element around each of the a- and b- cycles produces the *period matrix* of the Riemann surface

$$
A_{ij} = \oint_{a_j} \omega_i, \quad B_{ij} = \oint_{b_j} \omega_i.
$$

The *Riemann matrix* of the surface is then determined by

$$
\Omega = A^{-1}B \in \mathbb{C}^{g \times g}.
$$

The Riemann matrix is symmetric and its imaginary part is positive definite.


```
# for brevity, we only print the first four significant digits
import numpy
numpy.set_printoptions(precision=4, suppress=True)
```


```
Omega = X.riemann_matrix()
print Omega
```


```
symmetric_error = numpy.linalg.norm(Omega - Omega.T)
print symmetric_error
```


```
imag_part_evals = numpy.linalg.eigvals(Omega.imag)
print imag_part_evals
```

Demo - Places and Divisors
==========================

Places on a Riemann surface are given in terms of their Puiseux series. In this example, there is only one place on the surface "over" $x=0$ and three places "over" $x=2$. The places above two are "regular" in the sense that they lie above a non-discriminant point of the curve.

To declare places on a Riemann surface use the notation

```
places = RiemannSurface(xprojection)
```

where `xprojection` is the x-projection of the places you want to construct. That is, we compute all places above a given $x =$ `xprojection`.


```
places_above_zero = X(0)
for P in places_above_zero:
    print P
```


```
P.puiseux_series.extend(10)
```


```
print 'Places:'
places_above_two = X(2)
for P in places_above_two:
    print P
    
print '\nPuiseux series:'
series_at_two = puiseux(f,x,y,2)
for p in series_at_two:
    print p
```

**Divisors:** Once you have some places you can constuct divisors using the usual operations / notation.


```
P = places_above_zero[0]
Q = places_above_two[0]

D = 3*P + Q

print 'Divisor:', D
print 'Degree:', D.degree
```

Demo - The Abel Map
===================

The Abel map $A:X \to J(X)$ is a map from a Riemann surface to its Jacobian defined component-wise by

$$
A_j(P,Q) = \int_P^Q \omega_j.
$$

In `abelfunctions` this is computed using

    AbelMap(P,Q).

A single argument $P$ given to `AbelMap` will compute $A(P_0,P)$ where $P_0$ is a pre-determined "base place" of the Riemann surface. All paths on the Riemann surface (including the homology basis) begin at this base place.


```
J = Jacobian(X)
z1 = AbelMap(P)
z2 = AbelMap(Q)
z3 = AbelMap(P,Q)

print z1
print z2
print z3
```


```
print X.base_place()
```

Numerically verify that $A(P,Q) = A(P_0,Q) - A(P_0,P)$:


```
print numpy.linalg.norm( J((z2-z1) - z3) )
```

Demo - Localizing Differentials at Places
=========================================


```
nu = omega[0]
print nu
```


```
P
```


```
nu.localize(P)
```


```
P_oo = X('oo')[0]
P_oo.puiseux_series.extend(10)
print P_oo
```


```
nu.localize(P_oo)
```

Demo - Canonical Diviors
========================

Given an Abelian differential of the first kind you can compute its valuation divisor using the syntax

    Differential.valuation_divisor()


```
C0 = omega[0].valuation_divisor()

for place, multiplicity in C0:
    print '\nplace:       ', place
    print 'multiplicity:', multiplicity
    
print '\nDegree:', C0.degree
```


```
C1 = omega[1].valuation_divisor()

for place, multiplicity in C1:
    print '\nplace:       ', place
    print 'multiplicity:', multiplicity
    
print '\nDegree:', C1.degree
```


```
C2 = omega[2].valuation_divisor()

for place, multiplicity in C2:
    print '\nplace:       ', place
    print 'multiplicity:', multiplicity
    
print '\nDegree:', C0.degree
```


```
C3 = omega[3].valuation_divisor()

for place, multiplicity in C3:
    print '\nplace:       ', place
    print 'multiplicity:', multiplicity
    
print '\nDegree:', C3.degree
```

Demo - Riemann Constant Vector (RCV)
====================================

The function `RiemannConstantVector` is used to compute the Riemann constant vector on a Riemann surface. We make use of the fact

$$
    K(P) = K(P_0) + (g-1)A(P_0,P)
$$

to compute the RCV at various points. Most of the work is done in evaluating $K$ at the base place $P_0$ of the Riemann surface.


```
K = RiemannConstantVector # alias the function for brevity
```


```
P0 = X.base_place()
print K(P0)
```

**Theorem:** Given a canonical divisor $\mathcal{C}$

$$
2 K(P_0) + A(P_0, \mathcal{C}) \equiv 0.
$$

*Computational verification:*


```
C = omega[3].valuation_divisor()
print 'Divisor:'
print C

print '\nExpr:'
z = J(2*K(P0) + AbelMap(C))
print z
```

**Theorem:** Let $\theta : J(X) \times \mathfrak{h}_g$ denote the Riemann theta function. Then

$$
\theta(W,\Omega) = 0
$$

if and only if

$$
W = K(P_0) + A(P_0, \mathcal{D}).
$$

where $\mathcal{D}$ is a degree $g-1$ divisor.

*Computational verification:* we separate the "exponential" and "oscillatory" parts of the Riemann theta function (see DLMF)

$$
\theta(W,\Omega) = e^u v.
$$

Since $e^u$ never vanishes we only check for the vanishing of $v$ using the function `RiemannTheta.oscillatory_part()`.


```
W = K(P0)
v = RiemannTheta.oscillatory_part(W,Omega)
print abs(v)
```


```
D = sum(places_above_two)  # divisor of degree g-1 = 3
W = J(AbelMap(D) + K(P0))
v = RiemannTheta.oscillatory_part(W,Omega)
print abs(v)
```

**Verifying output:** we evaluate the oscillatory part using each half-lattice vector and plot the resulting magnitues in ascending order.


```
from abelfunctions.riemann_constant_vector import initialize_half_lattice_vectors

# list every half-lattice vector
h = initialize_half_lattice_vectors(X)

# evaluate the abel map at the valuation divisor C
AC = AbelMap(C)

# evaluate the oscillatory part of theta at
# each half-lattice vector minus the abel shift
W = numpy.array([J(hi.T - 0.5*AC).flatten() for hi in h])
v = RiemannTheta.oscillatory_part(W, Omega)
v = numpy.array(v, dtype=numpy.complex)

# compute the magnitues of the results and sort
magnitudes = sorted(abs(v))
print magnitudes[:3]
```


```
%matplotlib inline

import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.semilogy(magnitudes, marker='d', markersize=16, color='gray')
ax.axis([-32,256+32,1e-10,10])
fig.set_figwidth(16)
fig.set_figheight(8)
```


```

```

Comments / Suggestions for `abelfunctions`?
===========================================

Please send an email to Chris Swierczewski ([cswiercz@uw.edu](mailto:cswiercz@uw.edu)) or leave an "Issue" on the `abelfunctions` Github page.

> [http://github.com/cswiercz/abelfunctions](http://github.com/cswiercz/abelfunctions)


```

```

9a01865ea6785f958a3270dea6e7cfd600052f90
ipynb
Jupyter Notebook
notebooks/Riemann Constant Vector - Math Methods Talk.ipynb
RParini/abelfunctions
9263cd865e33cbd3ff4ffc1d59c14d77ff27f155
MIT
2016-03-22T14:08:44.000Z
2022-03-04T18:35:49.000Z
notebooks/Riemann Constant Vector - Math Methods Talk.ipynb
RParini/abelfunctions
9263cd865e33cbd3ff4ffc1d59c14d77ff27f155
MIT
2016-02-25T21:59:50.000Z
2021-06-26T11:05:29.000Z
notebooks/Riemann Constant Vector - Math Methods Talk.ipynb
RParini/abelfunctions
9263cd865e33cbd3ff4ffc1d59c14d77ff27f155
MIT
2016-03-22T13:27:41.000Z
2021-03-15T01:05:31.000Z
Qwen/Qwen-72B
1. YES
2. YES