sqrt mod p

[Steffen]

Hi, I need to find square roots in GF(prime). I did it like this:

y = sqrt(GF(prime)(ySquare))

So far I am not quite happy with the calculation periods and I would
like to know which algorithm is used. In my case is prime % 4 == 1,
which is the hardest case to find the square root mod p. I am sage
newbie and probably its again only a command that I cant find. I tried
sqrt? and similar things, but only got the information that sqrt is a
symbolic function. Could somebody tell me which algo is implemented or
better how to find the implemented algo.

Cheers, Steffen

[John Cremona]

Use the ?? operator to see the algorithm:

sage: a=GF(next_prime(10^6)).random_element()^2;
sage: a.sqrt??

Type: builtin_function_or_method
Base Class: <type 'builtin_function_or_method'>
String Form: <built-in method sqrt of
sage.rings.integer_mod.IntegerMod_int64 object at 0x99055a4>
Namespace: Interactive
Source:
def sqrt(self, extend=True, all=False):
r"""
Returns square root or square roots of self modulo n.

INPUT:
extend -- bool (default: True); if True, return a square
root in an extension ring, if necessary. Otherwise,
raise a ValueError if the square is not in the base ring.
all -- bool (default: False); if True, return *all* square
roots of self, instead of just one.

ALGORITHM: Calculates the square roots mod $p$ for each of the
primes $p$ dividing the order of the ring, then lifts them
p-adically and uses the CRT to find a square root mod $n$.

See also \code{square_root_mod_prime_power} and
\code{square_root_mod_prime} (in this module) for more
algorithmic details.

EXAMPLES:
sage: mod(-1, 17).sqrt()
4
sage: mod(5, 389).sqrt()
86
sage: mod(7, 18).sqrt()
5
sage: a = mod(14, 5^60).sqrt()
sage: a*a
14
sage: mod(15, 389).sqrt(extend=False)
Traceback (most recent call last):
...
ValueError: self must be a square
sage: Mod(1/9, next_prime(2^40)).sqrt()^(-2)
9
sage: Mod(1/25, next_prime(2^90)).sqrt()^(-2)
25

sage: a = Mod(3,5); a
3
sage: x = Mod(-1, 360)
sage: x.sqrt(extend=False)
Traceback (most recent call last):
...
ValueError: self must be a square
sage: y = x.sqrt(); y
sqrt359
sage: y.parent()
Univariate Quotient Polynomial Ring in sqrt359 over Ring
of integers modulo 360 with modulus x^2 + 1
sage: y^2
359

We compute all square roots in several cases:
sage: R = Integers(5*2^3*3^2); R
Ring of integers modulo 360
sage: R(40).sqrt(all=True)
[20, 160, 200, 340]
sage: [x for x in R if x^2 == 40] # Brute force verification
[20, 160, 200, 340]
sage: R(1).sqrt(all=True)
[1, 19, 71, 89, 91, 109, 161, 179, 181, 199, 251, 269,
271, 289, 341, 359]
sage: R(0).sqrt(all=True)
[0, 60, 120, 180, 240, 300]

sage: R = Integers(5*13^3*37); R
Ring of integers modulo 406445
sage: v = R(-1).sqrt(all=True); v
[78853, 111808, 160142, 193097, 213348, 246303, 294637, 327592]
sage: [x^2 for x in v]
[406444, 406444, 406444, 406444, 406444, 406444, 406444, 406444]
sage: v = R(169).sqrt(all=True); min(v), -max(v), len(v)
(13, 13, 104)
sage: all([x^2==169 for x in v])
True

Modulo a power of 2:
sage: R = Integers(2^7); R
Ring of integers modulo 128
sage: a = R(17)
sage: a.sqrt()
23
sage: a.sqrt(all=True)
[23, 41, 87, 105]
sage: [x for x in R if x^2==17]
[23, 41, 87, 105]

"""
if self.is_one():
if all:
return list(self.parent().square_roots_of_one())
else:
return self

