prng_4.py

from ecff.finitefield.finitefield import FiniteField
from ecff.elliptic import EllipticCurve, Point

"""
This will be the first prng using an elliptic curve

note that damn near all of this code assumes we're
using a prime-order finite field, not prime-power!
Thus, our finite field MUST essentially be Z mod p.
If not, all sorts of things will break here.

First EC candidate:
    F = FiniteField(104729,1)
    C = EllipticCurve(a=F(1),b=F(1))
    P = Point(C,F(32770),F(85186))

the subgroup generated by P has 211 elements (prime)


this subgroup takes forever to generate...
P = Point(C,F(4),F(8945))
this subgroup has 35026 elements!
however, the largest prime factor of 35026 is 211...
"""

"""
crap! we may need larger primes. 331337
is cool, let's try it out...

F = FiniteField(331337,1)
C = EllipticCurve(a=F(1),b=F(1))
P = Point(C, F(0), F(1))
P has order 331621 = 53 * 6257

"""

def is_quadratic_residue(n,p):
    """
    n is some integer, p is an odd prime.

    This is obviously only applicable 
    within a prime-order finite field,
    which thankfully is all we're using

    compute Legendre symbol
    """
    n = n%p
    if n == 0:
        return True #0**2 = 0, but BOOORING
    L = (n**((p-1)//2)) % p # integer division is fine if p is an odd prime
    if L == 1:
        return True
    else:
        return False

def tonelli_shanks(n,p):
    """
    n is some integer less than p
    p is an odd prime

    return integer x in Z mod p s.t. x**2 = n (mod p)

    INPUTS WHICH BREAK THIS:
    (2,41)
    (2,104729)

    TODO: return a tuple containing both solutions
    """
    if n == 0:
        return 0

    if is_quadratic_residue(n,p):
        #compute p-1 = Q * 2**S
        Q = p-1
        S = 0
        while Q % 2 == 0:
            #factor out powers of 2
            Q = Q//2
            S += 1
        #print("Q: ",Q," S: ",S,"(p-1,Q*2**S): ",(p-1,Q*2**S))
        if S == 1:
            #p = 3 (mod 4) thus p+1 % 4 = 0
            R = n**((p+1)//4) % p
            return (R,-R % p) # -R is also a solution

        z = 0 #scope
        for z in range(2,p):
            #lame algorithm to find a non-residue
            if not is_quadratic_residue(z,p):
                break
        #print("z: ",z)
        c = z**Q % p
        R = n**((Q+1)//2) % p #wiki says congruence
        t = n**Q % p #wiki says congruence
        M = S
        i = 1 #note: i=0 was screwing this up... how?
        while t % p != 1:
            #compute
            for i in range(1,M):
                #find i by repeated squaring (TODO: optimize me)
                #print("--i: ",i)
                if t**(2**i) % p == 1:
                    #print("break: ",t)
                    break
            b = c**(2**(M-i-1)) % p
            R = R*b % p
            t = t*(b**2) % p
            c = b**2 % p
            M = i
            #DEBUG
            #print("i: ",i)

        return (R, p-R) #2nd solution is p-R

    else:
        #not even a quadratic residue...
        raise Exception("Not a quadratic residue")

def find_point(C,x):
    """
    C: elliptic curve
    x: x coordinate (FiniteField object)

    attempts to find a point on the
    curve C with x coordinate x. Only
    one of the points is returned, as
    the other can be trivially found 
    with negation
    """
    F = C.a.field
    try:
        # y^2 = x^3 + ax + b
        val = x*x*x + C.a * x + C.b
        y = tonelli_shanks(val.n,F.p)[0]
        #print("point:",x,y)
        return Point(C,F(x),F(y))
    except Exception:
        #print("no possible point")
        return None #not a residue => no point on curve

def make_subgroup(P):
    """
    return the group generated by P.
    This is of course a subgroup of
    the group on P's elliptic curve

    this might be REALLY slow. consider
    if find_order is more appropriate
    """
    T=P
    I=0*P
    subgroup = [I]
    while T != I:
        subgroup.append(T)
        T+=P
    return subgroup

def find_order(P):
    """
    Given a point P on an elliptic curve,
    find its order.

    This lets us find the size of the subgroup
    generated by P without keeping track of
    ALL the elements in the subgroup
    """
    T = P
    I = 0*P
    i = 1
    while T != I:
        T+=P
        i+=1
    return i


class prng:
    def __init__(self,n=331337,seed=42):
        """
        create a new prng using a finite field
        of order n. n must be prime, as we basically
        use the integers mod some prime (Z mod p).
        Jeremy Kun's code lets us make this actual
        arithmetic code quite nice. We do take a
        performance hit, but not a *ton*
        """

        """
        TERRIBLE CHOICE:
        self.P = Point(C,F(32770),F(85186))
        Q = 5*P
        """

        """
        seeds with not-incredibly-terrible initial
        cycles: 2819,

        some seeds lead to idntical cycles, but
        there are at least 3 distinct cycles!
        identical cycles:
        4342,2819
        """
        #don't need to keep field or curve beyond constructor
        F = FiniteField(n,1)
        C = EllipticCurve(a=F(1),b=F(1)) #choice of curve assumes default
        self.state = F(seed).n
        #this gives a 'cycle' of 71 elements...
        self.Q = Point(C,F(153116),F(171795)) # <Q> has order 6257
        self.P = Point(C,F(285710),F(143307))

    def get_num(self):
        """
        produce a number, and update our
        internal state

        try to copy the EC_DRBG algorithm

        THIS IS SUPER BROKEN.
        """
        r = (self.state * self.P).x.n
        self.state = (r * self.P).x.n
        t = (r * self.Q).x.n
        #return t & (2**(19-4)-1) # throw away 4 leading bits, assuming our finite field is of order 331337
        return t & (0x7fff)

    def find_repeat(self):
        """
        get random numbers until we receive a 
        number which we've already seen.

        In this algorithm, we essentially are 
        generating output based on a coset of 
        the subgroup generated by P**2, so once
        we repeat an output, we've started to
        loop over again.
        """
        vals = {}
        output = self.get_num()
        while output not in vals:
            vals[output] = self.get_num()
            output = vals[output]
        return len(vals)

    def test_output(self):
        """
        Given the default finite field
        of order 104729, we almost always
        observe a cycle of 26182 elements.

        get this many random numbers, and see
        which are most frequent
        """
        vals = {}
        #found experimentally for default order
        for i in range(211):
            key = self.get_num()
            if key in vals:
                vals[key]+=1
            else:
                vals[key]=1
        print("uniques: ",len(vals))
        sorted_vals = sorted(vals.items(), key=lambda x:x[1], reverse=True)
        top_ten = sorted_vals[:10]
        print("top 10: ",top_ten)
        bottom_ten = sorted_vals[-10:]
        bottom_ten.reverse()
        print("bottom 10: ",bottom_ten)
        """
        gives WAAAY too much output if we use integers mod p
        without a mask

        for i in range(len(sorted_vals)-1):
            if sorted_vals[i+1] < sorted_vals[i]:
                print("break: ",i,", ",sorted_vals[i:i+2])
        """