Number Theory -- Sage

#Some basic operations on big numbers...

#gcd is part of SAGE's library. x = gcd(89765, 34513) x

#Here is our own code for the gcd. Help practice writing our own functions. def mygcd(m,n): if m*n == 0: return m + n else: ans = mygcd(n, m % n) return ans

""" The extended gcd is also part of SAGE's library. xgcd(a,b) returns 3 numbers: the gcd, and n,m such that na + mb = gc """ xgcd(89765, 34513)

(1, -3761, 9782)

#Checking -3761*89765 + 9782*34513

y = factorial(37) + 1 y

13763753091226345046315979581580902400000001

is_prime(y)

True

z = next_prime(y) z

13763753091226345046315979581580902400000311

""" Here is are three implementation of a^b mod c. This one uses a specific feature of SAGE: the Integer.binary method. """ def powmod(a,b,c): """ Returns a^b mod c. Assumes all inputs are positive integers. """ d = 1 for i in list(Integer.binary(b)): d = mod(d*d, c) if Integer(i) == 1: d = mod(d*a, c) return Integer(d)

""" python has the built-in function pow(a,b,c). It returns a^b mod c. And this implementation is the natural repeated squaring method that does not use the Integer.binary method. """ def powermod(a, b, c): x = a m = b temp = 1 while (m > 0): while ( (m % 2) == 0): x = (x * x) % c m = m // 2 temp = (temp * x) % c m = m - 1 return temp

# Test it first on known numbers. powmod(10, 3, 100000), pow(10,4,100000), powermod(10, 5, 100000000)

(1000, 10000, 100000)

#testing our function again x = powmod(2,10,100) x

#Recall, z is a prime number. Checking Fermat's theorem. powmod(564321, z-1, z)

y

13763753091226345046315979581580902400000001

""" Finding sqrt(a) mod p when p is a prime number. 1. Some facts about the finite field GF(p). a. GF(p) has primitive elements a is primitive if {1,a,a^2, a^3, ..., a^{p-2}} = {1,2,...,p-1} b. q is a quadratic residue mod p if q = x^2 mod p. c. (p-1)/2 numbers in GF(p) are quadratic residues. d. q is a quadratic residue mod p if and only if q^{(p-1)/2} = 1 mod p. 2. Given q, a quadratic residue mod p, we wish to find its square root, that is find x such that x^2 mod p = q. 3. Here are a couple of calculations with explanations. """

""" Does 111111 have a square root mod z? Be careful. Even though (z-1)/2 is an integer, SAGE casts it as rational, so you cannot execute powmod(111111, (z-1)/2, z). Remedy cast it as Integer. There are two ways to do it: Integer((z-1)/2) or (z-1)//2 """ powmod(111111, Integer((z-1)/2),z)

#O.K. 111111 is a quadratic residue mod z k = Integer((z-1)/2) k

6881876545613172523157989790790451200000155

#So we have 111111^k = 1 mod z or: 111111^{k+1} = 111111 mod z. And sqrt(111111) mod z is zr = powmod(111111, Integer((k+1)/2), z) zr

9731020214128031024190258558534610275772105

#Check zr^2 % z

#Great this was easy. It was easy because (z - 1)/2 is odd. So let us try a prime = 1 mod 4. w = next_prime(z) w

13763753091226345046315979581580902400000397

#We were lucky, w - 1 mod 4 = 1. This time we try to find the square root of 11111111 powmod(11111111, Integer((w-1)/2),w)

13763753091226345046315979581580902400000396

#Ooops, 11111111 is not a quadratic residue mod w. So try 1111111111 powmod(1111111111, Integer((w-1)/2),w)

#O.K. we have a quadratic residue. To see the intermediate exponents we shall display them. e1 = Integer((w-1)/2) e1

6881876545613172523157989790790451200000198

e2 = Integer(e1/2) powmod(1111111111, e2, w), e2

(13763753091226345046315979581580902400000396,
3440938272806586261578994895395225600000099)

