Fractals/mandel

< Fractals

Mandel[1] : software for real and complex dynamics by Wolf Jung. Here one can find unofficial docs about it.

Features

maps or functions

Families of functions (= maps ):

These famnilies are defined by classes

complex quadratic polynomial

The monic and centered form of complex quadratic polynomial :

" The Mandelbrot set is based on the one-parameter family of quadratic polynomials fc(z) = z2 + c. " ( from help )

// from  mndynamo.cpp  by Wolf Jung (C) 2007-2014.
void mndynamics::root(double x, double y, double &u, double &v)
{  v = sqrt(x*x + y*y);
   if (x > 0.0) { u = sqrt(0.5*(v + x)); v = 0.5*y/u; return; }
   if (x < 0.0)
   { v = sqrt(0.5*(v - x)); if (y < 0.0) v = -v; u = 0.5*y/v; return; }
   if (y >= 0.0) { u = sqrt(0.5*y); v = u; return; }
   u = sqrt(-0.5*y); v = -u;
}




// from mndlbrot.cpp  by Wolf Jung (C) 2007-2014. part of Mandel 5.12 
// input : c = a+b*i and z = x+y*i
// output : z = x+y*i

void mndlbrot::f(double a, double b, double &x, double &y) const
{  double u = x*x - y*y + a; // u is temporary variablke
   y = 2*x*y + b; 
   x = u; }

It can be used for :

// from mndlbrot.cpp  by Wolf Jung (C) 2007-2014. part of Mandel 5.12 
// input : c, z , mode
// c = A+B*i where A and B are global variables defined in mndynamo.h
// z = x+y*i
// 
// output : z = x+y*i
void mndlbrot::iterate(double &x, double &y, int mode) const // = 0
{  // Mapping f = forward itearation = f(z) = key f = mode 0
   if (!mode) f(A, B, x, y);
   if (mode > 0) { x = 0; y = 0; }
   if (mode >= 0 || mode < -2) return;

   // f^{-1}(z) = inverse with principal value
   if (A*A + B*B < 1e-20) 
   {  root(x - A, y - B, x, y); // 2-nd inverse function = key b 
      if (mode & 1) { x = -x; y = -y; } // 1-st inverse function = key a
      return;
   }

   //f^{-1}(z) =  inverse with argument adjusted
   // "When you write the real and imaginary parts in the formulas as complex numbers again,
   //    you see that it is sqrt( -c / |c|^2 )  *  sqrt( |c|^2 - conj(c)*z ) ,
   //  so this is just sqrt( z - c )  except for the overall sign:
   // the standard square-root takes values in the right halfplane,  but this is rotated by the squareroot of -c .
   // The new line between the two planes has half the argument of -c .
   // (It is not orthogonal to c ...  )" W Jung "
   double u, v, w = A*A + B*B; 
   root(-A/w, -B/w, u, v); 
   root(w - A*x - B*y, B*x - A*y, x, y);
   w = u*x - v*y; 
   y = u*y + v*x; 
   x = w;
   // 2-nd inverse function = key b 
   if (mode & 1) // mode = -1
       { x = -x; y = -y; } // 1-st inverse function = key a
   
}

It is called by : 

  // forward iteration
   fAct = new QAction(trUtf8("Mapping &f"), mapAG);
   fAct->setShortcut(QString("f"));

   // inverse iteration
   inv1Act = new QAction(trUtf8("&1st inverse"), mapAG);
   inv1Act->setShortcut(QString("a"));
   inv2Act = new QAction(trUtf8("&2nd inverse"), mapAG);
   inv2Act->setShortcut(QString("b"));

   fAct->setEnabled(ftype != 48);
   inv1Act->setEnabled(!ftype || (ftype == 28 && signtype*signtype == 4) || ftype == 38 || ftype == 48 || ftype == 58);
   inv2Act->setEnabled(!ftype || (ftype == 28 && signtype*signtype == 4) || ftype == 48 || ftype == 58);

void QmnShell::map(QAction *act)
{  double x, y; 
   dplane->getPoint(x, y);
   if (act == fAct) f->iterate(x, y); // mode = 0 ??
   if (act == inv1Act)  { f->iterate(x, y, -1); /*if (!ftype)*/ dplane->Move(12, f); }
   if (act == inv2Act)    { f->iterate(x, y, -2); /*if (!ftype)*/ dplane->Move(13, f); }
   dplane->setPoint(x, y);
}

Algorithms


Example programs which use code from Mandel

Misiurewicz point

"... the usual convention is to use the preperiod of the critical value. This has the advantage, that the angles of the critical value have the same preperiod under doubling as the point, and the same angles are found in the parameter plane." ( Wolf Jung )


Finding principal Misiurewicz points of the wake k/r of main cardioid


int_angle	prep,period	c							e_angles_of_the_wake 		                                                                                 e_engles_of_Misiurewicz	

1/2		2,1		c = -1.543689012692076  +0.000000000000000 i    	( 1/3  or  p01  and  2/3  or  p10 )							

1/3		3,1		c = -0.101096363845622  +0.956286510809142 i    	( 1/7  or  p001  and 2/7  or  p010 )										 9/56, 11/56 and 15/56

1/4		4,1		c = 0.366362983422764  +0.591533773261445 i    

1/5		5,1		c = 0.437924241359463  +0.341892084338116 i  

1/6		6,1		c = 0.424512719050040  +0.207530228166745 i    

1/7		7,1		c = 0.397391822296541  +0.133511204871878 i 

1/8		8,1		c = 0.372137705449577  +0.090398233158173 i 
	
1/9		9,1		c = 0.351423759052522  +0.063866559813293 i    		(1/511= p000000001 ; 2/511  = p000000010  )

1/10		10,1		c = 0.334957506651529  +0.046732666062027 i    		(1/1023  or  p0000000001  and 2/1023  or  p0000000010 )

--------------------------------------------------------------------------------------------------------------------------------------------

1/11		11,1		c = 0.321911396847221  +0.035204463294452 i    		(1/2047  or  p00000000001  and 2/2047  or  p00000000010 )

1/12		12,1		c = 0.311507660281508  +0.027173719501342 i    		(1/4095  or  p000000000001  and 2/4095  or  p000000000010 )

1/13		13,1		c = 0.303127979909719  +0.021411628038965 i    		(1/8191  or  p0000000000001  and 2/8191  or  p0000000000010 )   

1/14		14,1		c = 0.296304475349759  +0.017171379707062 i    		(1/16383  or  p00000000000001  and 2/16383  or  p00000000000010 )

1/15		15,1		c = 0.290687752431041  +0.013982147106557 i    		(1/32767  or  p000000000000001  and 2/32767  or  p000000000000010 )

1/16		16,1		c = 0.286016666695566  +0.011537401448441 i    		(1/65535  or  p0000000000000001  and 2/65535  or  p0000000000000010 )        

1/17		17,1		c = 0.282094678248954  +0.009631861589570 i    		(1/131071  or  p00000000000000001  and 2/131071  or  p00000000000000010 )

1/18		18,1		c = 0.278772459129383  +0.008124579648410 i    		( 1/262143  or  p000000000000000001  and 2/262143  or  p000000000000000010 )

1/19		19,1		c = 0.275935362441682  +0.006916613801737 i    		(1/524287  or  p0000000000000000001  and 2/524287  or  p0000000000000000010 )

1/20		20,1		c = 0.273494431676492  +0.005937124623590 i    		( 1/1048575  or  p00000000000000000001  and 2/1048575  or  p00000000000000000010 ) 
-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------

1/21		21,1		c = 0.271379927384804  +0.005134487480794 i   		( 1/2097151  or  p000000000000000000001  and 2/2097151  or  p000000000000000000010 )

1/22		22,1		c = 0.269536632623270  +0.004470475898698 i    		( 1/4194303  or  p0000000000000000000001  and 2/4194303  or  p0000000000000000000010 )

1/23		23,1		c = 0.267920417711850  +0.003916372840385 i    		( 1/8388607  or  p00000000000000000000001  and 2/8388607  or  p00000000000000000000010 )

1/24		24,1		c = 0.266495701963622  +0.003450320181976 i    		( 1/16777215  or  p000000000000000000000001  and 2/16777215  or  p000000000000000000000010 )

1/25		25,1		c = 0.265233559524147  +0.003055480072447 i    		( 1/33554431  or  p0000000000000000000000001  and 2/33554431  or  p0000000000000000000000010 )  

1/26		26,1		c = 0.264110291981537  +0.002718738820085 i    		( 1/67108863  or  p00000000000000000000000001  and 2/67108863  or  p00000000000000000000000010 )

1/27  		27,1		c = 0.263106342463072  +0.002429779655182 i    		( 1/134217727  or  p000000000000000000000000001  and 2/134217727  or  p000000000000000000000000010 )

