Mandelbrot boundary

Solution 1:

Edit: I answered the question of where the centers of the biggest bulbs are. I think the OP wants an approximation for the boundaries of the bulbs, which may require a more exact radius (size/2), and a better approximation of size; see my edit below.

All bulbs are periodic. If n is the period, than $f^{n}(0)=0$, where $f^{n}=(f^{n-1})^2+c$, and $f^{0}=0$. So the center of all of the bulbs (and cardiods), are all roots of algebraic equations. For n=2, you the two centers you mentioned: $c^2+c=0$, whose two roots are $c=0$ and $c=-1$. For all of these equations, you get the trivial root c=0, which can be trivially factored out. So, for n=3, after factoring out c=0, you get a cubic equation, $x^3+2x^2+x+1$, with one real root, and two complex roots. For n=4, you get a 8th order equation. Factoring out the trivial zeros, c=0 and c=-1, you get a sixth order equation, $c^6 + 3c^5 + 3c^4 + 3c^3 + 2c^2 + 1$, with two real roots, and two pairs of complex conjugate roots. This process can be continued ad nauseam, to find all values of c which are centers of bulbs. One way to find the approximate bulb size, is to say that if $f^{n}=0$, than $f^{2n}=0$ as well, but there is also another nearby zero for $f^{2n}=0$, which gives the approximate size of the bulb.

I implemented the algorithm above. Here are the results, for period(n)=2 through 19, for a size of 0.01. The size of the bulb is printed first, followed by the hyperbolic center for that bulb, which is a zero of the algebraic equation described above. The pari-gp code is posted below these results. There aren't that many big bulbs, and I didn't find any with a period>14! The algorithm uses the approximation that the size of a bulb is $1.6(z^n-z^{2n})$, where $z^{2n}$ is the nearest period 2n zero not equal to the $z^n$ zero. Edit $1.647(z^n-z^{2n})$ might be more accurate, based on the Feigenbaum constant, although that is a limiting bifurcation value, so maybe 1.61 would be better. I used an interesting routine, to estimate where the nearest 2n zero is with a cubic polynomial that factors into a quadratic, and then used Newton's method to get an exact location of the nearest 2n bulb. In the pari-gp code, if $f^n(c)$ is a zero, than "y=estim2nzero(c,n); y=invzero(y,n*2);" returns the nearest 2n zero.

You can run with smaller radius, as an optional second parameter in the periodnzero(n,l) routine. The periodnzero routine is exact, but I haven't run it for n>10 due to memory requirements. The randomperiodn routine matches the periodnzero routine results for n=5 through n=10, for a radius of 0.01, when used with 50 iterations. The randomperiodn routine only searches near the border of the main bug and the 2n bug, so it cannot find all of the smaller bulbs. For example, there is a period(6) mini-mandelbrot at -1.723, with a size of 0.0086, that the randomperiodn(6,50,0.005) doesn't find, that periodnzero(6,0.005) does find.

