Karim Belabas on Mon, 19 Jul 2010 15:06:42 +0200


[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

Re: precision issue


* John Cremona [2010-07-19 14:35]:
> Consider the following (using a recent svn on a 64-bit machine):
> 
> ? D = -37*4
> %1 = -148
> ? tau = (2+sqrt(D))/4
> %2 = 1/2 + 3.0413812651491098444998421226010335310*I
> ? ellj(tau)
> %3 = -199147904.00096667436353951991474362130 + 4.733165431326070833 E-30*I
> ? real(ellj(tau))
> %4 = -199147904.00096667436353951991474362130
> ? precision(real(ellj(tau)))
> %5 = 38
> ? precision(imag(ellj(tau)))
> %6 = 19
> ? real(tau)
> %7 = 1/2
> ? precision(real(tau))
> %8 = 9223372036854775807
> ? precision(imag(tau))
> %9 = 38
> 
> I don't like the way that j(tau)'s imaginary part only has half the
> precision of the real part.  Since tau is known to have real part
> *exactly* 1/2 we know that j(tau) is real, so could we not set the
> imaginary par of the returned value to exact 0?  (similarly when
> Re(tau)=0).  This must be quite a common special case.

1) Trying to answer your question, I just had a look at the code, and it
must be rewritten anyway. (Correct algorithm in terms of \eta(tau),
implemented in a highly inefficient way.)

2) This is analogous to 
  
  (14:44) precision((1.+1e-30) - 1.)
  %2 = 19


I'll see what I can do.

Cheers,

    K.B.

P.S:
> And secondly, why that bizarre value for precision(real(tau))?  I even get
> 
> ? precision(0)
> %12 = 9223372036854775807
> 
> while looks like a bug since I thought that exact objects had precision 0.

The documented behaviour is :

(14:38) gp > ??precision
precision(x,{n}):

   gives the precision in decimal digits of the PARI object x. If x is an exact
object, the largest single precision integer is returned.


--
Karim Belabas, IMB (UMR 5251)  Tel: (+33) (0)5 40 00 26 17
Universite Bordeaux 1          Fax: (+33) (0)5 40 00 69 50
351, cours de la Liberation    http://www.math.u-bordeaux1.fr/~belabas/
F-33405 Talence (France)       http://pari.math.u-bordeaux1.fr/  [PARI/GP]
`