# Sage Version 6.1.1, Release Date: 2014-02-04
# Numpy version 1.7.0
# Sympy version 0.7.4
# Abelfunctions version 0.1.0


### Packages needed.

import sympy
from sympy.abc import x,y
from abelfunctions import *
from numpy import dot, array

var('z1','z2','z3')

R.<z1,z2,z3> = PolynomialRing(CC, 3)

### Definition of even characteristics with binary code
### For example m17=[[0,0,1],[1,1,1]].
### Indeed 0*2^2+0*2^1+1*2^0=1
### and 1*2^2+1*2^1+1*2^0=7.

m11=array([[0,0,1],[0,0,1]])
m13=array([[0,0,1],[0,1,1]])
m15=array([[0,0,1],[1,0,1]])
m17=array([[0,0,1],[1,1,1]])
m22=array([[0,1,0],[0,1,0]])
m23=array([[0,1,0],[0,1,1]])
m27=array([[0,1,0],[1,1,1]])
m26=array([[0,1,0],[1,1,0]])
m31=array([[0,1,1],[0,0,1]])
m32=array([[0,1,1],[0,1,0]])
m35=array([[0,1,1],[1,0,1]])
m36=array([[0,1,1],[1,1,0]])
m44=array([[1,0,0],[1,0,0]])
m45=array([[1,0,0],[1,0,1]])
m46=array([[1,0,0],[1,1,0]])
m47=array([[1,0,0],[1,1,1]])
m51=array([[1,0,1],[0,0,1]])
m53=array([[1,0,1],[0,1,1]])
m54=array([[1,0,1],[1,0,0]])
m56=array([[1,0,1],[1,1,0]])
m62=array([[1,1,0],[0,1,0]])
m63=array([[1,1,0],[0,1,1]])
m64=array([[1,1,0],[1,0,0]])
m65=array([[1,1,0],[1,0,1]])
m71=array([[1,1,1],[0,0,1]])
m72=array([[1,1,1],[0,1,0]])
m74=array([[1,1,1],[1,0,0]])
m77=array([[1,1,1],[1,1,1]])


### Vector of the directions of the partial derivative.
v=array([[[1,0,0]],[[0,1,0]],[[0,0,1]]])


### Transformation factor from theta with null characteristic to theta with arbitrary characteristic m

def ee(m,t):
    return exp(pi*I*(dot(0.5*m[0],dot(t,0.5*m[0]))+2*dot(0.5*m[0],0.5*m[1])))

### Function which computes a Riemann Theta constant with characteristics m and Riemann matrix t.

def RT(m,t):
    return numerical_approx(ee(m,t)*RiemannTheta(dot(t,0.5*m[0])+0.5*m[1],t))

### Function which computes the partial derivative with respect z_j of the Riemann Theta functions with characteristic m, Riemann matrix t and evaluated at z=(0,0,0).

def djRT(m,t,j):
    return numerical_approx(ee(m,t)*(RiemannTheta(dot(t,0.5*m[0])+0.5*m[1],t,deriv=v[j])+pi*I*m[0,j]*RiemannTheta(dot(t,0.5*m[0])+0.5*m[1],t)))

### Function which computes the Jacobian determinant associated with the three odd characteristics m1, m2 and m3 and Riemann matrix om.

def jdet(m1,m2,m3,om):
    d1=djRT(m1,om,0)*djRT(m2,om,1)*djRT(m3,om,2)
    d2=djRT(m2,om,0)*djRT(m3,om,1)*djRT(m1,om,2)
    d3=djRT(m3,om,0)*djRT(m1,om,1)*djRT(m2,om,2)
    d4=djRT(m1,om,2)*djRT(m2,om,1)*djRT(m3,om,0)
    d5=djRT(m1,om,1)*djRT(m2,om,0)*djRT(m3,om,2)
    d6=djRT(m1,om,0)*djRT(m2,om,2)*djRT(m3,om,1)
    return d1 + d2 + d3 - d4 - d5 - d6


### Function which computes the bitangent associated with the odd characteristic bb and Riemann matrix om.

