On 9/12/07, John Voight <[EMAIL PROTECTED]> wrote:
> So I took the plunge and started trying to write optimized code in
> Cython for my enumeration problem. The code can be found at
> http://www.cems.uvm.edu/~voight/tr_data.spyx
> http://www.cems.uvm.edu/~voight/totallyreal.py
> and for the experts out there, I'd appreciate any feedback you can
> give, including tricks and such. I also have a bunch of questions,
> some of which are probably just due to my inexperience:
>
> (1) At a certain moment in the algorithm, I need to test if a
> polynomial with integer coefficients is squarefree (using exact
> arithmetic). Is it acceptable practice to do this in a brutal way,
> like:
>
> import sage.rings.all PolynomialRing, IntegerRing, gcd
> ZZx = PolynomialRing(IntegerRing(), 'x')
> [...]
> f = ZZx([...])
> df = ZZx([...])
> if gcd(f,df) <> 1:
>
> What exactly is happening internally there? And what is the overhead
> and effect on the timing?
It's building two polynomials over ZZ represented via ntl, then
calling the ntl integer GCD function.
> (1a) Shouldn't the integer polynomial gcd work using modular
> arithmetic to avoid coefficient blow-up? Or does it do this already?
> (If gcd(f,g) mod p is 1 for p >> 0 then gcd(f,g) = 1, etc.)
I have no idea -- I hope somebody more familiar with NTL internals
can just look and see. You could, of course, do the computation
modulo a prime right now if you wanted, but you'de have to choose
the prime.
> (2) I have my Cython code in a tr_data.spyx file and my Python code in
> a totallyreal.py file. How do I include the former into the latter--
> or do I have to include both at the prompt?
You can't trivially get at tr_data.spyx from totallyreal.py, but you can
from a .sage file, as illustrated in this example:
[EMAIL PROTECTED]:~/tmp$ more a.py
def foo(n):
return n*n
[EMAIL PROTECTED]:~/tmp$ more b.spyx
def bar(m):
return m+m
[EMAIL PROTECTED]:~/tmp$ more c.sage
load a.py
load b.spyx
print foo(10)
print bar(20)
[EMAIL PROTECTED]:~/tmp$ sage
----------------------------------------------------------------------
| SAGE Version 2.8.4.1, Release Date: 2007-09-09 |
| Type notebook() for the GUI, and license() for information. |
----------------------------------------------------------------------
sage: load c.sage
Compiling b.spyx...
100
40
That said, of course you really can get at b.spyx from any .py
file as follows:
[EMAIL PROTECTED]:~/tmp$ more d.py
from sage.all import pyrex
pyrex(open("b.spyx").read())
print bar(10)
[EMAIL PROTECTED]:~/tmp$ sage
----------------------------------------------------------------------
| SAGE Version 2.8.4.1, Release Date: 2007-09-09 |
| Type notebook() for the GUI, and license() for information. |
----------------------------------------------------------------------
sage: import d # pause while b.spyx is automatically compiled.
20
sage: d.bar(30)
60
> (3) I understand that there is significantly more overhead when a
By "more" I assume you mean "less".
> function is declared using
> cdef incr(self)
> instead of
> def incr(self).
> in tr_data.spyx. But if I do this, how do I then access incr in
> totallyreal.py?
You can't. It doesn't make sense, because
cdef incr(self)
is just making a C-only level function.
Any change you would make to make it possible would immediately
put back exactly the overhead you have with
def incr(self).
> (4) Sometimes I need quite verbose output, when debugging and testing,
> and othertimes I just want to turn it all off. Is there a way to set
> a global parameter which controls the verbosity for several routines,
> like the SetVerbose() in magma? Otherwise, I have to pass a "verbose"
> parameter around...
There are
set_verbose and get_verbose
commands in SAGE.
This might not work well for you though -- I wrote it right at the
beginning, SAGE-0.x days, and haven't changed it since, and I've since
learned
a lot more about Python. For example, there is a sophisticated logging
module:
sage: import logging
sage: help(logging)
That might be better for your application.
> (5) Is garbage collection automatic, or do I have to declare a __del__
> function and do a sage_free?
Garbage collection of what? Basically, if it's not obvious that garbage
collection isn't automatic then it is. For example if you put
cdef void* foo = malloc(100000)
randomly in the middle of some Cython code, that memory will
never be freed unless you explicitly free it.
> (6) Craig, did you ever get around to making the line
> K1._pari_().nfisisom(K2):
> for testing isomorphism between number fields any less ugly?
> Shouldn't it read
> is_isomorphic(K1, K2)?
> While you're at it, I need to be able to run an LLL-optimization
> routine to speed up these functions: the function is called polred().
It should be:
sage: k.<a> = NumberField(x^2 + 1)
sage: m.<b> = NumberField(x^2 + 4)
sage: k.is_isomorphic(m)
True
This is not implemented yet.
[... minutes pass ...]
This is now implemented. Do hg_sage.pull(); hg_sage.merge(); hg_sage.ci()
to get it.
> (7) What overhead happens in declaring a number field?
Basically nothing. The polynomial is checked for irreducibility and
that it is monic. If you do
sage: x = polygen(QQ,'x')
sage: %prun for i in range(100): k = NumberField(x**3 + x + 1393923489 + i,'a')
you'll see that the time is dominated by coercing numbers into the rationals.
It's definitey not calling pari's nfinit at all.
> Unless I am
> reading this incorrectly, it now dominates my computation now!
No, pari's nfinit is dominated your computation. That's maybe
because you might be using K1._pari_().nfisisom(K2) to check
isomorphism. Doing K1._pari_() creates the PARI number field
associated to K1, which is a very expensive operation:
sage: gp.nfinit?
nfinit(pol,{flag=0}): pol being a nonconstant irreducible polynomial, gives
the vector: [pol,[r1,r2],discf,index,[M,MC,T2,T,different] (see
manual),r1+r2
first roots, integral basis, matrix of power basis in terms of
integral basis,
multiplication table of basis].
^^^^^^^^^
With the new K1.is_isomorphic(K2) code I just posted this
is completely removed, since it just calls paris nfisisom code
with the two polynomials that define the fields as input.
> sage: prun enumerate_totallyreal_fields(5,10^5)
> 390120 function calls (390117 primitive calls) in 8.774 CPU
> seconds
>
> Ordered by: internal time
>
> ncalls tottime percall cumtime percall
> filename:lineno(function)
> 971 3.795 0.004 3.795 0.004 {method 'nfinit' of
> 'sage.libs.pari.gen.gen' objects}
> 1979 1.024 0.001 1.901 0.001 {method 'is_irreducible'
> of 'sage.rings.polynomial.polynomial_element.Polynomial' objects}
> 7600 0.403 0.000 0.710 0.000
> polynomial_element_generic.py:946(__init__)
> 21686 0.332 0.000 0.381 0.000 rational_field.py:
> 120(__call__)
> 971 0.230 0.000 0.431 0.000 {method '_pari_' of
> 'sage.rings.polynomial.polynomial_element.Polynomial' objects}
> 111586 0.223 0.000 0.223 0.000 {isinstance}
> 8571 0.211 0.000 1.318 0.000 polynomial_ring.py:
> 162(__call__)
> 7581 0.197 0.000 0.197 0.000
> {sage.libs.ntl.ntl.make_new_ZZX}
> 1996 0.169 0.000 0.419 0.000
> polynomial_element_generic.py:878(list)
> 1 0.164 0.164 8.774 8.774 totallyreal.py:
> 24(enumerate_totallyreal_fields)
> 22151 0.142 0.000 0.226 0.000
> polynomial_element_generic.py:1180(degree)
> 971 0.129 0.000 0.247 0.000
> polynomial_element_generic.py:518(__init__)
> 45 0.123 0.003 0.123 0.003 {method 'nfisisom' of
> 'sage.libs.pari.gen.gen' objects}
> 2120 0.119 0.000 0.444 0.000 {method 'incr' of
> 'tr_data_spyx_0.tr_data' objects}
> 90 0.112 0.001 0.130 0.001 {method '_pari_' of
> 'sage.structure.sage_object.SageObject' objects}
> 971 0.101 0.000 0.922 0.001 number_field.py:
> 157(NumberField)
> [...]
>
> It'd better not be the case that invariants of the field are
> automatically computed...!
They are -- see above. But that's I'm guessing in your isomorphism
code.
> (8) When printing a list of integers, why does it print vertically
> instead of horizontally, i.e.
> [1,
> 2,
> ...
> ]
> instead of
> [1, 2, ...]
> or is this a vagary of Python?
It's IPython's pretty printing, I think. Python itself won't do that:
[EMAIL PROTECTED]:~/d/sage$ python
Python 2.5.1 (r251:54863, May 2 2007, 16:27:44)
[GCC 4.1.2 (Ubuntu 4.1.2-0ubuntu4)] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> range(100)
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36,
37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52,
You can also -- I think -- turn it off by modifying
~/.sage/ipython/ipythonrc
Finally, if you do
sage: print range(100)
it all prints on 1 line.
If you make a SAGE sequence type it is also all on one line by
default, unless you pass cr=True as an option:
sage: Sequence(range(100))
[0, 1, 2, 3, 4, 5, 6, 7, 8,...
sage: Sequence(range(5), cr=True)
[
0,
1,
2,
3,
4
]
> Thanks,
Thanks for your good questions and using SAGE!
>
> John Voight
> Assistant Professor of Mathematics
> University of Vermont
> [EMAIL PROTECTED]
> [EMAIL PROTECTED]
> http://www.cems.uvm.edu/~voight/
>
>
> >
>
--
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/
-~----------~----~----~----~------~----~------~--~---
