{{{id=1|
##################################################
# RING-LWE ATTACK (Elias, Lauter, Ozman, Stange) #
##################################################

# General preparation of Sage:  Create a polynomial ring and import GaussianSampler, Timer
P.<y> = PolynomialRing(RationalField(), 'y')
from sage.stats.distributions.discrete_gaussian_lattice import DiscreteGaussianDistributionLatticeSampler
RP = RealField(300) # this sets the precision; if it is insufficient, the implementation won't be valid
from sage.doctest.util import Timer

# Give the Minkowski lattice for a given ring determined by a polynomial.  
# Also gives a key to which are real embeddings.
def cmatrix(): # returns a matrix, columns basis 1, x, x^2, x^3, ... given in the canonical embedding
    global N, a
    N.<a> = NumberField(f)
    fdeg = f.degree()
    key = [0 for i in range(fdeg)] # 0 = real, 1 = real part of complex emb, 2 = imaginary part
    embs = N.embeddings(CC)
    M = matrix(RP,fdeg,fdeg)
    print "Preparing an embedding matrix:  computing powers of the root."
    apows = [ a^j for j in range(n) ]
    print "Finished computing the powers of the root."
    i = 0
    while i < n:
        em = embs[i]
        if Mod(i,20)==Mod(0,20) or Mod(i,20)==Mod(1,20):
            print "Embedding matrix: ", i, " rows out of ", n, " complete."
        if em(a).imag() == 0:
            key[i] = 0
            for j in range(n):
                M[i,j] = em(apows[j]).real()
            i = i + 1    
        else:
            key[i] = 1
            key[i+1] = 2
            for j in range(n):
                M[i,j] = em(apows[j]).real()
                M[i+1,j] = (em(apows[j])*I).real()
            i = i + 2    
    return M, key

# Produce a random vector from (Z/qZ)^n
def random_vec(q, dim): 
    return vector([ZZ.random_element(0,q) for i in range(dim)])

# Useful function for real numbers modulo q
def modq(r,q):
    s = r/q
    t = r/q - floor(r/q)
    return t*q

# Call sampler
def call_sampler(): 
    e = sampler().change_ring(RP)
    return e

# Create samples using a lattice (given by latmat and its inverse), 
# a Gaussian sampler on that lattice, secret, prime
def get_sample(latmat, latmatinv, sec, qval, keyval): 
    e = call_sampler() # create error, in R^n
    dim = latmat.dimensions()[0] # detect dimension of lattice
    pre_a = random_vec(qval, dim) # create a uniformly randomly in terms of basis in cm
    a = latmat*pre_a # create a, in R^n
    b = vecmul_poly(a,sec,latmat,latmatinv) + e # create b, in R^n
    pre_b = latmatinv*b # move to basis in cm in order to reduce mod q
    pre_b_red = vector([modq(c,qval) for c in pre_b])
    b = latmat*pre_b_red
    return [a, b]

# Global choices: setup a field and prime, sampler.
# Set to dummy values that will be altered when an attack is run
q = 1
n = 1
sig = 1/sqrt(2*pi)
Zq = IntegerModRing(q)
R.<x> = PolynomialRing(Zq)
f = y + 1
N.<a> = NumberField(f)
S.<z> = R.quotient(f) # This is P_q
cm,key = cmatrix()
cmi = cm.inverse()
cm
cm53 = cm.change_ring(RealField(10))
cmqq = cm53.change_ring(QQ)
sampler = DiscreteGaussianDistributionLatticeSampler(cmqq.transpose(), sig)
matid = matrix.identity(n)

# Set the parameters for the attack
def setup_params(fval,qval,sval):
    global q,n,sig,f,S,x,z,Zq,matid
    f = fval
    n = f.degree()
    matid = matrix.identity(n)
    q = qval
    Zq = IntegerModRing(q)
    R.<x> = PolynomialRing(Zq)
    sig = sval/sqrt(2*pi)
    S.<z> = R.quotient(f)
    print "Setting up parameters, polynomial = ", f, " and prime = ", q, " and sigma = ", sig
    print "Verifying properties:  "
    print "Prime?", q.is_prime()
    print "Irreducible? ", f.is_irreducible()
    print "Value at 1 modulo q?", Mod(f.subs(y=1),q)
    return True

# Compute the lattices in Minkowski space
def prepare_matrices(polyonly):
    global cm, key, cmi, cmqq
    print "Preparing matrices."
    if polyonly:
        cm = matrix.identity(n)
    else:
        cm,key = cmatrix()
    print "Embedding matrix prepared."
    cmi = cm.inverse()
    print "Inverse matrix found."
    if polyonly:
        cmqq = cm.change_ring(ZZ)
    else:
        cm53 = cm.change_ring(RealField(10))
        cmqq = cm53.change_ring(QQ)
    print "All matrices prepared."
    return True

# Make a vector in R^n into a polynomial, given change of basis matrix and variable to use
def make_poly(a,matchange,var): 
    coeffs = matchange*a  #coefficients of the polynomial are given by the change of basis matrix
    pol = 0
    for i in range(n):
        pol = pol + ZZ(round(coeffs[i]))*var^i # var controls where it will live (what poly ring)
    return pol

# Make a polynomial into a vector in Minkowski space
def make_vec(fval,matchange):
    if fval == 0:
        coeffs = [0 for i in range(n)]
    else:
        coeffs = [0 for i in range(n)]
        colist = lift(fval).coefficients()
        for i in range(len(colist)):
            coeffs[i] = ZZ(colist[i])
    return matchange*vector(coeffs)

# Multiplication in the Minkowski space via moving to polynomial ring
def vecmul_poly(u,v,mat,matinv):
    poly_u = make_poly(u,matinv,z)
    poly_v = make_poly(v,matinv,z)
    poly_prod = poly_u*poly_v
    return make_vec(poly_prod,mat)

# Create the sampler on the lattice embedded in R^n or ZZ^n
def initiate_sampler():
    global sampler
    print "Initiating Sampler."
    sampler = DiscreteGaussianDistributionLatticeSampler(cmqq.transpose(), sig)
    print "Sampler initiated with sigma", RDF(sig)
    return True

# Produce error vectors, just a test to see how they look
def error_test(num):
    print "Testing the error vector production by producing ", num, " errors."
    errorlist = [sampler().norm().n() for _ in range(num)]
    meannorm = mean(errorlist) # average norm
    maxnorm = max(errorlist) # maximum norm
    print "The average error norm is ", RDF(meannorm/( sqrt(n)*sampler.sigma*sqrt(2*pi) )), " times sqrt(n)*s."
    maxratio = RDF(maxnorm/( sqrt(n)*sampler.sigma*sqrt(2*pi) ))
    print "The maximum error norm is ", maxratio, " times sqrt(n)*s."
    if maxratio > 1:
        print "~~~~~~~~~~~~~~~~~~~~~~~ ERROR ~~~~~~~~~~~~~~~~~~~~~~~~~"
        print "The errors do not satisfy a proven upper bound in norm."
    return True

# Create the secret
secret = 0
def create_secret():
    global secret
    secret = cm*random_vec(q,n)
    return True

# Produce samples
samps = []
numsamps = 1
def create_samples(numsampsval):
    global samps, numsamps
    samps = []
    print "Creating samples"
    for i in range(numsampsval):
        print "Creating sample number ", i
        samp = get_sample(cm, cmi, secret, q, key)
        samps.append(samp)
    numsamps = len(samps)
    print "Done creating ", numsamps, "samples."
    return True

# Function for going down to q
def go_to_q(a,matchange):
    pol = make_poly(a,matchange,x)
    #print "debug got pol:", pol
    pol_eval = pol.subs(x=1)
    #print "debug eval'd to:", pol_eval, " and then ", Zq((pol_eval))    
    return Zq(pol_eval)

# Check to make sure moving to q preserves product -- the last two lines should be equal
def sanity_check():
    print "Initiating sanity check"
    mat = cmi
    pvec1 = random_vec(q,n)
    vec1 = cm*pvec1
    pvec2 = random_vec(q,n)
    vec2 = cm*pvec2
    vprod2 = vecmul_poly(vec1,vec2,cm,cmi)
    first_thing = go_to_q(vprod2,mat)
    second_thing = go_to_q(vec1,mat)*go_to_q(vec2,mat)
    if first_thing == second_thing:
        print "Sanity confirmed."
    else:
        print "~~~~~~~~~~~~~~~~~~~~~~~ ERROR ~~~~~~~~~~~~~~~~~~~~~~~~~"
        print "Sanity problem:", first_thing, " is not equal to ", second_thing, "."
        print "Are you sure your ring has root 1 mod q?"
    return True

# Given a list of elements of Z/qZ, make a histogram and zero count
def histoq(data):
    hist = [0 for i in range(10)] # empty histogram
    zeroct=0 # count of zeroes mod q
    for datum in data:
        e = datum
        if e == 0:
            zeroct = zeroct+1
        histbit = floor(ZZ(e)*10/q)
        hist[histbit]=hist[histbit]+1
    return [hist, zeroct]

# Given a list of vectors in R^n, create a histogram of their 
# values in Z/qZ under make_poly, together with a zero count
def histo(data,cmi):
    return histoq([go_to_q(datum,cmi) for datum in data])
    
# Create a histogram of error vectors, transported to polynomial ring
def histogram_of_errors():
    print "Creating a histogram of errors mod q."
    errs = []
    for i in range(80):
        errs.append(sampler())
    hist = histo(errs,cmi)
    print "The number of error vectors that are zero:", hist[1]
    bar_chart(hist[0], width=1).show(figsize=2)
    return True
        
# Create a histogram of the a's in the samples, transported to polynomial ring
def histogram_of_as():
    print "Creating a histogram of sample a's mod q."
    a_vals = [samp[0] for samp in samps]
    hist = histo(a_vals,cmi)
    print "The number of a's that are zero:", hist[1]
    bar_chart(hist[0], width=1).show(figsize=2)
    return True
    
# Create a histogram of errors by correct guess
def histogram_of_errors_2():
    print "Creating a histogram of supposed errors if sample is guessed, mod q."
    hist = histoq([ lift(Zq(go_to_q(sample[1],cmi) - go_to_q(sample[0],cmi)*go_to_q(secret,cmi))) for sample in samps])
    print "The number of such that are zero:", hist[1]
    bar_chart(hist[0], width=1).show(figsize=2)
    return True
    
# Create the secret mod q
lift_s = 0
def secret_mod_q():
    global lift_s
    lift_s = go_to_q(secret,cmi)
    print "Storing the secret mod q.  The secret is ", secret, " which becomes ", lift_s
    return True

# Algorithm 2
# reportrate controls how often it updates the status of the loop; larger = less frequently
# quickflag = True will run only the secret and a few other values to give a quick idea if it works
def alg2(reportrate, quickflag = False): 
    print "Beginning algorithm 2."
    numsamps = len(samps)
    a = [ 0 for i in range(numsamps)]
    b = [ 0 for i in range(numsamps)]
    print "Moving samples to F_q."
    for i in range(numsamps):
        sample = samps[i]
        a[i] = go_to_q(sample[0],cmi)
        b[i] = go_to_q(sample[1],cmi)
    possibles = []
    winner = [[],0]
    print "Samples have been moved to F_q."
    for i in range(2):
        if i == 0:
            print "!!!!! ROUND 1: !!!!! First, checking how many samples the secret survives (peeking ahead)."
            iterat = [lift_s]
        if i == 1:
            print "!!!!! ROUND 2: !!!!! Now, running the attack naively."
            possibles = []
            if quickflag:
                print "We are doing it quickly (not a full test)."
                iterat = xrange(1000)
            else:
                iterat = xrange(q) 
        for g in iterat:
            if Mod(g,reportrate) == Mod(0,reportrate):
                print "Currently checking residue ", g
            g = Zq(g)
            potential = True
            ctr = 0
            while ctr < numsamps and potential:
                e = abs(lift(Zq(b[ctr]-g*a[ctr])))
                if e > q/4 and e < 3*q/4:
                    potential = False
                    if ctr == winner[1]:
                        winner[0].append(g)
                        print "We have a new tie for longest chain:", g, " has survived ", ctr, " rounds."
                    if ctr > winner[1]:
                        winner = [[g],ctr]
                        print "We have a new longest chain of samples survived:", g, " has survived ", ctr, " rounds."
                ctr = ctr + 1
            if potential == True:
                print "We found a potential secret: ", g
                possibles.append(g)
            if g == lift_s:
                if i == 0:
                    print "The real secret survived ", ctr, "samples."
                #break
    print "Full list of survivors of the ", numsamps, " samples:", possibles
    print "The real secret mod q was: ", lift_s
    if len(possibles) == 1 and possibles[0] == lift_s:
        print "Success!"
        return True
    else:
        print "Failure!"
        return False
    
# Run a simulation.
def shebang(fval,qval,sval,numsampsval,numtrials,quickflag=False,polyonly=False):
    global sig
    n = fval.degree()
    if polyonly:
        print "Welcome to the Poly-LWE Attack."
    else:
        print "Welcome to the Ring-LWE Attack."
        print "The attack should theoretically work if the following quantity is greater than 1."
        print "Quantity: ", RDF( qval/( 2*sqrt(2)*sval*n*(qval-1)^( (n-1)/2/n) ) )
    timer = Timer()
    timer2 = Timer()
    timer.start()
    print "********** PHASE 1: SETTING UP SYSTEM "
    setup_params(fval,qval,sval)
    prepare_matrices(polyonly)
    if not polyonly:
        print "Computing the adjustment factor for s."
        cembs = (n - len(N.embeddings(RR)))/2
        detscale = RP( ( 2^(-cembs)*sqrt(abs(f.discriminant())) )^(1/n) ) # adjust the sigma,s
        sval = sval*detscale
        sig = sig*detscale
        print "Adjusted s for use with this embedding, result is ", sval
    initiate_sampler()
    print "The sampler has been created with sigma = ", RDF(sampler.sigma)
    print "Sampled vectors will have expected norm ", RDF(sqrt(n)*sampler.sigma)
    error_test(5)
    print "Time for Phase 1: ", timer.stop()
    timer.start()
    count_successes = 0
    timer2.start()
    for trialnum in range(numtrials):
        print "*~*~*~*~*~*~*~*~*~*~*~*~* TRIAL NUMBER ", trialnum, "*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~"
        print "********** PHASE 2: CREATE SECRET AND SAMPLES"
        create_secret()
        create_samples(numsampsval)
        sanity_check()
        print "Time for Phase 2: ", timer.stop()
        timer.start()
        print "********** PHASE 3: HISTOGRAMS"
        histogram_of_errors()
        print "The histogram of errors (above) should be clustered at edges for success."
        histogram_of_as()
        print "The histogram of a's (above) should be fairly uniform."
        histogram_of_errors_2()
        print "The histogram of sample errors (above) should be clustered at edges for success."
        print "Time for Phase 3: ", timer.stop()
        timer.start()
        print "********** PHASE 4: ATTACK ALGORITHM"
        secret_mod_q()
        result = alg2(100000,quickflag)
        print "Result of Algorithm 2:", result
        print "Time for Phase 4: ", timer.stop()
        if result == True:
            count_successes = count_successes + 1
        print "*~*~*~*~*~*~*~*~*~*~*~*~* ", count_successes, " out of ", trialnum+1, " successes so far. *~*~*~*~*~*"
    totaltime = timer2.stop()
    print "Total time for ", trialnum+1, "trials was ", totaltime
    return count_successes
///
Preparing an embedding matrix:  computing powers of the root.
Finished computing the powers of the root.
Embedding matrix:  0  rows out of  1  complete.
}}}

{{{id=9|
n= 2^4
q = next_prime(2^10-2)
s = 3.192
shebang(y^n+q-1,q,s,40,1, quickflag=False, polyonly=True)
///
}}}

{{{id=10|

///
}}}

{{{id=11|

///
}}}

{{{id=3|

///
}}}

{{{id=5|

///
}}}

{{{id=7|

///
}}}

{{{id=8|

///
}}}