On Mon, Dec 22, 2008 at 8:40 AM, John Cremona <john.crem...@gmail.com> wrote: > > That looks very like the exponential integral you are computing. If > so, you can use Sage's function Ei() which calls scipy's > special.exp1().
Watch out, since scipy is double precision only. Pari has a real-only exponential integral that is arbitrary precision though. -- William > > John Cremona > > 2008/12/22 M. Yurko <myu...@gmail.com>: >> >> Alright, below is the original python implementation of the function: >> >> def python(x,bits): >> i = 1 >> sum_current = RealNumber(x,min_prec=bits) >> sum_last = 0 >> term = sum_current >> add_term = (ln(x)+euler_gamma).n(bits) >> while sum_current != sum_last: >> i+=1 >> term = term*(x)/(i) >> sum_last = sum_current >> sum_current += term/i >> return sum_current + add_term >> >> Then my original cython version at double precision: >> %cython >> cdef extern from "math.h": >> double log(double x) >> def cython_double(long double x): >> cdef int i = 1 >> cdef double sum_current = x >> cdef double sum_last = 0 >> cdef double term = sum_current >> cdef double add_term = log(x)+ 0.577215664901533 >> while sum_current != sum_last: >> i+=1 >> term = term*(x)/(i) >> sum_last = sum_current >> sum_current += term/i >> return sum_current + add_term >> >> And finally, the cython version using RealNumber: >> %cython >> from sage.rings.real_mpfr cimport RealNumber >> import sage.all >> from sage.all import log >> from sage.all import n >> def cython_arbitrary(x, int bits): >> cdef int i = 1 >> cdef RealNumber sum_current = sage.all.RealNumber(x,min_prec=bits) >> cdef RealNumber sum_last = sage.all.RealNumber(0, min_prec=bits) >> cdef RealNumber term = sum_current >> cdef RealNumber add_term = sage.all.RealNumber(log(x).n(bits) + >> 0.577215664901533, min_prec=bits) >> while sum_current != sum_last: >> i+=1 >> term = term*(x)/(i) >> sum_last = sum_current >> sum_current += term/i >> return sum_current + add_term >> >> When I timed these functions over 1 through 700 at 53 bits of >> precision, the python one took 5.46 sec., the double precision cython >> one took only .02 sec., and the arbitrary precision one took 6.77 sec. >> After looking at the .c file that cython generated, it seems to be >> doing a lot of conversions as simply initializing sum_current took >> almost 20 long lines of code. >> On Dec 22, 10:24 am, "Mike Hansen" <mhan...@gmail.com> wrote: >>> Hello, >>> >>> On Mon, Dec 22, 2008 at 6:10 AM, M. Yurko <myu...@gmail.com> wrote: >>> >>> > Thanks for your help. I tried your first and last suggestions, but >>> > they yielded code that was slower than the original python >>> > implementation. However, I'll take a look at sage.rings.real_mpfr and >>> > try to use mpfr directly. >>> >>> Well, If I were to guess, it's probably because of the way the Cython >>> code is written. Often when it is slower.that the Python, it means >>> that Cython is doing a lot of conversions behind the scene. If you >>> were to post the code somewhere, I'm sure someone could take a look at >>> it and let you know. Also the annotated version obtained by "cython >>> -a" is useful in tracking these things down. >>> >>> --Mike >> > >> > > > > -- William Stein Associate Professor of Mathematics University of Washington http://wstein.org --~--~---------~--~----~------------~-------~--~----~ To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support-unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-support URLs: http://www.sagemath.org -~----------~----~----~----~------~----~------~--~---