1/28		28,1		c = 0.262205461953178  +0.002180410330127 i  		( 1/268435455  or  p0000000000000000000000000001  and 2/268435455  or  p0000000000000000000000000010 )

1/29		29,1		c = 0.261394063652992  +0.001964069379345 i  		( 1/536870911  or  p00000000000000000000000000001  and 2/536870911  or  p00000000000000000000000000010 )

1/30		30,1		c = 0.260660718810682  +0.001775459345952 i    		( 1/1073741823  or  p000000000000000000000000000001  and 2/1073741823  or  p000000000000000000000000000010 )

----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

1/35		35,1		c = 0.257872123091544  +0.001121742345611 i    		( 1/34359738367  or  p00000000000000000000000000000000001  and 2/34359738367  or  p00000000000000000000000000000000010 )

------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
1/62 		62,1		c = 0.252537923584907  +0.000203841219482 i    		( 1/4611686018427387903 = p00000000000000000000000000000000000000000000000000000000000001 ;
                                                                                          2/4611686018427387903 = p00000000000000000000000000000000000000000000000000000000000010 )

1/63		63,1		c = 0.252458520819920  +0.000194332543868 i     	( 1/9223372036854775807 = p000000000000000000000000000000000000000000000000000000000000001  ;
                                                                                          2/9223372036854775807 = p000000000000000000000000000000000000000000000000000000000000010 )

1/64 		64,1		c = 0.252382784694904  +0.000185406363701 i    		( 1/18446744073709551615 = p0000000000000000000000000000000000000000000000000000000000000001 ; 
                                                                                          2/18446744073709551615 = p0000000000000000000000000000000000000000000000000000000000000010 )

/* 
from mndlbrot.cpp  part of Mandel 5.12 by Wolf Jung (C) 2007-2014.  
 the GNU General Public License as published by the Free Software Foundation; either
   version 3, or (at your option) any later version.

*/
int mndlbrot::find(int sg, uint preper, uint per, double &x, double &y) const
{  double a = A, b = B, fx, fy, px, py, w; uint i, j;
   for (i = 0; i < 30; i++)
   {  if (sg > 0) { a = x; b = y; }
      if (!preper)
      {  if (sg > 0) { fx = 0; fy = 0; px = 0; py = 0; }
         else { fx = -x; fy = -y; px = -1; py = 0; }
      }
      else
      {  fx = x; fy = y; px = 1.0; py = 0;
         for (j = 1; j < preper; j++)
         {  if (px*px + py*py > 1e100) return 1;
            w = 2*(fx*px - fy*py); py = 2*(fx*py + fy*px);
            px = w; if (sg > 0) px++;
            w = fx*fx - fy*fy + a; fy = 2*fx*fy + b; fx = w;
         }
      }
      double Fx = fx, Fy = fy, Px = px, Py = py;
      for (j = 0; j < per; j++)
      {  if (px*px + py*py > 1e100) return 2;
         w = 2*(fx*px - fy*py); py = 2*(fx*py + fy*px);
         px = w; if (sg > 0) px++;
         w = fx*fx - fy*fy + a; fy = 2*fx*fy + b; fx = w;
      }
      fx += Fx; fy += Fy; px += Px; py += Py;
      w = px*px + py*py; if (w < 1e-100) return -1;
      x -= (fx*px + fy*py)/w; y += (fx*py - fy*px)/w;
   }
   return 0;
}

critical orbit

If you want to draw critial orbit for more iterations then Mandel has as a standard

Functions :

algorithm 9 : zeros of qn(c)

Parameter plane colored according to sign of 5-th iteration. It uses algorithm 9, order n=5.
Dynamic plane colored according to sign of 8-th iteration. It uses algorithm 9, order q=8.

"It is meant to show the location of precritical points, i.e., the preimages of 0. So it shows in four different colors, in which quadrant f^{n-1}(z) is. This means that the real axis is pulled back n times, so it gives an approximation to the checkerboard when the parameter is real." ( Wolf Jung ) It can be used for making parabolic checkerboard ( chessboard)


This is radial 4th-decomposition of plane ( compare it with n-th decomposition of LSM/M )

4 colors are used because there are 4 quadrants :


"... when the four colors are meeting somewhere, you have a zero of q_n(c), i.e., a center of period dividing n. In addition, the light or dark color shows if c is in M." ( Wolf Jung ) See function quadrantqn from class mndlbrot in mndlbrot.cpp file [6]

// mndlbrot.cpp  by Wolf Jung (C) 2007-2014.  part of Mandel 5.12 
int mndlbrot::quadrant(double a, double b, double x, double y)
{  //exterior checked in Period (dist)
   int cl = 1, j;
   for (j = 1; j < subtype; j++)
   { f(a, b, x, y); if (x*x + y*y > 1e100) return 0; }
   if (x > 0) cl = 3; if (y < 0) cl++;
   if (Period) cl |= 8; return cl;
}//quadrant

int mndlbrot::quadrantqn(double a, double b)
{  //exterior checked in Period (dist)
   double x = a, y = b; int cl = 1, j;
   for (j = 1; j < subtype; j++)
   { f(a, b, x, y); if (x*x + y*y > 1e100) return 0; }
   if (x > 0) cl = 3; if (y < 0) cl++;
   if (Period) cl |= 8; return cl;
}//quadrantqn

int mndlbrot::quadrantqnp(double a, double b)
{  //exterior checked in Period (dist)
   int cl = 1, j;
   mndplex C(a, b), Q = C, P = 1.0;
   for (j = 1; j < subtype; j++)
   {  P = 2.0*Q*P + 1.0; Q = Q*Q + C;
      if (norm(Q) > 1e100 || norm (P) > 1e100) return 0;
   }
   if (real(P) > 0) cl = 3; if (imag(P) < 0) cl++;
   if (Period) cl |= 8; return cl;
} //quadrantqnp

