On Nov 9, 2007 11:19 PM, Steffen <[EMAIL PROTECTED]> 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
> 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" <[EMAIL PROTECTED]> wrote:
> > I see. In my example a was
> >
> > sage: type(a)
> > <type 'sage.rings.integer_mod.IntegerMod_int64'>
> >
> > John
> >
> > On 09/11/2007, Robert Bradshaw <[EMAIL PROTECTED]> wrote:
> >
> >
> >
> >
> >
> > > 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
> >
>
> > > > On 09/11/2007, Steffen <[EMAIL PROTECTED]> wrote:
> >
> > > >> 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
> >
> > > > --
> > > > John Cremona
> >
> > --
> > John Cremona
>
>
> >
>
--
William Stein
Associate Professor of Mathematics
University of Washington
http://wstein.org
--~--~---------~--~----~------------~-------~--~----~
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/
-~----------~----~----~----~------~----~------~--~---