periodnzero(2)
0.497124226138932614 -1.00000000000000000
periodnzero(3)
0.0288243794158896630 -1.75487766624669276
0.185709338183664371 -0.122561166876653620 + 0.744861766619744237*I
periodnzero(4)
0.113351748952365738 -1.31070264133683288
0.0129942224028632946 -0.156520166833755062 + 1.03224710892283180*I
0.0868986952447806464 0.282271390766913880 + 0.530060617578525299*I
periodnzero(5)
0.0782244856281390343 -0.504340175446244000 + 0.562765761452981964*I
0.0471208943401874649 0.379513588015923745 + 0.334932305597497587*I
periodnzero(6)
0.0511889369361296355 -1.13800066665096451 + 0.240332401262098302*I
0.0282242861365737314 0.389006840569771235 + 0.215850650870819108*I
0.0433219915118591584 -0.113418655949436572 + 0.860569472501573055*I
periodnzero(7)
0.0426479226891188707 -0.622436295041293588 + 0.424878436475629184*I
0.0181617636688071151 0.376008681846767560 + 0.144749371321632865*I
0.0345136362674653694 0.121192786105906486 + 0.610611692210754212*I
periodnzero(8)
0.0246366004359986754 -1.38154748443206147
0.0294627655353046600 -0.999442387206567375 + 0.265387532468407285*I
0.0318189851964300712 -0.359102390112449160 + 0.617353453398826859*I
0.0208607361162561978 0.324819701465459679 + 0.563815622140333520*I
0.0123381391559926738 0.359031062836614452 + 0.100934876864297563*I
periodnzero(9)
0.0265592163668188400 -0.672333174859118535 + 0.337714901737740193*I
0.0177260936911318770 0.328902947296482145 + 0.410912590892970139*I
0.0185193079849810855 -0.0315529748129236492 + 0.790783175406400223*I
0.0203944122420907310 -0.210705526779016425 + 0.804635638170098617*I
randomperiodn(10,50)
0.0191394909327384273 -0.919278545244102030 + 0.247048121273991288*I
0.0203674647559298651 -1.21039963741797970 + 0.152874772077754334*I
0.0182000934472523787 0.0502718348233550644 + 0.630468552293935962*I
0.0195932658053389602 -0.533089681450721463 + 0.602309758248200881*I
0.0115562613503280188 0.408518185116081373 + 0.340038064137432468*I
randomperiodn(11,50)
0.0179922313096332877 -0.697838195122424838 + 0.279304134101365780*I
0.0101876737961521636 0.376030292897131701 + 0.266467605276236516*I
0.0142926479408153391 0.172311886095790515 + 0.570759959055939233*I
0.0170058824618214608 -0.293902905530376044 + 0.632861961802271055*I
randomperiodn(12,50)
0.0116593892097745116 -1.15193862345878364 + 0.269129812348473540*I
0.0134209900466310253 -0.871266416358976791 + 0.222317914834670266*I
0.0155654577554219361 -0.562959829953243594 + 0.471465232041119193*I
0.0116746885689485088 -1.34419625941643939 + 0.0549457739045995046*I
randomperiodn(13,50)
0.0129307773546223782 -0.712577592671273148 + 0.237792568209519750*I
0.0110764626149245331 0.0109444639810614271 + 0.638228893362735870*I
0.0129409256506876911 -0.407104083085098332 + 0.584852842868102271*I
randomperiodn(14,50)
0.0109194116271260426 -1.22997158539723289 + 0.110671419951984272*I
0.0109394099532368402 -1.05674050872754759 + 0.248827389019942445*I
0.0104598592537462953 -0.257053859166882025 + 0.639339268776504382*I
0.0111324718829115060 -0.644083434218572230 + 0.440431179466495531*I
randomperiodn(15,50)
randomperiodn(16,50)
randomperiodn(17,50)
randomperiodn(18,50)
randomperiodn(19,50)

This is the pari-gp program that generated the above results.

default(format,"g0.18");
z=1.0;
precis=precision(z);
plim=10^(-precis/2); /* precision limit */
periodnzero(j,l) = {
  print("periodnzero("j")");
  C0=croot(j);
  for (n=1,length(C0),estimz2(C0[n],j,l));
}
/* return all the imag>0 roots of zn, iterating zn<=zn^2+x */
croot(m) = {
  local(zn,Cr,n,i,v);
  zn=x;
  for (n=1,m-1, zn=zn^2+x; );
  Cr=polroots(zn);
  i=0;
  v=vector(length(Cr));
  for (n=1,length(Cr),
    z=Cr[n];
    if (((real(z)<>0) || (imag(z)<>0)) && (imag(z)>=0),
      if (imag(z)==0, z=real(z));
      i++;
      v[i]=z;
    );
  );
  return(v);
}
estimz2(c,n,l) = {
  local(z,y,i);
  z=0;
  y=0;
  if (l==0,l=0.01);
  /* check to make sure that the z_n root is not repeating by a factor of n */
  for (i=1,n-1,z=z^2+c;if (abs(z)<plim,y=i));
  if (y==0,
    y=estim2nzero(c,n);  /* estimate for 2n zero */
    y=invzero(y,n*2); /* refinement for estimate */
    /* 1.6(zn-z2n) is the approximate size of the bug or bulb */
    if (abs(c-y)>l/1.6,
      print(1.6*abs(c-y)" "c );
      return(y);
    );
    return(0);
  ,
    return(0);
  );
}
/* use newton's method, centered at c */
/* iterating z=z^2+x+c with z only needing a0 and a1 terms */
invzero(c,n)={
  local(z,i,a1);
  z=1;
  while ((abs(z)>plim) && (abs(z)<2),
    z=0;
    i=0;
    while ((i<n) && (abs(z)<2), a1=2*a1*z+1;z=z^2+c;i++);
    c=c-z/a1;
  );
  if (abs(z)>plim, print("invzero covergence error "c" "c+z));
  return(c);
}
/* find nearest 2n zero, assuming that n is a zero */
estim2nzero(c,n) = {
  local(z1,z2,y,i,a0,a1,a2,a3);
  a0=0;a1=0;a2=0;a3=0;
  /* iterating 2n times, centered at c, */
  /* starting with z=0; z<=z^2+x+c, truncating to x^3 terms */
  for (i=1,2*n,
    a3=2*a3*a0+2*a1*a2;
    a2=2*a2*a0+a1^2;
    a1=2*a1*a0+1;
    a0=a0^2+c;
  );
  /* assume zn=0, so z2n=0, so a0=0, */
  /* divide by x, solve quadratic root closest to zero */
  a0=a1;
  a1=a2;
  a2=a3;
  y=a1^2-4*a0*a2;
  z1=(-a1+sqrt(y))/(2*a2);
  z2=(-a1-sqrt(y))/(2*a2);
  if (abs(z1)>abs(z2),z1=z2);
  return(c+z1);
}

