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
On 9 Nov., 18:24, "John Cremona" <[EMAIL PROTECTED]> wrote:
> 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
>
> On 09/11/2007, Steffen <[EMAIL PROTECTED]> wrote:
>
>
>
>
>
> > 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
--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/
-~----------~----~----~----~------~----~------~--~---
