misc.py

# Copyright 2013-2017 Edgar Costa 
# See LICENSE file for license details.

import itertools
from sage.all import Ideal, Partitions, PolynomialRing, QQ, Set, ZZ
from sage.all import binomial, vector 

# makes the vector v immutable and returns it
def immutable(v):
    v.set_immutable()
    return v
    
# returns vector obtained by decrementing entry i
# (i.e. differentiating in direction i)
def diff(v, i):
    return immutable(vector([v[j] - int(j == i) for j in range(len(v))]))    

# returns all tuples of degree d in n+1 variables    
def tuple_list(n,d):
    X=Set(range(n+1))
    
    if d==0:
        result=[immutable(vector(ZZ,n+1))]
    else:
        result=[];
        for variables in range(1, (n + 1) + 1):
            for subset in X.subsets(size=variables):
                for initialPartition in Partitions(d, length=variables):
                    for partition in Set(itertools.permutations(initialPartition)):
                        monomial=vector(ZZ,n+1);
                        for i in range(variables):
                            monomial[subset[i]]=partition[i]
                        result.append( immutable(monomial) )
    assert len(result)==binomial(d+n,n), "%s %s" %(len(result),binomial(d+n,n))
    return sorted(result);


# decrement first nonzero component of v
def tweak_step(v):
    for i,vi  in enumerate(v):
        if vi > 0:
            break
    return vector([v[j] - int(j == i) for j in range(len(v))])


# apply tweak_step r times
def tweak(v, r):
    return v if r == 0 else tweak(tweak_step(v), r - 1)

def poly_to_dic(fpoly):
    return dict([(immutable(vector(u)), c ) for (u, c) in fpoly.dict().iteritems()]);

def convertf(f):
    newf={}
    for i in f:
        newf[immutable(vector(i))]=f[i]
    return newf

def is_ND(f,ring=QQ):
    aMonomialOff = f.iterkeys().next()
    n = len(aMonomialOff)-1
    R=PolynomialRing(ring,n+1,"z")
    Z=R.gens()
    F = 0
    for i in f:
            mon = ring(f[i])
            for j in range(n+1):
                mon *= Z[j]**i[j]
            F+=mon
    Trivial = Ideal(Z)
    return Trivial == Ideal( F.derivative(Z[i]) for i in range(n+1) ).radical() and Trivial == Ideal( Z[i]*F.derivative(Z[i]) for i in range(n+1) ).radical()

def is_Smooth(f,ring=QQ):
    if f == {}:
        return False
    aMonomialOff = f.iterkeys().next()
    n = len(aMonomialOff)-1
    R=PolynomialRing(ring,n+1,"z")
    Z=R.gens()
    F = ring(0)
    for i in f:
            mon = ring(f[i])
            for j in range(n+1):
                mon *= Z[j]**i[j]
            F+=mon
    Trivial = Ideal(Z)
    return Trivial == Ideal( F.derivative(Z[i]) for i in range(n+1) ).radical()

def generateND(n, d, sparsity = 0, bound = 10, ring = QQ, verbose = False ):
    R = PolynomialRing(ring,n+1,"z")
    Z = R.gens()
    attempt = 0
    tl = tuple_list(n,d)
    while True:
        if sparsity == 0:
            f = dict([ ( u, ring(ZZ.random_element(-bound,bound) )) for u in tl ])
        else:
            f = {}
            for _ in range(n+sparsity+d):
                f[ tl[ ZZ.random_element( len(tl) )  ] ] = ring(ZZ.random_element(bound))
        F = 0
        for i in f:
            mon = f[i]
            for j in range(n+1):
                mon *= Z[j]**i[j]
            F+=mon
        Trivial = Ideal(Z)
        if Trivial == Ideal( F.derivative(Z[i]) for i in range(n+1) ).radical() and Trivial == Ideal( Z[i]*F.derivative(Z[i]) for i in range(n+1) ).radical():
            if verbose:
                print("Took %s attempts" % attempt)
            return f
        else:
            attempt += 1

def random_homogeneous(R, d, sparsity = 0):
    ring = R.base_ring();
    Z = R.gens()
    n = len(Z) - 1;
    tl = tuple_list(n,d)
    if sparsity == 0:
        f = dict([ ( u, ring.random_element() ) for u in tl ])
    else:
        f = {}
        for _ in range(n+sparsity+d):
            f[ tl[ ZZ.random_element( len(tl) )  ] ] = ring.random_element()
    F = 0
    for i in f:
        mon = f[i]
        for j in range(n+1):
            mon *= Z[j]**i[j]
        F+=mon
    return F;