shortperiodn(nt,nc,cnt,il) = {
  local(z);
  if (il==0,il=0.01,j);
  z=0;for (i=1,nt,z=truncs(z^2-0+x+nc,cnt+1));
  v=polroots(z);
  for (n=1,cnt,
    z=v[n]+nc;
    y=abs(gfunc(z,nt));
    if (y<0.04,
      z=invzero(z,nt);
      if (abs(imag(z))<1E-21, z=real(z));
      if (((real(z)<>0) || (imag(z)<>0)) && (imag(z)>=0),
        j=1;
        while (abs(oj[j]-z)>1E-21 && (oj[j]<>0), j++);
        if (oj[j]==0,
          if (estimz2(z,nt,il)<>0, oj[j]=z);
        );
      );
    );
  )
}
randomperiodn(nt,n,il)= {
  local(i,j);
  oj=vector(128);
  print("randomperiodn("nt","n")");
  for (i=1,n,
    shortperiodn(nt,cfromz(goldr^i),16,il); /* main cardiod */
    shortperiodn(nt,-1+0.25*(goldr^i),16,il); /* period2 bulb */
  );
}
/* truncate series */
truncs(z,n)={
  local(y,i);
  y=0;
  if (n==0,n=8);
  for (i=0,n-1,y=y+polcoeff(z,i)*x^i);
  return(y);
}
gfunc(c,n)={local(i,z);z=0;for (i=1,n,z=z^2+c);z}
cfromz(z) = {z=z/2-1/2;z=-(z^2-1/4);}
lfromc(c) = {1/2-sqrt(1/4-c);}
gold=2/(sqrt(5)+1);
goldr=exp(2*Pi*I*gold);

Solution 2:

In the discussion, the op asked wheather the bulbs are exactly circular, and if I could post the code to generate the Taylor series for the bulbs. The quick answer is that only the period2 bulb centered at -1 is exactly circular. My first answer gave the center and size of some of the bulbs and cardioids. The size is a rather crude approximation of a bulb or cardiod, since all but one of the bulbs isn't exactly circular, and since the cardioids also need a direction. I enhanced the original pari-gp code to also calculate the mapping of the bulbs/cardioid to the unit circle. The pari-gp code is below. islopetaylor(c,n) returns the taylor series for the period n cardiod or bulb near c.

All of these hyperbolic centers can have their bulbs and/or cardioids mapped to the unit circle, where the boundary of the unit circle corresponds to where the $|\text{slope}|=1$. The inside of the circle, is mapped to points inside the hyperbolic region, where the region has an attracting fixed point of period n.

Finding such a region involves getting an approximation for the nearest fixed point, in the neighborhood of the hyperbolic center. The fixed point is attracting, so we can iterate, starting with z=0 with period n, $z \mapsto f^{n}(z)$, where $f(z)=z^2+c$, and c is the point whose fixed point is desired, so that in the limit, the fixed point is $l=\lim_{k \to \infty} f^{kn}(0)$, $f^{n}(l)=l$.

At the fixed point l, calculate the derivitive. $\frac{dy}{dx}f^{n}(l+x)$, where l is the fixed point, described above. If we call this derivitive fixed point function g(c), then the desired function mapping the bulb/cariod to the unit circle is $g^{-1}(z)$. For a given bulb, $g^{-1}(0)$ is hyperbolic center of that bulb or cardiod. I generated the Taylor series by using a circular Cauchy integral numerical approximation at a radius of 0.5, using 50 sample points; which gives results accurate to approximately 20 decimal digits or so, with a 25 term Taylor series.

Here is an example of how to use the pari-gp code, below, to print out the taylor series for the period4 bulbs and cardioids. EDIT: I also included the results for all of the period3 and period5 bulbs and cardioids. Similar code would work for other periods, and sizes. The prtpoly routine is just for pretty printing, of the first few terms; The polynomial output of of islopetaylor can be used directly. For bulbs, the diameter is approximately twice the absolute value of the x coefficient; the Taylor series is obviously much more exact. My first answer estimate of 1.6x the distance to the nearest 2n bulb is also a good approximation. p=periodnzero(n,l); also works for n=6..10. There are a total of sixteen period6 bulbs and cardioids, most of which are pretty small. p=randomperiodn(n,l) works for the larger bulbs, for p>=5.

