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().
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 > > > --~--~---------~--~----~------------~-------~--~----~ 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 -~----------~----~----~----~------~----~------~--~---