if not self.is_square_c():
if extend:

and so on!

John

--
John Cremona

[Steffen]

Thx, I am wondering why I did not try the command a.sqrt?? on my own.

However, it seems as the implemented algorithm is not the most
efficient one. My result from a.sqrt?? from the latest release:

def sqrt(self, extend=True, all=False):
cdef int_fast32_t i, n = self.__modulus.int32
if n > 100:
moduli = self._parent.factored_order()
# Unless the modulus is tiny, test to see if we're in the
really
# easy case of n prime, n = 3 mod 4.
if n > 100 and n % 4 == 3 and len(moduli) == 1 and moduli[0]
[1] == 1:
if jacobi_int(self.ivalue, self.__modulus.int32) == 1:
# it's a non-zero square, sqrt(a) = a^(p+1)/4
i = mod_pow_int(self.ivalue, (self.__modulus.int32+1)/
4, n)
if i > n/2:
i = n-i
if all:
return [self._new_c(i), self._new_c(n-i)]
else:
return self._new_c(i)
elif self.ivalue == 0:
return self
elif not extend:
raise ValueError, "self must be a square"
# Now we use a heuristic to guess whether or not it will
# be faster to just brute-force search for squares in a c
loop...
# TODO: more tuning?
elif n <= 100 or n / (1 << len(moduli)) < 5000:
if all:
return [self._new_c(i) for i from 0 <= i < n if (i*i)
% n == self.ivalue]
else:
for i from 0 <= i <= n/2:
if (i*i) % n == self.ivalue:
return self._new_c(i)
if not extend:
raise ValueError, "self must be a square"
# Either it failed but extend was True, or the generic
algorithm is better
return IntegerMod_abstract.sqrt(self, extend=extend, all=all)

The easier %4 == 3 case seems to be implemented efficiently, but the
%4 == 1 not. The algo from Tonelli and Shanks might be a good solution
here. Any thoughts on other/better algorithm?

Steffen

[John Cremona]

You are right, that is a really stupid algorithm for p=1 (mod 4)!

I will suggest this as something easy to be fixed at Sage Days 6
(which is just starting). One could either use pari to do the sqrt
more efficiently, or implement something like Tonelli-Shanks as you
suggest.

John

--
John Cremona

[Robert Bradshaw]

On Nov 9, 2007, at 11:46 AM, Steffen wrote:

> The easier %4 == 3 case seems to be implemented efficiently, but the
> %4 == 1 not. The algo from Tonelli and Shanks might be a good solution
> here. Any thoughts on other/better algorithm?
>
> Steffen

I actually wrote this code. It eventually dispatches down to
square_root_mod_prime (for large enough p), which has some references
for the algorithm I used. It actually beats pari for large p (100+
digits)--after doing an initial pre-calculation.