p=periodnzero(3,0.001); /* rtn period3 bulbs size>0.001 */
/* taylor series for p3 bulbs, first 7 terms of gt poly  */
for (n=1,length(p),gt=islopetaylor(p[n]); prtpoly(gt,7));
p=periodnzero(4,0.001); /* rtn period4 bulbs size>0.001 */
for (n=1,length(p),gt=islopetaylor(p[n]); prtpoly(gt,7));
p=periodnzero(5,0.001); /* rtn period5 bulbs size>0.001 */

---- below is the output ----

(08:49) gp > p=periodnzero(3,0.001); /* rtn period3 bulbs size>0.001 */
periodnzero(3);
0.0288243794158896630 -1.75487766624669276
0.185709338183664371 -0.122561166876653620 + 0.744861766619744237*I
(09:06) gp > /* taylor series for p3 bulbs, first 7 terms of gt poly  */
(09:06) gp > for (n=1,length(p),gt=islopetaylor(p[n]); prtpoly(gt,7));
{z=   -1.75487766624669276
+x^ 1* 0.00951775795656622345
+x^ 2*-0.00439215019285891483
+x^ 3*-0.000255593673326239667
+x^ 4* 0.00000507529911656475223
+x^ 5* 0.00000242235932732248524
+x^ 6* 0.000000167831763783128979
}
{z=   (-0.122561166876653620 + 0.744861766619744237*I)
+x^ 1*(-0.00475887897828311173 - 0.0943369626738123763*I)
+x^ 2*( 0.00219607509642945741 - 0.00111573891029371424*I)
+x^ 3*( 0.000127796836663119834 + 0.0000966580666527584112*I)
+x^ 4*(-0.00000253764955828237611 + 0.0000131076055547205131*I)
+x^ 5*(-0.00000121117966366124262 + 0.000000329634003020475535*I)
+x^ 6*(-0.0000000839158818915644933 - 0.0000000951268060412508001*I)
}
(09:06) gp > p=periodnzero(4,0.001); /* rtn period4 bulbs size>0.001 */
periodnzero(4);
0.00157165406194145265 -1.94079980652948475
0.113351748952365738 -1.31070264133683288
0.0129942224028632946 -0.156520166833755062 + 1.03224710892283180*I
0.0868986952447806464 0.282271390766913880 + 0.530060617578525299*I
(09:06) gp > for (n=1,length(p),gt=islopetaylor(p[n]); prtpoly(gt,7));
{z=   -1.94079980652948475
+x^ 1* 0.000495813881701186015
+x^ 2*-0.000244518035475789364
+x^ 3*-0.00000233355263604428247
+x^ 4* 0.0000000534459067501123460
+x^ 5* 0.00000000181828170546195257
+x^ 6*-3.23710018614877151 E-12
}
{z=   -1.31070264133683288
+x^ 1* 0.0589801138156112997
+x^ 2* 0.00165047843656841052
+x^ 3* 0.0000693045614883186867
+x^ 4* 0.00000269896930057060465
+x^ 5* 0.0000000523314280419260434
+x^ 6*-0.00000000565936605551606439
}
{z=   (-0.156520166833755062 + 1.03224710892283180*I)
+x^ 1*( 0.00345889762902884997 - 0.00245474033761904195*I)
+x^ 2*(-0.00152469430548367362 + 0.00130179544716022133*I)
+x^ 3*(-0.000143828841409685183 - 0.0000378568944213399944*I)
+x^ 4*( 0.00000385828834515733905 - 0.00000994102399364532939*I)
+x^ 5*( 0.00000140472771445005478 + 0.000000673772055503541835*I)
+x^ 6*(-0.0000000523614237035883056 + 0.000000196431825642151357*I)
}
{z=   ( 0.282271390766913880 + 0.530060617578525299*I)
+x^ 1*(-0.0331968614776850928 - 0.0287926429442936276*I)
+x^ 2*( 0.000821714104937363043 - 0.00133634537333004157*I)
+x^ 3*( 0.000110343336983547981 + 0.0000572780596940444853*I)
+x^ 4*(-0.00000523449594881769755 + 0.0000117988699179895988*I)
+x^ 5*(-0.00000143180256932374878 - 0.000000528764025646146377*I)
+x^ 6*( 0.0000000551927252814394115 - 0.000000186571532024463594*I)
}
(09:06) gp > p=periodnzero(5,0.001); /* rtn period5 bulbs size>0.001 */
periodnzero(5);
0.00124142647819312318 -1.86078252220485487
0.00643004051683055908 -1.62541372512330374
0.0782244856281390343 -0.504340175446244000 + 0.562765761452981964*I
0.00750234541167943880 0.359259224758007439 + 0.642513737138542349*I
0.0471208943401874649 0.379513588015923745 + 0.334932305597497587*I
0.00540957192966561334 -0.0442123577040706231 + 0.986580976280892768*I
0.00109467916151353002 -0.198042099364253840 + 1.10026953729269853*I
0.00454737524002394192 -1.25636793006818076 + 0.380320963472722507*I
(09:06) gp > for (n=1,length(p),gt=islopetaylor(p[n]); prtpoly(gt,7));
{z=   -1.86078252220485487
+x^ 1* 0.000390175993253921407
+x^ 2*-0.000193606895721835075
+x^ 3*-0.00000111641882446931317
+x^ 4* 0.0000000945861637002427640
+x^ 5* 0.00000000182079818330302844
+x^ 6*-6.26858107084452022 E-11
}
{z=   -1.62541372512330374
+x^ 1* 0.00202943297897757605
+x^ 2*-0.00100562087600704118
+x^ 3*-0.0000100969021567375794
+x^ 4* 0.00000299007153447035234
+x^ 5* 0.0000000448549883194611938
+x^ 6*-0.0000000140393194397544057
}
{z=   (-0.504340175446244000 + 0.562765761452981964*I)
+x^ 1*( 0.0227804277538978753 - 0.0313962947316761485*I)
+x^ 2*(-0.000187493753264388258 + 0.000240573004638874982*I)
+x^ 3*(-0.0000173170154681859400 + 0.0000483439272385087499*I)
+x^ 4*( 0.00000186140591725642208 - 0.00000140612260499651276*I)
+x^ 5*(-0.0000000359687742865028775 - 0.000000238921821822236246*I)
+x^ 6*(-0.0000000140479390745627075 + 0.0000000154228924213440731*I)
}
{z=   ( 0.359259224758007439 + 0.642513737138542349*I)
+x^ 1*( 0.00125062477418773157 - 0.00205512908763713047*I)
+x^ 2*(-0.000473500958540455897 + 0.00108015520081310997*I)
+x^ 3*(-0.000108144537904337611 - 0.0000221136751964296838*I)
+x^ 4*( 0.00000365352385186322974 - 0.0000109205152171821385*I)
+x^ 5*( 0.00000143644262907468698 + 0.000000753636636313957017*I)
+x^ 6*(-0.000000133837667274211728 + 0.000000211504575842695306*I)
}
{z=   ( 0.379513588015923745 + 0.334932305597497587*I)
+x^ 1*(-0.0230612323455722479 - 0.00538785818574664379*I)
+x^ 2*( 0.000223715667692254576 - 0.00100096226615781152*I)
+x^ 3*( 0.0000924390752585046801 + 0.0000283948327212593745*I)
+x^ 4*(-0.00000444777220263711731 + 0.0000109991366640167731*I)
+x^ 5*(-0.00000147488363872001793 - 0.000000750149091456313290*I)
+x^ 6*( 0.000000131362683170454543 - 0.000000211363105809565828*I)
}
{z=   (-0.0442123577040706231 + 0.986580976280892768*I)
+x^ 1*(-0.00168631410697906826 + 0.000419905880430706866*I)
+x^ 2*( 0.000824122737737012684 - 0.000147452933665481287*I)
+x^ 3*( 0.0000148294436540590271 - 0.0000423012304388835166*I)
+x^ 4*(-0.00000179611219530982926 + 0.000000252886482712731910*I)
+x^ 5*( 0.000000132428561354388846 + 0.000000194421657850058569*I)
+x^ 6*( 0.0000000186904817675855266 - 0.0000000121128281556966309*I)
}
{z=   (-0.198042099364253840 + 1.10026953729269853*I)
+x^ 1*( 0.000338411530274354574 + 0.0000551258574846720436*I)
+x^ 2*(-0.000169328304150344563 - 0.0000234867261464734632*I)
+x^ 3*( 0.0000000298451123227361170 - 0.00000280280755567260578*I)
+x^ 4*( 0.0000000425917226359921085 + 0.0000000627926935964822686*I)
+x^ 5*(-0.00000000297733729350064710 + 0.00000000108485294983731985*I)
+x^ 6*( 6.65684728185643910 E-12 - 9.88920256175429103 E-11*I)
}
{z=   (-1.25636793006818076 + 0.380320963472722507*I)
+x^ 1*(-0.000846269215461327949 - 0.00116600657512957892*I)
+x^ 2*( 0.000389349729034723316 + 0.000591786929496725821*I)
+x^ 3*( 0.0000237849154393727229 - 0.00000452582161657443604*I)
+x^ 4*(-0.000000856099479577055100 - 0.00000108478239958016467*I)
+x^ 5*(-0.0000000783802180386627187 + 0.0000000644312584363213421*I)
+x^ 6*( 0.00000000487678960666910970 + 0.00000000500884893452725251*I)
}
(09:06) gp > p=periodnzero(5,0.001); /* rtn period3 bulbs size>0.001 */
periodnzero(5);
0.00124142647819312318 -1.86078252220485487
0.00643004051683055908 -1.62541372512330374
0.0782244856281390343 -0.504340175446244000 + 0.562765761452981964*I
0.00750234541167943880 0.359259224758007439 + 0.642513737138542349*I
0.0471208943401874649 0.379513588015923745 + 0.334932305597497587*I
0.00540957192966561334 -0.0442123577040706231 + 0.986580976280892768*I
0.00109467916151353002 -0.198042099364253840 + 1.10026953729269853*I
0.00454737524002394192 -1.25636793006818076 + 0.380320963472722507*I

