

# This file was *autogenerated* from the file sidon_cryptosystem.sage
from sage.all_cmdline import *   # import sage library

_sage_const_2 = Integer(2); _sage_const_0 = Integer(0); _sage_const_1 = Integer(1); _sage_const_100 = Integer(100); _sage_const_3 = Integer(3); _sage_const_4 = Integer(4)
def ConstructSidon(q, k, r): 
    ##Constructs a Sidon Space in Gq(rk, k) for r > 2 
    
    ##Define the field F_p^k and the polynomial ring over that field
    assert(r!=_sage_const_2 )
    F = GF(q**k, names=('z',)); (z,) = F._first_ngens(1)
    R = F['a']; (a,) = R._first_ngens(1)
    
    ##Construct another extension field of degree r over this ring using an irreducible polynomial of degree r over F_q^k
    irred_poly1 = R.irreducible_element(r)
    irred_poly2 = R.irreducible_element(r)
    F_r = F.extension(irred_poly1)['x']; (x,) = F_r._first_ngens(1)
    F_ = F.extension(irred_poly1)
    ##find a root of this in F_q^rk
    roots = irred_poly2(x).roots()

    y = roots[_sage_const_0 ][_sage_const_0 ]
    
    ##returns a tuple containing y, the subfield F_q^k and the field F_q^rk
    ##The Sidon space is defined as u + u^p*y for u in F 
    return y, F, F_

def ConstructSidon2k(q, k): 
    ##Constructs a Sidon Space in Gq(2k, k), q cannot equal 2 
    assert(q!=_sage_const_2 )
    
    ##Define the field F_p^k and the polynomial ring
    F = GF(q**k, names=('z',)); (z,) = F._first_ngens(1)
    R = F['a']; (a,) = R._first_ngens(1)
    
    p = F.characteristic()
    ##Find and irreducible polynomial of degree 2 with the constant term not in W_q-1 
    irred_poly = R.irreducible_element(_sage_const_2 )
    iterations = _sage_const_0 
   
    while irred_poly.list()[_sage_const_0 ]**((q**k - _sage_const_1 )/(q-_sage_const_1 )) == _sage_const_1  and iterations < _sage_const_100 : 
        irred_poly = R.irreducible_element(_sage_const_2 )
        iterations += _sage_const_1 
    
    assert(iterations != _sage_const_100 )
    
  
    ##Construct the second extension field as a polynomial ring over the first modded out by this irreducible element
    X  = F.extension(irred_poly)
    F_2k = F.extension(irred_poly)['x']; (x,) = F_2k._first_ngens(1)
    roots = irred_poly(x).roots() 
    
    ##Return a root of the polynomial as well as info about the two fields 
    y = roots[_sage_const_0 ][_sage_const_0 ]
    
    return y, F, X, irred_poly.list()[_sage_const_1 ], irred_poly.list()[_sage_const_0 ]
    
    
def factor(product, y, q, F, F_r): 
    ##Factoring algorithm for Sidon spaces with r > 2 
    r = len(list(F_r.modulus())) - _sage_const_1  
    p = F.characteristic()
    assert(r > _sage_const_2 )
    
    ##construct a basis from y 
    basis = [vector(y**n) for n in range(r)]
    mat = matrix(basis)
    
    ##calculate the change of basis matrix to go from the standard basis to our new basis, and compute 
    ##the representation of the product in that basis 
    cob_mat = mat.inverse() 
    product_representation = vector(product)*cob_mat
    
    ##Find the roots of the polynomial we derive from product_represntation, and then compute the original u and v
    assert(i == _sage_const_0  for i in product_representation[_sage_const_3 :])
    product_representation = product_representation[_sage_const_0 :_sage_const_3 ]
    F_ = F['x']; (x,) = F_._first_ngens(1)
    poly = F_(list(product_representation))
    roots = poly.roots()
    
    ##returns u and v up to multiplication from F_p, not the original things multiplied 
    if len(roots) == _sage_const_1 : 
        ans1 = (-_sage_const_1 *_sage_const_1 /roots[_sage_const_0 ][_sage_const_0 ]).nth_root(q - _sage_const_1 )
        ans2 = (-_sage_const_1 *_sage_const_1 /roots[_sage_const_0 ][_sage_const_0 ]).nth_root(q - _sage_const_1 )
        
    else: 
        ans1 = (-_sage_const_1 *_sage_const_1 /roots[_sage_const_0 ][_sage_const_0 ]).nth_root(F.characteristic() - _sage_const_1 )
        ans2 = (-_sage_const_1 *_sage_const_1 /roots[_sage_const_1 ][_sage_const_0 ]).nth_root(F.characteristic() - _sage_const_1 )
    
    if product/((ans1 + ans1**q*y)*(ans2 + ans2**q*y)) != _sage_const_1 : 
        ans1 = ans1*product/((ans1 + ans1**q*y)*(ans2 + ans2**q*y))
    return ans1, ans2

