sagemath / sage

Main repository of SageMath
https://www.sagemath.org
Other
1.45k stars 481 forks source link

Solving norm equation in a quadratic extension in Q_q (2|q) #24933

Open 742ccff4-512f-43e4-bdd0-5df91c799a52 opened 6 years ago

742ccff4-512f-43e4-bdd0-5df91c799a52 commented 6 years ago

Let b be a non square in Q_q (q=2n) and a an element of Q_q.

The problem is to find x,y in Q_q such that x2- b y22 =a i.e. solving the norm equation N(x + sqrt(b) y)=a. I have used it in order to find isotropic vectors for quadratic forms with coefficients in Q_q. But it is also interesting in order to compute the Hasse-Witt invariant of quadratic forms. I have a draft of code that I need to adapt to the coding standards of sage.

CC: @xcaruso @roed314

Component: padics

Keywords: padicIMA

Issue created by migration from https://trac.sagemath.org/ticket/24933

742ccff4-512f-43e4-bdd0-5df91c799a52 commented 6 years ago
comment:2

Here is my draft of code :


"""
Mod 4 Teichmuller lift
- input : x
- output : a Teichmuller lift of x mod 4
"""
def teichmuller_mod_four(x):
    R=parent(x)
    x0=x.residue()
    r=x0.parent().degree()
    return R(x0^(2^(r-2)))^4

"""
Mod 4 norm equation solver
- input : a, b
- output : alpha, beta such that alpha^2 - b*beta^2=a mod 4
Idea of the algorithm :
We reprents elements of Z_(2^r)/4 as Witt vectors i.e as couples (x_0, x_1) of elements of
F_(2^r) with laws given by :
 (x_0, x_1) + (y_0, y_1) = (x_0 + y_0,  x_1 + y_1 + x_0*y_0)
  (x_0, x_1) - (y_0, y_1) = (x_0 + y_0,  x_1 + y_1 + x_0*y_0 + y_0^2)
  (x_0, x_1) * (y_0, y_1) = (x_0*y_0,  x_0^2*y_1 + y_0^2*x_1)
The equation we want to solve can be written :
 (alpha,0)^2 - (b_0,b_1) * (beta,0)^2 = (a_0,a_1)
This gives :
    b_0*beta^2 = alpha^2 + a_0
    b_1*beta^4 = a_1 + a_0*alpha^2 + a_0^2
    and we find finally if b1!=0
    beta^4+a0*b0/b1*beta^2+a1/b1=0
    One can find beta by solving an Artin-Schreier equation. Then
    alpha = sqrt(b_0) beta +sqrt(a_0)
And we can compute a_0, a_1, b_0, b1 with :
  a_0 = a (mod 2)  ;  sqrt(a_1) = ([a] - a)/2  (mod 2)
  b_0 = b (mod 2)  ;  sqrt(b_1) = ([b] - b)/2  (mod 2)
where [.] is the Teichmuller representant
"""
def mod_four_norm_equation_solver(a,b):
    R=parent(a)
    r=a.residue().parent().degree()
    a0=a.residue()
    b0=b.residue()
    sqrta1=((teichmuller_mod_four(a)-a)/2).residue()
    sqrtb1=((teichmuller_mod_four(b)-b)/2).residue()
    a1=sqrta1^2
    b1=sqrtb1^2
    if b1!=0:
        c1=a0*b0/b1
        c2=a1/b1
        beta1,_=artin_schreier(c1,c2)
        if beta1 == None:
            return None, None
        else:
            beta = beta1^(2^(r-1))
            alpha=a0^(2^(r-1))*beta+a0^(2^(r-1))
        return R(alpha),R(beta)
    else:
        alpha=((a1+a0^2)/a0)^(2^(r-1))
        beta=((alpha^2+a0)/b0)^(2^(r-1))
        alpha1= R(alpha)  
        beta1= R(beta)
        return R(alpha), R(beta)

"""
Norm equation solver
- input a, b
- returns alpha, beta such that alpha^2 - b*beta^2=a 
"""
def norm_equation_solver(a,b):
    alpha0, beta0=mod_four_norm_equation_solver(a,b)
    if alpha0 == None:
        return None, None
    else:
        alpha0 = alpha0.lift_to_precision(3)
        beta0= beta0.lift_to_precision(3)
        R=((a+b*beta0^2-alpha0^2)/4).residue()
        c1=alpha0.residue()
        flag = True
        while flag: 
            beta1= c1.parent().random_element()
            c2=(b*beta0).residue()*beta1+b.residue()*beta1^2+R
            alpha1,_ = artin_schreier(c1,-c2)
            if alpha1 != None:
                flag = False
        R=alpha0.parent()
        prec=a.parent().precision_cap()
        beta=(beta0+2*R(beta1)).lift_to_precision(prec)
        alpha=(alpha0+2*R(alpha1)).lift_to_precision(prec)
        if alpha.valuation()==0:
            return my_sqrt(a+b*beta^2), beta
        else:
            return alpha, my_sqrt((alpha^2-a)/b)
xcaruso commented 6 years ago
comment:5

I would suggest to implement instead a function that solves the equation a x^2 + b y^2 = 1 which is related to the Hilbert symbol.

You could also have a look at arith/misc.py (line 3944-4045) and rings/number_field/number_field.py (line 8937-9165). There, the computation of global and local Hilbert symbols are implemented (using pari). Maybe you can rely on this code (though your implementation is probably more efficient since it avoids the creation of a number field).

xcaruso commented 6 years ago
comment:6

Moreover, the methods I mentioned only compute the Hilbert symbol but do not exhibit a solution to the quadratic equation a x^2 + b y^2 = 1. So they do not really solve the question you raised.

roed314 commented 6 years ago

Changed keywords from none to padicIMA