Here is the pari-gp program to generate this. This code also includes updates for the code in the original answer well. More updates, including the plotsetup routine, and maketheplot routine; see the picture below.

print ("periodnzero(n,l); /* calculate all periodn bulbs size>l, returns vector */");
print ("randomperiodn(n,l); /* n>4 approximate all periodn bulbs size>l rtn vect */");
print ("invzero(c,n); /* use Newton's method */");
print ("estim2nzero(c,n); /* find nearest 2n zero */");
print ("gt=islopetaylor(c,n); /* taylor series for c, n iterations */ ");
print ("prtpoly(gt,10); /* print taylor series for islopetaylor */ ");
print ("plotsetup(low,high,numb,radius,randn);   maketheplot;");
default(format,"g0.18");
z=1.0;
precis=precision(z);
plim=10^(-0.5*precis); /* precision limit */
/* plim=10^(-0.65*precis); */ /* precision limit */
periodnzero(n,l) = {
  local(z,i,pcur,otemp,Cp,oout);
  otemp=vector(2^(n-1));
  pcur=0;
  zn=0;
  print("periodnzero("n");");
  Cp=croot(n);
  for (i=1,length(Cp),
    z=estimz2(Cp[i],n,l);
    if (z<>0,pcur++;otemp[pcur]=Cp[i]);
  );
  oout=vector(pcur);
  i=1;while (otemp[i]<>0,oout[i]=otemp[i];i++);
  return(oout);
}
/* return all the imag>0 roots of zn, iterating zn<=zn^2+x */
croot(m) = {
  local(zn,Cr,n,i,v);
  zn=x;
  for (n=1,m-1, zn=zn^2+x; );
  Cr=polroots(zn);
  i=0;
  v=vector(length(Cr));
  for (n=1,length(Cr),
    z=Cr[n];
    if (((real(z)<>0) || (imag(z)<>0)) && (imag(z)>=0),
      if (imag(z)==0, z=real(z));
      i++;
      v[i]=z;
    );
  );
  return(v);
}
estimz2(C0,n,l) = {
  local(z,y,i);
  z=0;
  y=0;
  if (l==0,l=0.01);
  /* check to make sure that the z_n root is not repeating by a factor of n */
  for (i=1,n-1,z=z^2+C0;if (abs(z)<plim,y=i));
  if (y==0,
    y=estim2nzero(C0,n);  /* estimate for 2n zero */
    y=invzero(y,n*2); /* refinement for estimate */
    /* 1.647(zn-z2n) is approx size of the bug or bulb, based on Feignenbaum limit */
    if (1.6*abs(C0-y)>l,
      print(1.6*abs(C0-y)" "C0 );
      return(y);
    );
    return(0);
  ,
    return(0);
  );
}
sizez2(C0,n) = {
  1.6*abs(C0-invzero(estim2nzero(C0,n),n*2));
}
/* use newton's method, centered at C0 */
/* iterating z=z^2+x+C0 with z only needing a0 and a1 terms */
/* s<>0 is used to calculate inv(s) instead of inv(0) */
/* m<>0 is used to calculate Misiurewicz point near C0,n,m */
invzero(C0,n,s,m)={
  local(z,i,j,a1,sz,sa1);
  z=1;
  j=0;
  while ( ((abs(z)>plim)||(j<1)) && (abs(z)<2) && (j<50),
    j++;
    z=0;  a1=0;
    sz=0; sa1=0; /* Misiu */
    i=0;
    while ((i<n) && (abs(z)<2),
      a1=2*a1*z+1;
      z=z^2+C0;
      i++;
      if (i==m, sz=z; sa1=a1; ); /* Misiu */
    );
    z=z+sz-s;
    a1=a1+sa1;
    C0=C0-z/a1;
/*  if (j<3, C0=C0-z/(2*a1), C0=C0-z/a1); */
/*  C0=C0-(z-s)/a1; */
/*  if ((abs(z-s)>plim),print(z-s)); */
  );
  jm=0.001*a1^-2;
  if (abs(z)>plim, return(0));
  return(C0);
}
/* find nearest 2n zero, assuming that n is a zero */
estim2nzero(C0,n) = {
  local(z1,z2,y,i,a0,a1,a2,a3);
  if (C0==0,return(-1));
  a0=0;a1=0;a2=0;a3=0;
  /* iterating 2n times, centered at C0, */
  /* starting with z=0; z<=z^2+x+C0, truncating to x^3 terms */
  for (i=1,2*n,
    a3=2*a3*a0+2*a1*a2;
    a2=2*a2*a0+a1^2;
    a1=2*a1*a0+1;
    a0=a0^2+C0;
  );
  /* assume zn=0, so z2n=0, so a0=0, */
  /* divide by x, solve quadratic root closest to zero */
  a0=a1;
  a1=a2;
  a2=a3;
  y=a1^2-4*a0*a2;
  z1=(-a1+sqrt(y))/(2*a2);
  z2=(-a1-sqrt(y))/(2*a2);
  if (abs(z1)>abs(z2),z1=z2);
  return(C0+z1);
}