def factor2(product, y,q,F, F_r, b, c): 
    ##factoring algorithm for Sidon space with r = 2 
    p = F.characteristic()
    k = len(F.modulus().list()) - _sage_const_1 
    
    ##construct a basis for F_r that we will use to extract usefull info
    basis = [vector(y**n) for n in (_sage_const_0 , _sage_const_1 )]
    mat = matrix(basis)
    cob_mat = mat.inverse() 
    product_representation = vector(product)*cob_mat
    
    
    ##Figure out the linear transformation x - cx^q 
    basis_transf = []
    identity = matrix.identity(k)
    
    for i in range(k): 
        transformed = vector(F(list(identity[i])) - c*F(list(identity[i]))**q)
        basis_transf += [transformed]
        
    transformation = matrix(basis_transf).inverse() 
    
    ##invert this transformation and use it to calculate uv     
    uv = vector(product_representation[_sage_const_0 ])*transformation
    uv = F(uv)
 
    ##Calculate the last two terms and create a quadratic polynomial
    F_ = F['s']; (s,) = F_._first_ngens(1)
    second_term = product_representation[_sage_const_1 ] + b*(uv)**q
    third_term = uv**q
    poly = F_([vector(uv), vector(second_term), vector(third_term)])
    
    ##get the roots of the polynomial and extract u and v
    roots = poly.roots()
    
    ##returns u and v up to multiplication from F_q, not the original things multiplied 
    if len(roots) == _sage_const_1 : 
        ans1 = (-_sage_const_1 *_sage_const_1 /roots[_sage_const_0 ][_sage_const_0 ]).nth_root(q - _sage_const_1 )
        ans2 = (-_sage_const_1 *_sage_const_1 /roots[_sage_const_0 ][_sage_const_0 ]).nth_root(q - _sage_const_1 )
        
    else: 
        ans1 = (-_sage_const_1 *_sage_const_1 /roots[_sage_const_0 ][_sage_const_0 ]).nth_root(q - _sage_const_1 )
        ans2 = (-_sage_const_1 *_sage_const_1 /roots[_sage_const_1 ][_sage_const_0 ]).nth_root(q - _sage_const_1 )
    
    if product/((ans1 + ans1**q*y)*(ans2 + ans2**q*y)) != _sage_const_1 : 
        ans1 = ans1*product/((ans1 + ans1**q*y)*(ans2 + ans2**q*y))
    return ans1, ans2
    
def publicKey(y,q, F, F_r): 
    ##Generates the public key (and associated private key) using info from the given sidon space 
    p = F.characteristic()
    k = len(F.modulus().list()) - _sage_const_1 
    rk = (len(F_r.modulus().list()) - _sage_const_1 )*k
    basefield = GF(q)
    iterations = _sage_const_0 
    iterations2 = _sage_const_0 
    
    ##Construct bases for the sidon space as well as F_r
    sidonbasis = Matrix(basefield, k, lambda i,j: basefield.random_element())
    while sidonbasis.is_invertible() == False and iterations < _sage_const_100 :
        sidonbasis = Matrix(basefield, k, lambda i,j: basefield.random_element())
        iterations += _sage_const_1 
    F_r_basis = Matrix(basefield, rk, lambda i,j: basefield.random_element())
    while F_r_basis.is_invertible() == False and iterations2 < _sage_const_100 : 
         F_r_basis = Matrix(basefield, rk, lambda i,j: basefield.random_element())
         iterations2 += _sage_const_1 
    
    assert(iterations<_sage_const_100  and iterations2<_sage_const_100 )
    
    ##Get the multiplication table from the basis of the sidon space
    v = [F(list(sidonbasis[i])) for i in range(k)]
    origbasis = sidonbasis
    sidonbasis = [j + j**q*y for j in v]
    mult_table = vector(sidonbasis).column()*vector(sidonbasis).row()
    cob_matrix = F_r_basis.inverse()
    print(cob_matrix)
    vec_list = [[_sage_const_0  for i in range(k)] for j in range(k)]
    
    ##Generate the public key M(V,B)
    for i in range(k): 
        for j in range(k):
            element = mult_table[i][j]
            long_representation = convertToLong(element)
            vec_list[i][j] = list(long_representation*cob_matrix) 
    matrixlist = [_sage_const_0  for i in range(rk)]
    
    ##Return the public key 
    for i in range(rk): 
        matrixlist[i] = Matrix(basefield, k, lambda l,j: vec_list[l][j][i])
    return matrixlist, sidonbasis, mult_table, F_r_basis, origbasis
