Ticket #12561: 12561.patch

new file sage/rings/polynomial/padics/factor/__init__.py

@@ +1 @@
+# Initialization file for p-adic polynomial factoring

new file sage/rings/polynomial/padics/factor/associatedfactor.py

@@ +1 @@
+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

new file sage/rings/polynomial/padics/factor/factoring.py

@@ +1 @@
+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.<x> = 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)

new file sage/rings/polynomial/padics/factor/frame.py

@@ +1 @@
+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)

new file sage/rings/polynomial/padics/factor/frameelt.py

@@ +1 @@
+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<pi^'+repr(self._exponent)+'*'+repr(self._unit)+'>'
+        else:
+            return 'FET<'+repr(self._exponent)+','+repr(self._coefficient)+'>'

new file sage/rings/polynomial/padics/factor/segment.py

@@ +1 @@
+from sage.rings.polynomial.padics.factor.associatedfactor import AssociatedFactor
+from sage.rings.arith import gcd
+
+class Segment:
+
+    def __init__(self,frame,slope,vert,length):
+        self.frame = frame
+        self.vertex = vert
+        self.slope = slope
+        self.length = length
+        self.Eplus = self.e() / gcd(self.e(),int(self.frame.E))
+        if self.frame.E * self.frame.F * self.Eplus < self.frame.Phi.degree():
+            self.psi = self.frame.find_psi(self.slope*self.Eplus)
+            self._associate_polynomial = self.associate_polynomial(cached=False)
+            self.factors = [AssociatedFactor(self,afact[0],afact[1]) for afact in list(self._associate_polynomial.factor())]
+        else:
+            # Polynomial is irreducible
+            return
+
+    def h(self):
+        if isinstance(self.slope,int):
+            return self.slope.numerator
+        else:
+            return self.slope.numerator()
+    def e(self):
+        if isinstance(self.slope,int):
+            return self.slope.denominator
+        else:
+            return self.slope.denominator()
+
+    def vphi(self):
+        return self.slope
+
+    def associate_polynomial(self,cached=True):
+        if cached:
+            return self._associate_polynomial()
+        else:
+            a = self.frame.elt_phi_expansion()
+            vertx = [v[0] for v in self.vertex]
+            chiex = [int((v-vertx[0]) // self.Eplus) for v in vertx]
+            chi = [a[vertx[i]] * self.psi**chiex[i] for i in range(len(vertx))]
+            psitilde = self.frame.find_psi(chi[0].valuation())
+            Ahat = [(c/psitilde).reduce() for c in chi]
+            if self.frame.prev == None:
+                Az = sum([(Ahat[i].residue()) * self.frame.Rz.gen() ** chiex[i] for i in range(len(Ahat))])
+            else:
+                Az = sum([(Ahat[i].residue()) * self.frame.prev.FFz.gen() ** chiex[i] for i in range(len(Ahat))])
+            return Az
+
+    def __repr__(self):
+        return 'Segment of length '+repr(self.length)+' and slope '+repr(self.slope)

sage/rings/polynomial/padics/polynomial_padic_flat.py

@@ -96 +96 @@
             return "0"
         return s[1:]
 
+    def quo_rem(self,right):
+        """
+        Returns the quotient and remainder of division by right
+
+        EXAMPLES:
+            sage: Kx.<x> = PolynomialRing(Zp(7))
+            sage: (x^3+7*x+1).quo_rem(x-7)
+            ((1 + O(7^20))*x^2 + (7 + O(7^21))*x + (7 + 7^2 + O(7^21)), (O(7^20))*x^3 + (O(7^21))*x^2 + (O(7^21))*x + (1 + 7^2 + 7^3 + O(7^20)))
+        """
+        return self._quo_rem_naive(right)
+
+    def _quo_rem_naive(self, right):
+        """
+        Naive quotient with remainder operating on padic polynomials as their coefficient lists
+        """
+        if right == 0:
+            raise ZeroDivisionError, "cannot divide by a polynomial indistinguishable from 0"
+        P = self.parent()
+        F = list(self); G = list(right);
+        fdeg = self.degree()
+        gdeg = right.degree()
+        Q = [0 for i in range(fdeg-gdeg+1)]
+        j=1
+        while fdeg >= gdeg:
+            a = F[-j]
+            if a!=0:
+                for i in range(fdeg-gdeg,fdeg+1):
+                    F[i] -= a * G[i-(fdeg-gdeg)]
+                Q[fdeg-gdeg] += a
+            j+=1; fdeg-=1;
+        return (P(Q), P(F))
+
     def factor(self, absprec = None):
         r"""
         Returns the factorization of this Polynomial_padic_flat.

setup.py

@@ -975 +975 @@
                      'sage.rings.padics',
                      'sage.rings.polynomial',
                      'sage.rings.polynomial.padics',
+                     'sage.rings.polynomial.padics.factor',
                      'sage.rings.semirings',
                      
                      'sage.tests',