shortperiodn(n,nc,cnt,l) = {
  local(z,y,i,j);
  z=0;for (i=1,n,z=truncs(z^2-0+x+nc,cnt+1));
  v=polroots(z);
  for (i=1,cnt,
    z=v[i]+nc;
    y=0;
    for (j=1,n,y=y^2+z);
    y=abs(y);
    if (y<0.04,
      z=invzero(z,n);
      if (abs(imag(z))<1E-21, z=real(z));
      if (imag(z)<0,z=conj(z));
      if (((real(z)<>0) || (imag(z)<>0)) && (imag(z)>=0),
        j=1;
        while (abs(oj[j]-z)>1E-21 && (oj[j]<>0), j++);
        if (oj[j]==0,
          if (estimz2(z,n,l)<>0, oj[j]=z; if (j==15, oj[j]=v[i]+nc));
        );
      );
    );
  )
}
randomperiodn(n,l,j)= {
  local(i,y,oout);
  if (l==0,l=0.01);
  if (j==0, j=50);
  oj=vector(128);
  print("randomperiodn("n");");
  for (i=1,j,
    shortperiodn(n,cfromz(goldr^i),16,l);    /* main cardiod */
    shortperiodn(n,-1+0.25*(goldr^i),16,l);  /* 2n bulb */
    shortperiodn(n,-2+fwidth*((gold*i)%1),16,l); /* tip */
  );
  i=1;while (oj[i]<>0,i++);
  oout=vector(i-1);
  i=1;while (oj[i]<>0,oout[i]=oj[i];i++);
  return(oout);
}
/* truncate series */
truncs(z,n)={
  local(y,i);
  y=0;
  if (n==0,n=8);
  for (i=0,n-1,y=y+polcoeff(z,i)*x^i);
  return(y);
}
gfunc(C0,n)={local(i,z);z=0;for (i=1,n,z=z^2+C0);z}
cfromz(z) = subst(x/2-x^2/4,x,z);
lfromc(C0) = {1/2-sqrt(1/4-C0);}
gold=2/(sqrt(5)+1);
goldr=conj(exp(2*Pi*I*gold));
/* feigenbaum point+2, and feigenbaum ratio */
fwidth= 0.598844810907949399476173212106138707773702;
feign = 4.669201609102990671853203820466201617258186;

