""" This package provides a method of representing Jacobi forms and computing Jacobi theta series. (See page 1 of [EZ1985] for a definition of Jacobi theta series.) EXAMPLES: We can create the Jacobi-theta series for the E8 lattice with the vector (0,1,0,0,0,0,0,0) using the following code. sage: ME8 = [[2, 0, 0, 0, 0, 0, 0, 1],\ [0, 2, 1, 1, 1, 1, 1, 1],\ [0, 1, 2, 1, 1, 1, 1, 1],\ [0, 1, 1, 2, 1, 1, 1, 1],\ [0, 1, 1, 1, 2, 1, 1, 1],\ [0, 1, 1, 1, 1, 2, 1, 1],\ [0, 1, 1, 1, 1, 1, 2, 0],\ [1, 1, 1, 1, 1, 1, 0, 2]] sage: QE8 = QuadraticForm(matrix(QQ,8,8,ME8)) sage: J1=Jacobi_Form(QE8,[0,1,0,0,0,0,0,0],4) sage: J1 Jacobi Form of weight 4 and index 1 with Fourier expansion: 1 + (Z^-2 + 56*Z^-1 + 126 + 56*Z + Z^2)*q + (126*Z^-2 + 576*Z^-1 + 756 + 576*Z + 126*Z^2)*q^2 + (56*Z^-3 + 756*Z^-2 + 1512*Z^-1 + 2072 + 1512*Z + 756*Z^2 + 56*Z^3)*q^3 + (Z^-4 + 576*Z^-3 + 2072*Z^-2 + 4032*Z^-1 + 4158 + 4032*Z + 2072*Z^2 + 576*Z^3 + Z^4)*q^4 + O(q^5) We can create the Jacobi-theta series for the D16^+ lattice with the vector (0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0) using the following code. sage: MD16 = [[4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],\ [0, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],\ [0, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],\ [0, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],\ [0, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],\ [0, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],\ [0, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1],\ [0, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1],\ [0, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1],\ [0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1],\ [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1],\ [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1],\ [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1],\ [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1],\ [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 0],\ [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 2]] sage: QD16 = QuadraticForm(matrix(QQ,16,16,MD16)) sage: J2=Jacobi_Form(QD16,[0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0],4) sage: J2 Jacobi Form of weight 8 and index 1 with Fourier expansion: 1 + (Z^-2 + 56*Z^-1 + 366 + 56*Z + Z^2)*q + (366*Z^-2 + 14016*Z^-1 + 33156 + 14016*Z + 366*Z^2)*q^2 + (56*Z^-3 + 33156*Z^-2 + 260712*Z^-1 + 462392 + 260712*Z + 33156*Z^2 + 56*Z^3)*q^3 + (Z^-4 + 14016*Z^-3 + 462392*Z^-2 + 1987392*Z^-1 + 2998638 + 1987392*Z + 462392*Z^2 + 14016*Z^3 + Z^4)*q^4 + O(q^5) ALGORITHMS: Jacobi forms are created using the function compute_theta_series. This algorithm is described in Section 5.2 of [deQ2010]. It uses the method of Lagrange multiplies to find the vectors X in the lattice for which Q(X)<=the precision which we want for the Fourier expansion. ..WARNING:: If you do not input the correct types of inputs then the code will not throw an exception it will just not work. ..TODO:: Throw exceptions. Change matrix in compute_theta_series to an integral matrix. Consider other functions. Add documentation for all the functions. Copywrite statement. REFERENCES: .. [deQ2010] Victoria de Quehen. Jacobi Forms. Master's thesis, Department of Mathematics, McGill University, 2010. .. [EZ1985] Martin Eichler and Don Zagier. The Theory of Jacobi Forms, Progress in Mathematics, 55, Boston, MA: Birkhäuser Boston, 1985. AUTHORS: - Victoria de Quehen - Andrew Fiori - David Roe """ cdef extern from "math.h": double sqrt(double) double floor(double) long floorl(double) double ceil(double) long ceill(double) double round(double) include "../../ext/stdsage.pxi" #include "../ext/cdefs.pxi" #include "../ext/interrupt.pxi" from warnings import warn from copy import deepcopy from sage.matrix.matrix import Matrix from sage.matrix.constructor import Matrix from sage.modules.free_module_element import vector from sage.rings.integer_ring import IntegerRing, ZZ from sage.rings.rational_field import RationalField, QQ from sage.rings.ring import Ring from sage.rings.power_series_ring import PowerSeriesRing from sage.rings.laurent_series_ring import LaurentSeriesRing from sage.quadratic_forms.quadratic_form import QuadraticForm from sage.rings.integer import Integer from sage.rings.rational import Rational cdef compute_theta_series_precomp_nomaloc(self,MTX, MTXI, V, index, LRing, PRing, double* MTXfinalcolumn, double* Prec, double * MatricesItip, double** AdjustCentre): """ Secretly an internal function of jacobi forms for the purpose of using precomputed values. I would prefer to have a forward declaration and then put this code below the one inside to ease commenting the setup of this function, and how it interacts with compute_all_theta is very hackish. Functions the same as Jacobi_Form::compute_theta_series except that it takes as parameters values that are purely precomputed """ #### Initialize Various Values #### cdef int ZPrec = Integer(4*index*self.__prec).isqrt() + 1 Q = self.__Q cdef int dim = Q.dim() cdef int dim_minus_one = dim - 1 cdef double error=.01 cdef double precision=(self.__prec)+.01 cdef double orig_precision = precision #### Precomputations #### _BMTX = 2*MTX*V #### Convert Python Variables to C Variables #### cdef int i, k cdef int* BMTX = sage_malloc(dim*sizeof(int)) for i from 0 <= i < dim: BMTX[i] = _BMTX[i] cdef double* Boundary = sage_malloc(dim*sizeof(double)) for i from 0 <= i < dim: Boundary[i] = 0 Boundary[0]=sqrt(precision*(MTXI[0,0])) cdef int* X = sage_malloc(dim*sizeof(int)) for i from 0 <= i < dim: X[i] = 0 X[0] = ceill(-Boundary[0]-error) cdef int** Res = sage_malloc((self.__prec+1)*sizeof(int*)) for i from 0 <= i <= self.__prec: Res[i] = sage_malloc((ZPrec*2)*sizeof(int)) for k from 0 <= k < ZPrec*2: Res[i][k] = 0 cdef double* Cen = sage_malloc(dim*sizeof(double)) for i from 0 <= i < dim: Cen[i] = 0 #### Initialize #### cdef int lastspot=0 cdef int Bres = BMTX[0]*X[0] # Bres = , calculated iteratively cdef int Qres # Qres = Q(x), calculated iteratively cdef int B_X_res = (2*MTXfinalcolumn[0]*X[0]+error) cdef int top cdef double diff #### Main Iteration #### while 1: top = floorl(Boundary[lastspot] + Cen[lastspot]+error) while lastspot < dim - 1 and X[lastspot]<=top: lastspot=lastspot+1 ## Compute new boundary while we advance in hyperplanes precision=precision+Prec[lastspot]*(X[lastspot-1]-Cen[lastspot-1])**2 if -error(2*X[lastspot]*MTXfinalcolumn[lastspot]+error) top=floorl( Boundary[lastspot]+ Cen[lastspot]+error) ## Check the next vector if X[lastspot]<=top: ## In order to iteratively compute quadratic form we ## precompute its value at first integral point diff = Boundary[lastspot]-Cen[lastspot]+X[lastspot] Qres= (orig_precision + B_X_res*diff - diff**2 * MTXfinalcolumn[dim-1] + .5) # .5 so that it rounds to the nearest integer. We shouldn't need nearly this much. ## Check the first one before iterating Res[Qres][(Bres+ZPrec)] = Res[Qres][(Bres+ZPrec)]+ 1 while X[lastspot] < top: ## Tight iteration over final dimension X[lastspot]=X[lastspot]+1 Bres += BMTX[lastspot] ## Q(x+[0,..,1]) = Q(x)+B(x,[0,...,0,1])+Q([0,..,0,1]) Qres = (Qres + B_X_res + MTXfinalcolumn[dim-1]) ## B(X+1,1) = B(X,1) + 2Q(1) B_X_res += (2*MTXfinalcolumn[lastspot]+error) Res[Qres][(Bres+ZPrec)] = Res[Qres][(Bres+ZPrec)]+ 1 while X[lastspot] >=floorl(Boundary[lastspot]+Cen[lastspot]+error): ## Step back through to higher dimensional hyperplanes we have ## completed precision = precision - Prec[lastspot]*(X[lastspot-1]-Cen[lastspot-1])**2 Bres -= BMTX[lastspot]*X[lastspot] B_X_res -= (2*MTXfinalcolumn[lastspot]*X[lastspot]+error) lastspot=lastspot-1 if lastspot < 0: ## We are done now! Res_list = [] for i from 0 <= i <= self.__prec: Res_list.append([]) for k from 0 <= k < 2*ZPrec: Res_list[i].append(Res[i][k]) tmp = [ LRing(Res_list[i],-ZPrec) for i in range(self.__prec+1) ] returnresult = PRing(tmp, self.__prec+1) #### Clean up C variables #### sage_free(BMTX) sage_free(Boundary) sage_free(X) for i from 0 <= i <= self.__prec: sage_free(Res[i]) sage_free(Res) sage_free(Cen) return Jacobi_Form(Q, V, self.__prec,self.__weight, index, self.__level, returnresult) #### Advance the next iteration #### X[lastspot]=X[lastspot]+1 Bres += BMTX[lastspot] B_X_res += (2*MTXfinalcolumn[lastspot]+error) class Jacobi_Form(): """ TODO - add examples for all the functions. - add a description of what Jacobi forms are. """ def __reduce_basis(self, Q, V): """ Changes basis to try to minimize radius TODO - enable once basis of short vectors works (ie fix basis of short vectors) """ # M = Matrix(Q.basis_of_short_vectors()); # self.__Q = QuadraticForm(M*Q.matrix()*M.transpose()) # self.__V = M.transpose().inverse()*V self.__Q = Q self.__V = V return def __compute_radius(self): """ Computes the radius of a circle that encloses the ellipse our Quadratic Form defines """ from sage.functions.all import sqrt MinEigen = min( self.__Q.matrix().eigenvalues() ) radius = floorl(sqrt(self.__prec/MinEigen)) + 1 return radius def __compute_rectangle(self): """ Computes a bounding rectangle for our elipse """ from sage.functions.all import sqrt MTX=(self.__Q.matrix())/2 MTXI=MTX.inverse() Boundary=[floorl(sqrt(self.__prec*abs(MTXI[i,i]))) for i in range(self.__Q.dim())] def compute_theta_series(self): """ Compute the the Fourier coeffiecients of this Jacobi theta series. (See page 1 of [EZ1985] for a definition of Jacobi theta series.) INPUT: This internal function should is called on an object where .__Q is a positive-definite, unimodular quadratic form on a lattice and V is a vector in the lattice. OUTPUT: An object with the Jacobi form Fourier expansion of the Jacobi theta series associated to the quadratic form .__Q and vector .__V. ALGORITHM: This algorithm is described in Section 5.2 of [deQ2010]. It involves iterating over all the vectors X in the lattice for which Q(X)<= __prec (the precision which we want for the Fourier expansion). The vectors in the lattice which satisfy this property all lie within an n-dimensional ellipsoid, namely, <=Q(X), where n is the dimension of the lattice. Notice if we intersect the ellipsoid defined by this equation with a hyperplane we obtain an (n-1)-dimensional ellipsoid. If intersect the resulting ellipsoid with another hyperplane we obtain an (n-2)-ellipsoid. Repeating this process we eventually obtain a 0-dimensional point. If we iterate over all integral hyperplanes in one direction we get the (n-1)-ellipsoids that contain exactly those points (vectors X) such that Q(X)<=__prec. Following the recursive proceedure described in the last paragraph we can iterate over all of the integral points (vectors X) such that Q(X)<=__prec. In order to perform this proceedure efficiently we will perform the following steps: #. Precomputations which allow us to effectively compute the changing position of the centre and boundaries of the ellipsoids as we interate through hyperplanes. The boundaries are computed using the method of Lagrange multipliers. #. Next we convert the python variables to C variables. #. As we iterate through the hyperplane, we iteratively compute the the values that change with each iteration (for example, the center and boundaries of the ellipse, and the values of the quadratic and bilinear forms restricted to the hyperplane). The computation is iterative in that it adjusts the value from the previous iteration rather the completely recomputing it. #. With each new hyperplane and resulting ellipsoid we recurse to the lower dimensions. We also iteratively compute the changing values. #. Finally we convert the number of vectors with each value of #Q(X) and B(X,v) to coefficients for the Fourier expansion and #we return the object with the Fourier expansion. .. WARNING:: Does not throw exceptions if given the wrong type of data. .. TODO - use shortest basis and mention this in input - are the lattice condition in input the right ones? - general code cleanup - change tabs to space - change variable names (consistent capitalization) - Improve variable descriptions - Most Arithmetic could be made integral, this would make things a little faster (LCM of the det of the upper principle minors for some things, 2 for MTX stuff) - Properly document orig_prec, Prec, _BMTX and BMTX and Prec onward. - Change Centre to AdjustCen REFERENCES: .. [deQ2010] Victoria de Quehen. Jacobi Forms. Master's thesis, Department of Mathematics, McGill University, 2010. .. [EZ1985] Martin Eichler and Don Zagier. The Theory of Jacobi Forms, Progress in Mathematics, 55, Boston, MA: Birkhäuser Boston, 1985. AUTHORS: - Victoria de Quehen - Andrew Fiori - David Roe """ #### Variables that are precomputed # Q = self.__Q = The quadratic form whose Jacobi theta series we are # computing. # self.__V = The vector whose Jacobi theta series we are # computing. # dim = The dimension of the quadratic space. # self.__prec = The highest exponent for q whose coefficient # we will compute. # error = The assumed bound on rounding errors from taking square # roots. This is used to round back to values known to be integral. # orig_precision = self.__prec + epsilon to ensure even with # rounding errors we still find all points <= self.__prec. # ZPrec = The highest exponent on zeta whose Fourier # coefficient we will compute. (This is a result of the bound # self.__prec on q). # MXT = 1/2A, (where the quadratic form is Q(x)=(1/2)(X^T)AX). # MTXI = The inverse matrix of MTX. # BMTX = The vector 2 *MXT*(self.__V). This is used to compute # the bilinear form B(X,self.__V). # MTXfinalcolumn = The last column of MTX. This is used to # iteratively compute how the quadratic form changes as we # iterate. # Matrices[m] = Lower right (dim-m)x(dim-m)-submatrix of MTX. # MatriciesI[m] = The inverse matrix of Matrices[m]. # MatricesItip[m] = The upper left entry of MTXI[m]. This is # used to compute the boundary of iteration for mth-coordinate # of X. # AdjustCentre[m] = A vector calculated at the beginning that is # used to calculate how the center of the (dim-m)-dimensional # ellipsoid moves as we iterate. # Prec[m] = A value used to compute how the bound for the # (dim-m) dimensional ellipse changes as we iterate. #### TODO - change variable name of Prec to AdjustRadius #### ##### Variables changed during iterations # X = The current lattice point on which we are iterating. # precision = The distance to the boundary for the current hyperplane. #### TODO - rename precision to radius #### # Boundary[m] = The computed bounds for the current iteration # in the mth-dimension Cen = The centre of the current # ellipse. Res[n][r] = The number of vectors we have already # found with Q(X)=n and B(X,__V)=r. lastspot = The dimension # we are currently iterating in. Bres = The current value of # B(X,__V) (computed iteratively). Qres = The current value # of Q(X) (computed iteratively). B_X_res = A variable used # to iteratively compute Qres. top = The upper bound on # iteration for X in the current dimension. diff = Variable # to compute distance between boundary and first lattice point # on during an iteration. precision ??? #### Initialize Various Values #### cdef int ZPrec = Integer(4*self.__index*self.__prec).isqrt() + 1 Q = self.__Q cdef int dim = Q.dim() cdef int dim_minus_one = dim - 1 cdef double error=.0001 cdef double precision=(self.__prec)+error cdef double orig_precision = precision #### Precomputations #### MTX=(Q.matrix()) Matrices = [Q.matrix().copy()/2] MatricesI=[Matrices[0].inverse()] _AdjustCentre=[Matrix(QQ, [[0] for i in range(dim)])] _Prec=[0] for m in range(1, dim): Matrices[m-1].subdivide(1,1) M2=Matrices[m-1].subdivision(1,0) Matrices.append(Matrices[m-1].subdivision(1,1).copy()) MatricesI.append(Matrices[m].inverse()) _AdjustCentre.append(-MatricesI[m]*M2) _Prec.append(-(Matrices[m-1])[0,0]+(M2.transpose()*MatricesI[m]*M2)[0,0]) _BMTX = MTX*self.__V LRing = LaurentSeriesRing(ZZ, 'Z') PRing = PowerSeriesRing(LRing, 'q') #### Convert Python Variables to C Variables #### #### First Convert Matrix Variables### cdef int i, k cdef int* BMTX = sage_malloc(dim*sizeof(int)) for i from 0 <= i < dim: BMTX[i] = _BMTX[i] cdef int* MTX2finalcolumn = sage_malloc(dim*sizeof(int)) for i from 0 <= i < dim: MTX2finalcolumn[i] = (MTX[i, dim-1]) cdef int MTXEND = (MTX[dim-1,dim-1]/2) cdef double* MatricesItip = sage_malloc(dim*sizeof(double)) for i from 0 <= i < dim: MatricesItip[i] = MatricesI[i][0,0] #### Next Convert Variables Related to the Ellipsoid #### cdef double* Prec = sage_malloc(dim*sizeof(double)) for i from 0 <= i < dim: Prec[i] = _Prec[i] cdef double* Boundary = sage_malloc(dim*sizeof(double)) for i from 0 <= i < dim: Boundary[i] = 0 Boundary[0]=sqrt(precision*(MatricesItip[0])) cdef double** AdjustCentre = sage_malloc(dim*sizeof(double*)) for k from 0 <= k < dim: AdjustCentre[k] = sage_malloc(dim*sizeof(double)) for k from 0 <= k < dim: for i from 1 <= i <= k: AdjustCentre[i][k] = _AdjustCentre[i][k-i,0] cdef double* Cen = sage_malloc(dim*sizeof(double)) for i from 0 <= i < dim: Cen[i] = 0 cdef int* X = sage_malloc(dim*sizeof(int)) for i from 0 <= i < dim: X[i] = 0 X[0] = ceill(-Boundary[0]-error) cdef int** Res = sage_malloc((self.__prec+1)*sizeof(int*)) for i from 0 <= i <= self.__prec: Res[i] = sage_malloc((ZPrec*2)*sizeof(int)) for k from 0 <= k < ZPrec*2: Res[i][k] = 0 #### Initialize #### cdef int lastspot=0 cdef int Bres = BMTX[0]*X[0] # Bres = , calculated iteratively. cdef int Qres # Qres = Q(x), calculated iteratively. cdef int B_X_res = MTX2finalcolumn[0]*X[0] cdef int top cdef double diff #### Main Iteration #### while 1: top = floorl(Boundary[lastspot] + Cen[lastspot]+error) while lastspot < dim - 1 and X[lastspot]<=top: lastspot=lastspot+1 ## Compute the new boundary while we advance in the hyperplanes precision=precision+Prec[lastspot]*((X[lastspot-1]-Cen[lastspot-1])**2) if precision<0: # Correct any rounding errors Boundary[lastspot]=0 else: Boundary[lastspot]=sqrt(precision*MatricesItip[lastspot]) Cen[lastspot]=0 ## Compute the new centre for i from 1 <= i <= lastspot: Cen[lastspot]=Cen[lastspot]+AdjustCentre[i][lastspot]*(X[i-1]-Cen[i-1]) ## Initialize the next iteration X[lastspot] = ceill(-Boundary[lastspot]+Cen[lastspot]-error) Bres += BMTX[lastspot]*X[lastspot] B_X_res += X[lastspot]*MTX2finalcolumn[lastspot] top=floorl( Boundary[lastspot]+ Cen[lastspot] + error) ## Check the next vector if X[lastspot]<=top: ## In order to iteratively compute quadratic form we ## precompute its value at first integral point. diff = Boundary[lastspot]-Cen[lastspot]+X[lastspot] Qres= floorl(orig_precision + B_X_res*diff - (diff**2) * MTXEND + error) ## Check the first one before iterating Res[Qres][Bres+ZPrec] = Res[Qres][Bres+ZPrec]+ 1 while X[lastspot] < top: ## Tight iteration over final dimension X[lastspot]=X[lastspot]+1 Bres += BMTX[lastspot] ## Q(x+[0,..,1]) = Q(x)+B(x,[0,...,0,1])+Q([0,..,0,1]) Qres = (Qres + B_X_res + MTXEND) ## B(X+1,1) = B(X,1) + 2Q(1) B_X_res += (MTX2finalcolumn[lastspot]) Res[Qres][Bres+ZPrec] = Res[Qres][Bres+ZPrec]+ 1 while X[lastspot] >=floorl(Boundary[lastspot]+Cen[lastspot]+error): ## Step back through to higher dimensional hyperplanes as we ## complete th iteration of the lower dimensional ellipsoids. precision = precision - Prec[lastspot]*(X[lastspot-1]-Cen[lastspot-1])**2 Bres -= BMTX[lastspot]*X[lastspot] B_X_res -=( MTX2finalcolumn[lastspot]*X[lastspot]) lastspot=lastspot-1 if lastspot < 0: ## We we just need the information to a Fourier expansion. Res_list = [] for i from 0 <= i <= self.__prec: Res_list.append([]) for k from 0 <= k < 2*ZPrec: Res_list[i].append(Res[i][k]) tmp = [ LRing(Res_list[i],-ZPrec) for i in range(self.__prec+1) ] self.__Fourier_expansion = PRing(tmp, self.__prec+1) #### Clean up C variables #### sage_free(BMTX) sage_free(MTX2finalcolumn) sage_free(Prec) sage_free(MatricesItip) for k from 0 <= k < dim: sage_free(AdjustCentre[k]) sage_free(AdjustCentre) sage_free(Boundary) sage_free(X) for i from 0 <= i <= self.__prec: sage_free(Res[i]) sage_free(Res) sage_free(Cen) return #### Advance the next iteration #### X[lastspot]=X[lastspot]+1 Bres += BMTX[lastspot] B_X_res += MTX2finalcolumn[lastspot] def __init__(self, form=None, vec=None, Prec=10, weight=None, index=None, level=0, Fourier = None): """ Input: arg1 is the quadratic form, arg2 is the vector. ToDo-input, output, example """ if Fourier == None: #TODO check isinstance(form, QuadraticForm): #TODO - check ring of V/Q, check prec, fix level. V = vector(ZZ,vec) self.__prec = Prec self.__reduce_basis(form,V) self.__index = self.__Q(self.__V) self.__weight = self.__Q.dim()/2 self.__level = level self.compute_theta_series() return """ Creates a Jacobi form object from the data without needing a basis. Note that this will not remember a quadratic form nor vector as it may not be a theta series. Input: arg1 is the Fourier series, arg2 is the level. TODO - check the index, weight and level are the same - have minimum precision - check types """ if not (form == None): self.__Q = form if not (vec == None): self.__V = vec self.__prec =Prec self.__index=index self.__weight=weight self.__level=level self.__Fourier_expansion=Fourier return def __repr__(self): """ Outputs a representation TODO - level """ return "Jacobi Form of weight %i and index %i with Fourier expansion:\n %s"%(self.__weight,self.__index,self.__Fourier_expansion) def _latex_(self): """ Outputs latex rep TODO make pretty """ return "Jacobi Form of weight %i and index %i with Fourier expansion:\n %s"%(self.__weight,self.__index,self.__Fourier_expansion) def __getitem__(self, ij): """ Returns the c(i,j) coefficient TODO - Return 0 if 'j' is large/small - Return error if i > prec - do we want J[3] to return the coef of q^3??? """ i,j = ij i = int(i) j = int(j) # gets the coefficient of q^i L=self.__Fourier_expansion.padded_list(i+1)[i] # returns the coefficient of Z^j return L.list()[L.degree()+j] def __setitem__(self, ij, coeff): """ TODO - should maybe dissallow this... """ raise NotImplementedError, "Not allowed to set coefs of Jacobi forms" return def __eq__(self, right): """ check equality with right TODO - NYI... be sure that this accounds for precision. weight and index are implied by the above working (for non-constant) level... not sure how to care. """ if self.__weight != right.__weight: return False if self.__index != right.__index: return False return self.__Fourier_expansion == right.__Fourier_expansion def __add__(self, right): """ returns sum of jacobi forms ToDo check the congruence subgroup and level and create a new Jacobi form ToDo We might be able to do more things in the event of theta series if they are in the same space and vectors in diff irreducible lattices, this is a theta series """ if self.__index != right.__index : raise NotImplementedError, "Oops, the index should be the same." if self.__weight != right.__weight : raise NotImplementedError, "Oops, the weight should be the same." if self.__level != right.__level : raise NotImplementedError, "Currently we do not know how to add modular forms of different levels." return Jacobi_Form(Fourier = self.__Fourier_expansion+right.__Fourier_expansion,Prec = min(self.__prec,right.__prec) ,weight = self.__weight,index = self.__index) def __mul__(self, right): """ Multiply Not Yet implemented The product of two jacobi theta series is a jacobi theta series ToDo we know the Quadratic form and vector to use, could track this """ if self.__level != right.__level : raise NotImplementedError, "Currently we do not know how to multiply modular forms of different levels." return Jacobi_Form(Fourier = self.__Fourier_expansion*right.__Fourier_expansion,Prec = min(self.__prec,right.__prec) ,weight = self.__weight+right.__weight,index = self.__index+right.__index, level=self.__level) def __call__(self, n1, n2): """ Evaluate ... but on what? """ return self.__Fourier_expansion(n1)(n2) def weight(self): """ Returns the weight of the Jacobi form. """ return self.__weight def index(self): """ Returns the index of the Jacobi form. """ return self.__index def level(self): """ Returns the level of the Jacobi form. Not yet implemented. """ raise NotImplementedError, "NYI" return self.__level def Fourier_expansion(self): """ Returns the Fourier expansion of the Jacobi form. """ return self.__Fourier_expansion def compute_theta_series_precomp(self,MTX, MTXI, MatricesI,_AdjustCentre,_Prec,V,index): """ This function computes the Fourier coeffiecients of this Jacobi theta series. In order to do this we must iterate over all the lattice points such that Q(x) <= __prec. General Description of Algorithm: We iterate over all points of the lattice such that Q(x) <= __prec by observing that if we intersect the ellipse defined by this equation with a hyperplane we get a lower dimensional ellipse. Thus, if we iterate over all integral hyperplanes in one direction and follow a recursive proceedure we can iterate over all of the integral points. In order to perform this proceedure efficiently we will perform the following steps: 1. Precomputations which allow us to effectively compute the changing position of the centre and boundaries of the ellipses. 2. As we iterate through the hyperplane, we iteratively compute the the values that change with each iteration (for example the the center of the ellipse, the values of the quadratic and bilinear forms). The computation is iterative in that it adjusts the value from the previous iteration rather the completely recomputing it. 3. As we recurse to the lower dimensions, we also iteratively compute the changing values. TODO - general code cleanup - reduce boundary - change tabs to space - change variable names (consistend capitalization) - Improve variable descriptions """ #### Variables that are precomputed # Q = The Quadratic Form. # self.__V = The vector whose Jacobi theta series we are computing. # dim = The dimension of the quadratic space. # self.__prec = The highest exponent for q whose coefficient we will compute. # orig_precision = self.__prec + epsilon to ensure even with rounding errors we still find all points <= self.__prec. # ZPrec = The highest exponent on zeta whose Fourier coefficient we will compute (this is a result of the bound on q). # MXT = 1/2A (where the quadratic form is Q(x)=(1/2)(X^T)AX). # MTXI = The inverse matrix of MTX. # BMTX = The vector A*(self.__V) (used to compute the bilinear form B(X,__V)). # MTXfinalcolumn = The last column of MTX, used to iteratively compute how the quadratic form changes as we iterate. # Matrices[m] = Lower right (dim-m)x(dim-m) submatrix of MTX. # MatriciesI[m] = The inverse of Matrices[m] # MatricesItip[m] = The upper left entry of MTXI[m] which is used to compute the boundary of iteration for mth-coordinate of X. # AdjustCentre[m] = A vector calculated at the beginning which is used to calculate how the center of the (dim-m) dimensional ellipse moves as we iterate. #### TODO - change to AdjustCentre #### # Prec[m] = A value used to compute how the bound for the (dim-m) dimensional ellipse changes as we iterate. #### TODO - change variable name of Prec to AdjustRadius #### # error = Assumed bound on rounding errors from taking square roots, used to round back to values known to be integral. ##### Variables changed during iterations # X = The current lattice point on which we are iterating. # precision = The distance to the boundary for the current hyperplane. #### TODO - rename precision to radius #### # Boundary[m] = The computed bounds for the current iteration in the mth-dimension # Cen = The centre of the current ellipse. # Res[n][r] = The number of vectors we have already found with Q(X)=n and B(X,__V)=r. # lastspot = The dimension we are currently iterating in. # Bres = The current value of B(X,__V) (computed iteratively). # Qres = The current value of Q(X) (computed iteratively). # B_X_res = A variable used to iteratively compute Qres. # top = The upper bound on iteration for X in the current dimension. # diff = Variable to compute distance between boundary and first lattice point on during an iteration. #### Initialize Various Values #### cdef int ZPrec = Integer(4*index*self.__prec).isqrt() + 1 Q = self.__Q cdef int dim = Q.dim() cdef int dim_minus_one = dim - 1 cdef double error=.01 cdef double precision=(self.__prec)+.01 cdef double orig_precision = precision #### Precomputations #### _BMTX = 2*MTX*V LRing = LaurentSeriesRing(ZZ, 'Z') PRing = PowerSeriesRing(LRing, 'q') #### Convert Python Variables to C Variables #### cdef int i, k cdef int* BMTX = sage_malloc(dim*sizeof(int)) for i from 0 <= i < dim: BMTX[i] = _BMTX[i] cdef double* MTXfinalcolumn = sage_malloc(dim*sizeof(double)) for i from 0 <= i < dim: MTXfinalcolumn[i] = MTX[i, dim-1] cdef double* Prec = sage_malloc(dim*sizeof(double)) for i from 0 <= i < dim: Prec[i] = _Prec[i] cdef double* MatricesItip = sage_malloc(dim*sizeof(double)) for i from 0 <= i < dim: MatricesItip[i] = MatricesI[i][0,0] cdef double** AdjustCentre = sage_malloc(dim*sizeof(double*)) for k from 0 <= k < dim: AdjustCentre[k] = sage_malloc(dim*sizeof(double)) for k from 0 <= k < dim: for i from 1 <= i <= k: AdjustCentre[i][k] = _AdjustCentre[i][k-i,0] cdef double* Boundary = sage_malloc(dim*sizeof(double)) for i from 0 <= i < dim: Boundary[i] = 0 Boundary[0]=sqrt(precision*(MTXI[0,0])) cdef int* X = sage_malloc(dim*sizeof(int)) for i from 0 <= i < dim: X[i] = 0 X[0] = ceill(-Boundary[0]-error) cdef int** Res = sage_malloc((self.__prec+1)*sizeof(int*)) for i from 0 <= i <= self.__prec: Res[i] = sage_malloc((ZPrec*2)*sizeof(int)) for k from 0 <= k < ZPrec*2: Res[i][k] = 0 cdef double* Cen = sage_malloc(dim*sizeof(double)) for i from 0 <= i < dim: Cen[i] = 0 #### Initialize #### cdef int lastspot=0 cdef int Bres = BMTX[0]*X[0] # Bres = , calculated iteratively cdef int Qres # Qres = Q(x), calculated iteratively cdef int B_X_res = (2*MTXfinalcolumn[0]*X[0]+error) cdef int top cdef double diff #### Main Iteration #### while 1: top = floorl(Boundary[lastspot] + Cen[lastspot]+error) while lastspot < dim - 1 and X[lastspot]<=top: lastspot=lastspot+1 ## Compute new boundary while we advance in hyperplanes precision=precision+Prec[lastspot]*(X[lastspot-1]-Cen[lastspot-1])**2 if -error(2*X[lastspot]*MTXfinalcolumn[lastspot]+error) top=floorl( Boundary[lastspot]+ Cen[lastspot]+error) ## Check the next vector if X[lastspot]<=top: ## In order to iteratively compute quadratic form we ## precompute its value at first integral point diff = Boundary[lastspot]-Cen[lastspot]+X[lastspot] Qres= (orig_precision + B_X_res*diff - diff**2 * MTXfinalcolumn[dim-1] + .5) # .5 so that it rounds to the nearest integer. We shouldn't need nearly this much. ## Check the first one before iterating Res[Qres][(Bres+ZPrec)] = Res[Qres][(Bres+ZPrec)]+ 1 while X[lastspot] < top: ## Tight iteration over final dimension X[lastspot]=X[lastspot]+1 Bres += BMTX[lastspot] ## Q(x+[0,..,1]) = Q(x)+B(x,[0,...,0,1])+Q([0,..,0,1]) Qres = (Qres + B_X_res + MTXfinalcolumn[dim-1]) ## B(X+1,1) = B(X,1) + 2Q(1) B_X_res += (2*MTXfinalcolumn[lastspot]+error) Res[Qres][(Bres+ZPrec)] = Res[Qres][(Bres+ZPrec)]+ 1 while X[lastspot] >=floorl(Boundary[lastspot]+Cen[lastspot]+error): ## Step back through to higher dimensional hyperplanes we have ## completed precision = precision - Prec[lastspot]*(X[lastspot-1]-Cen[lastspot-1])**2 Bres -= BMTX[lastspot]*X[lastspot] B_X_res -= (2*MTXfinalcolumn[lastspot]*X[lastspot]+error) lastspot=lastspot-1 if lastspot < 0: ## We are done now! Res_list = [] for i from 0 <= i <= self.__prec: Res_list.append([]) for k from 0 <= k < 2*ZPrec: Res_list[i].append(Res[i][k]) tmp = [ LRing(Res_list[i],-ZPrec) for i in range(self.__prec+1) ] returnresult = PRing(tmp, self.__prec+1) #### Clean up C variables #### sage_free(BMTX) sage_free(MTXfinalcolumn) sage_free(Prec) sage_free(MatricesItip) for k from 0 <= k < dim: sage_free(AdjustCentre[k]) sage_free(AdjustCentre) sage_free(Boundary) sage_free(X) for i from 0 <= i <= self.__prec: sage_free(Res[i]) sage_free(Res) sage_free(Cen) return Jacobi_Form(Q, V, self.__prec,self.__weight, index, self.__level, returnresult) #### Advance the next iteration #### X[lastspot]=X[lastspot]+1 Bres += BMTX[lastspot] B_X_res += (2*MTXfinalcolumn[lastspot]+error) def compute_all_theta_series(self,coefs, extra=0): """ Code being messed with """ #### Variables that are precomputed # Q = The Quadratic Form. # self.__V = The vector whose Jacobi theta series we are computing. # dim = The dimension of the quadratic space. # self.__prec = The highest exponent for q whose coefficient we will compute. # orig_precision = self.__prec + epsilon to ensure even with rounding errors we still find all points <= self.__prec. # ZPrec = The highest exponent on zeta whose Fourier coefficient we will compute (this is a result of the bound on q). # MXT = 1/2A (where the quadratic form is Q(x)=(1/2)(X^T)AX). # MTXI = The inverse matrix of MTX. # BMTX = The vector A*(self.__V) (used to compute the bilinear form B(X,__V)). # MTXfinalcolumn = The last column of MTX, used to iteratively compute how the quadratic form changes as we iterate. # Matrices[m] = Lower right (dim-m)x(dim-m) submatrix of MTX. # MatriciesI[m] = The inverse of Matrices[m] # MatricesItip[m] = The upper left entry of MTXI[m] which is used to compute the boundary of iteration for mth-coordinate of X. # AdjustCentre[m] = A vector calculated at the beginning which is used to calculate how the center of the (dim-m) dimensional ellipse moves as we iterate. #### TODO - change to AdjustCentre #### # Prec[m] = A value used to compute how the bound for the (dim-m) dimensional ellipse changes as we iterate. #### TODO - change variable name of Prec to AdjustRadius #### # error = Assumed bound on rounding errors from taking square roots, used to round back to values known to be integral. ##### Variables changed during iterations # X = The current lattice point on which we are iterating. # precision = The distance to the boundary for the current hyperplane. #### TODO - rename precision to radius #### # Boundary[m] = The computed bounds for the current iteration in the mth-dimension # Cen = The centre of the current ellipse. # Res[n][r] = The number of vectors we have already found with Q(X)=n and B(X,__V)=r. # lastspot = The dimension we are currently iterating in. # Bres = The current value of B(X,__V) (computed iteratively). # Qres = The current value of Q(X) (computed iteratively). # B_X_res = A variable used to iteratively compute Qres. # top = The upper bound on iteration for X in the current dimension. # diff = Variable to compute distance between boundary and first lattice point on during an iteration. tmpprec = self.__prec #### Initialize Various Values #### cdef int ZPrec = Integer(4*self.__index*self.__prec).isqrt() + 1 Q = self.__Q cdef int dim = Q.dim() cdef int dim_minus_one = dim - 1 cdef double error=.01 cdef double precision=(self.__prec)+.01 cdef double orig_precision = precision #### Precomputations #### MTX=(Q.matrix())/2 MTXI=MTX.inverse() Matrices = [Q.matrix().copy()/2] MatricesI=[Matrices[0].inverse()] _AdjustCentre=[Matrix(QQ, [[0] for i in range(dim)])] _Prec=[0] for m in range(1, dim): Matrices[m-1].subdivide(1,1) M2=Matrices[m-1].subdivision(1,0) Matrices.append(Matrices[m-1].subdivision(1,1).copy()) MatricesI.append(Matrices[m].inverse()) _AdjustCentre.append(-MatricesI[m]*M2) _Prec.append(-(Matrices[m-1])[0,0]+(M2.transpose()*MatricesI[m]*M2)[0,0]) _BMTX = 2*MTX*self.__V LRing = LaurentSeriesRing(ZZ, 'Z') PRing = PowerSeriesRing(LRing, 'q') #### Convert Python Variables to C Variables #### cdef int i, k cdef int* BMTX = sage_malloc(dim*sizeof(int)) for i from 0 <= i < dim: BMTX[i] = _BMTX[i] cdef double* MTXfinalcolumn = sage_malloc(dim*sizeof(double)) for i from 0 <= i < dim: MTXfinalcolumn[i] = MTX[i, dim-1] cdef double* Prec = sage_malloc(dim*sizeof(double)) for i from 0 <= i < dim: Prec[i] = _Prec[i] cdef double* MatricesItip = sage_malloc(dim*sizeof(double)) for i from 0 <= i < dim: MatricesItip[i] = MatricesI[i][0,0] cdef double** AdjustCentre = sage_malloc(dim*sizeof(double*)) for k from 0 <= k < dim: AdjustCentre[k] = sage_malloc(dim*sizeof(double)) for k from 0 <= k < dim: for i from 1 <= i <= k: AdjustCentre[i][k] = _AdjustCentre[i][k-i,0] cdef double* Boundary = sage_malloc(dim*sizeof(double)) for i from 0 <= i < dim: Boundary[i] = 0 Boundary[0]=sqrt(precision*(MTXI[0,0])) cdef int* X = sage_malloc(dim*sizeof(int)) for i from 0 <= i < dim: X[i] = 0 X[0] = ceill(-Boundary[0]-error) cdef double* Cen = sage_malloc(dim*sizeof(double)) for i from 0 <= i < dim: Cen[i] = 0 #### Initialize #### cdef int lastspot=0 cdef int Bres = BMTX[0]*X[0] # Bres = , calculated iteratively cdef int Qres # Qres = Q(x), calculated iteratively cdef int B_X_res = (2*MTXfinalcolumn[0]*X[0]+error) cdef int top cdef double diff LIST = [] VLIST = [] COUNT=0 self.__prec = coefs #### Main Iteration #### while 1: top = floorl(Boundary[lastspot] + Cen[lastspot]+error) while lastspot < dim - 1 and X[lastspot]<=top: lastspot=lastspot+1 ## Compute new boundary while we advance in hyperplanes precision=precision+Prec[lastspot]*(X[lastspot-1]-Cen[lastspot-1])**2 if -error(2*X[lastspot]*MTXfinalcolumn[lastspot]+error) top=floorl( Boundary[lastspot]+ Cen[lastspot]+error) ## Check the next vector if X[lastspot]<=top: ## In order to iteratively compute quadratic form we ## precompute its value at first integral point diff = Boundary[lastspot]-Cen[lastspot]+X[lastspot] Qres= (orig_precision + B_X_res*diff - diff**2 * MTXfinalcolumn[dim-1] + .5) # .5 so that it rounds to the nearest integer. We shouldn't need nearly this much. ## Check the first one before iterating XX = vector([X[i] for i in range(dim)]) J3 = compute_theta_series_precomp_nomaloc(self, MTX, MTXI, XX, Qres, LRing, PRing, MTXfinalcolumn, Prec, MatricesItip, AdjustCentre) #J3 = self.compute_theta_series_precomp(MTX,MTXI,MatricesI,_AdjustCentre,_Prec,XX,Qres) #J3 = Jacobi_Form(Q,XX,coefs) COUNT = COUNT+1 if (COUNT % 100) == 0: print COUNT found = False for x in LIST: if J3 == x: found = True if not found: LIST.append(J3) VLIST.append(XX) while X[lastspot] < top: ## Tight iteration over final dimension X[lastspot]=X[lastspot]+1 Bres += BMTX[lastspot] ## Q(x+[0,..,1]) = Q(x)+B(x,[0,...,0,1])+Q([0,..,0,1]) Qres = (Qres + B_X_res + MTXfinalcolumn[dim-1]) ## B(X+1,1) = B(X,1) + 2Q(1) B_X_res += (2*MTXfinalcolumn[lastspot]+error) XX = vector([X[i] for i in range(dim)]) J3 = compute_theta_series_precomp_nomaloc(self, MTX, MTXI, XX, Qres, LRing, PRing, MTXfinalcolumn, Prec, MatricesItip, AdjustCentre) #J3 = self.compute_theta_series_precomp(MTX,MTXI,MatricesI,_AdjustCentre,_Prec,XX,Qres) #J3 = Jacobi_Form(Q,XX,coefs) COUNT = COUNT+1 if (COUNT % 100) == 0: print COUNT found = False for x in LIST: if J3 == x: found = True if not found: LIST.append(J3) VLIST.append(XX) while X[lastspot] >=floorl(Boundary[lastspot]+Cen[lastspot]+error): ## Step back through to higher dimensional hyperplanes we have ## completed precision = precision - Prec[lastspot]*(X[lastspot-1]-Cen[lastspot-1])**2 Bres -= BMTX[lastspot]*X[lastspot] B_X_res -= (2*MTXfinalcolumn[lastspot]*X[lastspot]+error) lastspot=lastspot-1 if lastspot < 0: ## We are done now! if extra > 0: LIST=[] self.__prec = self.__prec + extra for XX in VLIST: J3 = compute_theta_series_precomp_nomaloc(self, MTX, MTXI, XX, Q(XX), LRing, PRing, MTXfinalcolumn, Prec, MatricesItip, AdjustCentre) LIST.append(J3) #### Clean up C variables #### sage_free(BMTX) sage_free(MTXfinalcolumn) sage_free(Prec) sage_free(MatricesItip) for k from 0 <= k < dim: sage_free(AdjustCentre[k]) sage_free(AdjustCentre) sage_free(Boundary) sage_free(X) sage_free(Cen) self.__prec = tmpprec return LIST #### Advance the next iteration #### X[lastspot]=X[lastspot]+1 Bres += BMTX[lastspot] B_X_res += (2*MTXfinalcolumn[lastspot]+error)