# Ticket #12561: 12561.patch

File 12561.patch, 36.1 KB (added by bsinclai, 16 months ago)

# HG changeset patch
# User Brian Sinclair <basincla@uncg.edu>
# Date 1329985637 28800

new file mode 100644
 - # Initialization file for p-adic polynomial factoring

diff --git a/sage/rings/polynomial/padics/factor/associatedfactor.py b/sage/rings/polynomial/padics/factor/associatedfactor.py
new file mode 100644
 - from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing from sage.rings.finite_rings.constructor import GF from sage.matrix.constructor import Matrix from sage.rings.polynomial.padics.factor.frameelt import FrameElt from sage.rings.infinity import infinity class AssociatedFactor: def __init__(self,segment,rho,rhoexp): self.segment = segment self.rho = rho self.rhoexp = rhoexp self.Fplus = self.rho.degree() if self.segment.frame.is_first(): # In the first frame, so FFbase is the residue class field of O self.FFbase = self.segment.frame.R else: # Not the first frame self.FFbase = self.segment.frame.prev.FF if self.segment.frame.E * self.segment.frame.F * self.segment.Eplus * self.Fplus == self.segment.frame.Phi.degree(): # Polynomial is irreducible return if self.Fplus == 1: self.FF = self.FFbase self.FFz = PolynomialRing(self.FF,'z'+str(self.segment.frame.depth)) # rho is linear delta is the root of rho self.delta = self.rho.roots()[0][0] else: self.FF = GF(self.FFbase.order()**self.Fplus,'a'+str(self.segment.frame.depth)) self.FFz = PolynomialRing(self.FF,'z'+str(self.segment.frame.depth)) self.FFbase_gamma = (self.FFz(self.FFbase.modulus())).roots()[0][0] FFrho = self.FFz([self.FFbase_elt_to_FF(a) for a in list(rho)]) self.gamma = FFrho.roots()[0][0] basis = [(self.gamma**j*self.FFbase_gamma**i).polynomial() for j in range(0,self.Fplus) for i in range(0,self.FFbase.degree())] self.basis_trans_mat = Matrix([self.FF(b)._vector_() for b in basis]) def FF_elt_to_FFbase_vector(self,a): if self.segment.frame.is_first() and self.Fplus == 1: return a elif self.Fplus == 1: return self.segment.frame.prev.FF_elt_to_FFbase_vector(a) else: basedeg = self.FFbase.degree() avec = self.FF(a)._vector_() svector = self.basis_trans_mat.solve_left(Matrix(self.FF.prime_subfield(),avec)) s_list = svector.list() s_split = [ s_list[i*basedeg:(i+1)*basedeg] for i in range(0,self.Fplus)] s = [sum([ss[i]*self.FFbase.gen()**i for i in range(0,len(ss))]) for ss in s_split] return s def FFbase_elt_to_FF(self,b): if self.segment.frame.is_first(): return b elif self.Fplus == 1: return self.segment.frame.prev.FFbase_elt_to_FF(b) elif self.segment.frame.F == 1: return b * self.FFbase_gamma else: bvec = b._vector_() return sum([ bvec[i]*self.FFbase_gamma**i for i in range(len(bvec))]) def __repr__(self): return "AssociatedFactor of rho "+repr(self.rho) def lift(self,delta): # FrameElt representation of a lift of delta if self.segment.frame.F == 1: return FrameElt(self.segment.frame,self.segment.frame.Ox(delta)) elif self.segment.frame.prev.Fplus == 1: return FrameElt(self.segment.frame,self.segment.frame.prev.lift(delta),this_exp=0) else: dvec = self.segment.frame.prev.FF_elt_to_FFbase_vector(delta) return sum([self.segment.frame.prev.gamma_frameelt**i*FrameElt(self.segment.frame,self.segment.frame.prev.lift(dvec[i]),this_exp=0) for i in range(len(dvec)) if dvec[i] != 0]) def next_frame(self,length=infinity): from frame import Frame if self.Fplus == 1 and self.segment.Eplus == 1: next = Frame(self.segment.frame.Phi,self.segment.frame.prev,self.segment.frame.iteration) else: next = Frame(self.segment.frame.Phi,self,self.segment.frame.iteration) self.next = next self.gamma_frameelt = FrameElt(next,self.segment.psi**-1,self.segment.Eplus) if self.Fplus == 1 and self.segment.frame.F == 1: next_phi = self.segment.frame.phi**self.segment.Eplus-(self.segment.psi.polynomial()*self.segment.frame.Ox(self.delta)) self.reduce_elt = FrameElt(next,self.segment.psi*self.lift(self.delta),0) next.seed(next_phi,length=length) elif self.Fplus == 1 and self.segment.Eplus == 1: delta_elt = self.lift(self.delta) next_phi_tail = self.segment.psi*delta_elt.reduce() next_phi = self.segment.frame.phi-next_phi_tail.polynomial() self.reduce_elt = FrameElt(next,next_phi_tail,0) next.seed(next_phi,length=length) else: lifted_rho_coeffs = [self.lift(r) for r in list(self.rho)] lifted_rho_coeffs_with_psi = [FrameElt(next,(self.segment.psi**(self.Fplus-i)*lifted_rho_coeffs[i]).reduce(),0) for i in range(len(lifted_rho_coeffs))] phi_elt = FrameElt(next,self.segment.frame.Ox(1),1) next_phi_tail = sum([phi_elt**(self.segment.Eplus*i)*lifted_rho_coeffs_with_psi[i] for i in range(len(lifted_rho_coeffs_with_psi)-1)]) next_phi = (phi_elt**(self.segment.Eplus*self.Fplus)+next_phi_tail).polynomial() self.reduce_elt = FrameElt(next)+(-next_phi_tail) # that is -next_phi_tail next.seed(next_phi,length=length) return next