/* slope at the fixed point of C0, using n iterations */
slopeh(C0,n)={
  local(l,l0,l1,i,j,co,s);
  s=1;
  co=0;
  while (abs(s)>plim, /* one past plim */
    l0=co;l1=1; /* l1=slope */
    for (i=1,n,l0=l0^2+C0;l1=l1*2*l0);
    s=l0-co; /* as co approaches fixed point, (l0-co) approaches zero */
    co=co-(l0-co)/(l1-1); /* co=updated approximation for fixed point */
  );
  l0=co;l1=1; /* one more iteration to improve l1 slope approximation */
  for (i=1,n,l0=l0^2+C0;l1=l1*2*l0);
  return(l1);
}
invslopeh(z,C0,n,est) = {
  local (y,s,slop,lastyz,curyz,lest,ly,pgoal);
  lastyz=100;
  y=slopeh(C0+est,n);
  curyz=abs(y-z);
  lest=est+curyz*jm;
  ly=slopeh(C0+lest,n);
  pgoal=10^(-precis/1.3);
  /* generate the fixed point for pentation by iteration slog */
  s=1;
  while ((curyz>pgoal) && ((curyz<lastyz) || (s<3)),
    est=precision(est,precis);
    y=precision(y,precis);
    slop=(y-ly)/(est-lest);
    lest=est;
    ly=y;
    est=est+(z-y)/slop;
    lastyz=curyz;
    y=slopeh(C0+est,n);
    curyz=abs(y-z);
    s++;
  );
  if (curyz>0.1, print (curyz " bad result, need better initial est"));
  return(est);
}
/* default use r=0.5, samples=50 */
islopetaylor(C0,n,r,samples) = {
  local(rinv,s,t,x1,y,z,tot,t_est,tcrc,halfsamples,wtaylor,terms);
  if (n==0,
    n=1;
    z=C0;
    while (abs(z)>plim,z=z^2+C0;n++);
  );
  if (r==0,r=0.5);
  C0=invzero(C0,n);
  if (samples==0, samples=50); 
  halfsamples=samples/2;
  terms = floor(samples*0.51);
  t_est    = vector (samples,i,0);
  tcrc     = vector (samples,i,0);
  if (r==0,r=1);
  rinv = 1/r;
  wtaylor=C0;
  for(s=1, samples, x1=-1/(samples)+(s/halfsamples); tcrc[s]=exp(Pi*I*x1); );

  for (t=1,samples, t_est[t] = invslopeh(r*tcrc[t],C0,n); );

  for (s=0,terms-1,
    tot=0;
    for (t=1,samples,
      tot=tot+t_est[t];
      t_est[t]=t_est[t]*conj(tcrc[t]);
    );
    tot=tot/samples;
    if (s>=1, tot=tot*(rinv)^s);
    wtaylor=wtaylor+tot*x^s;
  );
  wtaylor=precision(wtaylor,precis);
  if (imag(C0)==0, wtaylor=real(wtaylor));
  return(wtaylor);
}
prtpoly(wtaylor,t) = {
  local(s,z,iprt);
  if (t==0,t=7);
  z=polcoeff(wtaylor,0);
  if (imag(z)<>0,iprt=1,iprt=0);
  print1 ("{z=   ");
  if (iprt,print1("("));
  if (real(z)<0, print1(z), print1(" " z));
  if (iprt,print(")"),print() );
  for (s=1,t-1,
    z=polcoeff(wtaylor,s);
    if (s>9, print1("+x^" s), print1("+x^ " s));
    if (iprt,print1("*("),print1("*") );
    if (real(z)<0, print1(z), print1(" " z));
    if (iprt,print(")"),print() );
  );
  print("}");
}

