Mercurial > btquotients
changeset 13:ca02f8d93048 tip
Some additions made in functionality.
| author | Marc Masdeu <masdeu@math.columbia.edu> |
|---|---|
| date | Thu, 08 Sep 2011 14:01:26 -0400 |
| parents | 6cf8fe6afa3b |
| children | |
| files | btquotients.sage |
| diffstat | 1 files changed, 201 insertions(+), 115 deletions(-) [+] |
line wrap: on
line diff
--- a/btquotients.sage Wed Jul 27 09:21:09 2011 +0200 +++ b/btquotients.sage Thu Sep 08 14:01:26 2011 -0400 @@ -46,7 +46,7 @@ extrapow: gets added to the power attribute, and it is used for the Hecke action. """ -class EdgeDecomposition(SageObject): +class DoubleCosetReduction(SageObject): def __init__(self,Y,x,extrapow=0): e1=Y._BT.edge(x) try: @@ -78,14 +78,23 @@ def igamma(self,iota=None): Y=self._parent if iota is None: + # The default behavior: use the precision of the bruhat tits quotient (usually too low) iota=Y.iota - if(self._igamma_prec>=Y._prec): - return self._cached_igamma - self._igamma_prec=Y._prec + prec=Y._prec + else: + try: + # The user wants higher precision + prec=ZZ(iota) + iota=Y.get_iota(prec) + print 'prec=',prec + print 'iota=',iota + except TypeError: + # The user knows what she is doing, so let it go + return iota(self.gamma) + if(self._igamma_prec<prec): + self._igamma_prec=prec self._cached_igamma=iota(self.gamma) - return self._cached_igamma - else: - return iota(self.gamma) + return self._cached_igamma def t(self): Y=self._parent @@ -106,7 +115,7 @@ class BruhatTitsTree: def __init__(self,p): self._p=p - self._MQp=None + self._MQp=MatrixSpace(Qp(self._p),2,2) self._Mat_22=MatrixSpace(ZZ,2,2) self._prec=-1 @@ -167,12 +176,12 @@ try: return self._edges_leaving_origin except: p=self._p - self._edges_leaving_origin=[[self._Mat_22([0,-1,p,0]),self.target([0,-1,p,0])]] - self._edges_leaving_origin.extend([[self._Mat_22([p,i,0,1]),self.target([p,i,0,1])] for i in range(p)]) + self._edges_leaving_origin=[self._Mat_22([0,-1,p,0])] + self._edges_leaving_origin.extend([self._Mat_22([p,i,0,1]) for i in range(p)]) return self._edges_leaving_origin def leaving_edges(self,M): - return [[self.edge(M*A[0]),self.vertex(M*A[1])] for A in self.edges_leaving_origin()] + return [self.edge(M*A) for A in self.edges_leaving_origin()] def opposite(self,M): x=copy(M) @@ -181,7 +190,7 @@ return self.edge(x) def entering_edges(self,M): - return [[self.opposite(e[0]),M] for e in self.leaving_edges(M)] + return [self.opposite(e) for e in self.leaving_edges(M)] r""" Given a list of edges (which we think of finite union of open balls) @@ -201,12 +210,11 @@ edge=self._Mat_22(edge) origin=self.origin(edge) newE=self.leaving_edges(self.target(edge)) - newEgood.extend([e[0] for e in newE if e[1]!=origin]) + newEgood.extend([e for e in newE if self.target(e)!=origin]) return self.subdivide(newEgood,level-1) def get_ball(self,center=1,level=1): - newE=self.leaving_edges(center) - return self.subdivide([e[0] for e in newE],level) + return self.subdivide(self.leaving_edges(center),level) def whereis(self,z): #Assume z belongs to some extension of QQp. @@ -275,10 +283,10 @@ class BTQuotient(SageObject,UniqueRepresentation): @staticmethod - def __classcall__(cls,p,Nminus,Nplus=1,backend='magma'): + def __classcall__(cls,p,Nminus,Nplus=1,backend='sage'): return super(BTQuotient,cls).__classcall__(cls,p,Nminus,Nplus,backend) - def __init__(self,p,Nminus,Nplus=1,backend='magma'): + def __init__(self,p,Nminus,Nplus=1,backend='sage'): Nminus=Integer(Nminus) Nplus=Integer(Nplus) p=Integer(p) @@ -300,7 +308,7 @@ if(backend=='magma' or not self._Nminus.is_prime() or self._Nplus!=1 or self._p==2): try: self._magma=magma - magmap=magma.p + magmap=self._magma(p) except RuntimeError: raise NotImplementedError,'Sage does not know yet how to work with the kind of orders that you are trying to use. Try installing Magma first and set it up so that Sage can use it.' self._backend='magma' @@ -314,6 +322,7 @@ self._cached_paths=dict() self._cached_decomps=dict() self._cached_equivalent=dict() + self._CM_points=dict() self._V=(QQ**4).ambient_module().change_ring(ZZ) self._Mat_44=MatrixSpace(ZZ,4,4) @@ -336,11 +345,17 @@ self._S=Graph(ndict['_S'],multiedges=True,weighted=True) def _repr_(self): - return "Quotient of the Bruhat Tits tree of GL_2(QQ_%s) with discriminant %s and level %s at the prime %s"%(self.Nminus().factor(),self.Nplus().factor(),self.prime()) + return "Quotient of the Bruhat Tits tree of GL_2(QQ_%s) with discriminant %s and level %s"%(self.prime(),self.Nminus().factor(),self.Nplus().factor()) def _latex_(self): return "X(%s,%s)\\otimes_{\\mathbb{Z}} \\mathbb{F}_{%s}"%(latex(self.level().factor()),latex(self.Nplus().factor()),latex(self.prime())) + def get_boundary(self): + try: return self._boundary + except AttributeError: + self._compute_quotient() + return self._boundary + def get_vertex_list(self): try: return self._vertex_list except AttributeError: @@ -357,7 +372,8 @@ def genus(self): Nplus=self._Nplus lev=self._p*self._Nminus - #Compute the genus using the genus formulas + # Compute the genus using the genus formulas + # Source: Rotger, V. "Non-elliptic Shimura curves of genus one", Proposition 5.2 e4=1 e3=1 mu=Nplus @@ -414,7 +430,7 @@ disc=D if(D%4!=1): disc*=4 - ff=self._level.factor() + ff=self.level().factor() for f in ff: if kronecker_symbol(disc,f[0])!=-1: return False @@ -527,12 +543,21 @@ def iota(self,g,prec=None): return Matrix(self._R,2,2,self.get_Iota(prec)*g) + def get_iota(self,prec=None): + return lambda g:Matrix(self._R,2,2,self.get_Iota(prec)*g) + def get_edge_stabs(self): try: return self._edge_stabs except AttributeError: - self._edge_stabs=[self._edge_stabilizer(e.rep) for e in self.get_edge_list()] + self._edge_stabs=[self._stabilizer(e.rep,as_edge=True) for e in self.get_edge_list()] return self._edge_stabs + def get_vertex_stabs(self): + try: return self._vertex_stabs + except AttributeError: + self._vertex_stabs=[self._stabilizer(e.rep,as_edge=False) for e in self.get_vertex_list()] + return self._vertex_stabs + def get_quaternion_algebra(self): try: return self._A except AttributeError: pass @@ -581,18 +606,15 @@ O_units.append(vec) return O_units - @cached_method - def get_CM_points(self,O,prec=None,ignore_units=False): + def get_CM_points(self,O,ignore_units=False): p=self._p - if prec is None: - prec=self._prec - if(not isinstance(O,AbsoluteOrder)): + try: disc=Integer(O) R = QQ['x'] f = R([-disc, 0, 1]) K=NumberField(f, 'sq', check=False) O=K.maximal_order() - else: + except TypeError: disc=O.discriminant() K=O.fraction_field() try: all_elts=self._CM_points[disc] @@ -648,6 +670,7 @@ Kp=Qq(p**2,prec=prec,names='g') W=[] for m in all_elts_split: + m.set_immutable() v=self._fixed_points(m,Kp) W.append(v[0]) return W @@ -657,13 +680,13 @@ and each entry consists of the corresponding data for the matrix [p,a,0,1] where a varies from 0 to p-1. The data is a tuple (acter,edge_images), with edge images being of type - EdgeDecomposition + DoubleCosetReduction """ @cached_method def _get_Up_data(self): E=self.get_edge_list() vec_a=self._BT.subdivide([1],1) - return [[alpha.inverse(),[EdgeDecomposition(self,e.rep*alpha) for e in E]+[EdgeDecomposition(self,e.opposite.rep*alpha) for e in E]] for alpha in vec_a] + return [[alpha.inverse(),[DoubleCosetReduction(self,e.rep*alpha) for e in E]+[DoubleCosetReduction(self,e.opposite.rep*alpha) for e in E]] for alpha in vec_a] @cached_method def _get_atkin_lehner_data(self,q): @@ -682,7 +705,7 @@ try: x=self.iota(beta1) nn=ceil(x.determinant().valuation()) - T=[beta1,[EdgeDecomposition(self,x.adjoint()*e.rep,extrapow=nn) for e in E]] + T=[beta1,[DoubleCosetReduction(self,x.adjoint()*e.rep,extrapow=nn) for e in E]] success=True except PrecisionError: self.increase_precision(10) @@ -732,7 +755,7 @@ try: x=self.iota(v1)*alphamat nn=ceil(x.determinant().valuation()) - T.append([v1,[EdgeDecomposition(self,x.adjoint()*e.rep,extrapow=nn) for e in E]]) + T.append([v1,[DoubleCosetReduction(self,x.adjoint()*e.rep,extrapow=nn) for e in E]]) success=True except PrecisionError: self.increase_precision(10) @@ -790,19 +813,20 @@ m=self._mat_p001 new_v,h=vertex_with_hash(v) chain=[] + boundary=self.get_boundary() while h[1]!=0 or h[0]<0: - if self._boundary.has_key(h): - return chain,self._boundary[h] + if boundary.has_key(h): + return chain,boundary[h] chain.append(new_v) new_v,h=vertex_with_hash(new_v*m) - if self._boundary.has_key(h): - return chain,self._boundary[h] + if boundary.has_key(h): + return chain,boundary[h] chain.append(new_v) a=h[0] while(a>0): a-=1 - if self._boundary.has_key((a,0)): - return chain,self._boundary[(a,0)] + if boundary.has_key((a,0)): + return chain,boundary[(a,0)] new_v=self._Mat_22([new_v[0,0]/self._p,0,0,1]) new_v.set_immutable() chain.append(new_v) @@ -825,11 +849,11 @@ OM=self.get_eichler_order_quadmatrix() return E,E*OM*E.transpose() - def _edge_stabilizer(self,e): + def _stabilizer(self,e,as_edge=True): p=self._p m=valuation(e.determinant(),p) twom=2*m - E,A = self._find_lattice(e,e,True,twom) + E,A = self._find_lattice(e,e,as_edge,twom) n_units=len(self.get_units_of_order()) ## Using PARI to get the shortest vector in the lattice (via LLL) try: @@ -840,9 +864,9 @@ mat=v[2].python().transpose() n_vecs=mat.nrows() - B=self.get_eichler_order_basis() if(n_vecs<=1): - return 1 + return [[self._B_one,0,False]] + #return 1 stabs=[] for jj in range(n_vecs): vec=Matrix(ZZ,1,4,list(mat.row(jj))) @@ -854,7 +878,6 @@ wt=w.transpose() wt.set_immutable() stabs.append([wt,m,x!=p**m]) - #self._elements.add(wt) return stabs def _are_equivalent(self,v1,v2,as_edges=False,twom=None,comp=False): @@ -880,8 +903,8 @@ if(nrd.valuation(p)==twom and nrd.is_power_of(p**2)): A=(Matrix(ZZ,1,4,[v[2][ii,0].python() for ii in range(4)])*E).transpose() A.set_immutable() - if(comp): - self._elements.add(A) + if comp and not as_edges: + self._elements.add(A) self._cached_equivalent[(v1,v2,as_edges)]=(A,m) return (A,m) else: @@ -949,20 +972,20 @@ boundary=dict({(0,0):v0}) self._num_edges=0 num_verts+=1 - B_one=self._are_equivalent(v0.rep,v0.rep,True,0,comp=True) + self._B_one=self._are_equivalent(v0.rep,v0.rep,True,0,comp=True) while len(V)>0: v=V.pop(0) E=self._BT.leaving_edges(v.rep) #print 'V = %s, E = %s, G = %s (target = %s)'%(num_verts,self._num_edges,1+self._num_edges-num_verts,self.genus()) for e in E: new_edge=True - edge_det=e[0].determinant() + edge_det=e.determinant() edge_valuation=edge_det.valuation(p) - g,e1=self._find_equivalent_edge(e[0],v.leaving_edges,comp=True,valuation=edge_valuation) + g,e1=self._find_equivalent_edge(e,v.leaving_edges,comp=True,valuation=edge_valuation) if(g): e1.links.append(g) else: - target=e[1] + target=self._BT.target(e) target.set_immutable() new_det=target.determinant() new_valuation=new_det.valuation(p) @@ -980,11 +1003,11 @@ boundary[self._BT.vertex_with_hash(v1.rep)[1]]=v1 # Add the edge to the list - new_e=Edge(self,self._num_edges,e[0],v,v1,[],[],edge_det,edge_valuation) - new_e.links.append(B_one) + new_e=Edge(self,self._num_edges,e,v,v1,[],[],edge_det,edge_valuation) + new_e.links.append(self._B_one) # Find the opposite edge - opp=self._BT.opposite(e[0]) + opp=self._BT.opposite(e) # Add the edge to the graph S.add_edge(v.label,v1.label,self._num_edges) @@ -1024,7 +1047,10 @@ @cached_method def _fixed_points(self,x,K): - a=x[0,0],b=x[0,1],c=x[1,0],d=x[1,1] + a=x[0,0] + b=x[0,1] + c=x[1,0] + d=x[1,1] t=x.trace() n=x.determinant() A=K(a-d) @@ -1090,11 +1116,18 @@ import sage.modular.hecke.hecke_operator class HarmonicCocycle(HeckeModuleElement): - def __init__(self,_parent,vec=None): - HeckeModuleElement.__init__(self,_parent,vec) + def __init__(self,_parent,vec=None,from_values=False): + HeckeModuleElement.__init__(self,_parent,None) self._parent=_parent + if(from_values): + self._R=_parent._U._R + self._wt=_parent._k + self._nE=len(_parent._E) + self._F=copy(vec) + return + if(isinstance(vec,HarmonicCocycle)): #The first argument is just a modular form that we copy as is self._R=vec._R @@ -1103,8 +1136,8 @@ self._F=[_parent._U.element_class(_parent._U,vec._F[ii],quick=True) for ii in range(self._nE)] return - try: vec=Sequence(vec) - except TypeError: pass + try: vec=vec.list() + except AttributeError: pass # When vec contains coordinates for the basis if(isinstance(vec,(list,tuple)) and vec[0].parent() is self._parent._R): try: @@ -1128,22 +1161,6 @@ E=_parent._E self._F=[_parent._U.element_class(_parent._U,vec._F[ii]).l_act_by(E[ii].rep) for ii in range(self._nE)] return - elif(isinstance(vec,dict)): - assert(vec['p'])==_parent._X._p - assert(vec['t'])==_parent._U._depth - assert(vec['level'])==_parent._X._level - assert(vec['k']==_parent._U.weight()+2) - self._R=_parent._U._R - self._wt=_parent._k - self._nE=len(_parent._E) - self._F=[_parent._U(c) for c in mydata['value']] - return - elif(isinstance(vec,(list,tuple)) and vec[0].parent() is _parent._U): - self._R=_parent._U._R - self._wt=_parent._k - self._nE=len(_parent._E) - self._F=copy(vec) - return else: try: veczp=_parent._U._R(vec) @@ -1153,21 +1170,24 @@ self._F=[_parent._U(veczp) for ii in range(self._nE)] return except: + print vec.parent() print vec raise ValueError def _add_(self,g): #Should ensure that self and g are modular forms of the same weight and on the same curve - vec=[self._F[e]+g._F[e] for e in range(self._nE)] - return HarmonicCocycle(self._parent,vec) + #vec=[self._F[e]+g._F[e] for e in range(self._nE)] + return HarmonicCocycle(self._parent,self.element()+g.element()) def _sub_(self,g): #Should ensure that self and g are modular forms of the same weight and on the same curve - vec=[self._F[e]-g._F[e] for e in range(self._nE)] - return HarmonicCocycle(self._parent,vec) + #vec=[self._F[e]-g._F[e] for e in range(self._nE)] + return HarmonicCocycle(self._parent,self.element()+g.element()) def _rmul_(self,a): #Should ensure that 'a' is a scalar - return HarmonicCocycle(self._parent,[a*self._F[e] for e in range(self._nE)]) + #return HarmonicCocycle(self._parent,[a*self._F[e] for e in range(self._nE)]) + return HarmonicCocycle(self._parent,a*self.element()) + def _repr_(self): tmp='Harmonic cocycle on '+str(self._parent)+':\n' @@ -1190,19 +1210,19 @@ else: return min([self._F[e].valuation() for e in range(self._nE)]) - def list(self): + def _compute_element(self): R=self._R A=self._parent.basis_matrix().transpose() B=Matrix(R,self._nE*(self._parent._k-1),1,[self._F[e]._val[ii,0] for e in range(self._nE) for ii in range(self._parent._k-1) ]) res=(A.solve_right(B)).transpose() - return res.list() + return res #In HarmonicCocycle def evaluate(self,e1): X=self._parent._X p=X._p - u=EdgeDecomposition(X,e1) - return (u.sign()*self._F[u.label]).l_act_by(u.igamma(self.iota)*(p**(-u.power))) + u=DoubleCosetReduction(X,e1) + return (u.sign()*self._F[u.label]).l_act_by(u.igamma(self._parent.iota)*(p**(-u.power))) #In HarmonicCocycle def riemann_sum(self,f,center=1,level=0,E=None): @@ -1254,6 +1274,9 @@ def __init__(self,X,k,prec=None,basis_matrix=None): self._k=k self._X=X + self._E=self._X.get_edge_list() + self._V=self._X.get_vertex_list() + if(prec is None): self._is_exact=True self._prec=None @@ -1281,8 +1304,6 @@ self.__matrix.set_immutable() self.__rank=self.__matrix.nrows() AmbientHeckeModule.__init__(self, self._R, self.__rank, self._X.prime()*self._X.Nplus()*self._X.Nminus(), weight=self._k) - self._E=self._X.get_edge_list() - self._V=self._X.get_vertex_list() self._populate_coercion_lists_() def rank(self): @@ -1296,15 +1317,14 @@ return self.rank()==1 def _repr_(self): - s='Space of harmonic cocycles of weight '+str(self._k)+' on '+str(self._X) - return s + return 'Space of harmonic cocycles of weight %s on %s'%(self._k,self._X) def _latex_(self): s='\\text{Space of harmonic cocycles of weight }'+latex(self._k)+'\\text{ on }'+latex(self._X) return s def _an_element_(self): - return HarmonicCocycle(self,1) + return self.basis()[0] def _coerce_map_from_(self, S): """ @@ -1335,7 +1355,7 @@ #Code how to coherce x into the space #Admissible values of x? if isinstance(x,(HarmonicCocycle,pAutomorphicForm)): - return HarmonicCocycle(self,x) + return HarmonicCocycle(self,x.element()) return HarmonicCocycle(self,x) def is_exact(self): @@ -1369,7 +1389,7 @@ S=self._X.get_edge_stabs() p=self._X._p d=self._k-1 - for e in filter(lambda e:S[e.label]!=1,self._E): + for e in self._E: try: g=filter(lambda g:g[2],S[e.label])[0] C=self._U.l_matrix_representation(self.iota(g[0])) @@ -1415,7 +1435,7 @@ for jj in range(len(self._E)): t=d1[jj] tmp[jj]+=(t.sign()*f._F[t.label]).l_act_by(p**(-t.power)*mga*t.igamma(self.iota)) - return HarmonicCocycle(self,tmp) + return HarmonicCocycle(self,tmp,from_values=True) def __apply_hecke_operator(self,l,f): R=self._R @@ -1434,7 +1454,7 @@ for jj in range(len(self._E)): t=d1[jj] tmp[jj]+=(t.sign()*f._F[t.label]).l_act_by(p**(-t.power)*mga*t.igamma(self.iota)) - return HarmonicCocycle(self,[factor*x for x in tmp]) + return HarmonicCocycle(self,[factor*x for x in tmp],from_values=True) def _compute_atkin_lehner_matrix(self,d): res=self.__compute_operator_matrix(lambda f:self.__apply_atkin_lehner(d,f)) @@ -1452,9 +1472,14 @@ for rr in range(len(basis)): g=T(basis[rr]) B.set_block(0,rr,Matrix(R,len(self._E)*(self._k-1),1,[g._F[e]._val[ii,0] for e in range(len(self._E)) for ii in range(self._k-1) ])) - res=(A.solve_right(B)).transpose() - res.set_immutable() - return res + try: + res=(A.solve_right(B)).transpose() + res.set_immutable() + return res + except ValueError: + print A + print B + raise ValueError class HarmonicCocyclesSubmodule(sage.modular.hecke.submodule.HeckeSubmodule,HarmonicCocycles): """ @@ -1482,7 +1507,6 @@ class pAutomorphicForm(ModuleElement): def __init__(self,_parent,vec,quick=False): - ModuleElement.__init__(self,_parent) self._parent=_parent self._nE=2*len(_parent._E) # We record the values at the opposite edges @@ -1502,10 +1526,17 @@ assert(2*len(vec._F)==self._nE) assert(isinstance(_parent._U,OCVn)) E=self._parent._E + + + MMM=vec._parent._U.element_class + tmp=[] for ii in range(len(vec._F)): - self._F.append(_parent._U(self._MElement(vec._parent._U,vec._F[ii]).l_act_by(E[ii].rep.inverse()))) + newtmp=MMM(vec._parent._U,vec._F[ii]).l_act_by(E[ii].rep.inverse()) + tmp.append(newtmp) + self._F.append(_parent._U(newtmp)) + A=Matrix(QQ,2,2,[0,-1/_parent._X._p,-1,0]) for ii in range(len(vec._F)): - self._F.append(_parent._U((-1*(self._MElement(vec._parent._U,self._F[ii])).r_act_by(Matrix(QQ,2,2,[0,-1/_parent._X._p,-1,0]))))) + self._F.append(_parent._U(-1*tmp[ii].r_act_by(A))) self._make_invariant() elif(isinstance(vec,list) and len(vec)==self._nE): @@ -1535,11 +1566,12 @@ M=[e.rep for e in E]+[e.opposite.rep for e in E] lS=len(S) assert(2*len(S)==self._nE) + MMM=self._parent._U.element_class for ii in range(self._nE): Si=S[ii%lS] - if(Si!=1): - x=self._MElement(self._F[ii].parent(),self._F[ii]._val,quick=True) - self._F[ii]=self._MElement(self._F[ii].parent(),0) + if(any([v[2] for v in Si])): + x=MMM(self._F[ii].parent(),self._F[ii]._val,quick=True) + self._F[ii]=MMM(self._F[ii].parent(),0) s=QQ(0) m=M[ii] for v in Si: @@ -1554,7 +1586,6 @@ mydata=dict([]) mydata['p']=self._parent._X._p mydata['t']=self._parent._U._depth - mydata['level']=self._parent._X._level mydata['weight']=self._parent._U.weight() mydata['value']=[c._val for c in self._F] return mydata @@ -1563,7 +1594,6 @@ mydata=load(s) assert(mydata['p'])==self._parent._X._p assert(mydata['t'])==self._parent._U._depth - assert(mydata['level'])==self._parent._X._level assert(mydata['weight']==self._parent._U.weight()) self._F=[self._parent._U(c) for c in mydata['value']] @@ -1581,7 +1611,7 @@ def evaluate(self,e1): X=self._parent._X - u=EdgeDecomposition(X,e1) + u=DoubleCosetReduction(X,e1) return (self._F[u.label+u.parity*self._nE/2]).r_act_by((u.t())*X._p**(u.power)) def _rmul_(self,a): @@ -1594,7 +1624,7 @@ p=Y._p Tf=[] for e in E: - u=EdgeDecomposition(Y,alpha*e.rep) + u=DoubleCosetReduction(Y,alpha*e.rep) r=u.t()*p**(-(u.power)) if(u.parity==0): tmp=self._F[u.label].r_act_by(r) @@ -1602,7 +1632,7 @@ tmp=self._F[u.label+len(E)].r_act_by(r) Tf.append(tmp) for e in E: - u=EdgeDecomposition(Y,alpha*e.opposite.rep) + u=DoubleCosetReduction(Y,alpha*e.opposite.rep) r=u.t()*gg*p**(-(u.power)) if(u.parity==0): tmp=self._F[u.label].r_act_by(r) @@ -1631,17 +1661,18 @@ def improve(self): MMM=self._parent - h2=MMM.__apply_Up_operator(self,True) + h2=MMM._apply_Up_operator(self,True) ii=0 current_val=0 - old_val=-1 - while(current_val<MMM._prec and current_val>old_val): + old_val=-Infinity + init_val=self.valuation() + while(current_val>old_val): old_val=current_val ii+=1 + self._F=[self._parent._U(c) for c in h2._F] + h2=MMM._apply_Up_operator(self,True) + current_val=(h2-self).valuation()-init_val print ii,"th iteration, precision=",current_val - self._F=[self._parent._U(c) for c in h2._F] - h2=MMM.__apply_Up_operator(self,True) - current_val=(h2-self).valuation() self._F=[self._parent._U(c) for c in h2._F] def integrate(self,f,center=1,level=0,method='moments'): @@ -1839,7 +1870,7 @@ def iota(self,g): return self._X.iota(g,self._prec) - def __apply_Up_operator(self,f,scale=False): + def _apply_Up_operator(self,f,scale=False): HeckeData=self._X._get_Up_data() if(scale==False): factor=(self._X._p)**(self._U.weight()/2) @@ -1910,7 +1941,8 @@ if(B is None): B=self._parent.basis() A=Matrix(self._parent._R,self._parent.dimension(),self._parent.dimension(),[[b._val[ii,0] for b in B] for ii in range(self._depth)]) - return A.solve_right(self._val) + tmp=A.solve_right(self._val) + return tmp def _add_(self,y): val=self._val+y._val @@ -1937,7 +1969,7 @@ factor=R.prime()**(-t) try: x=self._parent._powers[(factor*a,factor*b,factor*c,factor*d)] - return self.__class__(self._parent,extrafactor*factor**(-self._n)*x*self._val,quick=True) + return self.__class__(self._parent,(extrafactor*factor**(-self._n))*(x*self._val),quick=True) except KeyError: tmp=self._parent._get_powers_and_mult(factor*a,factor*b,factor*c,factor*d,extrafactor*factor**(-self._n),self._val) return self.__class__(self._parent,tmp,quick=True) @@ -1989,7 +2021,7 @@ def valuation(self,l=None): if not self._parent.is_exact(): - if(l!=self._parent._R.prime()): + if(not l is None and l!=self._parent._R.prime()): raise ValueError, "This function can only be called with the base prime" return min([self._val[ii,0].valuation() for ii in range(self._depth)]) else: @@ -2091,6 +2123,7 @@ def weight(self): return self._n + def l_matrix_representation_list(self,v,B=None): if B is None: B=self.basis() @@ -2214,3 +2247,56 @@ n.append(0) wd=prod([v[x] for x in n]) yield wd + +# # Assume properly formatted +# # [(x,xinv),(y,yinv),(z,zinv)] +# def enum_freegp_words(v): +# mydict=dict() +# L=2*len(v) +# n=[] +# while(True): +# add_new=True +# for jj in range(len(n)): +# newnj=(n[jj][0]+1)%L +# if n[jj][2]==0: +# n[jj]=[newnj,n[jj][1],1] +# else: +# n[jj]=[newnj,(n[jj][1]+1)%len(v),0] +# if n[jj][0]!=0: +# add_new=False +# break +# if add_new: +# n.append([0,0,0]) +# +# is_legal=len(n)==1 or not any([n[jj][1]==n[jj+1][1] and n[jj][2]!=n[jj+1][2] for jj in range(len(n)-1)]) +# while not is_legal: +# add_new=True +# for jj in range(len(n)): +# newnj=(n[jj][0]+1)%L +# if n[jj][2]==0: +# n[jj]=[newnj,n[jj][1],1] +# else: +# n[jj]=[newnj,(n[jj][1]+1)%len(v),0] +# if n[jj][0]!=0: +# add_new=False +# break +# if add_new: +# n.append([0,0,0]) +# is_legal=len(n)==1 or not any([n[jj][1]==n[jj+1][1] and n[jj][2]!=n[jj+1][2] for jj in range(len(n)-1)]) +# +# #wd=''.join([v[floor(x/2)][x%2] for x in n]) +# try: +# a=mydict[tuple([n[ii][0] for ii in range(len(n)-1)])] +# except: +# a=1 +# x=n[len(n)-1] +# wd=a*v[x[1]][x[2]] +# +# mydict[tuple([n[ii][0] for ii in range(len(n))])]=wd +# #wd=prod([v[x[1]][x[2]] for x in n]) +# yield (wd,len(n),n) +# +# #I=enum_freegp_words([('u','v'),('p','d')]) +# I=enum_freegp_words([(132,1/132),(411,1/411)]) +# I10=itertools.islice(I,100000) +# # time Wnew=[x for x in I10]