2008/12/22 William Stein <wst...@gmail.com>:
>
> 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
>
Using the precision= parameter of eint1 works here since -10 has been
converted to pari format as exact quantity, so its eint can be
computed to any precision. (Here the precision is measured in bits.)
It is perhaps safer to give pari a value which is already at the
desired precision, like this:
sage: pari(RealField(150)(-10)).eint1()
-2492.22897624188
Now that looks as if it has low precision, but in fact it is only
printing a feault small number of digits. (William's
pari.set_real_precision() nly has the effect of changing the number of
digits output). The other thing about the output of the previous line
is that its type is still in pari-land:
sage: type(pari(RealField(150)(-10)).eint1())
<type 'sage.libs.pari.gen.gen'>
sage: pari(RealField(150)(-10)).eint1().type()
't_REAL'
and to do more with it in Sage you should convert back to Sage/python:
sage: pari(RealField(150)(-10)).eint1().python()
-2492.2289762418777591384401439985248489896471013
where now all the relevant digits are displayed. If you had wanted
500 bits, that's just
sage: pari(RealField(500)(-10)).eint1().python()
-2492.22897624187775913844014399852484898964710143094234538818526713774122742888744417794599665663156560488342454657568480015672868779475213684965774390402
I'm not sure why Sage is not using pari's eint1 by default here.
pari's function does not handle complex arguments, but seems better
for real arguments than scipy.
John
> 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
-~----------~----~----~----~------~----~------~--~---