It looks like for the range you're dealing with, pari is 3-7 times
faster. However, as much as 25% or so of that is overhead (Sage can
handle sqrts of composite numbers, and return all square roots. The
conversion sage->pari->pari sqrt->sage is slower than doing the
operation simply using Sage.

If someone thinks Tonelli and Shanks will be faster and wants to
implement it, that could be good.

- Robert

sage: p = next_prime(11^8); print p, p%16
214358891 11
sage: K = GF(p)
sage: a = K.random_element()^2
sage: time a.sqrt()
69738242
Time: CPU 0.00 s, Wall: 0.00 s
sage: a = K.random_element()^2
sage: time L = [a.sqrt() for _ in range(10000)]
CPU time: 0.22 s, Wall time: 0.22 s
sage: b = pari(a)
sage: time L = [b.sqrt() for _ in range(10000)]
CPU time: 0.07 s, Wall time: 0.07 s
sage: time L = [K(pari(a).sqrt()) for _ in range(10000)]
CPU time: 0.73 s, Wall time: 0.73 s
sage: from sage.rings.integer_mod import square_root_mod_prime
sage: time L = [square_root_mod_prime(a,p) for _ in range(10000)]
CPU time: 0.15 s, Wall time: 0.15 s

[John Cremona]

[Hello Robert -- are you in Bristol yet?]

I'm puzzled now. My comment on inefficiency and Steffen's were based
on the code

> for i from 0 <= i <= n/2:
> if (i*i) % n == self.ivalue:
> return self._new_c(i)
>

but now when I do a.sqrt?? I do not see that code. Steffen, are you
using 2.8.12?

John

--
John Cremona

[Robert Bradshaw]

On Nov 9, 2007, at 1:43 PM, John Cremona wrote:

> [Hello Robert -- are you in Bristol yet?]

Just got in.

> I'm puzzled now. My comment on inefficiency and Steffen's were based
> on the code
>> for i from 0 <= i <= n/2:
>> if (i*i) % n == self.ivalue:
>> return self._new_c(i)
>>
>
> but now when I do a.sqrt?? I do not see that code. Steffen, are you
> using 2.8.12?

This depends on the size of the modulus. There are three types--
IntegerMod_int32, IntegerMod_int64, and IntegerMod_gmp. The 32-bit
one has that code, and only for small p. If your modulus is beg
enough, a.sqrt?? won't give you that at all.

[John Cremona]

I see. In my example a was

sage: type(a)
<type 'sage.rings.integer_mod.IntegerMod_int64'>

John

--
John Cremona

[Steffen]

Hi, I have implemented the Tonelli and Shanks algorithm which has a
complexity of O(ln^4(p)) and appears to be amazing fast. So far its a
quick and dirty implementation directly in my Sage file and without
error handling. I will ask Martin (when meeting him next time in the
office) how to integrate the implentation properly in Sage. Of course,
if somebody else does it, I am even more happy :-) Here is the current
code which returns x such that x^2 == a mod p:

def step3(b,p,r,x):
# Step 3: Find exponent
if GF(p)(b) == GF(p)(1):
return b,r,x,0
m = 0
while GF(p)(b**(2**m)) != 1:
m = m + 1
if m == r:
return b,r,0,0
return b,r,x,m

def s_root(a,p):
# Step 0: Determine q:
q = 0
e = 0
while q % 2 != 1:
e = e+1
q = (p-1) / 2**e
# Step 1: Find generator
n = ZZ.random_element()
while kronecker(n,p) != -1:
n = ZZ.random_element()
n = GF(p)(n)
z = GF(p)(n**q)
# Step 2: Initialize
y = z
r = e
a = GF(p)(a)
x = GF(p)(a**((q-1)/2))
b = GF(p)(a*(x**2))
x = GF(p)(a*x)
# Step 3:
b,r,x,m = step3(b,p,r,x)
# Step 4: Reduce exponent
while ZZ(m) != ZZ(0):
t = GF(p)(y**(2**(r-m-1)))
y = GF(p)(t**2)
r = GF(p)(m)
x = GF(p)(x*t)
b = GF(p)(b*y)
b,r,x,m = step3(b,p,r,x)
return x

a = GF(17)(13)
print s_root(a,17)

Steffen

On 9 Nov., 21:50, "John Cremona" <john.crem...@gmail.com> wrote:
> I see. In my example a was
>
> sage: type(a)
> <type 'sage.rings.integer_mod.IntegerMod_int64'>
>
> John
>

[William Stein]

On Nov 9, 2007 11:19 PM, Steffen <sre...@gmx.de> wrote:
> Hi, I have implemented the Tonelli and Shanks algorithm which has a
> complexity of O(ln^4(p)) and appears to be amazing fast. So far its a

Could you post some cool impressive benchmarks? :-)

william

--
William Stein
Associate Professor of Mathematics
University of Washington
http://wstein.org

[Steffen]

Hi,