diff --git a/sage/rings/polynomial/padics/factor/factoring.py b/sage/rings/polynomial/padics/factor/factoring.py
new file mode 100644
 - from sage.rings.padics.factory import ZpFM from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing from sage.rings.polynomial.padics.factor.frame import Frame def pfactortree(Phi): r""" Return a factorization of Phi This is accomplished by constructing a tree of Frames of approximations of factors INPUT:: - Phi -- squarefree, monic padic polynomial with fixed precision coefficients OUTPUT:: - list -- polynomial factors of Phi EXAMPLES:: Factoring polynomials over Zp(2)[x]:: sage: from sage.rings.polynomial.padics.factor.factoring import jorder4,pfactortree sage: f = ZpFM(2,24,'terse')['x']( (x^32+16)*(x^32+16+2^16*x^2)+2^34 ) sage: pfactortree(f) # long (3.8s) [(1 + O(2^24))*x^64 + (65536 + O(2^24))*x^34 + (32 + O(2^24))*x^32 + (1048576 + O(2^24))*x^2 + (256 + O(2^24))] See the irreducibility of x^32+16 in Zp(2)[x]:: sage: pfactortree(ZpFM(2)['x'](x^32+16)) [(1 + O(2^20))*x^32 + (2^4 + O(2^20))] Test the irreducibility of test polynomial jorder4 for Zp(3):: sage: len(pfactortree(jorder4(3))) == 1 True Factor jorder4 for Zp(5) and Zp(7) and check that the products return the original:: sage: ff = pfactortree(jorder4(5)) sage: prod(ff) == jorder4(5) True sage: ff = pfactortree(jorder4(7)) sage: prod(ff) == jorder4(7) True AUTHORS:: - Brian Sinclair (2012-02-22): initial version """ from sage.misc.flatten import flatten def followsegment(next,Phi): """ Returns next if it corresponds to an irreducible factor of $\Phi$ and follows every branch if not """ # Handle the unlikely event that our approximation is actually a factor if next.phi_divides_Phi(): return next # With E potentially increased, Check to see if E*F == deg(Phi) (and thus Phi is irreducible) if next.E * next.F * next.polygon[0].Eplus == Phi.degree(): return next # With F potentially increaded, Check to see if E*F == deg(Phi) (and thus Phi is irreducible) if next.E * next.F * next.polygon[0].Eplus * next.polygon[0].factors[0].Fplus == Phi.degree(): return next # Check if we should begin Single Factor Lifting if sum([seg.length for seg in next.polygon]) == 1: return next return [[followsegment(fact.next_frame(fact.rhoexp+1),Phi) for fact in seg.factors] for seg in next.polygon] # Handle the situation that x is a factor of $\Phi(x)$ if Phi.constant_coefficient() == 0: x_divides = True Phi = Phi >> 1 else: x_divides = False # Construct and initialize the first frame (phi = x) next = Frame(Phi) next.first() # With E potentially increased, Check to see if E*F == deg(Phi) (and thus Phi is irreducible) if next.E * next.F * next.polygon[0].Eplus == Phi.degree(): return ([Phi.parent().gen()] if x_divides else []) + [Phi] # With F potentially increaded, Check to see if E*F == deg(Phi) (and thus Phi is irreducible) if next.E * next.F * next.polygon[0].Eplus * next.polygon[0].factors[0].Fplus == Phi.degree(): return ([Phi.parent().gen()] if x_divides else []) + [Phi] # Construct the next level of the tree by following every factor of the # residual polynomials of every Newton polygon segment in our frame tree = [[followsegment(fact.next_frame(fact.rhoexp+1),Phi) for fact in seg.factors] for seg in next.polygon] # tree contains the leaves of the tree of frames and each leaf corresponds # to an irreducible factor of Phi, so we flatten the list and start lifting flat = flatten(tree) # If we only have one leaf, Phi is irreducible, so we do not lift it. if len(flat) == 1: return ([Phi.parent().gen()] if x_divides else []) + [Phi] # quo_rem is faster than Improve, so Phi = f*g are specially handled if len(flat) == 2: if flat[0].phi.degree() < flat[1].phi.degree(): fact = Improve(flat[0]) return ([Phi.parent().gen()] if x_divides else []) + [fact,Phi.quo_rem(fact)[0]] else: fact = Improve(flat[1]) return ([Phi.parent().gen()] if x_divides else []) + [fact,Phi.quo_rem(fact)[0]] # Phi has more than two factors, so we lift them all return ([Phi.parent().gen()] if x_divides else []) + [Improve(frame) for frame in flatten(tree)] def jorder4(p): r""" Produce a particularly complicated example of polynomials for factorization over Zp(p) INPUT:: - p a prime number EXAMPLES:: sage: from sage.rings.polynomial.padics.factor.factoring import jorder4 sage: jorder4(3) (1 + O(3^20))*x^24 + (24 + O(3^20))*x^23 + (276 + O(3^20))*x^22 + (2048 + O(3^20))*x^21 + (11310 + O(3^20))*x^20 + (51180 + O(3^20))*x^19 + (201652 + O(3^20))*x^18 + (709092 + O(3^20))*x^17 + (2228787 + O(3^20))*x^16 + (6232484 + O(3^20))*x^15 + (15469950 + O(3^20))*x^14 + (34143276 + O(3^20))*x^13 + (67323664 + O(3^20))*x^12 + (119268300 + O(3^20))*x^11 + (190652502 + O(3^20))*x^10 + (275456480 + O(3^20))*x^9 + (359189415 + O(3^20))*x^8 + (420635664 + O(3^20))*x^7 + (438402286 + O(3^20))*x^6 + (400618284 + O(3^20))*x^5 + (313569267 + O(3^20))*x^4 + (203945072 + O(3^20))*x^3 + (105227142 + O(3^20))*x^2 + (38341248 + O(3^20))*x + (10912597 + O(3^20)) Input must be prime:: sage: jorder4(4) Traceback (most recent call last): ... ValueError: p must be prime """ K = ZpFM(p,20,print_mode='terse') Kx = PolynomialRing(K,names='x') x = Kx.gen() f1 = (x+1)**3+p; f2 = f1**2+p**2*(x+1); f3 = f2**2+4*p**3*f1*(x+1)**2; f4 = f3**2+20*p**2*f3*f2*(x+1)**2+64*p**9*f1; return f4; def Improve(frame,prec=2**20): r""" Lift a frame to a factor EXAMPLES:: sage: from sage.rings.polynomial.padics.factor.factoring import Improve sage: from sage.rings.polynomial.padics.factor.frame import Frame sage: Kx. = PolynomialRing(ZpFM(3,20,'terse')) sage: f = (x**3+x-1)*(x**2+1) sage: fr = Frame(f) sage: fr.first() sage: fr = fr.polygon[0].factors[0].next_frame() sage: fact = Improve(fr) ; fact (1 + O(3^20))*x + (904752403 + O(3^20)) sage: f % fact 0 """ def _reduce(poly,phi,d): """ returns poly mod phi and simplifies the denominator of poly """ poly = phi.parent()(poly) % phi if d != 1: g = min([d] + [p.valuation() for p in poly]) if g > 0: poly = poly.parent( [p >> g for p in poly] ) d = d - g return poly,d if frame.phi_divides_Phi(): return frame.phi LiftRing = ZpFM(frame.O.uniformizer(),2*frame.O.precision_cap()) Lifty = PolynomialRing(LiftRing,names='y') Phi = Lifty(frame.Phi) phi = Lifty(frame.phi) a0,a1 = frame.elt_phi_expansion()[0:2] Psi = frame.find_psi(-a1.valuation()) A0 = Psi * a0 A1 = Psi * a1 Psi,Psiden = Psi.polynomial(True) Psi = Lifty(Psi) C1inv = frame.polygon[0].factors[0].lift(1/(A1.residue())) C1inv,C1invden = C1inv.polynomial(True) C1inv,C1invden = _reduce(C1inv,phi,C1invden) A0,A0den = A0.polynomial(True) A0,A0den = _reduce(Lifty(A0),phi,A0den) C,Cden = _reduce(frame.Ox(A0*C1inv),phi,A0den+C1invden) phi = (phi + C) h = 2 oldphi = None while h < prec and phi != oldphi: oldphi = phi C1, C0 = Phi.quo_rem(phi) C0,C0den = _reduce((Psi*C0),phi,Psiden) C1,C1den = _reduce((Psi*C1),phi,Psiden) x1,x1den = _reduce((LiftRing(2)<<(C1den+C1invden))-C1*C1inv,phi,C1den+C1invden) C1inv,C1invden = _reduce(C1inv*x1,phi,C1invden+x1den) C,Cden = _reduce((C0*C1inv),phi,C0den+C1invden) phi = (phi + C) h = 2 * h return frame.Ox(phi)

