2008/12/22 William Stein <wst...@gmail.com>: > > On Mon, Dec 22, 2008 at 9:57 AM, John Cremona <john.crem...@gmail.com> wrote: >> >> 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 would bet money the user (or sage) is somehow mis-using pari, and > that there isn't a bug in pari itself. It's easy to get confused > about how to do arbitrary precision using pari.
That's one reason I asked for an example. I wanted to make sure that precision was not being lost in the pari/Sage interface. But the original poster might have been using pari (or gp) outside Sage. John > > William > >> >> 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 >>> > >>> >> >> > >> > > > > -- > 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 -~----------~----~----~----~------~----~------~--~---