> Could you post some cool impressive benchmarks? :-)
>
> william

ok :-)

Here the result for finding one root in GF(p) where p is 128, 256 and
512 bit prime:

Results for 128 bit prime:
Original:
Time: CPU 0.04 s, Wall: 0.04 s
Number of roots found:
1
Tonelli:

Time: CPU 0.00 s, Wall: 0.00 s

Number of roots found:
1
Results for 256 bit prime:
Original:
Time: CPU 0.96 s, Wall: 0.96 s
Number of roots found:
0
Tonelli:

Time: CPU 0.00 s, Wall: 0.00 s

Number of roots found:
0
Results for 512 bit prime:
Original:
Time: CPU 9.94 s, Wall: 10.09 s
Number of roots found:
0
Tonelli:

Time: CPU 0.00 s, Wall: 0.00 s

Number of roots found:
0

This shows only that the old implementation was inefficient for large
p. The results for the original implementation for 256 and 512 bit are
not contained in the following results, since sqrt(a) returned an
error for the 256 bit and 512 bit case. Here are the results for 100
random values for a. Find x: x^2==a mod p

Results for 128 bit prime:
Original:
Time: CPU 5.36 s, Wall: 5.41 s
Number of roots found:
47
Tonelli:
Time: CPU 0.07 s, Wall: 0.08 s
Number of roots found:
47
Results for 256 bit prime:
Tonelli:
Time: CPU 0.10 s, Wall: 0.10 s
Number of roots found:
48
Results for 512 bit prime:
Tonelli:
Time: CPU 0.30 s, Wall: 0.33 s
Number of roots found:
50

Here is the complete code I used:

A = {}

def go1(p):
X = {}
for i in range(0,len(A)):
X[i] = sqrt(GF(p)(A[i]))
return X

def go2(p):
X ={}
for i in range(0,len(A)):
X[i] = s_root(A[i], p)
return X

def analyseResult(A,X,p):
numberRoots = 0
for i in range(0,len(A)):
x = X[i]
a = A[i]
if x in GF(p):
if (GF(p)(x) != GF(p)(0)) & (GF(p)(x**2) != GF(p)(a)):
print 'Error in implemenation'
elif GF(p)(x) != GF(p)(0):
numberRoots = numberRoots + 1
print 'Number of roots found:'
print numberRoots

n = 100

print 'Results for 128 bit prime:'
p = next_prime(2^(128))
while p % 4 != 1:
p = next_prime(p+1)

for i in range(0,n):
A[i] = GF(p).random_element()

print 'Original:'
time X = go1(p)
analyseResult(A,X,p)

print 'Tonelli:'
time X = go2(p)
analyseResult(A,X,p)

print 'Results for 256 bit prime:'
p = next_prime(2^(256))
while p % 4 != 1:
p = next_prime(p+1)

for i in range(0,n):
A[i] = GF(p).random_element()

#print 'Original:'
#time X = go1(p)
#analyseResult(A,X,p)

print 'Tonelli:'
time X = go2(p)
analyseResult(A,X,p)

print 'Results for 512 bit prime:'
p = next_prime(2^(512))
while p % 4 != 1:
p = next_prime(p+1)

for i in range(0,n):
A[i] = GF(p).random_element()

#print 'Original:'
#time X = go1(p)
#analyseResult(A,X,p)

print 'Tonelli:'
time X = go2(p)
analyseResult(A,X,p)

-- Steffen

[William Stein]

On Nov 10, 2007 11:53 AM, Steffen <sre...@gmx.de> wrote:

> Original:
> Time: CPU 9.94 s, Wall: 10.09 s
> Number of roots found:
> 0
> Tonelli:
> Time: CPU 0.00 s, Wall: 0.00 s
> Number of roots found:
> 0

Wow, that's an improvement!

I've made this trac #1138

http://trac.sagemath.org/sage_trac/ticket/1138

--