Open Copy - A (Programmed) Course in Arithmetic.sagews

Ring = PolynomialRing(QQ,"X",3)
X = Ring.gens()
INF = maxima('infinity').sage()

def epsilon(n): #Recibe un primo impar, devuelve 0 si -1Rn y 1 si -1Nn.
    eps = ((n-1)/2) % 2
    return eps

def omega(n): #Recibe un primo impar, devuelve 0 si nR2 y 1 si nN2.
    ome = ((n^2 - 1)/8) % 2
    return ome

def Legendre(n,p): #Recibe un entero n y un primo p. Devuelve 1 si nRp y -1 si nNp.
    n = (n % p)
    if n == 0:
        Lgn = 0
    elif n == 1:
        Lgn = 1
    elif n == (p-1):
        Lgn = (-1)^epsilon(p)
    elif n == 2:
        Lgn = (-1)^omega(p)
    elif n in Primes():
        Lgn = Legendre(p,n)*(-1)^(epsilon(n)*epsilon(p))
    else:
        fact = prime_factors(squarefree_part(n))
        Lgn = 1
        for L in fact:
            Lgn *= Legendre(L,p)
#    else:
#        Lgn = power_mod(n,(p-1)//2,p)
#        if Lgn == (p-1):
#            Lgn = -1
    return Lgn

def valuation(n,p): #Recibe un número racional y un primo. Devuelve el exponente de p en la factorización de n.
    if n in Integers():
        if n == 0:
            v = INF
        elif abs(n) == 1:
            v = 0
        else:
            v = 0
            while (n % p == 0):
                v+= 1
                n /= p
    else:
        v = valuation(n.numerator(),p) - valuation(n.denominator(),p)
    return v

def Hilbert(a,b,p): #Recibe dos racionales y un primo. Devuelve 1 si la forma z^2 - ax^2 - by^2 tiene algún cero no trivial en Q_p; si no, devuelve -1
    alpha = valuation(a,p)
    beta = valuation(b,p)
    u = a//(p^alpha)
    v = b//(p^beta)
    if p == 2:
        if ((epsilon(u)*epsilon(v) + alpha*omega(v) + beta*omega(u)) % 2) == 0:
            H = 1
        else:
            H = -1
    else:
        U = (u % p)
        V = (v % p)
        H = ((-1)^(alpha*beta*epsilon(p)))*(Legendre(U,p)^beta)*(Legendre(V,p)^alpha)
    return H

def invariants(P): #Recibe una forma cuadrática ternaria. Devuelve el discriminante y la signatura de esta forma.
    c = P.coefficients()
    c = [1,squarefree_part(c[1]/c[0]),squarefree_part(c[2]/c[0])]
    d = 1
    r = 0
    for i in c:
        d *= i
        if i > 0:
            r += 1
    s = 3 - r
    return [d,r,s]

def has_rational_zero(P): #Recibe una forma cuadrática ternaria. Devuelve True si tiene algún cero racional no trivial y Fase en caso contrario.
    c = P.coefficients()
    c = [1,squarefree_part(c[1]/c[0]),squarefree_part(c[2]/c[0])]
    [d,r,s] = invariants(P)
    if (r*s) != 0:
        has = True
        div = prime_factors(2*d)
        for p in div:
            if Hilbert(-c[1],-c[2],p) == -1:
                has = False
    else:
        has = False
    return has

def TonelliShanks(n,p): #Recibe un entero y un primo de manera que nRp. Devuelve el menor R positivo tal que R^2 = n (mod p).
    n = (n % p)
    if n == 0:
        R = 0
    elif n == 1:
        R = 1
    elif ((p+1) % 4) == 0:
        R = power_mod(n,(p+1)//4,p)
    else:
        S = valuation(p-1,2)
        Q = (p-1)//(2^S)
        z = 1
        while Legendre(z,p) == 1:
            z += 1
        c = power_mod(z,Q,p)
        R = power_mod(n,(Q+1)//2,p)
        t = power_mod(n,Q,p)
        M = S
        while t != 1:
            temp_t = t
            i = 0
            while temp_t != 1:
                temp_t = power_mod(temp_t,2,p)
                i += 1
            if M > (i+1):
                for j in range(M-i-1):
                    c = power_mod(c,2,p)
            b = c
            R = (R*b) % p
            c = power_mod(b,2,p)
            t = (t*c) % p
            M = i
    return R

def extended_crt(N,P): #Recibe una L-tupla de enteros y una L-tupla de primos. Resuelve el sistema x = N[i] (mod P[i]).
    L = len(N)
    if L == 1:
        X = N[0]
    elif L == 2:
        X = crt(N[0],N[1],P[0],P[1])
    else:
        X = N[0]
        B = P[0]
        for i in range(L-1):
            X = crt(X,N[i+1],B,P[i+1])
            B *= P[i+1]
    return X

def mod_sqrt(A,B): #Recibe un entero y un entero libre de cuadrados. Devuelve el menor entero positivo R tal que R^2 = A (mod B).
    B = abs(B)
    P = prime_factors(B)
    N = []
    for p in P:
        N.append(TonelliShanks(A,p))
    R = extended_crt(N,P)
    if 2*R > B:
        R = B - R
    return R

def find_solution(A,B): #Recibe dos racionales. Encuentra un cero no trivial de X^2 - AY^2 - BZ^2.
    sqfA = squarefree_part(A)
    sqfB = squarefree_part(B)
    if (A != sqfA) or (B != sqfB):
        factor = [1,sqrt(sqfA/A),sqrt(sqfB/B)]
        temp_sol = find_solution(sqfA,sqfB)
        sol = [factor[i]*temp_sol[i] for i in range(3)]
    elif A == 1:
            sol = [1,1,0]
    elif B == 1:
            sol = [1,0,1]
    elif abs(B) < abs(A):
        temp_sol = find_solution(B,A)
        sol = [temp_sol[0],temp_sol[2],temp_sol[1]]
    else:
        T = mod_sqrt(A,B)
        B1 = (T^2 - A)//B
        B2 = squarefree_part(B1)
        U = sqrt(B1/B2)
        temp_sol = find_solution(A,B2)
        sol = [T*temp_sol[0] + A*temp_sol[1],temp_sol[0] + T*temp_sol[1], B2*U*temp_sol[2]]
    return sol

def ternary_solve(P,check=False): #Recibe una forma cuadrática ternaria. Avisa de que no tiene ceros racionales si es el caso y si no, encuentra uno.
    if has_rational_zero(P):
        c = P.coefficients()
        c= [1,c[1]/c[0],c[2]/c[0]]
        sol = find_solution(-c[1],-c[2])
        if check:
            checksum = 0
            for i in range(3):
                checksum += (c[i]*sol[i]^2)
            if checksum == 0:
                print "Correct"
            else:
                print "Something went wrong"
    else:
        sol = "This quadratic form has no nontrivial rational zeros"
    return sol

#P = X[0]^2 + X[1]^2 - X[2]^2
#P = 3*X[0]^2 + (1/2)*X[1]^2 - (1/3)*X[2]^2
#P = 13*X[0]^2 + (23/2)*X[1]^2 - (9/59)*X[2]^2
#P = 21*X[0]^2 - (91/17)*X[1]^2 + (14/109)*X[2]^2
#P = 653896783*X[0]^2 - X[1]^2 + (7093/681677)*X[2]^2
#P = 37*X[0]^2 + 2017*X[1]^2 - X[2]^2
P = X[0]^2 -13*X[1]^2 -61*X[2]^2
ternary_solve(P)

[69, 18, 3]