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