[sage-devel] Re: number_field_element coercion
On Nov 8, 2007 9:52 PM, mabshoff
[EMAIL PROTECTED] wrote:
[...]
Woah! Can someone explain to me the various calls above? I'd think
this should take epsilon time to coerce the elements of the sequence.
Or perhaps is there another better way to coerce into Z_F (or,
equivalently for me, F)?
There is without a doubt something fishy going on with coercion. See
^
also malb's report with polynomial rings at
http://www.sagetrac.org/sage_trac/ticket/1046
I have some doubt that John Voight's observation above has to do with
Malb's speed regression report.I think it's just that a particular way
of constructing elements in an order (coercing from a list) hasn't
been optimized
one speck since when we implement orders a month ago. And code that
has had zero optimization tends to be slow. The sort answer is that *right now*
it's vastly faster to construct the element of the order via doing arithmetic
instead of explicitly coercing in a list, since we've optimized arithmetic more.
See the timings and examples in the worksheet below.
I've made this
http://trac.sagemath.org/sage_trac/ticket/1134
coerce speed question from john voight
system:sage
{{{id=0|
def stupid_function(n):
Z_F = NumberField(x^2-x-1, 't').maximal_order()
for i in range(n):
Z_F([5,1])
}}}
{{{id=1|
time stupid_function(10^4)
///
CPU time: 7.88 s, Wall time: 9.31 s
}}}
{{{id=10|
def stupid_function(n):
Z_F = NumberField(x^2-x-1, 't').maximal_order()
a,b = Z_F.gens()
for i in range(n):
w = a + 5*b
}}}
{{{id=11|
time stupid_function(10^4)
///
CPU time: 0.05 s, Wall time: 0.05 s
}}}
{{{id=2|
def stupid_function(n):
K = NumberField(x^2-x-1, 't')
for i in range(n):
K([5,1])
}}}
{{{id=3|
time stupid_function(10^4)
///
CPU time: 4.81 s, Wall time: 4.88 s
}}}
{{{id=4|
def stupid_function(n):
K = NumberField(x^2-x-1, 't')
v = [5,1]
for i in range(n):
K(v)
}}}
{{{id=5|
time stupid_function(10^4)
///
CPU time: 4.78 s, Wall time: 4.81 s
}}}
{{{id=6|
def stupid_function(n):
K = NumberField(x^2-x-1, 't')
one = K(1); t = K.gen(); five = K(5)
for i in range(n):
w = five*t + one
}}}
{{{id=7|
time stupid_function(10^4)
///
CPU time: 0.04 s, Wall time: 0.04 s
}}}
{{{id=8|
def stupid_function(n):
K = NumberField(x^2-x-1, 't')
t = K.gen()
for i in range(n):
w = 5*t + 1
}}}
{{{id=9|
time stupid_function(10^4)
///
CPU time: 0.38 s, Wall time: 0.38 s
}}}
{{{id=12|
}}}
--~--~-~--~~~---~--~~
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/
-~--~~~~--~~--~--~---
[sage-devel] sqrt mod p
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 --~--~-~--~~~---~--~~ 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/ -~--~~~~--~~--~--~---
[sage-devel] sqrt mod p
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 --~--~-~--~~~---~--~~ 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/ -~--~~~~--~~--~--~---
[sage-devel] Re: random polynomial generation
On 11/7/07, Martin Albrecht [EMAIL PROTECTED] wrote: Hi everybody, I've attached a 'random_monomial.py' to http://trac.sagemath.org/sage_trac/ticket/980 which implements Steffen's and my proposal. Hey guys, I've attached a patch for this at http://trac.sagemath.org/sage_trac/ticket/980 . Let me know if it works for you. didier Thoughts? Martin -- name: Martin Albrecht _pgp: http://pgp.mit.edu:11371/pks/lookup?op=getsearch=0x8EF0DC99 _www: http://www.informatik.uni-bremen.de/~malb _jab: [EMAIL PROTECTED] --~--~-~--~~~---~--~~ 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/ -~--~~~~--~~--~--~---
[sage-devel] Re: sqrt mod p
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/
-~--~~~~--~~--~--~---
[sage-devel] Re: sqrt mod p
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
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,
[sage-devel] Re: number_field_element coercion
On Fri, 09 Nov 2007 15:07:03 -, John Voight [EMAIL PROTECTED] wrote: Thanks, somehow I knew this was going to become a trac ticket. It is also my suspicion that it is an optimization issue with number fields. It seems really bizarre that it should be calling a polynomial ring constructor! Why? Elements of number fields are stored as polynomial ring elements. (The cost right now is absolutely killing me right now. I've started enumerating relative extensions of number fields, and doing actual computation right now is taking negligible time in comparison to actual stuff! I've already spent a couple of hours trying to figure out exactly where in my code this is happening, but because of the many ways in which elements are being created--and needing to use pari in between for some functionality which has not been implemented, e.g. finding an LLL-reduced basis of the ring of integers--even after my attempts, I had only managed to slow things down!) Oh no. But it's really not surprising given how little time has been spent optimizing number fields. Keep in mind, e.g., that just over a month ago Sage didn't even have rings of integers. -- William --~--~-~--~~~---~--~~ 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/ -~--~~~~--~~--~--~---
[sage-devel] Re: sqrt mod p
[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
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
[sage-devel] Re: sqrt mod p
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
[sage-devel] Re: sqrt mod p
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 =
[sage-devel] Re: sqrt mod p
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 [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:
[sage-devel] Re: Fwd: [sage-support] sage-2.8.12 build report
I tried to sage -upgrade my 2.8.9 this evening on an Ubuntu x686, and I got: Building sage/matrix/misc.c because it depends on sage/matrix/ misc.pyx. touch sage/matrix/misc.pyx; cython --embed-positions --incref-local- binop -I/home/kostadm/sage/devel/sage-main -o sage/matrix/misc.c sage/ matrix/misc.pyx Error converting Pyrex file to C: ... 0, 0, 0) cdef mpz_t* L_row cdef mod_int* A_row for i from 0 = i A._nrows: L_row = L._matrix[i] A_row = A._matrix[i] ^ /home/kostadm/sage/devel/sage-main/sage/matrix/misc.pyx:54:25: Cannot assign type 'sage.matrix.matrix_modn_dense.mod_int *' to 'sage.matrix.misc.mod_int *' sage: Error running cython. sage: There was an error installing modified sage library code. JV --~--~-~--~~~---~--~~ 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/ -~--~~~~--~~--~--~---
[sage-devel] Re: sqrt mod p
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
[sage-devel] Re: Fwd: [sage-support] sage-2.8.12 build report
I just downloaded a new tarball from sage.math. JV --~--~-~--~~~---~--~~ 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/ -~--~~~~--~~--~--~---
[sage-devel] My talk with martin Albrecht
Hi everybody, The slides of Martin and my talk are available at http://sage.math.washington.edu/tmp/talk/. An accompanying SAGE worksheet can be found at http://sage.math.washington.edu/home/malb/SAGE_Demo.sws . Feedback is very welcome :-) William -- 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/ -~--~~~~--~~--~--~---
[sage-devel] Re: ?Matrix
Type sage: Matrix? And the file line tells you that its defined in sage/matrix/ constructor.py - Robert On Nov 9, 2007, at 4:22 PM, John Voight wrote: How do I use matrix (or Matrix?) inside python files which are compiled? Suppose I have a file: def my_matrix(d): return Matrix(ZZ,d,d,[ [i+j for i in range(d)] for j in range(d)]) I include it in my favorite directory, and compile SAGE, and run: sage: my_matrix(3) -- - type 'exceptions.NameError' Traceback (most recent call last) /home/jvoight/sage/devel/sage-main/sage/rings/number_field/ipython console in module() /home/jvoight/sage/local/lib/python2.5/site-packages/sage/rings/ number_field/hot.py in my_matrix(d) 1 def my_matrix(d): 2 return Matrix(ZZ,d,d,[ [i+j for i in range(d)] for j in range(d)]) type 'exceptions.NameError': global name 'Matrix' is not defined In other situations, I type ?foo to see where I need to import foo from. I can't figure that out with Matrix! I tried various combinations, but I thought I would just ask: What is the officially correct way to do this? Thanks, JV --~--~-~--~~~---~--~~ 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/ -~--~~~~--~~--~--~---
[sage-devel] Re: My talk with martin Albrecht
A typo: on page 11, replace permutations groups with permutation groups. On page 38, you say #sage-devel. You might also want to give the irc server's address. Looks like a great intro to Sage. Good luck with the talk. -Bobby On Nov 9, 2007 4:17 PM, William Stein [EMAIL PROTECTED] wrote: Hi everybody, The slides of Martin and my talk are available at http://sage.math.washington.edu/tmp/talk/. An accompanying SAGE worksheet can be found at http://sage.math.washington.edu/home/malb/SAGE_Demo.sws . Feedback is very welcome :-) William -- William Stein Associate Professor of Mathematics University of Washington http://wstein.org -- Bobby Moretti [EMAIL PROTECTED] --~--~-~--~~~---~--~~ 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/ -~--~~~~--~~--~--~---