diff --git a/sage/rings/polynomial/padics/factor/frame.py b/sage/rings/polynomial/padics/factor/frame.py
new file mode 100644
 - from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing from sage.rings.polynomial.padics.factor.frameelt import FrameElt, FrameEltTerm from sage.rings.polynomial.padics.factor.segment import Segment from sage.rings.infinity import infinity from sage.misc.functional import denominator class Frame: """Everything for one iteration of local field polynomial factorization""" def __init__(self,Phi,pre=None,iteration_count=0): self.prev = pre self.Phi = Phi self.O = Phi.base_ring() self.Ox = Phi.parent() self.x = self.Ox.gen() self.R = self.O.residue_class_field() self.Rz = PolynomialRing(self.R,names='z') self.phi = None self.iteration = iteration_count + 1 if self.is_first(): # that is self.prev == None self.E = 1 self.F = 1 self.depth = 0 else: self.E = self.prev_frame().E * self.prev.segment.Eplus self.F = self.prev_frame().F * self.prev.Fplus self.depth = self.prev_frame().depth + 1 def first(self): """ Computes all of the frame data based on being the root of the tree of frames """ self._phidividesPhi = False self.phi = self.x self._expansion_of_Phi() self.polygon = self._newton_polygon([e.valuation() for e in self._phiexpelt]) # list of segments def is_first(self): return self.prev == None def seed(self,phi,length=infinity): """ Computes all of the frame data with the new phi """ self.phi = phi self._phidividesPhi = False self._expansion_of_Phi(length=length) # If phi divides Phi, we return to the main factorization loop, where # this property should be checked for if self._phiexppoly[0] == 0: return else: self.polygon = self._newton_polygon([e.valuation() for e in self._phiexpelt]) # list of segments def find_psi(self,val): psielt = FrameElt(self) if self.prev == None: psielt.rep = [FrameEltTerm(psielt,self.O(1),val)] else: vphi = self.prev.segment.slope d = self.prev_frame().E vprime = val*d e = vphi * d psimod = denominator(e) s = 0 if not psimod == 1: s = vprime / e if denominator(s) == 1: s = s % psimod else: s = int(s % psimod) val = val - s * vphi psielt.rep = [FrameEltTerm(psielt,self.prev_frame().find_psi(val),s)] return psielt def _expansion_of_Phi(self,length=infinity): self._phiexppoly = [] if self.phi.degree() > self.Phi.degree(): self._phiexppoly = [self.Phi] q, r = self.Phi.quo_rem(self.phi) self._phiexppoly.append(r) while q != 0 and length > len(self._phiexppoly): q, r = q.quo_rem(self.phi) self._phiexppoly.append(r) self._phiexpelt = [FrameElt(self,a) for a in self._phiexppoly] def _newton_polygon(self,a,xoffset=0): if len(a) == 0: raise ValueError, "Cannot compute Newton polygon from empty list" slope = (a[0]-a[len(a)-1]) / (len(a)-1) verts = [(xoffset,a[0])] length = 0 for i in range(1,len(a)): y = a[0] - i*slope if a[i] == y: verts.append((xoffset+i,y)) length = i elif a[i] < y: verts = [(xoffset,a[0]),(xoffset+i,a[i])] slope = (a[0]-a[i]) / i length = i elif y < 0: if len(a[verts[len(verts)-1][0]-xoffset:]) == 0: return [Segment(self,slope,verts,length)] else: return [Segment(self,slope,verts,length)]+self._newton_polygon(a[verts[len(verts)-1][0]-xoffset:],verts[len(verts)-1][0]) return [Segment(self,slope,verts,length)] # Data Access Methods def prev_frame(self): if self.prev == None: return None else: return self.prev.segment.frame def elt_phi_expansion(self): return self._phiexpelt def phi_divides_Phi(self): """ Returns whether or not phi divides Phi Checks to see if the constant term of the phi-expansion of Phi is 0 """ return self._phiexppoly[0] == 0 def __repr__(self): if self.phi: return 'Frame with phi '+repr(self.phi) else: return 'Unseeded Frame regarding '+repr(self.Phi)

