Working on Cython Code in Sage
Sage Days 24
William Stein
{{{id=101|
///
}}}
Trial Division -- Sage Integers and Preparsing
{{{id=79|
def trial_division_python(n):
if n == 1: return 1
for p in [2, 3, 5]:
if n%p == 0: return p
# Algorithm: only trial divide by numbers that
# are congruent to 1,7,11,13,17,29,23,29 mod 30=2*3*5.
dif = [6, 4, 2, 4, 2, 4, 6, 2]
m = 7; i = 1
while m*m <= n:
if n%m == 0: return m
m += dif[i%8]
i += 1
return n
///
}}}
{{{id=82|
n = 2011*10000000019; n
///
20110000038209
}}}
{{{id=81|
timeit('trial_division_python(20110000038209)')
///
625 loops, best of 3: 686 µs per loop
}}}
{{{id=112|
///
}}}
Trial Division in Pure Python
{{{id=100|
%python
def trial_division_python(n):
if n == 1: return 1
for p in [2, 3, 5]:
if n%p == 0: return p
# Algorithm: only trial divide by numbers that
# are congruent to 1,7,11,13,17,29,23,29 mod 30=2*3*5.
dif = [6, 4, 2, 4, 2, 4, 6, 2]
m = 7; i = 1
while m*m <= n:
if n%m == 0: return m
m += dif[i%8]
i += 1
return n
///
}}}
{{{id=99|
timeit('trial_division_python(20110000038209r)')
///
625 loops, best of 3: 279 µs per loop
}}}
{{{id=105|
///
}}}
Using the Python/C API
{{{id=111|
nm = 'trial%s'%randint(1,10^8)
open(nm+'.c','w').write("""
/************** edit below here ************/
#include
static PyObject * trial_division(PyObject *self, PyObject *args) {
unsigned long n, m=7, i=1, dif[8]={6,4,2,4,2,4,6,2};
if (!PyArg_ParseTuple(args, "k", &n)) return NULL;
if (n==1) return PyInt_FromLong(1);
if (n%2==0) return PyInt_FromLong(2);
if (n%3==0) return PyInt_FromLong(3);
if (n%5==0) return PyInt_FromLong(5);
// Algorithm: only trial divide by numbers that
// are congruent to 1,7,11,13,17,29,23,29 mod 30=2*3*5.
while (m*m <= n) {
if (n%m == 0) return PyInt_FromLong(m);
m += dif[i%8]; i += 1;
}
return PyInt_FromLong(m);
}
static PyMethodDef TrialMethods[] = {
{"trial_division", trial_division, METH_VARARGS, "Trial division."},
{NULL, NULL, 0, NULL} /* Sentinel */
};
/************** edit above here ************/
PyMODINIT_FUNC initMODULE_NAME(void) {
(void) Py_InitModule("MODULE_NAME", TrialMethods);
}
""".replace("MODULE_NAME",nm))
open('setup.py','w').write("""
from distutils.core import setup, Extension
setup (name = 'TrialDivision', version = '1.0',
description = 'This is a C extension',
ext_modules = [Extension('%s', sources = ['%s.c'])])
"""%(nm,nm))
os.system("python setup.py install")
trial = __import__(nm)
///
running install
running build
running build_ext
building 'trial49278231' extension
creating build
creating build/temp.macosx-10.6-i386-2.6
gcc -fno-strict-aliasing -g -O2 -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -I/Users/wstein/sage/build/sage-4.5/local/include/python2.6 -c trial49278231.c -o build/temp.macosx-10.6-i386-2.6/trial49278231.o
creating build/lib.macosx-10.6-i386-2.6
gcc -L/Users/was/build/sage-4.5.rc1/local/lib -bundle -undefined dynamic_lookup build/temp.macosx-10.6-i386-2.6/trial49278231.o -o build/lib.macosx-10.6-i386-2.6/trial49278231.so
ld: warning: directory '/Users/was/build/sage-4.5.rc1/local/lib' following -L not found
running install_lib
copying build/lib.macosx-10.6-i386-2.6/trial49278231.so -> /Users/wstein/sage/build/sage-4.5/local/lib/python2.6/site-packages
running install_egg_info
Removing /Users/wstein/sage/build/sage-4.5/local/lib/python2.6/site-packages/TrialDivision-1.0-py2.6.egg-info
Writing /Users/wstein/sage/build/sage-4.5/local/lib/python2.6/site-packages/TrialDivision-1.0-py2.6.egg-info
}}}
{{{id=109|
timeit('trial.trial_division(20110000038209r)')
///
625 loops, best of 3: 6.16 µs per loop
}}}
{{{id=108|
686/6.1
///
112.459016393443
}}}
{{{id=104|
trial.trial_division(20110000038209)
///
Traceback (most recent call last):
File "", line 1, in
File "_sage_input_137.py", line 10, in
exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("dHJpYWwudHJpYWxfZGl2aXNpb24oMjAxMTAwMDAwMzgyMDkp"),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))' + '\n', '', 'single')
File "", line 1, in
File "/private/var/folders/FE/FEo498bGEOeewp0B4AIIP++++TM/-Tmp-/tmpZG_siq/___code___.py", line 3, in
exec compile(u'trial.trial_division(_sage_const_20110000038209 )' + '\n', '', 'single')
File "", line 1, in
TypeError: argument 1 must be integer, not sage.rings.integer.Integer
}}}
{{{id=113|
///
}}}
Implementing Trial Division using Cython
{{{id=83|
%cython
def trial_division_cython1(n):
if n == 1: return 1
for p in [2, 3, 5]:
if n%p == 0: return p
# Algorithm: only trial divide by numbers that
# are congruent to 1,7,11,13,17,29,23,29 mod 30=2*3*5.
dif = [6, 4, 2, 4, 2, 4, 6, 2]
m = 7; i = 1
while m*m <= n:
if n%m == 0: return m
m += dif[i%8]
i += 1
return n
///
}}}
{{{id=84|
timeit('trial_division_cython1(20110000038209)')
///
125 loops, best of 3: 1.78 ms per loop
}}}
{{{id=85|
%cython
def trial_division_cython2(unsigned long n):
if n == 1: return 1
cdef unsigned long p
for p in [2, 3, 5]:
if n%p == 0: return p
# Algorithm: only trial divide by numbers that
# are congruent to 1,7,11,13,17,29,23,29 mod 30=2*3*5.
dif = [6, 4, 2, 4, 2, 4, 6, 2]
cdef unsigned long m = 7, i = 1
while m*m <= n:
if n%m == 0: return m
m += dif[i%8]
i += 1
return n
///
}}}
{{{id=86|
timeit('trial_division_cython2(20110000038209)')
///
625 loops, best of 3: 43.3 µs per loop
}}}
{{{id=88|
%cython
def trial_division_cython3(unsigned long n):
if n == 1: return 1
cdef unsigned long p
for p in [2, 3, 5]:
if n%p == 0: return p
# Algorithm: only trial divide by numbers that
# are congruent to 1,7,11,13,17,29,23,29 mod 30=2*3*5.
cdef unsigned long dif[8]
dif[0]=6;dif[1]=4;dif[2]=2;dif[3]=4;dif[4]=2;dif[5]=4;dif[6]=6;dif[7]=2
cdef unsigned long m = 7, i = 1
while m*m <= n:
if n%m == 0: return m
m += dif[i%8]
i += 1
return n
///
}}}
{{{id=87|
timeit('trial_division_cython3(20110000038209)')
///
625 loops, best of 3: 6.54 µs per loop
}}}
{{{id=89|
trial_division_cython3(20110000038209)
///
2011L
}}}
{{{id=90|
%cython
def trial_division_cython4(unsigned long n):
if n == 1: return 1
cdef unsigned long p
if n%2==0: return 2
if n%3==0: return 3
if n%5==0: return 5
# Algorithm: only trial divide by numbers that
# are congruent to 1,7,11,13,17,29,23,29 mod 30=2*3*5.
cdef unsigned long dif[8]
dif[0]=6;dif[1]=4;dif[2]=2;dif[3]=4;dif[4]=2;dif[5]=4;dif[6]=6;dif[7]=2
cdef unsigned long m = 7, i = 1
while m*m <= n:
if n%m == 0: return m
m += dif[i%8]
i += 1
return n
///
}}}
{{{id=91|
trial_division_cython4(20110000038209)
///
2011L
}}}
{{{id=92|
timeit('trial_division_cython4(20110000038209)')
///
625 loops, best of 3: 6.33 µs per loop
}}}
{{{id=93|
///
}}}
{{{id=95|
timeit('sage.rings.arith.trial_division(20110000038209)')
///
125 loops, best of 3: 2.75 ms per loop
}}}
OUCH!
{{{id=115|
///
}}}
{{{id=98|
///
}}}
{{{id=78|
///
}}}
other stuff (irrelevant for now...)
{{{id=77|
///
}}}
{{{id=1|
timeit('f*g')
///
625 loops, best of 3: 4.67 µs per loop
}}}
{{{id=40|
R. = ZZ[]
time f = R.random_element(degree=10^6)
///
Time: CPU 1.90 s, Wall: 1.92 s
}}}
{{{id=16|
R. = ZZ[]
f = R([1,2,5,-9]); f
///
-9*x^3 + 5*x^2 + 2*x + 1
}}}
{{{id=6|
g = R([1,2,3,4]); g
///
4*x^3 + 3*x^2 + 2*x + 1
}}}
{{{id=52|
%time
for i in range(10^6):
a = f*g
///
CPU time: 1.95 s, Wall time: 1.96 s
}}}
{{{id=53|
1.95/0.05
///
39.0000000000000
}}}
{{{id=15|
f*g
///
-9*x^3 + 5*x^2 + 2*x + 1
}}}
{{{id=7|
timeit('f*g')
///
625 loops, best of 3: 126 ns per loop
}}}
{{{id=38|
%time
for i in range(10^6):
a = 5r*8r
///
CPU time: 0.21 s, Wall time: 0.21 s
}}}
{{{id=18|
%time
for i in range(10^6):
a = f*g
///
CPU time: 0.30 s, Wall time: 0.30 s
}}}
{{{id=14|
time f.mulbench(g,10^6)
///
-36*x^6 - 7*x^5 + 5*x^4 + 11*x^3 + 12*x^2 + 4*x + 1
Time: CPU 1.32 s, Wall: 1.32 s
}}}
{{{id=25|
time f.mulbench2(g,10^6)
///
Time: CPU 0.17 s, Wall: 0.17 s
}}}
{{{id=44|
time f.mulbench3(g,10^6)
///
Time: CPU 0.06 s, Wall: 0.06 s
}}}
{{{id=33|
time f.mulbench3(g,10^6)
///
Time: CPU 0.06 s, Wall: 0.06 s
}}}
{{{id=34|
time f.mulbench4(g,10^6)
///
Time: CPU 0.06 s, Wall: 0.06 s
}}}
{{{id=35|
1.32 + .17 + .06 + .3
///
1.85000000000000
}}}
{{{id=36|
2.14 - 1.85
///
0.290000000000000
}}}
{{{id=37|
timeit('2*3')
///
625 loops, best of 3: 512 ns per loop
}}}
{{{id=29|
%time
c=39r;d=394r
for i in range(10^6):
a = c*d
///
CPU time: 0.27 s, Wall time: 0.27 s
}}}
{{{id=30|
.32 + .17 + 1.57
///
2.06000000000000
}}}
{{{id=31|
%time
for i in range(10^6):
a = f.degree()
///
CPU time: 0.38 s, Wall time: 0.38 s
}}}
{{{id=46|
int(-2^31)
///
}}}
{{{id=45|
%cython
def f(int n):
print n
///
}}}
{{{id=32|
f(-2^31)
///
-2147483648
}}}
{{{id=47|
f(-2^31-1)
///
Traceback (most recent call last):
File "", line 1, in
File "_sage_input_28.py", line 10, in
exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("ZigtMl4zMS0xKQ=="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))' + '\n', '', 'single')
File "", line 1, in
File "/private/var/folders/FE/FEo498bGEOeewp0B4AIIP++++TM/-Tmp-/tmpFbWhUY/___code___.py", line 3, in
exec compile(u'f(-_sage_const_2 **_sage_const_31 -_sage_const_1 )' + '\n', '', 'single')
File "", line 1, in
File "_Users_wstein__sage_sage_notebook_sagenb_home_admin_478_code_sage22_spyx_0.pyx", line 6, in _Users_wstein__sage_sage_notebook_sagenb_home_admin_478_code_sage22_spyx_0.f (_Users_wstein__sage_sage_notebook_sagenb_home_admin_478_code_sage22_spyx_0.c:416)
def f(int n):
OverflowError: value too large to convert to int
}}}
{{{id=48|
f(2^31-1)
///
2147483647
}}}
{{{id=8|
%magma
R := PolynomialRing(Integers());
f := R![1,2,5,-9];
g := R![1,2,3,4];
t := Cputime();
for i in [1..10^6] do
a := f*g;
end for;
print Cputime(t);
///
1.050
}}}
{{{id=21|
f = 9292; g = 2949
///
}}}
{{{id=20|
%time
for i in range(10^6):
a = f*g
///
CPU time: 0.33 s, Wall time: 0.33 s
}}}
{{{id=9|
R. = Integers(2^12)[]
f = R([1,2,5,-9]); print f
g = R([1,2,3,4]); print g
///
4087*x^3 + 5*x^2 + 2*x + 1
4*x^3 + 3*x^2 + 2*x + 1
}}}
{{{id=73|
f*g
///
4060*x^6 + 4089*x^5 + 5*x^4 + 11*x^3 + 12*x^2 + 4*x + 1
}}}
{{{id=74|
timeit??
///
}}}
{{{id=67|
f = R.random_element(degree=350)
///
hi from modulus
hi from modulus
hi from modulus
hi from modulus
}}}
{{{id=70|
a = R.an_element()
///
}}}
{{{id=71|
timeit('a*a')
///
625 loops, best of 3: 483 ns per loop
}}}
{{{id=69|
timeit('f*f')
///
625 loops, best of 3: 24.7 µs per loop
}}}
{{{id=68|
timeit('-x')
///
625 loops, best of 3: 636 ns per loop
}}}
{{{id=12|
R. = Integers(2^6)[]
f = R([1,2,5,-9]); print f
g = R([1,2,3,4]); print g
///
hi from modulus
hi from modulus
hi from modulus
hi from modulus
hi from modulus
hi from modulus
hi from modulus
hi from modulus
55*x^3 + 5*x^2 + 2*x + 1
hi from modulus
hi from modulus
hi from modulus
hi from modulus
hi from modulus
hi from modulus
hi from modulus
hi from modulus
4*x^3 + 3*x^2 + 2*x + 1
}}}
{{{id=72|
f*g
///
x^5 + x^4 + x^3 + 1
}}}
{{{id=63|
%time
for i in range(10^6):
a = f*g
///
CPU time: 0.97 s, Wall time: 0.98 s
}}}
{{{id=64|
.7/.05
///
14.0000000000000
}}}
{{{id=66|
type(f)
///
}}}
{{{id=65|
///
}}}
{{{id=55|
type(f)
///
}}}
{{{id=61|
f*g
///
4060*x^6 + 4089*x^5 + 5*x^4 + 11*x^3 + 12*x^2 + 4*x + 1
}}}
{{{id=60|
timeit('f*g')
///
625 loops, best of 3: 21.7 µs per loop
}}}
{{{id=56|
timeit('f._mul_zn_poly(g)')
///
625 loops, best of 3: 795 µs per loop
}}}
{{{id=58|
%time
for i in range(10^6):
a = f*g
///
CPU time: 26.12 s, Wall time: 26.91 s
}}}
{{{id=57|
type(f)
///
}}}
{{{id=59|
///
}}}