Thank you to Adrien-Marie Legendre for his symbol

Okay. Let’s say you want to determine the square root of 16. Well, that’s easy … as the answer is 4. But now, what is the square root of…

Thank you to Adrien-Marie Legendre for his symbol

Okay. Let’s say you want to determine the square root of 16. Well, that’s easy … as the answer is 4. But now, what is the square root of this:

6131408466595397710012614540363234062359700887928125190599453844504986420213585686778819217466610878266563405415556

Ans: 2476168101441297080746512578325117519920374855425678540834L

In a world of public key cryptography, we often perform our operations within a finite field and use a prime number to define is finite field. Our operations are then undertaken with a (mod p) operation. Within elliptic curve methods we define the curve as:

y² = x³ + ax + b (mod p)

and where we have our (mod p) operations. But if we calculate a point x on this curve, how do we find out if we have an integer value squared and (mod p)?

If we have the form of =a (mod p), we must find a value of y which results in a value of a (mod p). It is actually a difficult problem to solve. If a solution exists, the value of a is a quadradic residue (mod p). In modular arithmetic this operation is equivalent to a square root of a number (and where x is the modular square root of a modulo p).

For this we turn to Adrien-Marie Legendre who, in 1798, defined the Legendre symbol. In the following we will try and solve for the value of x, and also generate the Legendre symbol value [link]:

Solve x²=12 (mod 13)

with his method we can determine that the answer is 8, as 64 (mod 13) is 12. Some sample code is:

import sys

def modular_sqrt(a, p):
""" Find a quadratic residue (mod p) of 'a'. p
must be an odd prime.

Solve the congruence of the form:
x^2 = a (mod p)
And returns x. Note that p - x is also a root.

0 is returned is no square root exists for
these a and p.

The Tonelli-Shanks algorithm is used (except
for some simple cases in which the solution
is known from an identity). This algorithm
runs in polynomial time (unless the
generalized Riemann hypothesis is false).
"""
# Simple cases
#
if legendre_symbol(a, p) != 1:
return 0
elif a == 0:
return 0
elif p == 2:
return p
elif p % 4 == 3:
return pow(a, (p + 1) / 4, p)

# Partition p-1 to s * 2^e for an odd s (i.e.
# reduce all the powers of 2 from p-1)
#
s = p - 1
e = 0
while s % 2 == 0:
s /= 2
e += 1

# Find some 'n' with a legendre symbol n|p = -1.
# Shouldn't take long.
#
n = 2
while legendre_symbol(n, p) != -1:
n += 1

# Here be dragons!
# Read the paper "Square roots from 1; 24, 51,
# 10 to Dan Shanks" by Ezra Brown for more
# information
#

# x is a guess of the square root that gets better
# with each iteration.
# b is the "fudge factor" - by how much we're off
# with the guess. The invariant x^2 = ab (mod p)
# is maintained throughout the loop.
# g is used for successive powers of n to update
# both a and b
# r is the exponent - decreases with each update
#
x = pow(a, (s + 1) / 2, p)
b = pow(a, s, p)
g = pow(n, s, p)
r = e

while True:
t = b
m = 0
for m in xrange(r):
if t == 1:
break
t = pow(t, 2, p)

if m == 0:
return x

gs = pow(g, 2 ** (r - m - 1), p)
g = (gs * gs) % p
x = (x * gs) % p
b = (b * g) % p
r = m


def legendre_symbol(a, p):
""" Compute the Legendre symbol a|p using
Euler's criterion. p is a prime, a is
relatively prime to p (if p divides
a, then a|p = 0)

Returns 1 if a has a square root modulo
p, -1 otherwise.
"""
ls = pow(a, (p - 1) / 2, p)
return -1 if ls == p - 1 else ls

While we have a trivial example here, we can use it for more complex ones, such has finding a point on the elliptic curve [here]:

A sample run is:

Elliptic curve is:		P-192
Finding elliptic point closest to: 1
Prime number: 6277101735386680763835789423207666416083908700390324961279
a,b -3 2455155546008943817740293915197451784769108058161191238065
(2, 1126956676167578795924565825825899020268914906345645360775L)
(3, 2476168101441297080746512578325117519920374855425678540834L)
(5, 936760824408609109609580731987662341845728162027345586443L)
(6, 61374494507529673497365598443020935064779457192199494327L)
(8, 1539168359597512271047259505090133446672063593980132990812L)
(12, 3464753203279792192409824182683870253677262339932562461307L)
(13, 3288234558942609973454802567986887155175778959720199156770L)
(15, 4548834217212027584647316131131523554591911664904227806291L)
(17, 2148916484007672061843886225501299518817815267521173400039L)
(18, 1600977792967480259538850281480651298625682822208237361467L)
(22, 1682016893107185458056834822961338463540516386180178478778L)

In this case we calculate the y² value for x³+ax + b and then use the Legendre symbol to determine if there is an integer solution for y² (mod p). After this we can perform an modular square root:

def modular_sqrt(a, p):
""" Find a quadratic residue (mod p) of 'a'. p
must be an odd prime.

Solve the congruence of the form:
x^2 = a (mod p)
And returns x. Note that p - x is also a root.

0 is returned is no square root exists for
these a and p.

The Tonelli-Shanks algorithm is used (except
for some simple cases in which the solution
is known from an identity). This algorithm
runs in polynomial time (unless the
generalized Riemann hypothesis is false).
"""
# Simple cases
#
if legendre_symbol(a, p) != 1:
return 0
elif a == 0:
return 0
elif p == 2:
return p
elif p % 4 == 3:
return pow(a, (p + 1) / 4, p)

# Partition p-1 to s * 2^e for an odd s (i.e.
# reduce all the powers of 2 from p-1)
#
s = p - 1
e = 0
while s % 2 == 0:
s /= 2
e += 1

# Find some 'n' with a legendre symbol n|p = -1.
# Shouldn't take long.
#
n = 2
while legendre_symbol(n, p) != -1:
n += 1

# Here be dragons!
# Read the paper "Square roots from 1; 24, 51,
# 10 to Dan Shanks" by Ezra Brown for more
# information
#

# x is a guess of the square root that gets better
# with each iteration.
# b is the "fudge factor" - by how much we're off
# with the guess. The invariant x^2 = ab (mod p)
# is maintained throughout the loop.
# g is used for successive powers of n to update
# both a and b
# r is the exponent - decreases with each update
#
x = pow(a, (s + 1) / 2, p)
b = pow(a, s, p)
g = pow(n, s, p)
r = e

while True:
t = b
m = 0
for m in xrange(r):
if t == 1:
break
t = pow(t, 2, p)

if m == 0:
return x

gs = pow(g, 2 ** (r - m - 1), p)
g = (gs * gs) % p
x = (x * gs) % p
b = (b * g) % p
r = m

Conclusions

While cryptography seems so modern, the core methods it uses are often grounded in the work of individuals on the past.