Re: Bug#1291: Complex AGM

[Karim Belabas <Karim.Belabas@math.u-bordeaux1.fr>]

* Karim Belabas [2012-02-16 02:37]:
> Package: pari
> Version: since f59d9e0c [2005]
> 
> * Bill Allombert [2012-02-16 01:11]:
> > On Wed, Feb 15, 2012 at 03:24:21PM +0000, John Cremona wrote:
> > > For a definition of what "optimal" means and why it matters for
> > > elliptic curve period computations, see http://arxiv.org/abs/1011.0914
> > > 
> > > I am pretty sure that taking the principal square root will never give
> > > a sequence converging to zero.  Using the optimal branch always gives
> > > the largest limit (and hence the smallest periods), though there is al
> > > ittle more to the question than that.
> > 
> > In theory, I agree that should converge, but in practice (git-9749657)
> > 
> > parisize = 8000000, primelimit = 500509
> > ? agm(.1+I/1000,1)
> >   ***   at top-level: agm(.1+I/1000,1)
> >   ***                 ^----------------
> >   *** agm: the PARI stack overflows !
> 
> That's quite unrelated, the bug is in precision(), precrealexact() to be
> precise.
> 
> (01:45) gp > agm(.1+I/1000.,1)
>                           ^------ important !
> time = 0 ms.
> %1 = 0.42504345251657375080360809311316525277 + 0.0011373632973048871855491442682160760332*I
> 
> Looking at that function precrealexact() [ which I wrote about 7 years ago ],
> I have no idea why it should ever have been useful or necessary. As I see it
> now, 
> 
>   precision(t_COMPLEX with a *small, exact* component [= t_FRAC]) 
> 
> is just broken (larger than it should be).

In case the above was too vague, what commit f59d9e0c did in principle was
[essentially] to change the formula for precision(t_COMPLEX) from 

  precision(x + I*y) := min ( precision(x), precision(y) )

to

  precision(x + I*y) := precision( |x + I*y| ) ~ precision( |x| + |y| )

The new "precision" is in general larger than the old one since
|x| + |y| may have a larger precision than one of x and y. E.g.

  (09:15) gp > precision(1.)
  %1 = 38
  (09:16) gp > precision(1. + 10^100)
  %2 = 154

I was wondering why the formula implemented in precrealexact() [ to be called
when exactly one of x or y has infinite precision ] does not match the above:

  (09:16) gp > precision(.1+I/1000)
  %3 = 57

  (09:16) gp > precision(.1 + 1/1000)
  %4 = 38

Cheers,

    K.B.
--
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]
`