On Mon, Dec 22, 2008 at 1:51 PM, M. Yurko <myu...@gmail.com> wrote:
>
> First, about my issues with PARI's precision. I just tested the
> following:
> -pari(-10).eint1().n(150)
> and got 2492.2289762418777541164160993503173813223839 which is
> inaccurate after the 14th decimal place. As Stein said, Its quite
> likely that I didn't call PARI correctly, as this is the first time
> that I have directly used that interface.
Yes, you're using the interface incorrectly.
sage: pari.set_real_precision(150)
sage: -pari(-10).eint1(precision=150)
2492.228976241877759138440143998524848989647101430942345358
sage: timeit('-pari(-10).eint1(precision=150)')
625 loops, best of 3: 41.8 µs per loop
Or
sage: gp.set_real_precision(150)
154
sage: gp('-eint1(-10)')
2492.22897624187775913844014399852484898964710143094234538818526713774122742888744417794599665663156560488342454657568480015672868779475213684965774390
> To Bradshaw: Thanks for your
> suggestion as that was pretty much what I was looking for from the
> beginning. The reason that I divide term by i again is that the
> denominator of each term should be i times i factorial, so by dividing
> by I again, but not altering term, I save the need to multiply by i-1
> in the next iteration. About the algorithm, I did some testing for the
> numbers 1 through 10,000 with a step of .1 and I found that
> occasionally it would be one bit off, but it was never more than that,
> so I should probably add one or two bits of precision to the
> requested, but overall, I didn't see a single major deviation
> (although that obviously isn't concrete proof).
> On Dec 22, 2:17 pm, Robert Bradshaw <rober...@math.washington.edu>
> wrote:
>> On Dec 22, 2008, at 8:28 AM, M. Yurko wrote:
>>
>> > 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.
>>
>> I'm not sure why the Python function is faster, this is something to
>> look at. However, I notice that you're not really using the fact that
>> these are real numbers (i.e. not calling any cdef'd methods/using any
>> cdef'd values, etc.) so it's doing a bunch of extra typechecking on
>> each assignment. However, even the pure Python version is slightly
>> slower when compiled (not sure why, shouldn't be the case, but I've
>> madehttp://trac.cython.org/cython_trac/ticket/174for later
>> investigation).
>>
>> To really take advantage of cython, one wants to call mpfr directly.
>> E.g.
>>
>> %cython
>>
>> from sage.all import RealNumber as makeRealNumber, euler_gamma, ln,
>> Integer
>> from sage.rings.real_mpfr cimport RealNumber
>> from sage.libs.mpfr cimport *
>>
>> def cython_arbitrary_mpfr(xx,bits):
>> cdef int i = 1
>> cdef RealNumber x = makeRealNumber(xx,min_prec=bits)
>> cdef RealNumber sum_current = x + 0 # need a copy
>> cdef RealNumber sum_last = makeRealNumber(0)
>> cdef RealNumber term = x + 0 # need a copy
>> cdef RealNumber tmp = makeRealNumber(0, min_prec=bits)
>> cdef add_term = (ln(x)+euler_gamma).n(bits)
>> while mpfr_cmp(sum_current.value, sum_last.value) != 0:
>> i+=1
>> mpfr_mul(term.value, term.value, x.value, GMP_RNDN)
>> mpfr_div_ui(term.value, term.value, i, GMP_RNDN)
>> mpfr_set(sum_last.value, sum_current.value, GMP_RNDN)
>> mpfr_div_ui(tmp.value, term.value, i, GMP_RNDN)
>> mpfr_add(sum_current.value, sum_current.value, tmp.value,
>> GMP_RNDN)
>> return sum_current + add_term
>>
>> which is 12.5x faster than the pure Python version, and the advantage
>> will grow as the required precision goes up (the part outside the
>> loop is stilly pretty slow). I transcribed your algorithm directly
>> for a fair comparison, but I notice you're dividing term by i twice
>> (once at the bottom of the loop, once at the top) which could cut
>> your time by another 35%. Note here I'm using mpfr calls directly,
>> but creating RealNumbers to store the results so I don't have to
>> worry about memory management. The loop is pure C.
>>
>> All that being said, I agree it's worth looking into how we're using
>> pari, and if it's been fixed in the meantime. I'm not confident that
>> the algorithm above will always give the result to the desired
>> precision.
>>
>> - Robert
> >
>
--
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
-~----------~----~----~----~------~----~------~--~---