def print_map(map):
    for k in sorted(map.keys()):
        print("%s: %s" %(k, map[k]))



def str_vector(v):
    out = "["
    for i in range(len(v)):
        out += str(v[i])
        if i != len(v) - 1:
            out+=" "
    out += "]"
    return out


def str_f(f):
    out = ""
    out += str(len(f));
    out += "\n"
    for (u,y) in f.iteritems():
        out += str_vector(u)
        out += "\n"
        out += str(y)
        out += "\n"
    return out

def str_mapmatrix(f):
    out = ""
    out += str(len(f));
    out += "\n"
    for (u,y) in f.iteritems():
        out += str_vector(u)
        out += "\n"
        out += "["+y.str()+"\n]"
        out += "\n"
    return out


def str_matrix(M):
    out = ""
    out += "["
    for row in M.rows():
        out += str_vector(row)+"\n"

    out += "]"
    return out;

def poly_to_vector(f):
    d = f.degree();
    n = len(f.variables());
    tuples_list = sorted(tuple_list(n-1,d));
    tuples_dict = dict(zip(tuples_list, range(len(tuples_list))))
    result =[0]*len(tuples_list);
    R = PolynomialRing(f.base_ring(), n, f.variables());
    for mon, coef in R(f).dict().iteritems():
        result[ tuples_dict[immutable(vector(mon))] ] = coef;
        
    return result;


    

precision.py

# Copyright 2013-2017 Edgar Costa 
# See LICENSE file for license details.

# -*- coding: utf-8 -*-
from misc import tuple_list, poly_to_vector, str_vector
from sage.all import Infinity, Integer, Matrix, O, PolynomialRing, Qp, ZZ
from sage.all import binomial, ceil, floor, log, vector 

# -*- coding: utf-8 -*-