def bit(bb,om):
     return djRT(bb,om,0)*z1 + djRT(bb,om,1)*z2 + djRT(bb,om,2)*z3

### Definitions of the constants A, B and C, as in Section 6.

def A(om): return jdet(m77,m46,m51,om)*jdet(m31,m13,m26,om)*jdet(m77,m54,m26,om)/(jdet(m77,m31,m26,om)*jdet(m77,m23,m51,om)*jdet(m54,m13,m26,om))
def B(om): return jdet(m77,m46,m51,om)*jdet(m77,m62,m26,om)*jdet(m31,m13,m26,om)/(jdet(m77,m31,m26,om)*jdet(m62,m13,m26,om)*jdet(m77,m15,m51,om))
def C(om): return jdet(m77,m46,m51,om)*jdet(m31,m13,m26,om)*jdet(m77,m45,m26,om)/(jdet(m77,m31,m26,om)*jdet(m77,m32,m51,om)*jdet(m45,m13,m26,om))

### Function which computes the 8x8 bitangent matrix L associated with the Riemann matrix om.

def L(om):
    l=matrix(R,8)
    l[0,1]=jdet(m31,m13,m26,om)*bit(m77,om)
    l[0,2]=jdet(m22,m13,m35,om)*bit(m64,om)
    l[0,3]=jdet(m77,m64,m46,om)*bit(m51,om)
    l[0,4]=jdet(m77,m64,m51,om)*bit(m46,om)
    l[0,5]=jdet(m77,m64,m51,om)*bit(m23,om)
    l[0,6]=jdet(m77,m64,m51,om)*bit(m15,om)
    l[0,7]=jdet(m77,m64,m51,om)*bit(m32,om)    
    l[1,2]=(jdet(m22,m13,m35,om)*jdet(m77,m31,m26,om)/jdet(m77,m46,m51,om))*bit(m13,om)
    l[1,3]=jdet(m77,m13,m31,om)*bit(m26,om)
    l[1,4]=jdet(m77,m13,m26,om)*bit(m31,om)
    l[1,5]=(jdet(m31,m13,m26,om)*jdet(m77,m13,m26,om)/(A(om)*jdet(m54,m13,m26,om)))*bit(m54,om)
    l[1,6]=(jdet(m31,m13,m26,om)*jdet(m77,m13,m26,om)/(B(om)*jdet(m62,m13,m26,om)))*bit(m62,om)
    l[1,7]=(jdet(m31,m13,m26,om)*jdet(m77,m13,m26,om)/(C(om)*jdet(m45,m13,m26,om)))*bit(m45,om)
    l[2,3]=jdet(m64,m13,m22,om)*bit(m35,om)
    l[2,4]=jdet(m64,m13,m35,om)*bit(m22,om)
    l[2,5]=(jdet(m22,m13,m35,om)*jdet(m64,m13,m35,om)/(A(om)*jdet(m47,m13,m35,om)))*bit(m47,om)
    l[2,6]=(jdet(m22,m13,m35,om)*jdet(m64,m13,m35,om)/(B(om)*jdet(m71,m13,m35,om)))*bit(m71,om)
    l[2,7]=(jdet(m22,m13,m35,om)*jdet(m64,m13,m35,om)/(C(om)*jdet(m56,m13,m35,om)))*bit(m56,om)
    l[3,4]=(jdet(m77,m64,m46,om)*jdet(m51,m26,m35,om)/jdet(m17,m26,m35,om))*bit(m17,om)
    l[3,5]=(jdet(m77,m64,m46,om)*jdet(m51,m26,m35,om)/(A(om)*jdet(m72,m26,m35,om)))*bit(m72,om)
    l[3,6]=(jdet(m77,m64,m46,om)*jdet(m51,m26,m35,om)/(B(om)*jdet(m44,m26,m35,om)))*bit(m44,om)
    l[3,7]=(jdet(m77,m64,m46,om)*jdet(m51,m26,m35,om)/(C(om)*jdet(m63,m26,m35,om)))*bit(m63,om)
    l[4,5]=(A(om)*jdet(m77,m23,m51,om)-jdet(m77,m46,m51,om))*(jdet(m22,m31,m17,om)*jdet(m64,m13,m35,om)/(A(om)*jdet(m65,m31,m17,om)*jdet(m22,m13,m35,om)))*bit(m65,om)
    l[4,6]=(B(om)*jdet(m77,m15,m51,om)-jdet(m77,m46,m51,om))*(jdet(m22,m31,m17,om)*jdet(m64,m13,m35,om)/(B(om)*jdet(m53,m31,m17,om)*jdet(m22,m13,m35,om)))*bit(m53,om)
    l[4,7]=(1-C(om)*(jdet(m77,m64,m32,om)/jdet(m77,m64,m46,om)))*(jdet(m77,m64,m51,om)*jdet(m46,m31,m22,om)/(C(om)*jdet(m74,m31,m22,om)))*bit(m74,om)
    l[5,6]=(1-(jdet(m54,m13,m26,om)*jdet(m77,m23,m51,om)*jdet(m15,m64,m51,om)*jdet(m77,m62,m26,om)/jdet(m23,m64,m51,om)*jdet(m77,m54,m26,om)*jdet(m77,m15,m51,om)*jdet(m62,m13,m26,om)))*(jdet(m77,m64,m51,om)*jdet(m23,m47,m72,om)/(B(om)*jdet(m36,m47,m72,om)))*bit(m36,om)
    l[5,7]=(1/C(om)-(jdet(m77,m64,m32,om)/(A(om)*jdet(m77,m64,m23,om))))*(jdet(m23,m54,m47,om)*jdet(m77,m64,m51,om)/(jdet(m11,m54,m47,om)))*bit(m11,om)
    l[6,7]=(1/C(om)-(jdet(m77,m64,m32,om)/(B(om)*jdet(m77,m64,m15,om))))*(jdet(m15,m62,m71,om)*jdet(m77,m64,m51,om)/(jdet(m27,m62,m71,om)))*bit(m27,om)
    return l+l.transpose()


