#= This is a Julia code to recover the algebraic curve from its Riemann matrix via the Dubrovin threefold of the curve, based on a paper: Numerical reconstruction of curves from theirJacobians by Daniele Agostini, Türkü Özlüm Çelik and Demir Eken. =# using LinearAlgebra using Theta using DynamicPolynomials import Base.Cartesian using Base.Cartesian using HomotopyContinuation using Pkg using MultivariatePolynomials using Random: shuffle function D2Theta(doubletau::RiemannMatrix,caratt::Array) #= ---------- doubletau : gxg Riemann Matrix A Riemann matrix of a genus g Riemann surface. caratt : array of characteristic Returns ------- gxg Matrix Hessian Matrix of the theta function with the specific characteristics i.e. Q[eps] in Lemma 2.1 =# g = doubletau.g; z = zeros(Int64,g) H = zeros(Complex{Float64},g,g); E = Array[zeros(Int64,g)] for i=1:g-1 append!(E,Array[zeros(Int64,g)]) end for i=1:g E[i][i]+=1 end for i=1:g for j=i:g H[i,j] = theta(z,doubletau,derivs=[E[i],E[j]],char=caratt); H[j,i] = H[i,j]; end end return H; end function D4Theta(doubletau::RiemannMatrix,caratt::Array) #= ---------- doubletau : gxg Riemann Matrix A Riemann matrix of a genus g Riemann surface. caratt : array of characteristic Returns ------- gxgxgxg tensor Matrix of theta constants of the fourth derivatives evaluated at the specific characteristics (all the theta constants of fourth order needed for Lemma 2.1) =# g = doubletau.g; z = zeros(Int64,g) F = zeros(Complex{Float64},g,g,g,g); E = Array[zeros(Int64,g)] for i=1:g-1 append!(E,Array[zeros(Int64,g)]) end for i=1:g E[i][i]+=1 end for i=1:g for j=i:g for k=j:g for l=k:g F[i,j,k,l]=theta(z,doubletau,derivs=[E[i],E[j],E[k],E[l]],char=caratt); F[i,j,l,k]=F[i,j,k,l]; F[i,k,j,l]=F[i,j,k,l]; F[i,l,j,k]=F[i,j,k,l]; F[i,k,l,j]=F[i,j,k,l]; F[i,l,k,j]=F[i,j,k,l]; F[j,i,k,l]=F[i,j,k,l]; F[j,i,l,k]=F[i,j,k,l]; F[j,k,i,l]=F[i,j,k,l]; F[j,l,i,k]=F[i,j,k,l] F[j,k,l,i]=F[i,j,k,l]; F[j,l,k,i]=F[i,j,k,l]; F[k,j,i,l]=F[i,j,k,l]; F[k,j,l,i]=F[i,j,k,l]; F[k,i,j,l]=F[i,j,k,l]; F[k,l,j,i]=F[i,j,k,l]; F[k,i,l,j]=F[i,j,k,l]; F[k,l,i,j]=F[i,j,k,l]; F[l,j,k,i]=F[i,j,k,l]; F[l,j,i,k]=F[i,j,k,l]; F[l,k,j,i]=F[i,j,k,l]; F[l,i,j,k]=F[i,j,k,l]; F[l,k,i,j]=F[i,j,k,l]; F[l,i,k,j]=F[i,j,k,l]; end end end end return F; end function U4D(doubletau::RiemannMatrix,caratt::Array,u::Array) #= Parameters ---------- doubletau : gxg Riemann Matrix A Riemann matrix of a genus g Riemann surface. caratt : epsilon value, array consisting of characteristic u: Array of variables. For instance for a g 4 curve this input will be u[1:4]. Returns ------- Polynomial whose coefficients are theta constants with specific characteristics, =# g = doubletau.g; F = D4Theta(doubletau,caratt) z = zeros(Int64,g) pol = 0 for i=1:g for j=1:g for k=1:g for l=1:g pol+=u[i]*u[j]*u[k]*u[l]*F[i,j,k,l] end end end end return pol; end function generateBit(size::Int64,combination::Array) #=Supportive function to generate bits of n size that will be used in lambdaValues function. Parameters ---------- size : Integer Genus of curve combation: Empty array A array that will be used for recursion purposes Returns ------- Array String array consisting of 0,1 values =# if size== 0 return combination else if length(combination) ==0 generateBit(size - 1, ["0", "1"]) else recursionArray=[] for i in combination addZero=string(i , "0") append!(recursionArray,[addZero]) end for i in combination addOne=string(i , "1") append!(recursionArray,[addOne]) end return generateBit(size-1,recursionArray) end end end function getEpsilon(n::Int64,value::Int64) #=Supportive function that will be used in lambdaValues function. Parameters ---------- n : Integer Genus of curve value: Integer Value that corresponds to the specific epsilon.(For instance for a genus 4 curve there will be 2^4 epsilon values. This value corresponds to the each of those epsilon values) Returns ------- Array Array consisting of 0,1 values that corresponds to the epsilon value. =# bitArray=generateBit(n,[]) epsVal=[parse(Int, x) for x in bitArray[value]] return [epsVal, zeros(Int64,n)] end function lambdaValues(doubletau::RiemannMatrix) #=Supportive function for findQuartics. Compute lambda_epsilon values that will make the dubrovin quartic lambda*F(epsilon) independent of the two unknowns c and d Parameters ---------- doubletau : Riemann Matrix Riemann Matrix that corresponds to curve Returns ------- 2^g x (2^g- g*(g+1)/2 -1) Array where g is the genus of the curve 2^g part corresponds to the lambda_epsilon values of each quartic and (2^g- g*(g+1)/2 -1) corresponds to different quartics. =# g = doubletau.g; z = zeros(Int64,g); indp= Int((g*(g+1)/2)+1); bMatrix=zeros(Complex{Float64},indp,2^g); for k=1:2^g dummyA=[] hessian=D2Theta(doubletau,getEpsilon(g,k)) for i=1:g for j=i:g complexvalue=hessian[i,j] append!(dummyA,[complexvalue]) end end for i=1:indp-1 bMatrix[i,k]+=dummyA[i] end thetaat= theta(z,doubletau,char=getEpsilon(g,k)) bMatrix[indp,k] = thetaat end return nullspace(bMatrix) end function findQuartics(doubletau::RiemannMatrix) #=Compute the quartics which defines the curve of the Riemann matrix as a set, i.e. the quartics in Lemma 2.1 of the paper Parameters ---------- doubletau : Riemann Matrix Riemann Matrix that corresponds to curve Returns ------- Array that has lenght (2^g- g*(g+1)/2 -1) where g is the genus of the curve Array of quartics corresponding in u[1:g] coming from the given Riemann Matrix =# g = doubletau.g; z = zeros(Int64,g) lambdas = lambdaValues(doubletau) FD4Matrix=zeros(Polynomial{true,Complex{Float64}},2^g,1) @polyvar u[1:g] finalQuartic=[] for k=1:2^g derValue=U4D(doubletau,getEpsilon(g,k),u) FD4Matrix[k]=derValue end for i=1:size(lambdas,2) finalQuartic = append!(finalQuartic,[dot(conj(lambdas[:,i]),FD4Matrix)]) end return finalQuartic end #Example 3.2 of the paper #We define the Riemann Matrix using the period matrix #to be able to go back to the desired coordinates for a verification. (cf. Equation (18) in paper.) PMatrix=[[ -1.3096323621833203800780650010e-31 + 0.023488474383265381159603412007im 0.031549149334681191609436229192 + 9.8607613152626475676466070660e-32im 0.031549149334681191609436229192 + 1.1170393677445967947724672067e-31im -0.031549149334681191609436229192 + 0.023488474383265381159603412007im -2.0800043399382147213004561780e-32 + 2.3207455829866192029324534208e-31im 9.7837241174871582620800145415e-32 + 0.023488474383265381159603412007im]; [ 4.0059342843254505743564341206e-32 - 0.013840159414824083576257113403im 6.9333477997940490710015205933e-32 0.024982524765593877234425949619 + 5.3926038442842603885567382392e-32im 0.024982524765593877234425949619 - 0.013840159414824083576257113403im 3.7988968153038227201529164917e-32 - 2.7348205210298749113394886785e-32im -4.7763062620803446536377880492e-32 + 0.013840159414824083576257113403im]; [ -7.8963127719876669975295095646e-33 + 0.013840159414824083576257113403im -0.024982524765593877234425949619 - 0.027680318829648167152514226806im 2.6963019221421301942783691196e-33 + 1.5869662741750823429181258247e-31im -3.4281553010092798184396407378e-32 + 0.013840159414824083576257113403im 1.5600032549536610409753421335e-32 - 0.027680318829648167152514226806im 2.7489812021255453595631052306e-32 + 0.013840159414824083576257113403im]]; g=size(PMatrix)[1] piA=PMatrix[1:g,1:g] piB=PMatrix[1:g,g+1:2*g] RMMatrix=RiemannMatrix(2*inv(piA)*piB,siegel=false); #find the Dubrovin Quartics that is corresponding to RMMatrix quartics=findQuartics(RMMatrix); #### This function is to test the resulting quartics ## whether cutting out the canonical model. For this #we need to pass to the starting coordinate system. # (cf. Equation (18)). function changeOfBasis(periodMatrix::Array,quartic::Polynomial{true,Complex{Float64}}) #=Compute the Dubrovin quartics of any given Riemann Matrix Parameters ---------- periodMatrix : Array of length g*2g Period Matrix that corresponds to curve quartic: Polynomial Dubrovin Quartics that corresponds to the curve which can be found with findQuartics method Returns ------- Polynomial Dubrovin quartics that the basis has been changed to the normalized basis of differential =# varBasis=[] g=size(periodMatrix)[1] n=length(variables(quartic)) for i in variables(quartic) append!(varBasis,[i]) end varBasis=reshape(varBasis,n,1) newBasis = inv(periodMatrix[1:g,1:g])*varBasis varBasis=reshape(varBasis,n) newBasis=reshape(newBasis,n) varBasis=convert(Array{PolyVar{true},1}, varBasis) newBasis=convert(Array{Polynomial{true,Complex{Float64}},1}, newBasis) quartic =subs(quartic,varBasis=>newBasis); return quartic end #The test for Example 3.2. One estimates the equation of the Trott curve after the basis change. changeOfBasis(PMatrix,quartics[1])*81/coefficients(changeOfBasis(PMatrix,quartics[1]))[1]