# Listings included in the paper 
# "A bootstrap for the number of F_{q^r}-rational points on a curve over F_q"
# by S. Molina, N. Sayols and S. Xambó. arXiv ...

from PyM import *

# Function XN(q,X,k)
# An improved procedure to get the number X_k of points of a curve C/F_q
# in arbitrary extensions F_{q^k} taking X=[X_1,...,X_g] as input, g the genus of C.
# It outputs [X_1,...,X_g,...,X_k]

# We assume that Q_ is an implementation of the rational numbers and that
# x >> Q_ is the inclusion of the integer x in Q_
def XN(q,X,k):
    g = len(X)
    if k<=g: return X[:k]
    X = [0]+X           # trick so that X[j] refers to F_{q^j}
    X = [x>>Q_ for x in X]
    S = [q**(j)+1-X[j] for j in range(1,g+1)] # Newton sums
    S = [0]+S             # similar trick
    # Computation of c1,...,cg; set c0=1 
    c = [1>>Q_]           
    for j in range(1,g+1):
        cj = S[j]
        for i in range(1,j):
            cj += c[i]*S[j-i]
        c += [-cj/j]
    # Add c_{g+i}, for i=1,...,g
    for i in range(1,g+1):
        c += [q**i*c[g-i]]
    # Find Sj for j = g+1,...,k
    for j in range(g+1,k+1):
        if j>2*g:
            Sj=0
        else: 
            Sj = j*c[j]
        for i in range(1,j):
            if i>2*g: break
            Sj += c[i]*S[j-i]
        S += [-Sj]
    # Find X[i] for i = g+1,...,k
    for i in range(g+1,k+1):
        X += [q**i+1-S[i]]
    return create_vector(X[1:])
    
# This function implements the Deuring's algorithm to list the possible cardinals #E(F_q)
# of the elliptic curves E/F_q. Our main reference here has been Waterhouse-1969. 
# We have split the computation in two parts: the function Deuring_offsets(q) 
# (which computes the list of integers t in the segment [-m,m], m=\floor{2\sqrt{q}},
# such that #E(F_q)=q+1-t for some E), and Deuring_set(q), that outputs the list in question.    

def Deuring_offsets(q):
    P = prime_factors(q)  # prime_factors(12) => [2, 2, 3]
    p = P[0]; n = len(P)
    m = int(2*sqrt(q))
    D = [t for t in range(-m,m+1) if gcd(p,t)==1]
    if n%2==0: 
        r = p**(n//2)
        D += [-2*r,2*r] 
        if p%3 != 1: 
            D += [-r,r]
    if n%2 and (p==2 or p==3):
        r = p**((n+1)//2)
        D += [-r,r]
    if n%2 or (n%2==0 and p%4!=1):
        D += [0]
    return sorted([t for t in D])

def Deuring_set(q):
    D =Deuring_offsets(q)
    return [t+q+1 for t in D]

def Serre(q, g=1):
    D = ifactor(q) # ifactor(12) => {2: 2, 3: 1}
    if len(D)>1:
        return 'Serre: {} is not a prime power'.format(q)
    p = list(D)[0]
    e = D[p] # q = p^e
    m = int(2*sqrt(q)) # g
m is the HWS bound
    if g==1:
        if e%2 and e>=3 and m%p==0: return q+m
        else: return q+m+1
    if g==2:
        if q==4: return 10
        if q==9: return 20
        if e%2==0: return q+1+2*m
        def special(s):
            if m%p==0 or is_square(s-1) or \
                         is_square(4*s-3) or \
                         is_square(4*s-7): 
                return True
            else: return False
        if special(q):
            if 2*sqrt(q)-m > (sqrt(5)-1)/2: return q+2*m
            else: return q+2*m-1
        return q+1+2*m
    if g==3:
        P = [2,3,4,5,7,8,9]
        T={2:7,3:10,4:14,5:16,7:20,8:24,9:28}
        if P.count(q): return T[q]
    return "Serre: I do not know the value of N({},{})".format(q,g)

show('Deuring examples')
Q = [2,3,4,5,7, 8,9,11,13,16, 17,19]
for q in Q:
    m = int(2*sqrt(q))
    D = Deuring_set(q)
    show([q,2*m+1-len(D)],[q+1-m,q+1+m], D)

show('Serre examples')
for q in Q:
    print("q =",q,": ",Serre(q), Serre(q,2), Serre(q,3))