[seqfan] Extension of A097486, Calculating Pi in the Mandelbrot Set

[Robert Munafo]

Kerry Mitchell pointed out that they usually initialize z to 0, and Jason
Kimberly helped me figure out that I really should have been using 2.0 as
the max amount of the modulus, which means I should have been using the
following PARI/GP code:


A097486(n)=local(a,c,z);c=0.1^n*I-0.75;z=0;a=0;while(abs(z)<2.0,{z=z^2+c;a=a+1});a

Confusing the issue was the fact that the terms I reported earlier were
computed partly with that PARI code and partly with a faster C version that
was using the 2.0 limit.

The difference first shows up in the a(8) value, which should be 314159266,
although my erroneous PARI/GP program was giving 314159267.

If the corrected PARI formula is used, I get a(8)=314159266.

So now I am claiming the extension of A097486 should be:

3, 33, 315, 3143, 31417, 314160, 3141593, 31415927, 314159266, 3141592655,
31415926537

My PARI output with corrected formula:

?
A097486(n)=local(a,c,z);c=0.1^n*I-0.75;z=0;a=0;while(abs(z)<2.0,{z=z^2+c;a=a+1});a
? A097486(0)
%1 = 3
? A097486(1)
%2 = 33
? A097486(2)
%3 = 315
? A097486(3)
%4 = 3143
? A097486(4)
%5 = 31417
? A097486(5)
%6 = 314160
? A097486(6)
%7 = 3141593
? A097486(7)
%8 = 31415927
? A097486(8)
%9 = 314159266


On Fri, Jan 8, 2010 at 21:38, Hans Havermann <pxp at rogers.com> wrote:

> Robert Munafo:
>
> > I discovered that the last term in the current entry, A(7)=31415928,
> > is off
> > by one...
>
> > Can anyone with Mathematica or Maple or Matlab verify?
>
> > First 11 terms: 3, 33, 315, 3143, 31417, 314160, 3141593, 31415927,
> > 314159266, 3141592655, 31415926537
>
> It took me a while to figure out how to squeeze the extra precision
> out of the coding in Mathematica but I think/hope that the following,
> which utilizes 128-digit accuracy, does the job:
>
> $MinPrecision=128; Do[c=SetPrecision[.1^n*I-.75,128]; z=c; a=0;
> While[Abs[z]<4, z=z^2+c; a++]; Print[a], {n,0,8}]
>
> 3
> 33
> 315
> 3143
> 31417
> 314160
> 3141593
> 31415927
> 314159267
>
> So, yes and no.
>

-- 
 Robert Munafo  --  mrob.com