""" Note: we have here the following equality: 1111111111^e2 = -1 mod w. So if we find a non-quadratic residue say m we'll have m^e1 = -1 and m^e1*1111111111^e2 = 1 or m^e1*1111111111^(e2+1) = 1111111111 and the sqrt(1111111111) = m^(e1/2)*1111111111^((e2+1)/2). """ powmod(7, e1, w)

13763753091226345046315979581580902400000396

# rot will be the sqrt of 1111111111 mod z. rot = (powmod(7, e2, w)*powmod(1111111111, Integer((e2+1)/2), w)) % w

#Checking rot^2 % w, rot

(1111111111, 5378803641663964606359448815127374402828073)

#And the square root of 1111111111 mod w is: 5378803641663964606359448815127374402828073

# Testing our own gcd mygcd(factorial(90), 7^30)

96889010407

factor(96889010407)

7^13

""" Finding sqrt(m) mod p, p a prime number. Let p-1 = (2k+1)2^m. 1. If qr is a quadratic residue then powermod(qr, (2k+1)*2^(m-1), p) = 1. 2. Claim: we can find an integer s such that powmod(qr, 2k+1, p)powmod(s, 2n, p) % p = 1 3. powmod(qr, k+1, p)*powmod(s, n, p) % p = qr. 4. O.K. lets get started """ #p = 3630793, qr = 1234567. Lets test them. p = 3630793 qr = 1234567 is_prime(p), p % 8, powmod(1234567, (p-1)//2, p)

(True, 1, 1)

#Great. All checked. So first lets find an integer s1 such that powmod(s1, (p-1)//2, p) = -1 s1 = 3 while powmod(s1, (p-1)//2, p) == 1: s1 = s1+1 s1

# Now we evaluate qr^(p-1)/2^j and find s. powmod(qr, (p-1)//4, p)

#Great, no adjustment needed. powmod(qr, (p-1)//8, p)

#So we have qr^(p-1)/8 = qr^453849 = -1 mod p. On the other hand 5^(p-1)/2 = -1 mod p. # This means that: qr^453850*5^(p-1)/2 mod p = qr mod p. # Or sqrt(qr) = qr^(453850/2)*5^(p-1)/4 mod p. sqrtqr = powmod(qr, ((p-1)//8 + 1)//2, p)*powmod(5, (p-1)//4, p) % p sqrtqr

#Verifying sqrtqr^2 % p, (p-1)//8

(1234567, 453849)

""" Square roots of a quadratic residue mod pq. If x^2 mod pq = K then K has four different square roots: 1. Since p and q are distinct primes, gcd(p,q) = 1. 2. Let a, b, be integers such that ap + bq = 1 (use the extended gcd to find a and b). 3. +/-apx +/-bqx will be square roots of K mod pq. """

#Finding all four square roots of k mod pq. p = next_prime(23^25) q = next_prime (29^31) key = p*q n = 11111^2 key, n

(2385249619862761231582193966103068754396726923268011122228540812860\
7467299484479, 123454321)

""" We have a large key and a small quadratic residue: 123454321. One of its square roots is 11111 A second square root mod key is key - 11111. We shall find two more square roots mod key, but first a simple verification. """ (key - 11111), (key - 11111)^2 % key

(2385249619862761231582193966103068754396726923268011122228540812860\
7467299473368, 123454321)

#Evaluation verified the second root. The third root. xgcd(p,q)

(1, -688956198424922698242292894275475302861634032,
3524110986023803162673307008147975)

rt3 = (688956198424922698242292894275475302861634032*p + 3524110986023803162673307008147975*q)*11111 % key rt3, rt3^2 % key

(2018632647615889960160069080060061350859538693487589410370743157533\
6142465468408, 123454321)

#Once again, we found another square root and the fourth is key - rt3, (key-rt3)^2 % key

(3666169722468712714221248860430074035371882297804217118577976553271\
324834016071, 123454321)