// mndlbrot.cpp  by Wolf Jung (C) 2007-2014.  part of Mandel 5.12,
uint mndlbrot::pixcolor(double x, double y)
{  uint j, cl = 0; double a, b;
  
   if (sign > 0) { a = x; b = y; } //z = c, not z = 0 !
     else { a = A; b = B; }
   if (drawmode == 9 && sign < 0) return quadrant(a, b, x, y);
   if (drawmode == 9)             return quadrantqn(a, b);
   if (drawmode == 10)            return quadrantqnp(a, b);


The differences from standard loop are :

Compare :

Here is fragment of c code for computing 8-bit color for point c = cx + cy * i  :

/* initial value of orbit = critical point Z= 0 */
                       Zx=0.0;
                       Zy=0.0;
                       Zx2=Zx*Zx;
                       Zy2=Zy*Zy;
                       /* the same number of iterations for all points !!! */
                       for (Iteration=0;Iteration<IterationMax; Iteration++)
                       {
                           Zy=2*Zx*Zy + Cy;
                           Zx=Zx2-Zy2 +Cx;
                           Zx2=Zx*Zx;
                           Zy2=Zy*Zy;
                       };
                       /* compute  pixel color (8 bit = 1 byte) */
                       if ((Zx>0) && (Zy>0)) color[0]=0;   /* 1-st quadrant */
                       if ((Zx<0) && (Zy>0)) color[0]=10; /* 2-nd quadrant */
                       if ((Zx<0) && (Zy<0)) color[0]=20; /* 3-rd quadrant */
                       if ((Zx>0) && (Zy<0)) color[0]=30; /* 4-th quadrant */

Tiling of parabolic Julia set interior

Tiling of parabolic Julia set interior ( made with plain C program )

It works for interior of parabolic Julia sets, only for cases when :

It makes parabolic checkerboard

// from mndlbrot.cpp  by Wolf Jung (C) 2007-2015.

uint mndlbrot::parabolic(double x, double y)
{  uint j; double u, v;


   for (j = 1; j <= maxiter; j++)
   {  u = x*x; v = y*y;
      if (u + v <= bailout)
      { y = 2*x*y; x = u - v + A; }
      else return j;
      if (A > 0 && x >= 0 && x <= 0.5 && (y > 0 ? y : -y) <= 0.5 - x)
         return (y > 0 ? 65281u : 65289u);
      if (A < 0 && x >= -0.5 && x <= 0 && (y > 0 ? y : -y) <= 0.3 + 0.6*x)
      {  if (j & 1) return (y > 0 ? 65282u : 65290u);
         else return (y > 0 ? 65281u : 65289u);
      }
      
      // c para 1/3 ; not working now 
      // if (x > -0.25 && y > 0 && x + y < 0.183012701892219)
      // {  j %= 3; if (!j) return 65290;
      //      else if (j & 1) return 65291; else return 65289;} 
      //
   }

   return 65293u;
}

External rays

Methods :

// qmnshell.cpp  by Wolf Jung (C) 2007-2014. Mandel 5.12
void QmnShell::setRay(QAction *act)
{  int error = 0;
   QString enterAngleString = trUtf8(
   "<b></b>Enter the angle &theta; "
   "with 0 &le; &theta; &lt; 1. Notation :<br>"
   "\"1/7\" or \"p001\" for 1/7 = "
   ".<nobr style=\"text-decoration:overline\">001</nobr> "
   "(periodic angles),<br>"
   "\"9/56\" or \"001p010\" for 9/56 = "
   ".001<nobr style=\"text-decoration:overline\">010</nobr> "
   "(preperiodic angles),<br>"
   "\"1/4\" or \"01\" for 1/4 = .01 (dyadic angles).");
   QString pullBackString = trUtf8(
   "<br>Then hit Return repeatedly to perform the iteration<br>"
   "of the Thurston Algorithm. A connecting path between<br>"
   "point configurations is used instead of spider legs.<br>"
   "Use Home or enter 0 to quit.");
   if (act == greenAct)
   {  double x, y; theplane->getPoint(x, y); y = f->green(signtype, x, y);
      QmnDoubleDialog *dialog = new QmnDoubleDialog(trUtf8(
         "To draw an equipotential line, enter the\n"
         "potential  log|\xce\xa6(%1)| > 0 :\n"
         "The suggested value is the potential of the\n"
         "current point. Just hit Return to accept it."
         ).arg((signtype > 0 ? QChar('c') : QChar('z'))), &y, 0, 3, this);
      connect(dialog, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
      if (!dialog->exec() || y <= 0 || y > 10.0) return;
      if (signtype > 0) { drawLater = 0; dplane->stop(); }
      else pplane->stop();
      theplane->green(f, signtype, y, 10);
   }
   if (act == rayAct)
   {  uint method, q = 0; double x, y; f->getParameter(x, y);
      qulonglong N1 = 1LL, N2 = 2LL;
      if (signtype < 0) { N1 = 0LL; if (dplane->getType() < 0) N2 = 1LL; }
      if (ftype == 1) { N1 = 2LL; N2 = 2LL; }
      if (ftype == 28) { N1 = 1LL; N2 = 1LL; }
      QmnRayDialog *dialog = new QmnRayDialog(enterAngleString, trUtf8(
         "Adjust the quality, 2...10. It is the number of intermediate\n"
         "points, or it concerns the distance to the starting point."),
         &N1, &N2, &method, &q, 3, this);
      connect(dialog, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
      if (!dialog->exec()) return; if (ftype == 28) method = 5;
      if (!method && dplane->backRay(N1, N2, x, y, q, Qt::white, 1))
         method = 1;
      if (method & 1)
         theplane->newtonRay(signtype, N1, N2, x, y, q, Qt::white, method);
      if (method & 2)
      {  if (signtype > 0) { drawLater = 0; dplane->stop(); }
         else pplane->stop();
         theplane->traceRay(signtype, double(N1)/double(N2), f, x, y, q);
      }
   }
   if (act == rayPointAct)
   {  double x, y; f->getParameter(x, y);
      qulonglong N1 = 0LL, N2;
      QmnCombiDialog *dialog = new QmnCombiDialog(enterAngleString,
         &N1, &N2, 3, this);
      connect(dialog, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
      if (!dialog->exec()) return;
      int k, p = mndAngle::normalize(N1, N2, k), q = 0; qulonglong n1 = N1;
      if (!p)
      {  QmnUIntDialog *dialog1 = new QmnUIntDialog(trUtf8(
            "The angle  %1/%2\n"
            "has  preperiod + period > 64.").arg(N1).arg(N2),
            0, 0, 0u, 3, this, 0);
         connect(dialog1,SIGNAL(needHelp(int)), helptabs,SLOT(showPage(int)));
         dialog1->exec(); return;
      }
      QString bin; QmnCombiDialog::numbersToBinary(N1, N2, bin); double l;
      if (p == 1) l = mndAngle::lambda(N1, N2, 1.0e-10, 1000);
      else
      {  qulonglong L = 1ULL; L <<= 60;
         l = ((double)(L))*((double)(N1))/((double)(N2));
         l = -mndAngle::lambda(((qulonglong)(l)), L, 1.0e-10, 1000);
      }
      QString text = trUtf8(
         "The angle  %1/%2  or  %3\n"
         "has  preperiod = %4  and  period = %5.\n"
         ).arg(N1).arg(N2).arg(bin).arg(k).arg(p);
      if (l > 0.0 && signtype > 0) text += trUtf8(
        "Entropy: e^h = 2^B = \xce\xbb = %1\n").arg(
           QString::number(l, 'f', 8));
      if (k && signtype < 0) text += trUtf8(
         "The corresponding dynamic ray is landing\n"
         "at a preperiodic point of preperiod %1 and\n"
         "period dividing %2.\n"
         "Do you want to draw the ray and to shift z\n"
         "to the landing point?").arg(k).arg(p);
      if (k && signtype > 0) text += trUtf8(
         "The corresponding parameter ray is landing\n"
         "at a Misiurewicz point of preperiod %1 and\n"
         "period dividing %2.\n"
         "Do you want to draw the ray and to shift c\n"
         "to the landing point?").arg(k).arg(p);
      if (!k && signtype < 0) text += trUtf8(
         "The dynamic ray is landing at a repelling\n"
         "or parabolic point of period dividing %1.\n"
         "Do you want to draw the ray and to shift\n"
         "z to the landing point?").arg(p);
      if (!k && signtype > 0)
      {  q = mndAngle::conjugate(n1, N2);
         QmnCombiDialog::numbersToBinary(n1, N2, bin); text += QString(trUtf8(
         "The conjugate angle is  %1/%2  or  %3 .\n"
         )).arg(n1).arg(N2).arg(bin);
         if (q < p) bin = trUtf8("satellite"); else bin = trUtf8("primitive");
         QString t1, t2; mndCombi c; c.fromAngle(N1, N2); qulonglong b;
         c.getKneading(b); QmnCombiDialog::codeToKneading(b, t1);
         c.getAddress(b); QmnCombiDialog::codeToAddress(b, t2);
         text += trUtf8(
            "The kneading sequence is  %1  and\n"
            "the internal address is  %2 .\n").arg(t1).arg(t2);
         text += trUtf8(
            "The corresponding parameter rays are landing\n"
            "at the root of a %1 component of period %2.\n").arg(bin).arg(p);
         if (q && q < p) text += trUtf8(
            "It is bifurcating from period %1.\n").arg(q);
         text += trUtf8(
            "Do you want to draw the rays and to shift c\n"
            "to the corresponding center?");
      }
      QmnUIntDialog *dialog1 = new QmnUIntDialog(text, 0, 0, 0u, 3, this);
      connect(dialog1, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
      if (!dialog1->exec()) return;
      if (ftype == 58)
      {  x = ((double)(N1))/((double)(N2)); if (l < 0.0) l = -l;
         uint n; pplane->getNmax(n);
         pplane->setPoint(x, l*0.01*n); return;
      }
      if (q) pplane->newtonRay(1, n1, N2, x, y, 5, Qt::white, 1);
      if (signtype < 0) q = dplane->backRay(N1, N2, x, y, 5, Qt::white, 2);
      else q = pplane->newtonRay(1, N1, N2, x, y, 5, Qt::white, 2);
      if (q <= 1 && !f->find(signtype, k, p, x, y)) theplane->setPoint(x, y);
   }
   if (act == portraitAct)
   {  double x, y; f->getParameter(x, y);
      qulonglong n1, N1 = 0ULL, N2;
      QmnCombiDialog *dialog = new QmnCombiDialog(trUtf8(
         "Enter the characteristic angle \xce\xb8 with 0 < \xce\xb8 < 1,\n"
         "2 \xe2\x89\xa4 period \xe2\x89\xa4 64. "
         "Notation : \"2/7\" or \"p010\" ."),
         &N1, &N2, 3, this);
      connect(dialog, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
      if (!dialog->exec()) return;
      int k, p = mndAngle::normalize(N1, N2, k); if (k || p <= 1) return;
      n1 = N1; int q = mndAngle::conjugate(n1, N2);
      if (signtype < 0)
      {  for (k = 0; k < p; k++)
         { dplane->backRay(N1, N2, x, y); mndAngle::twice(N1, N2); }
         if (q == p) for (k = 0; k < p; k++)
         { dplane->backRay(n1, N2, x, y); mndAngle::twice(n1, N2); }
         return;
      }
      bool sphere = (dplane->getType() < 0); if (sphere) dplane->setType(0);
      dplane->setPlane(0, 0, 1.2, 0); dplane->draw(f, 0, themode);
      dplane->drawEllipse(0, 0, 1.0, 1.0, Qt::green); x = double(N2);
      for (k = 0; k < p; k++)
      {  dplane->drawOrtho(double(N1)/x, double(n1)/x);
         mndAngle::twice(N1, N2); mndAngle::twice(n1, N2);
      }
      updateRegion = true; drawLater = 0; dplane->setPoint(0, 0);
      if (sphere) dplane->setType(-1); else dplane->setPlane(0, 0, 2.0, 0);
   }
   if (act == laminatAct)
   {  double x, y, u, v; dplane->getPoint(x, y);
      qulonglong N = 0ULL, N1 = 0ULL, N2;
      QmnCombiDialog *dialog = new QmnCombiDialog(enterAngleString,
         &N1, &N2, 3, this);
      connect(dialog, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
      if (!dialog->exec() || !N1) return;
      uint n; int k, p = mndAngle::normalize(N1, N2, k); u = (double)(N2);
      if (!k && p)
      { N = N1; mndAngle::conjugate(N, N2); v = ((double)(N))/u; }
      u = ((double)(N1))/u; dplane->setPoint(0.5*u, 0.0);
      if (N && (N1 & 1)) dplane->setPoint(0.5*v, 0.0);
      //problem 1/7 avoided, why? new problem 41/127, need rational angles?
      pplane->getNmax(n); if (n > 20u) n = 12u;
      if (signtype < 0) //no crash when n + k + p ~ 64
      {  dplane->setPoint(x, y); f->getParameter(x, y); qulonglong K, D = 1ULL;
	 for (k = 0; k <= n; k++)
         {  for (K = 0ULL; K < D; K++)
	    {  dplane->backRay(N1 + K*N2, D*N2, x, y);
	       if (N) dplane->backRay(N + K*N2, D*N2, x, y);
	    }
	    D <<= 1;
	 }
         return;
      }
      bool sphere = (dplane->getType() < 0); if (sphere) dplane->setType(0);
      dplane->setPlane(0, 0, 1.2, 0); dplane->draw(f, 0, themode);
      dplane->drawEllipse(0, 0, 1.0, 1.0, Qt::green);
      if (N)
      {  dplane->drawOrtho(u, v);
	 dplane->drawLamination(0.5*u, 0.5*v + 0.5, n);
	 dplane->drawLamination(0.5*v, 0.5*u + 0.5, n);
      }
      else if (N2 == 56ULL && (N1 == 9ULL || N1 == 11ULL || N1 == 15ULL))
      {  dplane->drawOrtho(9.0/56.0, 11.0/56.0);
         dplane->drawOrtho(11.0/56.0, 15.0/56.0);
         dplane->drawOrtho(9.0/56.0, 15.0/56.0);
	 dplane->drawLamination(9.0/112.0, 11.0/112.0, n);
	 dplane->drawLamination(11.0/112.0, 15.0/112.0, n);
	 dplane->drawLamination(15.0/112.0, 65.0/112.0, n);
	 dplane->drawLamination(65.0/112.0, 67.0/112.0, n);
	 dplane->drawLamination(67.0/112.0, 71.0/112.0, n);
	 dplane->drawLamination(71.0/112.0, 9.0/112.0, n);
      }
      else dplane->drawLamination(0.5*u, 0.5*u + 0.5, n);
      updateRegion = true; drawLater = 0; dplane->setPoint(0, 0);
      if (sphere) dplane->setType(-1); else dplane->setPlane(0, 0, 2.0, 0);
   }
   if (act == wakeAct)
   {  uint k = 1, r = 2;
      QmnUIntDialog *dialog = new QmnUIntDialog(trUtf8(
         "Determine the wake of a limb at the main cardioid.\n"
         "Enter a fraction  k/r  for the rotation number,\n"
         "in lowest terms, with  1 \xe2\x89\xa4 k < r \xe2\x89\xa4 64 :"),
         &k, &r, 65001u, 3, this);
      connect(dialog, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
      if (!dialog->exec()) return;
      qulonglong n, d = mndAngle::wake(((int) (k)), ((int) (r)), n);
      if (!d) return;
      QString b1, b2; QmnCombiDialog::numbersToBinary(n, d, b1);
      QmnCombiDialog::numbersToBinary(n + 1LL, d, b2);
      QmnUIntDialog *dialog1 = new QmnUIntDialog(trUtf8(
         "The %1/%2-wake of the main cardioid is\n"
         "bounded by the parameter rays with the angles\n"
         "%3/%4  or  %5  and\n"
         "%6/%4  or  %7 .\n"
         "Do you want to draw the rays and to shift c\n"
         "to the center of the satellite component?").arg(k).arg(r).arg(
         n).arg(d).arg(b1).arg(n + 1LL).arg(b2), 0, 0, 0u, 3, this);
      connect(dialog1, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
      if (!dialog1->exec()) return;
      if (ftype == 58)
      {  uint nn; pplane->getNmax(nn);
         pplane->setPoint(((double)(n))/((double)(d)), 0.01*nn); return;
      }
      double x, y; pplane->newtonRay(1, n, d, x, y, 5, Qt::white, 1); n++;
      if (pplane->newtonRay(1, n, d, x, y, 5, Qt::white, 2) <= 1
          && !f->find(1, 0, r, x, y)) pplane->setPoint(x, y);
   }
   if (act == kneadAct)
   {  qulonglong N1 = 1LL, N2, n1, n2, d;
      QmnCombiDialog *dialog = new QmnCombiDialog(trUtf8(
         "Enter a *-periodic kneading sequence,\n"
         "e.g, \"ABAAB*\" or \"10110*\".\n"
         "Or enter an internal address,\n"
         "e.g., \"1-2-4-5-6\".\n"
         "(The periods \xe2\x89\xa4 64 are increasing.)"), &N1, &N2, 3, this);
      connect(dialog, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
      if (!dialog->exec()) return;
      int p, q; QString t1, t2, text; mndCombi c;
      if (!N2) { p = c.setKneading(N1); c.getAddress(N2); }
      else { p = c.setAddress(N2); c.getKneading(N1); }
      QmnCombiDialog::codeToKneading(N1, t1);
      QmnCombiDialog::codeToAddress(N2, t2);
      text = trUtf8(
         "The kneading sequence  %1  corresponds\n"
         "to the internal address  %2 .\n"
         "The period is %3.\n").arg(t1).arg(t2).arg(p);
      q = c.renorm();
      if (q <= 1) text += trUtf8(
         "It is not simply renormalizable.\n");
      else text += trUtf8(
         "It is simply renormalizable with lowest period %1.\n").arg(q);
      q = c.failsBS();
      if (q) text += trUtf8(
         "This combinatorics is not realized by quadratic polynomials,\n"
         "since it fails the Bruin-Schleicher admissibility condition:\n"
         "the abstract Hubbard tree has an evil branch point of period %1."
         ).arg(q);
      else
      {  q = c.count();
         if (q == 1)
         {  text += trUtf8(
            "This combinatorics is realized once on the real axis.\n");
            t1 = trUtf8("external");
         }
         else
         {  text += trUtf8(
            "This combinatorics is realized for %1 complex parameters.\n"
            ).arg(q);
            t1 = trUtf8("smallest");
         }
         q = c.toAngles(n1, n2, d);
         if (q) text += QString("Angles not computed: Error %1").arg(q);
         else text += trUtf8(
         "The %4 angles are %1/%3 and %2/%3 .\n"
         "Do you want to draw the rays and to shift c\n"
         "to the corresponding center?").arg(n1).arg(n2).arg(d).arg(t1);
      }
      QmnUIntDialog *dialog1
         = new QmnUIntDialog(text, 0, 0, 0u, 3, this, (q ? 0 : 1) );
      connect(dialog1, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
      if (!dialog1->exec() || q) return;
      if (ftype == 58)
      {  double l; qulonglong L = 1ULL; L <<= 60;
         l = ((double)(L))*((double)(n1))/((double)(d));
         l = mndAngle::lambda(((qulonglong)(l)), L, 1.0e-12, 500);
         uint n; pplane->getNmax(n);
         pplane->setPoint(((double)(n1))/((double)(d)), l*0.01*n); return;
      }
      q = 0; double x, y; pplane->newtonRay(1, n1, d, x, y, 5, Qt::white, 1);
      if (pplane->newtonRay(1, n2, d, x, y, 5, Qt::white, 2) <= 1
          && !f->find(1, 0, p, x, y)) pplane->setPoint(x, y);
      while (1)
      {  q++; N2 -= 1LL << (p - 1); p = c.setAddress(N2); if (p <= 1) break;
         c.toAngles(n1, n2, d);
         pplane->newtonRay(1, n1, d, x, y, 5,
            (q & 1 ? Qt::yellow : Qt::white), 1);
         pplane->newtonRay(1, n2, d, x, y, 5,
            (q & 1 ? Qt::yellow : Qt::white), 1);
      }
      pplane->newtonRay(1, 0LL, 1LL, x, y, 5,
         (q & 1 ? Qt::yellow : Qt::white), 1);
   }
   if (act == spiderAct || act == slowspiderAct)
   {  if (act == spiderAct) pullBackString = trUtf8(
      "<br>Then hit Return repeatedly to perform the iteration of<br>"
      "the Spider Algorithm. This tentative implementation is not<br>"
      "refining the discretization when spider legs get twisted.<br>"
      "Use Home or enter 0 to quit.");
      qulonglong N1 = 0LL, N2;
      QmnCombiDialog *dialog = new QmnCombiDialog(
         enterAngleString + pullBackString, &N1, &N2, 3, this);
      connect(dialog, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
      if (!dialog->exec()) return;
      int k, p = mndAngle::normalize(N1, N2, k);
      if (!k && p == 1)
      { if (gamma < 0.0) setRay(0); else pplane->setPoint(0.0, 0.0); return; }
      if (!p)
      {  QmnUIntDialog *dialg = new QmnUIntDialog(trUtf8(
            "The angle  %1/%2\n"
            "has  preperiod + period > 64.").arg(N1).arg(N2),
            0, 0, 0u, 3, this, 0);
         connect(dialg, SIGNAL(needHelp(int)), helptabs,SLOT(showPage(int)));
         dialg->exec(); return;
      }
      if (act == spiderAct)
      {  path->spiderInit(N1, N2); double t = mndAngle::radians(N1, N2);
         dplane->setPoint(cos(t), sin(t));
      }
      else path->angleInit(N1, N2);
      error = 1;
   }
   if (act == pathAct)
   {  updateRegion = true; lsppp = 0; uint n, m = 0;
      if (dplane->getType()) resizeWin(sphereAct);
      double x, y; pplane->getPoint(x, y); n = f->period(x, y);
      if (n < 3 || n > 64)
      { n = 3; if (gamma < -1.0) { n = (uint)(-gamma); gamma = 0.0; } }
      pplane->stop(); if (gamma < 0.0) setRay(0); else pMoved();
      gamma = 0.0;
      QmnUIntDialog *dialog = new QmnUIntDialog(trUtf8(
         "<b></b>A polynomial may be deformed such that the resulting<br>"
         "postcritically finite branched covering is equivalent<br>"
         "to a polynomial again:"
         "<ul><li>"
         "You may shift the n-periodic critical value along a<br>"
         "closed path back to itself, moving around and between<br>"
         "the other points of the critcal orbit. The resulting<br>"
         "branched covering is equivalent to a polynomial with<br>"
         "critical period n again."
         "</li><li>"
         "If the path is enclosing a single postcritcal point,<br>"
         "this gives a Dehn twist about that point and the<br>"
         "critical value. Note that a right-handed path gives a<br>"
         "left-handed twist. You may turn around several times<br>"
         "to define a Dehn twist with higher winding number."
         "</li><li>"
         "Or shift the critcal value c to another point z<sub>1</sub> ,<br>"
         "which is mapped to c in n iterations. (You can find<br>"
         "it with the keys a and b, or approximately from your<br>"
         "intuition of the dynamics.) E.g., if the internal<br>"
         "address 1-...-k-n is realized in the 1/2-sublimb of<br>"
         "period k, you may take the center with 1-...-k as the<br>"
         "initial parameter c and choose z<sub>1</sub> by iterating<br>"
         "0 backwards according to the kneading sequence.<br>"
         "For a direct path, the branched covering will be<br>"
         "equivalent to the polynomial realizing 1-...-k-n.<br>"
         "But recapture surgery is defined in more general<br>"
         "cases: the initial c may be any parameter except 0."
         "</ul>"
         "Examples of period n = 3 are available with the key Ctrl+d.<br>"
         "Now enter the period 3 &le; n &le; 64 and draw the path by<br>"
         "dragging the mouse. (You may cancel it by releasing the<br>"
         "left button outside of the image. To zoom in, you may<br>"
         "turn the path off temporarily with the context menu.)")
         + pullBackString, &n, &m, 64u, 3, this, 1);
      connect(dialog, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
      if (!dialog->exec()) return;
      if (n < 3) { setRay(0); return; }
      gamma = -((double)(n)); if (signtype > 0) paraDyn();
      dplane->setPoint(x, y); dplane->drawOrbit(f, x, y, 0, 4000);
      dplane->move(9); return;
   }
   if (act == redrawAct) //from dMoved(), from user path
   {  if (dplane->hasPathLength() <= 0 || gamma > -2.0) return;
      ldouble *X = new ldouble[101], *Y = new ldouble[101];
      dplane->getUserPath(100, X, Y); int M = (int)(X[0]);
      double x, y; f->getParameter(x, y);
      X[0] = (ldouble)(x); Y[0] = (ldouble)(y); dplane->move(10); paraDyn();
      path->shiftInit((int)(-gamma), M, X, Y); delete[] X; delete[] Y;
      error = 1; act = pathAct;
   }
   if (act == twistAct)
   {  const int M = 100; int m, p = 0;
      ldouble a, b, A, B, u, v; double w, z;
      a = -0.122561166876654L; b = 0.744861766619744L; //rabbit
      //a = -1.266423235901739L; b = 0.357398971748218L; //6347/16383
      pplane->getPoint(w, z); A = (ldouble)(w); B = (ldouble)(z); w = 1.0;
      if (dplane->getType()) dplane->setType(0); if (gamma < 0.0) setRay(0);
      updateRegion = true; lsppp = 0; pplane->stop(); pplane->setPoint(a, b);
      QmnDoubleDialog *dialog = new QmnDoubleDialog(trUtf8(
         "<b></b>To perform a Dehn twist around the Rabbit's ears,<br>"
         "enter the winding number: 1 ... 10 is left handed<br>"
         "and -1 ... -10 is right-handed.<br>"
         "Or enter 400 or 500 to see examples of recapture<br>"
         "surgery: the critical value is shifted to a<br>"
         "preimage along an arc.") + pullBackString, &w, 0, 3, this);
      connect(dialog, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
      if (!dialog->exec() || w < -10.0 || w > 500.0) p = 1;
      else { m = (int)(w); if(!m || (m > 10 && m != 400 && m != 500)) p = 1; }
      if (p) { pplane->setPoint(A, B); return; }
      dplane->setPoint(a, b);
      ldouble *X = new ldouble[M+1], *Y = new ldouble[M+1];
      if (m >= 400) //shift to preimage along straight line segment
      {  if (m == 400)
         { p = 4; A = -0.069858045513714L; B = 0.978233344895305L; }
         else
         { p = 5; A = 0.027871676743202L; B = 0.902736237344649L; }
         A -= a; B -= b; A /= M; B /= M;
         for (m = 0; m <= M; m++) { X[m] = a + m*A; Y[m] = b + m*B; }
      }
      else //shift to itself along circle around  z = c^2 + c
      {  w = (-2.0*m)*PI/M; //left-handed twist: exterior left, interior right
         p = 3; A = -0.662358978622373L; B = 0.562279512062301L;//rabbit
         //p = 14; A = -1.253762663411539L; B = 0.368547951368429L;//6347/16383
         A = 0.7*A + 0.3*a; B = 0.7*B + 0.3*b; a -= A; b -= B;
         for (m = 0; m <= M; m++)
         {  v = w*m; u= cos(v); v = 0.3L*sin(v);
            X[m] = A + u*a - v*b; Y[m] = B + u*b + v*a;
         }
      }
      path->shiftInit(p, M, X, Y); delete[] X; delete[] Y; error = 1;
   }
   if (error)
   {  gamma = -1.0; drawLater = 0; lsppp = 0;
      pplane->setCursorPix(spiderPix); dplane->setCursorPix(spiderPix);
      updateRegion = true; if (dplane->getType()) dplane->setType(0);
      dplane->setNmax(0); updateActions();
      dplane->setPlane(0.0, 0.0, 2.0, 0.0); dplane->draw(f, 0, themode);
      if (act == slowspiderAct) act = stepAct;
      else dplane->drawPathSegment(path);
   }
   error = 0;
   if (act == stepAct && gamma == -1.0)
   {  ldouble x, y; error = path->pathStep(x, y);
      dplane->setPoint(x, y); pplane->setPoint(x, y);
      f->setParameter(x, y); dplane->drawPathSegment(path);
   }
   /*if (error > 0)
   {  QmnUIntDialog *dialog = new QmnUIntDialog(trUtf8(
         "The homotopy type has changed because the\n"
         "discretization is too coarse.\n"
         "The algorithm may converge to a wrong value."),
         0, 0, 0u, 3, this, 0);
      connect(dialog, SIGNAL(needHelp(int)), helptabs,SLOT(showPage(int)));
      dialog->exec(); return;
   } //optional quit?*/
   if (error < 0)
   {  QmnUIntDialog *dialog = new QmnUIntDialog(trUtf8(
         "Further iteration is pointless due\n"
         "to floating-point cancellations."),
         0, 0, 0u, 3, this, 0);
      connect(dialog, SIGNAL(needHelp(int)), helptabs,SLOT(showPage(int)));
      dialog->exec(); act = 0;
   }
   if (act == 0 && gamma < 0.0)
   {  gamma = 0.0; path->deactivate(); dplane->move(10);
      pplane->setCursorPix(0); dplane->setCursorPix(0);
      updateActions(); pMoved(); //twice when called from homeAct
   }
} //setRay

Algorithms

Spider algorithm

"The Thurston Algorithm is iterating the critical orbit backwards to compute the parameter c corresponding to a given external angle. To choose the correct preimages, a connecting path between point configurations is used by Chéritat for slow mating. The Spider Algorithm employs rays to infinity. Use s or Ctrl+s to illustrate these algorithms. With d you may move the critical value along some path to define a Dehn twist or a recapture. Special examples are available with Ctrl+d. Hit Return repeatedly to observe the iteration. Use Home or enter 0 to quit." ( from help )

Wake
Parameter plane ( left with wake /limb of Mandelbrot set) and dynamic plane ( right with Julia set and external rays) for combinatorial rotation number = 13/34

Examples : 

/*
 
This is not official program by W Jung,
but it usess his code ( I hope in a good way)
   These functions are part of Mandel 5.9 by Wolf Jung (C) 2007-2013,
   which is free software; you can
   redistribute and / or modify them under the terms of the GNU General
   Public License as published by the Free Software Foundation; either
   version 3, or (at your option) any later version. In short: there is
   no warranty of any kind; you must redistribute the source code as well.
   
   http://www.mndynamics.com/indexp.html
   
   to compile :
   g++ w.cpp -Wall -lm
   ./a.out
   
   
*/

#include <iostream> // std::cout
#include <complex> // std::complex, std::norm

#define  PI      3.1415926535897932385L //from mndynamo.h 

// type qulonglong  = unsigned long long int 
// n is a numerator of external angle that land on root point of the wake k/r
// d is a denominator 
// funcion mndAngle::wake from mndcombi.cpp  by Wolf Jung (C) 2007-2013
unsigned long long int wake(int k, int r, unsigned long long int  &n)
{  
   if (k <= 0 || k >= r || r > 64) return 0LL; // 
   unsigned long long int  d = 1LL; int j, s = 0; n = 1LL;
   
   for (j = 1; j < r; j++)
   {  s -= k; 
      if (s < 0) s += r; 
      if (!s) return 0LL;
      if (s > r - k) n += d << j;
   }
   d <<= (r - 1); 
   d--; 
   d <<= 1; 
   d++; //2^r - 1 for r <= 64
   
   return d;
}

// from mndynamo.cpp  by Wolf Jung (C) 2007-2013
void root(double x, double y, double &u, double &v)
{  v = sqrt(x*x + y*y);
   if (x > 0.0) { u = sqrt(0.5*(v + x)); v = 0.5*y/u; return; }
   if (x < 0.0)
   { v = sqrt(0.5*(v - x)); if (y < 0.0) v = -v; u = 0.5*y/v; return; }
   if (y >= 0.0) { u = sqrt(0.5*y); v = u; return; }
   u = sqrt(-0.5*y); v = -u;
}

// int mndlbrot::bifurcate from mndlbrot.cpp  by Wolf Jung (C) 2007-2013
// type mndplex = complex 
int bifurcate(double t, double &a, double &b)
{  int per = 1; if (a < -0.75) per = 2;
   if (a < -1.5 || b > sqrt(27/64.0L) || b < -sqrt(27/64.0L) ) per = 3;
   if (t <= -1.0) return per;
   t *= (2*PI);
   if (per == 1)
   { a = 0.5*cos(t) - 0.25*cos(2*t); b = 0.5*sin(t) - 0.25*sin(2*t); }
   if (per == 2) { a = 0.25*cos(t) - 1.0; b = 0.25*sin(t); }
   if (per <= 2) return per;
   std::complex<double> u, c, c1, l = std::complex<double>(cos(t), sin(t));
   if (a < -1.54) c1 = -1.754877666;
   else
   { c1 = std::complex<double>(-.122561167, .744861767); if (b < 0) c1 = conj(c1); }
   c = 81.0*l*l-528.0*l+4416.0; root(real(c), imag(c), a, b);
   u = pow(-13.5*l*l+144.0*l-800.0 + (-1.5*l+12.0)*std::complex<double>(a, b), 1/3.0);
   c = (1/3.0)*(.25*u + (1.5*l+4.0)/u - 2.0);
   if (norm(c - c1) > .25)
   { u *= std::complex<double>(-0.5, sqrt(0.75)); c = (1/3.0)*(.25*u + (1.5*l+4.0)/u - 2.0); }
   if (norm(c - c1) > .25)
   { u *= std::complex<double>(-0.5, sqrt(0.75)); c = (1/3.0)*(.25*u + (1.5*l+4.0)/u - 2.0); }
   a = real(c); b = imag(c); return 3;
} //bifurcate



// mndlbrot::find from mndlbrot.cpp  by Wolf Jung (C) 2007-2013
// sign 		int. Defined in mndynamo.h  positive is parameter plane, negative is dynamic plane."
int find(int sg, unsigned int preper, unsigned int per, double cx, double cy, double &x, double &y) 
{  double a = cx, b = cy, fx, fy, px, py, w; 
   unsigned int i, j;
   for (i = 0; i < 30; i++)
   {  if (sg > 0) { a = x; b = y; }
      if (!preper)
      {  if (sg > 0) { fx = 0; fy = 0; px = 0; py = 0; }
         else { fx = -x; fy = -y; px = -1; py = 0; }
      }
      else
      {  fx = x; fy = y; px = 1.0; py = 0;
         for (j = 1; j < preper; j++)
         {  if (px*px + py*py > 1e100) return 1;
            w = 2*(fx*px - fy*py); py = 2*(fx*py + fy*px);
            px = w; if (sg > 0) px++;
            w = fx*fx - fy*fy + a; fy = 2*fx*fy + b; fx = w;
         }
      }
      double Fx = fx, Fy = fy, Px = px, Py = py;
      for (j = 0; j < per; j++)
      {  if (px*px + py*py > 1e100) return 2;
         w = 2*(fx*px - fy*py); py = 2*(fx*py + fy*px);
         px = w; if (sg > 0) px++;
         w = fx*fx - fy*fy + a; fy = 2*fx*fy + b; fx = w;
      }
      fx += Fx; fy += Fy; px += Px; py += Py;
      w = px*px + py*py; if (w < 1e-100) return -1;
      x -= (fx*px + fy*py)/w; y += (fx*py - fy*px)/w;
   }
   return 0;
}



// =====================================================================================================================
// ====================================================================================================================
// ========================================================================================================================

int main()
{
  unsigned long long int p, q;
  unsigned long long int num, denom;
  double cx,cy;
  double zx,zy;
  double t;
  
  // --------------------------------------------------------------------------------------------------------------------
  // --------------------  initial value ------------------------------------------------------------------------------
  // ------------------------------------------------------------------------------------------------------------------
   std::cout << "Determine the wake of a limb at the main cardioid of Mandelbrot set. " << "\n";
   std::cout << "Enter a fraction  k/r  for the rotation number, in lowest terms, with  1 ≤ k < r ≤ 64 " << "\n";    
  while (1)
    {
        std::cout << " Enter numerator of the rotation number, it is integer  1 ≤ k < 64 :  " << "\n";
        std::cin >> p ;

        if (std::cin.fail()) // no extraction took place
        {
            std::cin.clear(); // reset the state bits back to goodbit so we can use ignore()
            std::cin.ignore(1000, '\n'); // clear out the bad input from the stream
            continue; // try again
        }

        std::cin.ignore(1000, '\n'); // clear out any additional input from the stream

        if (p >= 0) break; // if input value is in good range then exit 
            

	
    }


  
  while (1)
    {
        std::cout << "Enter the denominator of the rotation number, it is integer 1 ≤ r < 64 :  " << "\n";
        std::cin >> q ;

        if (std::cin.fail()) // no extraction took place
        {
            std::cin.clear(); // reset the state bits back to goodbit so we can use ignore()
            std::cin.ignore(1000, '\n'); // clear out the bad input from the stream
            continue; // try again
        }

        std::cin.ignore(1000, '\n'); // clear out any additional input from the stream

        if (q > 0) break; // if input value is in good range then exit 
            

	
    }
    
  std::cout.precision( 15 );  
  std::cout << "The rotation number is " << p << "/" << q << "\n";
  denom = wake(p,q,num);  
    
  if (denom!=0LL)
  {
    
    std::cout << "The "<< p << "/" <<q <<" wake of the main cardioid is bounded by the parameter rays with the angles :\n";
    std::cout << num << "/" << denom << " and "<< num+1 << "/" << denom << "\n";
  }
  else std::cout << "(k <= 0 || k >= r || r > 64) \n";
  
  t=(double)p/((double)q);
  bifurcate(t, cx,cy);
  std::cout << "The root point of wake is c = "<< cx << " ; " << cy << ":\n";
  
  find(-1,0,1,cx,cy,zx,zy);
  std::cout << "The fixed point alfa is z = "<< zx << " ; " << zy << ":\n";

   return 0;
}

Methods

Backwards iteration 
// qmnplane.cpp by Wolf Jung (C) 2007-2014. Mandel 5.12, 
//rewrite double array allocation as in mndPath

int QmnPlane::backRay(qulonglong num, qulonglong denom, double &a, double &b,
  const int quality, QColor color, int mode) // = 5, white, 1
{  //Draw a dynamic ray with angle  num/denom  by backwards iteration with
   //quality points per step.  At present only for quadratic polynomials.
   //mode : 0 all rays, 1 one ray, 2 return endpoint in a, b.
   if (type > 0) return 2; //works on sphere and on plane
   if (quality <= 1) return 3;
   int i, j, k, preper = 0, pppp = 0; //pppp = preper + per
   mndAngle t; pppp = t.setAngle(num, denom, preper); if (!pppp) return 4;
   //pppp += preper; double c, s, X[quality][], Y[quality][];
   //X = new double[quality][pppp + 1]; Y = new double[quality][pppp + 1];
   pppp += preper; double c, s,
   X[quality][pppp + 1], Y[quality][pppp + 1];

   //The second index corresponds to the i-th iterate of the angle. Initial
   //radii between 2^12 and 2^24 : 2^(24*2^(-k/quality))
   s = log(2); c = 24*s; s = exp(-s/double(quality));
   for (k = 0; k < quality; k++)
   { c *= s; X[k][pppp] = exp(c); Y[k][pppp] = 0.5/X[k][pppp]; }
   //Using the approximation Phi_c^{-1}(w) ~ w - 0.5*c/w :
   for (i = 0; i < pppp;i++)
   {  s = t.radians(); c = cos(s); s = sin(s);
      for (k = 0; k < quality; k++)
      {  X[k][i] = c*X[k][pppp] - (b*s + a*c)*Y[k][pppp];
         Y[k][i] = s*X[k][pppp] + (a*s - b*c)*Y[k][pppp];
      }
      t.twice();
   }
   for (k = 0; k < quality; k++)
   { X[k][pppp] = X[k][preper]; Y[k][pppp] = Y[k][preper]; }
   stop(); QPainter *p = new QPainter(buffer); p->setPen(color);
   //2 more for bailout |z| = 4 :
   for (j = nmax + 2; j; j--) for (k = 0; k < quality; k++)
   {  for (i = 0; i < pppp; i++)
      {  c = X[(k ? k - 1 : quality - 1)][i];
         s = Y[(k ? k - 1 : quality - 1)][i];
         mndynamics::root(X[k][i + 1] - a, Y[k][i + 1] - b, X[k][i], Y[k][i]);
         if (c*X[k][i] + s*Y[k][i] < 0) //choose closest preimage
         { X[k][i] = -X[k][i]; Y[k][i] = -Y[k][i]; }
         //if (k & 1) p->setPen(Qt::red); else p->setPen(Qt::blue);
         int i1, k1, i2, k2;
         if ( (!i || !mode) && pointToPos(c, s, i1, k1) <= 1
            && pointToPos(X[k][i], Y[k][i], i2, k2) <= 1)
            p->drawLine(i1, k1, i2, k2);
      }
      X[k][pppp] = X[k][preper]; Y[k][pppp] = Y[k][preper];
   }
   if (mode == 2) { a = X[quality - 1][0]; b = Y[quality - 1][0]; }
   delete p; /*delete[] Y; delete[] X;*/ update(); return 0;
} //backRay
Newton method 
// qmnplane.cpp by Wolf Jung (C) 2007-2014. Mandel 5.12, 
//Time ~ nmax^2 , therefore limited  nmax .
int QmnPlane::newtonRay(int signtype, qulonglong N, qulonglong D,
   double &a, double &b, int q, QColor color, int mode) //5, white, 1
{  uint n; int j, i, k, i0, k0; mndAngle t; t.setAngle(N, D, j);
   double logR = 14.0, x, y, u;
   u = exp(0.5*logR); y = t.radians();
   x = u*cos(y); y = u*sin(y);
   if (pointToPos(x, y, i0, k0) > 1) i0 = -10;
   stop(); QPainter *p = new QPainter(buffer); p->setPen(color);
   for (n = 1; n <= (nmax > 5000u ? 5000u : nmax + 2); n++)
   {  t.twice();
      for (j = 1; j <= q; j++)
      {  if (mode & 4 ? tricornNewton(signtype, n, a, b, x, y,
                exp(-j*0.69315/q)*logR, t.radians())
           : rayNewton(signtype, n, a, b, x, y,
                exp(-j*0.69315/q)*logR, t.radians()) )
         { n = (n <= 20 ? 65020u : 65010u); break; }
         if (pointToPos(x, y, i, k) > 1) i = -10;
         if (i0 > -10)
         {  if (i > -10) { if (!(mode & 8)) p->drawLine(i0, k0, i, k); }
            else { n = 65020u; break; }
         }
         i0 = i; k0 = k;
      }
   }
   //if rayNewton fails after > 20 iterations, endpoint may be accepted
   delete p; update(); if (n >= 65020u) return 2;
   if (mode & 2) { a = x; b = y; }
   if (n >= 65010u) return 1; else return 0;
} //newtonRay

int QmnPlane::rayNewton(int signtype, uint n, double a, double b,
   double &x, double &y, double rlog, double ilog)
{  uint k, l; double fx, fy, px, py, u, v = 0.0, d = 1.0 + x*x + y*y, t0, t1;
   for (k = 1; k <= 60; k++)
   {  if (signtype > 0) { a = x; b = y; }
      fx = cos(ilog); fy = sin(ilog);
      t0 = exp(rlog)*fx - 0.5*exp(-rlog)*(a*fx + b*fy);
      t1 = exp(rlog)*fy + 0.5*exp(-rlog)*(a*fy - b*fx);
      fx = x; fy = y; px = 1.0; py = 0.0;
      for (l = 1; l <= n; l++)
      {  u = 2.0*(fx*px - fy*py); py = 2.0*(fx*py + fy*px);
         px = u; if (signtype > 0) px++;
         u = fx*fx; v = fy*fy; fy = 2.0*fx*fy + b; fx = u - v + a;
         u += v; v = px*px + py*py; if (u + v > 1.0e100) return 1;
      }
      fx -= t0; fy -= t1; if (v < 1.0e-50) return 2;
      u = (fx*px + fy*py)/v; v = (fx*py - fy*px)/v;
      px = u*u + v*v; if (px > 9.0*d) return -1;
      x -= u; y += v; d = px; if (px < 1.0e-28 && k >= 5) break;
   }
   return 0;
} //rayNewton
Argument tracing 
//  mndlbrot.cpp  by Wolf Jung (C) 2007-2014. Mandel 5.12, 

int mndlbrot::turnsign(double x, double y)
{/*Calculate the turn, i.e., the argument of Phi. Return +-1 by comparing
   temp[0] and the turn, 0 for failure or far from the ray. Using two
   tricks to reduce the ambiguity from the multi-valued argument:
   First, the argument should jump on the Julia set instead of the
   line [0, c]. Approximate the Julia set by the lines [0, alpha] and
   [alpha, c] and change the argument accordingly within the triangle.
   Second, keep track of the arguments in temp[j] to detect jumps.
   Before searching a starting point for drawing the external ray,
   calling prepare(201) sets temp[0] = theta and temp[j] = 0,
   and it disables comparison by setting drawmode = 0. Later on, before
   tracing the ray, prepare(202) enables comparison by drawmode = 1.
   */
   double a = x, b = y; if (sign < 0) { a = A; b = B; }
   int j; double s = 1.0, dr = 0.5, theta, u, v, r, eps = 0.004;
   if (x*x + y*y < 1e-12) return 0; //prevent atan2(0, 0) if disconnected
   theta = atan2(y, x); root(.25 - a, -b, u, v); double X = .5 - u, Y = -v;
   for (j = 1; j <= 63; j++)
   {  s *= dr; u = x*x; v = y*y; r = u + v; if (r < 1e-12) return 0;
      u -= v; v = 2*x*y; x = u + a; y = v + b;
      //theta += s*u; First adjust argument in triangle:
      u = atan2(u*y - v*x, u*x + v*y);
      if ( (y*a - x*b)*(Y*a - X*b) > 0
        && (y*X - x*Y)*(b*X - a*Y) > 0
        && ((b-y)*(a-X) - (a-x)*(b-Y))*(a*Y - b*X) > 0)
      { if (u < 0) u += 2*PI; else u -= 2*PI; }
      //Second compare and shift.  3.6 is ok for initial value 0:
      if (drawmode)
      {  if (u > temp[j] + 3.6) u -= 2*PI;
         if (u < temp[j] - 3.6) u += 2*PI;
         temp[j] = u;
      }
      theta += s*u; if (r > 1e18*s) break;
   }
   //Problem: j larger is inaccuarte, but thus ray ends at esctime 64:
   if (r < 100) return 0; //prevent strong inaccuracy. Or r < 1e10*s ?
   theta *= (.5/PI);
   theta -= temp[0]; theta -= floor(theta); // 0 <= theta < 1
   if (theta < eps) return 1; if (1 - eps < theta) return -1;
   return 0;
} //turnsign



// qmnplane.cpp by Wolf Jung (C) 2007-2014. Mandel 5.12, 
//rewrite with posToPoint?
int QmnPlane::traceRay(int signtype, double t, mndynamics *f,
  double &x0, double &y0, int quality, QColor color) //5, white
{ //Draw the external ray for the turn t with the argument tracing method.
  //Return 0 if the endpoint was found, 1 no startpoint, 2 or 3 no endpoint.
  //f->turnsign() checks if the turn is in a small intervall around t,
  //returns +-1 for the side and 0 otherwise, may adjust jumps.
  if (type) return 4;
  int i, i0, i1, i2, k, k0, k1, k2, j, sign, draw = 0, iold, kold, u, v;
  //Set signtype and t in mndynamics,  disable adjusting of jumps:
  uint m = 201; f->prepare(signtype, 0, m, &t);
  //First,  search the startpoint on the boundary enlarged by  quality >= 1:
  i = quality*imax; k = -quality*kmax - 1; i2 = -2; k2 = 0; //go left on top
  j = f->turnsign(hmid + pw*i + ph*k, vmid + ph*i - pw*k);
  while (1)
  {  i += i2; k += k2;
     sign = f->turnsign(hmid + pw*i + ph*k, vmid + ph*i - pw*k);
     if (j < 0 && sign > 0)
     {  iold = i2/2; kold = k2/2;
        i1 = i; k1 = k;
        i0 = i - iold; k0 = k - kold;
        if (f->turnsign(hmid + pw*i0 + ph*k0, vmid + ph*i0 - pw*k0) > 0)
        { i1 = i0; k1 = k0; i0 -= iold; k0 -= kold; }
        i2 = i0 + kold; k2 = k0 - iold; u = i1; v = k1;
        break; //found startpoint
     }
     j = sign;
     if (i < -quality*imax && i2 < 0) { i2 = 0; k2 = 2; } //down on left line
     if (k >= quality*kmax && k2 > 0) { i2 = 2; k2 = 0; } //right on bottom
     if (i >= quality*imax && i2 > 0) { i2 = 0; k2 = -2; } //up on right line
     if (k < -quality*kmax && k2 < 0) return 1;
  }
  //Second, trace the ray by triangles with sign(z0) < 0, sign(z1) > 0:
  stop(); QPainter *p = new QPainter(buffer); p->setPen(color);
  m = 202; f->prepare(signtype, 0, m, &t); //adjusting jumps enabled
  for (j = 1; j < quality*3*(imax + kmax); j++)
  {  sign = f->turnsign(hmid + pw*i2 + ph*k2, vmid + ph*i2 - pw*k2);
     if (sign > 0) { i = i1; k = k1; i1 = i2; k1 = k2; }
     else { i = i0; k = k0; i0 = i2; k0 = k2; }
     //New reflected point z2:
     i2 = i0 + i1 - i; k2 = k0 + k1 - k;
     //If not reflected at diagonal edge, take maximal distance to (u, v) :
     if (i0 == i1) { if (mabs(k1-v) > mabs(k0-v)) k2 = k1; else k2 = k0; }
     if (k0 == k1) { if (mabs(i1-u) > mabs(i0-u)) i2 = i1; else i2 = i0; }
     u = i; v = k;
     //Check viewport and draw at z1 (t = 0 | 1/2 too low with z0; z2 rough):
     if (-imax <= i1 && i1 < imax && -kmax <= k1 && k1 < kmax)
     {  if (draw)
        {  // >= 8 is less smooth but not going in circles at 399/992:
           if ((i1-iold)*(i1-iold)+(k1-kold)*(k1-kold) >= 5)
           {  p->drawLine(imax+iold, kmax+kold, imax+i1, kmax+k1);
              iold = i1; kold = k1;
           }
           if (!sign)
           {  x0 = hmid + pw*i1 + ph*k1; y0 = vmid + ph*i1 - pw*k1;
              delete p; update(); return 0;
           }
        }
        else
        { draw = 1; iold = i1; kold = k1; }
     }
     else if (draw) //has left the viewport, no endpoint
     { delete p; update(); return 2; }
  }
  delete p; update(); return 3;
} //traceRay

Colors

Code

classes

Classes types :

Base classes :

mndynamics

mndynamics is an abstract base class (escape to infinity)

Classes derived from mndynamics


Most classes describe a one-parameter family of complex functions  :

where :

Global variables

and the same angles are found in the parameter plane.

viewport

Wolf Jung uses center and width for describing rectangle viewport of complex plane: 

/* 
   from mndlbrot.cpp  by Wolf Jung (C) 201
   These classes are part of Mandel 5.7, which is free software; you can
   redistribute and / or modify them under the terms of the GNU General
   Public License as published by the Free Software Foundation; either
   version 3, or (at your option) any later version. In short: there is
   no warranty of any kind; you must redistribute the source code as well
*/
void mndlbrot::startPlane(int sg, double &xmid, double &rewidth) const
{  if (sg > 0) 
     { xmid = -0.5; rewidth = 1.6; } // parameter plane
     else { xmid = 0; rewidth = 2.0; } // dynamic plane

So initial values of corners are :

Images

There is a category on commons : Category:Fractals created with Mandel

Demos

See :

Demo 2 : periodic points and bifurcations

page 7

Golden Mean Quadratic Siegel Disk ( with some orbits inside) ~ initial image

Finding c parameters (on the boundary of main cardioid) with Siegel Disk :


If you choose any ( ~ random) number from <0.0; 1.0 ) then you will have Siegel disk. Here is the explanation

//
t= .61803398874989484820L;
for (i=0; i<10; ++i) { 
  t *= (2*PI); // from turns to radians
  cx = 0.5*cos(t) - 0.25*cos(2*t); 
  cy = 0.5*sin(t) - 0.25*sin(2*t); 
  // increase t
  t += dt; // t = t + dt
}


Here is the original code : 

//  qmndemos.cpp  by Wolf Jung (C) 2007-2013. part of Mandel 5.12,
// 
   pplane->draw(pf, 1, &mode); double x = 0, y = 0;
   pf->bifurcate(.61803398874989484820L, x, y); pplane->setPoint(x, y);

void QmnDemoPB::go()
{
// .....
if (page == 7)
   {  double x, y; pplane->getPoint(x, y); // 2alpha :
      mndynamics::root(1 - 4*x, -4*y, x, y); // 
      x = 1 - x; 
      y = -y;
      y = atan2(y, x) + 0.01;
      x = 0.5*cos(y) - 0.25*cos(2*y); 
      y = 0.5*sin(y) - 0.25*sin(2*y);
      pplane->setPoint(x, y);}
}

opening images

Eample image which was checked and modified with Mandel

One can open png image in Mandel and check it.

Image must be :

Steps :

Then you can :

You can't:

Videos

videos by Wolf Jung produced with FFmpeg, combining images made with Mandel[9]

Wishlist

Papers with images made with Mandel

References

  1. ↑ Mandel by W Jung
  2. ↑ Mandelbrot set in C++
  3. ↑ Mandelbrot set in C++ and SDL
  4. ↑ ArgPhi
  5. ↑ symbolic dynamics
  6. ↑ Multiplatform cpp program Mandel by Wolf Jung
  7. ↑ Formation of Escape-Time Fractals By Christopher Olah
  8. ↑ wikipedia : Golden_ratio_conjugate
  9. ↑ Mandelbrot set and quadratic polynomials - videos by Wolf Jung
This article is issued from Wikibooks. The text is licensed under Creative Commons - Attribution - Sharealike. Additional terms may apply for the media files.