2008/12/22 M. Yurko <myu...@gmail.com>: > > 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
You should really report this as a bug to pari (though before you do, please check if there is still a problem with the latest version (2.4.2) of pari. The version currently used in Sage is older). Can you give an example showing how pari is wrong? I have a multiprecision version of this function in eclib (in C++, using NTL's RR type) but it is not currently wrapped in Sage. John Cremona > 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 -~----------~----~----~----~------~----~------~--~---
