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| /// }}}