class pAdicLocalAnalyticFunction(object):

    def __init__(self,p):
        assert(is_prime(p))
        self._p=p
        self._residue_classes=set([])
        self._center=[0..p-1]
        self._power_series=[0..p-1]

    def __repr__(self):
        return "p-adic local analytic function with p=%s"%(self._p)

    def define(self,a,zeta,f):
        self._residue_classes.add(a)
        self._center[a]=zeta
        self._power_series[a]=f

    def eval(self,z):
        a=mod(z,self._p)
        assert a in self._residue_classes,"function not defined on residue class" 
        f=self._power_series[a]
        K=f.base_ring()
        return f(K(z)-self._center[a])

    def newton_pol(self,a):
        f=self._power_series[a]
        vc=[c.valuation() for c in f.padded_list()]
        return newton_pol(vc),vc

    def zeroes(self):
        Z=[]
        for a in self._residue_classes:
            Za=padic_zeroes(self._power_series[a])
            Z+=[z+self._center[a] for z in Za]
        return Z


#--------------------------------------------------------------------------


def padic_zeroes(f):

    R=f.parent()
    K=R.base_ring()
    p=K.prime()
    t=R.gen()
    Z=[]
    
    g=(f(p*t)).truncate(15)
    FL=g.factor()
    for F in FL:
        F=F[0]
        if F.degree()==1 and F[0].valuation()>=F[1].valuation():
            z=-p*F[0]/F[1]
            z=newton_approx(f,z)
            Z+=[z]
    return Z


def padic_common_zeroes(F,prec):

    n=len(F)
    Z=F[0].zeroes()
    for i in [1..n-1]:
        Z1=[]
        for z in Z:
            if F[i].eval(z).valuation()>=prec:
                Z1+=[z]
        Z=Z1
    return Z

                
        

#-----------------------------------------------------------------------------


def subs_pol(P,F):

    p=F[0]._p
    n=len(F)

    m=P.parent().ngens()
    print p,n,m

    residue_classes=set([0..p-1])
    for i in range(n):
    	residue_classes.intersection_update(F[i]._residue_classes)
    
    g=pAdicLocalAnalyticFunction(p)
    for a in residue_classes:
        zeta=F[0]._center[a]
        f=[F[i]._power_series[a] for i in range(n)]
        for i in [n..m-1]: f.append(1)
        ga=P(f)
        g.define(a,zeta,ga)
        
    return g

    
#------------------------------------------------------------------------------



def next_kink(a,i):

    N=len(a)-1
    if a[i]=='infty':
        j=i+1
        while j<N and a[j]=='infty':
            j=j+1
        return(j,'infty')

    mu=0
    k=i
    j=i
    while j<N:
        j=j+1
        if a[j]!='infty':
            l=(a[i]-a[j])/(j-i)
            if l>=mu:
                mu=l
                k=j
    return(k,mu)


def newton_pol(a):

    N=len(a)-1
    i=0
    np=[]
    while i<N:
        k,l=next_kink(a,i)
        if k==i:
            i=N
        else:    
            np.append([k-i,l])
            i=k
    return np

    


#---------------------------------------------------------------

# Compute the zeroe of a function by Newton approximation


def newton_approx(f,x0):

    f1=f.derivative()
    x=x0
    for k in range(0,20):
        x=x-f(x)/f1(x)
    return(x)