def bound_f(m, n, p, bound = 20, A0 = []):
    """
    an equivalent of the algorithm 3.4.10
    returns A with at least m+1 elements
    where A[j] >= f(p * j + i) for 1 <= i <= p
    note: A[j] >= f( p * (j + 1))
    """

    def f0(m, n, p):
        """
        returns f_0 (m) as defined in 3.4.10@AKR
        """
        return sum(floor( log( max(m - i, 1) , p)) for i in range(1,n + 1))

    def g(m, i, p):
        """
        returns g(m,i) as defined in 3.4.8@AKR
        """
        bin = binomial(-m,i)
        return val(bin,p)

    def val(m, p):
        """
        returns \nu_p(m)
        """
        if m == 0:
            return 0
        if m % p == 0:
            return 1 + val(ZZ(m)/p,p)
        else:
            return 0


    # A0[j] is an upper bound for f(p * (j + 1)) = f(p * j + i) for 1 <= i <= p
    # A0[j] >= f(p * j + i) for 1 <= i <= p
    A = A0[:]
    for i in range(len(A0),m+1):
        # we want to be sure that we have A[m], hence if no initial bound is given we use f0 to get such bound
        # so we set A[i] = f0( p*i + 1) >= f(p*i + 1)
        A.append(f0(p * i + 1, n , p))
    for iteration in range(bound):
        Ai = A[:]
        j = 1
        while j < len(A):
            # why the loop stops
            # for j big enough such that p | j
            # we will have l1 = 1 and A[j//p] = f0( p * (j//p) + 1)
            # hence n - 1 + A[ j//p ] >= n - 1 + n*log(j + 1 -n)
            # and m_log = n + n * log( j+ l) -l
            # thus n - 1 + A[ j//p ] - m_log >= n * log( (j + 1 -n)/(j + l) ) + l - 1 >= 0 for l>=2 and j >> l
            # same argument works to show that m  >= f0_jl - l - g_jl for l = 1, since g_jl > 0 
            
            # each iteration tries to improve the upper bound on f(p * (j+1)) = f(p * j + i )  1 <= i <= p
    
            m = n - 1 + A[j//p] 
    
            # for l >= l1 the fucntion n*log( ((j + 1) + l)*p) - l is decreasing 
            # and is an upper bound for f( ((j + 1) + l)*p ) - l
            l1 = max(n - j, 1)
    
            l = 1
            m_log = n * log( (j + 1 + l) * p, p) - l
            
            while not( l >= l1 and m_log <= m ):
                # every loop m_log decreases (for l big enough) and m never decreases
    
                if not f0( p * (j + l) + 1, n, p) - l -  g(j+1,l, p) <= m:
                # maybe A[j+l] is a shaper upper bound for f(p*(j + l + 1)) = f(p * (j + l) + 1)
                    if j + l >= len(A):
                        for i in range(len(A),j+l+1):
                            A.append(f0(p * i + 1, n , p))
    
                    m = max(m, A[j + l] - l - g(j + 1, l, p) )
                l += 1
                m_log = n * log( (j + 1 + l) * p , p) - l
            
            A[j] = min(A[j], m)
            if A[j] < A[j-1]:
                for i in range(j-1,0,-1):
                    if A[i] > A[i+1]:
                        A[i] = A[i+1]
                    else:
                        break
            j += 1
        if Ai == A:
            # we got to a fixed point
            break
    return A

def precision_s(r, m, n, p, bound = 20, A = []):
    """
    algorithm 3.5.0
    returns s such that
        m - 1 + j - f(p * (m +j)) >= r \forall j >= s - m + 1
    """
    s = r
    j = s - m + 1
    while not( j > 0 and n * log(p * (m + j - 1), p) <= m - 1 + j - r):
        # n * log(p * (m + j - 1), p) is an upper bound for f( p * (m + j)) = f( p * (m + j - 1) + 1)
        if m - 1 + j >= len(A):
            A = bound_f( m + j - 1, n, p, bound, A)
        if m - 1 + j - r < A[ m - 1 + j]: # >= f(p * (m + j)
            s += 1
            j = s - m + 1
        else:
            j += 1
    return (s, A)



def Hodge_numbers(n,d):
    ZZpoly = PolynomialRing(ZZ, "t")
    p2 = ZZpoly( [ Integer(1) ] * ( d -1 ) )
    Hilbert_J = ( p2**( n + 1) ).list()
    return [ Hilbert_J[ d*m - (n+1) ] if  d*m - (n+1) >= 0 and  d*m - (n+1) < len(Hilbert_J) else 0 for m in range(1, n + 1)  ]

def RH_bounds(n,d, p = None):
    n = ZZ(n)
    d = ZZ(d)
    degree = sum(Hodge_numbers(n,d));
    if p is not None:
        return [ 2*binomial(degree, k) * p**((n-1)/2 * (degree-k)) for k in range(degree + 1) ]
    else:
        return [ ((n-1)/2 * (degree-k),2* binomial(degree,k)) for k in range(degree+1) ]


def improved_bounds(n,d, p = None):
    n = ZZ(n)
    d = ZZ(d)
    degree = sum(Hodge_numbers(n,d));
    if p is not None:
        return [2 * degree/(degree-k) * p**((n-1)/2 * (degree-k)) if k != degree else 1  for k in range(degree + 1) ];
    else:
        return [ ((n-1)/2 * (degree-k), 2*degree/(degree-k) if (k != degree and k!=0) else 1 ) for k in range(degree+1) ];


    
# extremely badly written...
def Mazur_Inequalities(n,d):
    """
    returns a tuple
    (basis_precision, MI, MI_prec, MI_prec_without_symmetry)
    basis_precision - to compute the characteristic polynomial coefficients  with r-digits of precision we need to compute F(x^\beta \Omega /f^i) modulo p^(n-m+r-basis_precision[i])
    MI - valuation of characteristic polynomial coefficients
    MI_prec - same as basis_precision but now a tuple for each coefficient of the characteristic polynomial
    MI_prec_without_symmetry - same as above, but whitout using the symmetry given by RH
    """

    H = Hodge_numbers(n,d)
    degree = sum(H)
    val = []
    
    for k in range(degree + 1):
        val.append({})
        for part in tuple_list(n-1,k):
            
            if all( H[i] >= part[i] for i in range(n) ):
                v = sum(i*xi for i,xi in enumerate(part))
                
                if v in val[-1]:
                    val[-1][v].append(part)
                else:
                    val[-1][v]=[part]
    val.reverse()

    # prec to compute each coefficient of det(T - F) with r-digits of precision we need to compute F(x^\beta \Omega /f^(m+1)) modulo p^(n-m+r-basis_precision[m])
    prec = []
    for x in val:
        local = [+Infinity] * n
        for i,ti in enumerate(sorted(x.keys())):
            sum_of_xti = sum(x[ti])
            # local[j] = min(i,local[j]) if for some k in  x[ti] k[j] > 0 else local[j]
            # if term j doesnt contribute to i then we can always use less precision to compute this coefficient
            local = [ min(i,local[j]) if sum_of_xti[j]>0 else local[j] for j in range(n)]
        # since the beginning we had the partitions inverted
        local.reverse()
        prec.append(local)
    MI_prec_without_symmetry = prec[:];  
    for k in range(degree+1):
        prec[k] = [ max(prec[k][j], prec[degree - k][j]) for j in range(n) ]
    
        
    basis_precision = [ min( prec[k][j] for k in range(degree+1) ) for j in range(n) ]
    MI = [ min( val[k].keys()) for k in range(degree+1) ]
    MI_prec = [ prec[k] for k in range(degree+1) ]
    return basis_precision, MI, MI_prec, MI_prec_without_symmetry                

        
def precision_for_zeta(n, d, p , bound = 20, A=[]):
    """
        returns an array A of tuples where
        A[r-1] = (r, s, r_vector, N_vector, ...)
        r - the characteristic polynomial can be computed with at least r-digits of precision 
        s - we need to perform the calculations modulo p^s
        r_vector - (r1, r2, ... ,rn) F(x^\beta \Omega /f^i) will be computed with modulo p^(n-i+ri), also r = max(ri)
        N_vector - (N1, ..., Nn) to obtain F(x^\beta \Omega /f^i) modulo p^(n-i+ri) we only add the first with Ni terms of the series
        relative_precision - let det(T - F) = \sum ai T^i. If zeta_precision[i]==0 the ai term can be computed, otherwise it can be compute with zeta_precision[i] digits of precision
        absolute_precision - same as above but now for absolute precision, if 0 means that you can lift the coeff to ZZ and use the bound to decide the sign
        relative_precision_without_symmetry - by working modulo p^s and using the N_vector the characteristic polynomial will have initially these digits of precision
        absolute_precision_without_symmetry - same as above but for absolute precision  
    """

    basis_precision, MI , MI_prec, MI_prec_without_symmetry = Mazur_Inequalities(n,d)
    RH = improved_bounds(n,d); #RH_bounds(n,d)
    degree = len(MI)-1

    # precision needed to compute each coefficient
    relative_precision_needed = [ ceil(RH[k][0] - MI[k] + log(RH[k][1],p)) for k in range(degree+1)]
    # applying RH
    # whenever we use the improved bounds we cannot use the symmetry
    # relative_precision_needed = [ min(relative_precision_needed[k],relative_precision_needed[degree-k]) for k in range(degree+1) ]
    max_r = max(relative_precision_needed)
    s, working_precision, N = precision_tuple(max_r,  n, p, bound, A)    
    result = []

    for r in range(1, max_r + 1):
        relative_precision = [ 0 if relative_precision_needed[k] <= r else r for k in range(degree+1) ]
        absolute_precision = [ 0 if relative_precision[k] == 0 else relative_precision[k] + MI[k] for k in range(degree+1)]
        
        # compute the relative precision need from each column
        precision_needed_basis = [ vector( [ max(0, min(relative_precision_needed[k],r) - y) for y in MI_prec[k]]) for k in range(1,degree) ]
        
        # get to which relative precision we will compute column
        r_vector = [max( pn[m] for pn in precision_needed_basis) for m in range(n) ]
        s_vector = [ s[r_vector[m]][m] for m in range(n) ]
        
        wp_vector = [ working_precision[r_vector[m]][m] for m in range(n) ]
        N_vector = [ max(s_vector[m] - m,0) for m in range(n) ]
        #by working modulo p^s and using the N_vector the characteristic polynomial will have initially these digits of precision
        relative_precision_without_symmetry = [ min( [r_vector[i] +  MI_prec_without_symmetry[k][i] for i in range(n)] ) for k in range(degree)]
        relative_precision_without_symmetry.append(+Infinity)
        # same for absolute precision
        absolute_precision_without_symmetry = [ relative_precision_without_symmetry[k] + MI[k] for k in range(degree) ]
        absolute_precision_without_symmetry.append(+Infinity)
        result.append((r,max(wp_vector), r_vector, N_vector,relative_precision, absolute_precision, relative_precision_without_symmetry, absolute_precision_without_symmetry))

        
    return result
    
    
def precision_tuple(max_r,  n, p, bound = 20, A = []):
    """
        returns a tuple (s, working_precision, N)
        where s[r][m], working_precision[r][m] and N[r][m]
        are the values required to compute F(x^\beta \Omega /f^m) modulo p^(n-m+r)
    """
    def factorial_val(m):
        return (m-sum(Qp(p)(m).expansion()))/ZZ(p-1)

    s = [[0]*n]
    working_precision = [ [0]*n ]
    N = [ [0]*n ]
    A0 = A
    for r in range(1, max_r + 1):
        s_r = []
        wp_r = []
        N_r = []
        for m in range(1,n+1):
            s0, A0 = precision_s(r, m, n, p, bound, A0)
            s_r.append(s0)
            N_r.append(s0 - m + 1)
            # p*(m + N - 1) = p*(m + s0 - m + 1 - 1) = p*s0
            wp_r.append(int(r + factorial_val(p * s0 - 1) - m + 1))
        s.append(s_r)
        N.append(N_r)
        working_precision.append(wp_r)

    return (s, working_precision, N)


def compute_zeta(n, d, F, p, pfz):
    __, __, rvector, __, __, __, __, z_apwrh = pfz;
    n = ZZ(n)
    
    Hodge = Hodge_numbers(n,d)
    RH_bound = RH_bounds(n,d,p)
    degree = sum(Hodge);
    Orow = []
    for i, hi in enumerate(Hodge):
        Orow += [O(p**(n-i+rvector[i])) for _ in range(hi) ]
    
    assert len(Orow) == degree
    Fcorrect = F + Matrix([Orow for _ in range(degree)])
    FZZ = Matrix(ZZ,Fcorrect)
    ZZX = PolynomialRing(ZZ, 1, "T")
    T = ZZX.gens()
    R = FZZ.charpoly(T)
    
    Rcorrect = []
    
    for i,x in enumerate(R.list()):
        if z_apwrh[i] == +Infinity:
            Rcorrect.append(x)
        else:   
            x = x + O(p**(z_apwrh[i]))
            if 2*RH_bound[i] <= p**z_apwrh[i]:
                if x.lift() > RH_bound[i]:
                    Rcorrect.append(-(-x).lift())
                elif x.lift() <= RH_bound[i]:
                    Rcorrect.append(x.lift())
            else:
                Rcorrect.append(x)
    print(Rc)
    if Rcorrect[1] == - Rcorrect[degree-1] * p**((degree-2*1)*(n-1)/2):
        signal = -1
    else:
        signal = 1
        assert Rcorrect[1] == Rcorrect[degree-1] * p**((degree-2*1)*(n-1)/2)

    print(Rcorrect)
    for i,x in enumerate(Rcorrect):
        if hasattr(x,"lift"):
            if hasattr(Rcorrect[degree-1],"lift"):
                print("Cannot pinpoint the "+str(i)+"th coefficient.")
            else:
                Rcorrect[i] = signal * p**((degree-2*i)*(n-1)/2) * Rcorrect[degree -i]
        else:
            assert Rcorrect[i] == signal * p**((degree-2*i)*(n-1)/2) * Rcorrect[degree -i]
    
    Q = Rcorrect[:]
    Q.reverse()
    return Q
    
    
    
def generate_input(f, p):
    d = f.degree();
    n = len(f.variables()) - 1;
    if n == 3 and  d == 4 and p > 42:
        pfz = (2, 4, [1, 2, 2], [3, 3, 2], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0], [1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,2, 2, 2, 2, 2, +Infinity], [22, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12,11, 10, 9, 8, 7, 6, 5, 4, 3, 2, +Infinity])
    if n == 2 and d == 4 and p > 16:
        pfz = (2, 3, [1, 2], [2, 2], [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0], [1, 1, 1, 2, 2, 2, +Infinity], [4, 3, 2, 2, 2, 2, +Infinity])
    else:
        pfz = precision_for_zeta(n, d, p, bound = 5 if p == 3 else 10)[-1];
    input = "";
    input = input + "%s\n" % (p,)
    input = input + "%s\n" % (pfz[1],)
    input = input + "%s\n" % (n,)
    input = input + "%s\n" % (d,)
    input = input + "%s\n" % (str_vector(poly_to_vector(f)),)
    input = input + "%s\n" % (str_vector(pfz[3]),)
    abp = pfz[-1];
    abp[-1] = abp[-2];
    input = input + "%s\n" % (str_vector(abp),)
    return input

# deprecated
def produce_input(f, p):
    return generate_input(f, p)