It is the exponential integral, but I want greater than double precision. I tried PARI, and it reports the number in arbitrary precision, but it is only seems to be accurate to double precision as all digits after are completely wrong. Also, I'm trying to compare a few a acceleration ideas for the defining sequence, and this is just the baseline to which I want to compare the other versions.
On Dec 22, 11:43 am, "William Stein" <wst...@gmail.com> wrote: > 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 Washingtonhttp://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 -~----------~----~----~----~------~----~------~--~---