plotsetup(low,high,numb,radius,randn) = {
local(n,i,zeros,zeron);
  zp=0;
  if (numb==0,numb=85);
  zeros=vector(numb);zeron=vector(numb);zerov=vector(numb);
  if (low==0,low=1);
  if (high==0, high=14);
  for (n=low,high,
    if (n>7,C=randomperiodn(n,radius,randn),C=periodnzero(n,radius));
    for (i=1,length(C),
      zp++;
      zeros[zp]=C[i];
      zeron[zp]=n;
    );
  );
  ip=0;
  print(zp);
  for (i=1,zp,
    ip++;
    zerov[ip]=truncs(islopetaylor(zeros[i],zeron[i]),7);
    if (imag(zeros[i])<>0,
      ip++;
      zerov[ip]=conj(zerov[ip-1]);
    );
  );
}
evaln(z)={local(y,s,i);
  y=vector(ip*2);
  for (i=1,ip,
    s=subst(zerov[i],x,z);
    y[i*2-1]=real(s);
    y[i*2]=imag(s);
  );
  return(y)
}
maketheplot(w) = local(t,z); { ploth(t=0,2*Pi,z=exp(t*I);evaln(z),1); }

Mandelbrot plot using bulbs/cardioids from period 1..20, bigger than 0.004. I truncated the Taylor series to seven terms, and plotted 215 different bulbs and cardioids; 112 are unique if not counting complex conjugates twice. For the randomperiodn, I used (n,0.004,350). It took 350 random tries to get both of the period(15) sub-bulbs on the period(5) bulb; some bulbs are still missing. Directly searching for bulbs near $\frac{k\pi i}{n}$ and other algorithms would probably work better; I'm pretty new at this.

Solution 3:

The algebraic solutions one can find : John Stephenson : "Formulae for cycles in the Mandelbrot set", Physica A 177, 416-420 (1991); "Formulae for cycles in the Mandelbrot set II", Physica A 190, 104-116 (1992); "Formulae for cycles in the Mandelbrot set III", Physica A 190, 117-129 (1992)

Example image using this technique is here :