def convertToLong(element):
    ##Converts an element from a vector over F to a vector over F_q
    long_representation = []
   
    for l in range(len(list(element))): 
        long_representation += list(vector(list(element)[l]))
    long_representation = vector(long_representation)
    long_representation
    return long_representation

def convertFromLong(element, F, F_r):
    ##Converts an element from F_r represented as a vector over F_q to a vector over F
    
    k = len(F.modulus().list()) - _sage_const_1 
    r = len(F_r.modulus().list()) - _sage_const_1 
    final_list = []
    
    for i in range(r): 
        final_list += [F(list(element)[i*k:(i+_sage_const_1 )*k])]
    return(F_r(final_list))
    
def getIndices(y,x, k, r): 
    y = int(y)
    x = int(x)
    rk  = r*k
    t1 = x//rk**_sage_const_3 
    rem = x% rk**_sage_const_3  
    a1 = rem//rk**_sage_const_2 
    rem = rem%rk**_sage_const_2 
    t2 = rem//rk
    rem = rem%rk
    a2= rem 
    
    
    d = y//((k**_sage_const_2 -k)*(k**_sage_const_2  - k)/_sage_const_4 )
    rem = y%((k**_sage_const_2 -k)*(k**_sage_const_2  - k)/_sage_const_4 )
    l = rem//int(((k**_sage_const_2  - k)/_sage_const_2 )) + _sage_const_1 
    rem = rem%((k**_sage_const_2  - k)/_sage_const_2 )
    j = rem + _sage_const_1  
    
   
    lx  = _sage_const_1 
    jx = _sage_const_1 
    ly = _sage_const_0 
    jy = _sage_const_0 
    continuel = True
    continuej = True
    startl = k*(k-_sage_const_1 )/_sage_const_2 
    startj = k*(k-_sage_const_1 )/_sage_const_2 
    for i in range(k-_sage_const_1 ):
        if continuel == True:
            startl -= k-_sage_const_1  - i
            if l> startl:
                lx = k-_sage_const_1  - i 
                ly = l - startl - _sage_const_1 
                continuel = False
        if continuej == True:
            startj -= k-_sage_const_1  - i
            if j> startj:
                jx = k-_sage_const_1  - i 
                jy = j - startj - _sage_const_1 
                continuej = False
    return int(t1), int(a1), int(t2), int(a2), int(d), int(lx), int(ly), int(jx), int(jy)
            
            
        
def getLinearizedCoeff(y, x, k, r, matrixList, c_matrix):
    t1,a1,t2,a2,d,lx,ly,jx,jy = getIndices(y,x,k,r)

    if t1*r*k + a1 > t2*r*k + a2: 
      
        return _sage_const_0 
    c = c_matrix[a1][a2][d]
    coeff = c*(matrixList[t1][jy][ly]*matrixList[t2][jx][lx] - matrixList[t1][jy][lx]*matrixList[t2][jx][ly])
    coeff2 = _sage_const_0 
    if t1*r*k + a1 != t2*r*k + a2:
        
        temp = t1
        t1 = t2
        t2 = temp

        temp = a1
        a1 = a2
        a2 = temp

        c2 = c_matrix[a1][a2][d]
        coeff2 = c2*(matrixList[t1][jy][ly]*matrixList[t2][jx][lx] - matrixList[t1][jy][lx]*matrixList[t2][jx][ly])
    return coeff + coeff2
    
    