### Function which computes the 4x4 matrix Q associated with the Riemann matrix om.

def Q(om):
    q=matrix(R,4)
    q[0,1]=(jdet(m31,m13,m26,om)/jdet(m77,m31,m26,om))*bit(m77,om)
    q[0,2]=(jdet(m22,m13,m35,om)/jdet(m77,m31,m26,om))*bit(m64,om)
    q[0,3]=(jdet(m77,m64,m46,om)/jdet(m77,m31,m26,om))*bit(m51,om)
    q[1,2]=(jdet(m22,m13,m35,om)/jdet(m77,m46,m51,om))*bit(m13,om)
    q[1,3]=(jdet(m77,m13,m31,om)/jdet(m77,m31,m26,om))*bit(m26,om)
    q[2,3]=(jdet(m64,m13,m22,om)/jdet(m77,m31,m26,om))*bit(m35,om)
    return q+q.transpose()

### A Riemann sruface of genus 3 of which we want to compute the bitangent matrix.
f=y**4+x**4+1 


### The Riemann surface X associated with f
X = RiemannSurface(f,x,y)

print 'The Riemann surface X is:\n', f,'= 0'

### Compute the genus of X.
print 'The genus of the curve is:\n', X.genus()

### Compute the Riemann matrix of X.
Omega = X.riemann_matrix()

print 'The Riemann matrix is:\n', Omega


### The bitangents matrix L evaluated at the period matrix Omega
### It takes few seconds

print 'Computing the bitangents matrix...'
ll=L(Omega)
print 'The bitangents matrix L evaluated at the period matrix Omega is:\n', ll

### The matrix Q evaluated at the period matrix Omega
### It takes few seconds

qq=Q(Omega)
print 'The matrix Q evaluated at the period matrix Omega is:\n', qq


### det Q=0 is a quartic in P^2 projectively equivalent to the canonical model of the original Riemann surface f=y**4+x**4+1

print 'A quartic in P2 projectively equivalent to the canonical model of the original Riemann surface f=y**4+x**4+1 is:\n', qq.det(),'= 0'
#