diff --git a/sage/rings/polynomial/padics/factor/frameelt.py b/sage/rings/polynomial/padics/factor/frameelt.py
new file mode 100644
 - class FrameElt: """Polynomials recursively represented by powers of the approximations $\phi_t(x)$ to a factor of $\Phi(x)$.""" def __init__(self,frame,a=None,this_exp=None): # deg(a) < deg(frame.phi)*frame.Eplus*frame.Fplus self.frame = frame if this_exp != None: # initializes a FrameElt of phi_t ^ this_exp * a self.rep = [FrameEltTerm(self,a,this_exp)] elif a == None: self.rep = [] elif frame.is_first(): self.rep = [FrameEltTerm(self,a,0)] else: b = _phi_expansion(a,self.frame.prev_frame().phi) self.rep = [FrameEltTerm(self,b[i],i) for i in range(len(b)) if b[i].is_zero() == False] def is_single_term(self): if len(self.rep) == 1: if self.frame.prev == None: return True else: return self.rep[0]._coefficient.is_single_term() return False def valuation(self,purge=True): if not purge: return min([a.valuation() for a in self.rep]) else: v = min([a.valuation() for a in self.rep]) self.rep = [a for a in self.rep if a.valuation() == v] return v def residue(self): if not self.is_reduced(): self = self.reduce() if self.frame.is_first(): if self.rep[0]._exponent == 0: # unable to coerce in Zq return self.frame.R(self.rep[0]._unit) else: return self.frame.R(0) Eplus = int(self.frame.prev.segment.Eplus) if self.frame.prev.Fplus > 1: psi = self.frame.prev.segment.psi gamma = self.frame.prev.gamma residuelist = [gamma**(a._exponent/Eplus)*self.frame.prev.FFbase_elt_to_FF((psi**(a._exponent/Eplus)*a.value()).residue()) for a in self.rep if int(a._exponent) % Eplus == 0] else: residuelist = [a.value().residue() for a in self.rep if int(a._exponent) == 0] if len(residuelist) > 0: return sum(residuelist) else: return self.frame.R(0) def reduce(self): if self.frame.is_first(): return self Eplus = self.frame.prev.segment.Eplus Fplus = self.frame.prev.Fplus reduced_elt = FrameElt(self.frame) # zero FrameElt for a in self.rep: if a._exponent >= Eplus * Fplus: b = FrameElt(self.frame) b.rep = [FrameEltTerm(b,a.value(),a._exponent - Eplus * Fplus)] b *= self.frame.prev.reduce_elt reduced_elt += b.reduce() elif a._exponent < 0: b = FrameElt(self.frame) b.rep = [FrameEltTerm(b,a.value(),a._exponent + Eplus * Fplus)] b *= (self.frame.prev.reduce_elt)**(-1) reduced_elt += b.reduce() else: summand = FrameElt(self.frame) summand.rep = [FrameEltTerm(reduced_elt,a.value().reduce(),a._exponent)] reduced_elt += summand return reduced_elt def sumcollapse(self): sumt = FrameElt(self.frame) sumc = FrameElt(self.frame) for i in range(len(self.rep)): sumc.rep = [self.rep[i]] sumt += sumc return sumt def is_reduced(self): return all([a.is_reduced() for a in self.rep]) def find_denominator(self): """ Find the lowest power (possibly negative) of the uniformizer in a FrameElt """ if self.frame.is_first(): return self.rep[0]._exponent else: return min([fet._coefficient.find_denominator() for fet in self.rep]) def polynomial(self,denominator=False): if denominator: piexp = self.find_denominator() if piexp < 0: return (self * FrameElt(self.frame,self.frame.Ox(self.frame.O.uniformizer()**(-piexp)))).polynomial(),-piexp else: return self.polynomial(),0 else: if self.frame.is_first(): return self.frame.O.uniformizer()**int(self.rep[0]._exponent)*self.rep[0]._unit else: return sum([self.frame.prev_frame().phi**int(a._exponent)*a._coefficient.polynomial() for a in self.rep]) def __neg__(self): if self.frame.is_first(): return FrameElt(self.frame,-self.polynomial()) else: neg = FrameElt(self.frame) neg.rep = [-r for r in self.rep] return neg def __radd__(self,l): return self.__add__(l) def __add__(self,r): if isinstance(r,int) and r == 0: return self if self.frame.phi != r.frame.phi: raise ValueError, "Cannot add FrameElts with different values of phi" if len(self.rep) == 0: return r sum = FrameElt(self.frame) if self.frame.prev == None: v = min(self.rep[0]._exponent,r.rep[0]._exponent) pi = self.frame.O.uniformizer() usum = self.rep[0]._unit * pi ** (self.rep[0]._exponent - v) usum = usum + r.rep[0]._unit * pi ** (r.rep[0]._exponent - v) sum.rep = [FrameEltTerm(sum,usum,v)] else: if self.rep == []: for k in range(len(r.rep)): sum.rep.append(FrameEltTerm(sum,r.rep[k].value(),r.rep[k]._exponent)) elif r.rep == []: for k in range(len(self.rep)): sum.rep.append(FrameEltTerm(sum,self.rep[k].value(),self.rep[k]._exponent)) else: i = 0 ; j = 0 while i < len(self.rep) and j < len(r.rep): if self.rep[i]._exponent < r.rep[j]._exponent: sum.rep.append(FrameEltTerm(sum,self.rep[i].value(),self.rep[i]._exponent)) i = i + 1 elif self.rep[i]._exponent > r.rep[j]._exponent: sum.rep.append(FrameEltTerm(sum,r.rep[j].value(),r.rep[j]._exponent)) j = j + 1 elif self.rep[i]._exponent == r.rep[j]._exponent: sum.rep.append(FrameEltTerm(sum,self.rep[i].value()+r.rep[j].value(),self.rep[i]._exponent)) i = i + 1; j = j + 1 if j < len(r.rep): for k in range(j,len(r.rep)): sum.rep.append(FrameEltTerm(sum,r.rep[k].value(),r.rep[k]._exponent)) elif i < len(self.rep): for k in range(i,len(self.rep)): sum.rep.append(FrameEltTerm(sum,self.rep[k].value(),self.rep[k]._exponent)) return sum def __rmul__(self,l): return self.__mul__(l) def __mul__(self,r): if isinstance(r,int) and r == 0: return self if self.frame.depth != r.frame.depth: raise ValueError, "Cannot multiply FrameElts with different frame depths" product = FrameElt(self.frame) if self.frame.prev == None: v = self.rep[0]._exponent + r.rep[0]._exponent uprod = self.rep[0]._unit uprod = uprod * r.rep[0]._unit product.rep = [FrameEltTerm(product,uprod,v)] else: for i in range(len(self.rep)): summand = FrameElt(self.frame) for j in range(len(r.rep)): summand.rep.append(FrameEltTerm(product,self.rep[i].value()*r.rep[j].value(),self.rep[i]._exponent+r.rep[j]._exponent)) product = product + summand return product def __pow__(self,n): if not self.is_single_term(): raise NotImplementedError, "Cannot take a power of a non-single term FrameElt" else: product = FrameElt(self.frame) product.rep = [self.rep[0]**n] return product def __div__(self,right): if not right.is_single_term(): raise NotImplementedError, "Cannot divide by a non-single term FrameElt" else: quotient = FrameElt(self.frame) quotient.rep = [a / right.rep[0] for a in self.rep] return quotient def __repr__(self): return repr(self.rep) def _phi_expansion(P,phi): """ Compute the phi-expansion of P """ if phi.degree() > P.degree(): return [P] expansion = [] q, r = P.quo_rem(phi) expansion.append(r) while q != 0: q, r = q.quo_rem(phi) expansion.append(r) return expansion class FrameEltTerm: def __init__(self,frelt,a,e): self.frameelt = frelt self._scalar_flag = (self.frameelt.frame.prev == None) self._exponent = e self._cached_valuation = None self._zero_flag = False if a in self.frameelt.frame.Ox or a in self.frameelt.frame.O: self._zero_flag = a.is_zero() if self._scalar_flag: a = self.frameelt.frame.O(a) if self._zero_flag: self._unit = a else: self._unit = a.unit_part() if a.valuation() > 0 and a.valuation(): self._exponent += a.valuation() else: self._coefficient = FrameElt(self.frameelt.frame.prev_frame(),a) else: self._coefficient = a a._zero_flag = False def valuation(self): if self._cached_valuation == None: if self.frameelt.frame.prev == None: self._cached_valuation = self._exponent else: self._cached_valuation = self.frameelt.frame.prev.segment.slope * self._exponent + self._coefficient.valuation() return self._cached_valuation def reduce(self): if self.frameelt.frame.prev == None: return Eplus = self.frameelt.frame.prev.segment.Eplus Fplus = self.frameelt.frame.prev.Fplus if self._exponent >= Eplus * Fplus: q,r = int(self._exponent).quo_rem(int(Eplus)) self._exponent = r self._coefficient *= (self.frameelt.frame.prev.segment.psi ** (q*Fplus)) self._coefficient.reduce() return def is_reduced(self): if self.frameelt.frame.prev == None: return True if self._exponent < self.frameelt.frame.prev.segment.Eplus: return self._coefficient.is_reduced() return False def is_scalar(self): return self._scalar_flag def is_zero(self): return self._zero_flag def is_single_term(self): if self._scalar_flag: return True else: return self._coefficient.is_single_term() def value(self): if self._scalar_flag: return self._unit else: return self._coefficient #def __add__(self,right): #def __mul__(self,right): #    We don't add or multiply on FrameEltTerms directly -- the parent FrameElt does this def __pow__(self,n): if not self.is_single_term(): raise NotImplementedError, "Cannot take a power of a non-single term FrameEltTerm" else: n = int(n) if self._scalar_flag: return FrameEltTerm(self.frameelt, self._unit**n, self._exponent*n) else: return FrameEltTerm(self.frameelt, self._coefficient**n, self._exponent*n) def __div__(self,right): if not right.is_single_term(): raise NotImplementedError, "Cannot divide by a non-single term FrameEltTerm" else: if self._scalar_flag: return FrameEltTerm(self.frameelt, self._unit/right._unit, self._exponent-right._exponent) else: return FrameEltTerm(self.frameelt, self._coefficient/right._coefficient, self._exponent-right._exponent) def __neg__(self): if self.frameelt.frame.is_first(): return FrameEltTerm(self.frameelt,-self._unit,self._exponent) else: return FrameEltTerm(self.frameelt,-self._coefficient,self._exponent) def __repr__(self): if self._zero_flag: return 'FET<0>' if self._scalar_flag: return 'FET' else: return 'FET<'+repr(self._exponent)+','+repr(self._coefficient)+'>'
diff --git a/sage/rings/polynomial/padics/factor/segment.py b/sage/rings/polynomial/padics/factor/segment.py
new file mode 100644
diff --git a/sage/rings/polynomial/padics/polynomial_padic_flat.py b/sage/rings/polynomial/padics/polynomial_padic_flat.py
diff --git a/setup.py b/setup.py