def make_Li(p,S,n,prec_p):

    global Lip,zetap,logp

    K=Qp(p,prec_p)
    R.<v>=PowerSeriesRing(K,default_prec=p*prec_p)
    T.<t>=PowerSeriesRing(K,default_prec=100)
    
    if S=='all':
        A=[2..p-1]
    else:
        A=[]
        for  s in S:
            A.append(K(s).residue())
            A.append(K(1-s).residue())

    P.<x>=K[]
    f=x^(p-1)-1
    l=f.roots(multiplicities=False)
    zeta=[0..p-1]
    for z in l:
        a=z.residue().lift()
        zeta[a]=z

    g=[0..n]
    g[0]=v-1-(v-1)^p/(v^p-(v-1)^p)
    for k in [1..n]:
        g[k]=(-g[k-1].shift(-1)/(1-v)).integral()

    c=[0..n]
    F=[0..n]
    F[0]=[0..p-1]
    Li=[0..n]

    for k in [1..n]:
        c[k]=[0..p-1]
        F[k]=[0..p-1]
        Li[k]=pAdicLocalAnalyticFunction(p)
            
    for a in A:
        F[0][a]=1/(1-zeta[a]-t)
        c[1][a]=g[1].polynomial()(1/(1-zeta[a]))*p/(p-1)
        F[1][a]=c[1][a]+F[0][a].integral()
        Li[1].define(a,zeta[a],F[1][a])
        
    for k in [2..n]:
        for a in A:
            c[k][a]=g[k].polynomial()(1/(1-zeta[a]))*p^k/(p^k-1)
            F[k][a]=c[k][a]+(F[k-1][a]/(zeta[a]+t)).integral()
            Li[k].define(a,zeta[a],F[k][a])
    
    Li[0]=pAdicLocalAnalyticFunction(p)
    for a in A:
        Li[0].define(a,zeta[a],-F[1][p+1-a].polynomial()(1-zeta[a]-zeta[p+1-a]-t))
   
    Lip= lambda k,z: Li[k].eval(z)

    
    zetap= lambda k: 2^(k-1)/(1-2^(k-1))*Lip(k,-1)

    logp = lambda z: Lip(0,z)

    return Li