Fractals/mandel

From Wikibooks, open books for an open world
Jump to navigation Jump to search

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

install[edit | edit source]

from sources[edit | edit source]

  • download
  • extract
  • qmake
  • make
  • ./mandel

Features[edit | edit source]

  • multiplatform ( the platform-independent )
  • written in C++ with source code on GNU GPL licence
  • libraries :
    • GUI-toolkit Qt
  • many algorithms
  • many maps
  • numerical precision
    • floating point limited to double
    • integer :


angle[edit | edit source]

Program mandel uses own notation for binary fractions:

Enter the angle  
with  
Notation:
  decimal 1/7 or binary p001 for  (binary periodic angles)
  decimal 9/56 or binary 001p010 for  (binary preperiodic angles)
  decimal 1/4 or binary 01 for 1/4 = 0.01 (binary dyadic angles with both finite and equal infinite preperiodic binary expansion)

maps or functions[edit | edit source]

Families of functions (= maps ):

  • Mandelbrot set z^2+c ( key Ctrl+0 ). It is default map. Here map degree q = 2
  • Multibrot sets z^q+c ( key Ctrl+1 ). Use key q to change map degree
  • Branner-Fagella cz(1+z/q)^q ( key Ctrl+2 )
  • cubic polynomials
  • quartic polynomials
  • rational mappings
  • Newton mappings
  • transcendental
  • Non-analytic mappings

These famnilies are defined by classes

complex quadratic polynomial[edit | edit source]

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 :

  • forward iteration , key f = mode 0
  • 2 inverse ( or backward ) iterations ( multivalued with argument adjusted ) :
    • 1-st inverse function = key a =
    • 2-nd inverse function = key b =
   "Use the keys a and b for the inverse mapping. (The two branches are approximately mapping to the parts A and B of the itinerary.)



See also example c code

// 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 ...  )" 
    ...
    "the argument adjusting in the inverse branch has nothing to do with computing external arguments.  It is related to itineraries and kneading sequences,  ...
    Kneading sequences are explained in demo 4 or 5, in my slides on the stripping algorithm, and in several papers by Bruin and Schleicher.
 
    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);
}

internal address[edit | edit source]

Internal addres

  • the numbers (period) must be increasing by definition

Algorithms[edit | edit source]

  • drawing
    • plane ( drawmode or algoActs)
      • LSM Escape time ( algorithm 2)
      • DEM (algorithm 5)
    • curves
    • points
      • orbits ( iteration)
        • critical orbit
        • forward and inverse orbit
  • computing :

Example programs which use code from Mandel

  • Mandelbrot set in C++ [2]
  • Mandelbrot set in C++ and SDL[3]
  • ArgPhi [4]
  • symbolic dynamics [5]
// qmnshell.cpp  by Wolf Jung (C) 2007-2019.
QActionGroup *algoAG = new QActionGroup(this);
   algoActs[0] = new QAction(QString("&0"), algoAG);
   algoActs[0]->setShortcut(QString("0"));
   algoActs[1] = new QAction(trUtf8("&1: Comparing neighbors"), algoAG);
   algoActs[1]->setShortcut(QString("1"));
   algoActs[2] = new QAction(trUtf8("&2: Escape time"), algoAG);
   algoActs[2]->setShortcut(QString("2"));
   algoActs[3] = new QAction(trUtf8("&3: Escape time, rainbow colors"),algoAG);
   algoActs[3]->setShortcut(QString("3"));
   algoActs[4] = new QAction(trUtf8("&4: Marty normality"), algoAG);
   algoActs[4]->setShortcut(QString("4"));
   algoActs[5] = new QAction(trUtf8("&5: Distance estimate"), algoAG);
   algoActs[5]->setShortcut(QString("5"));
   algoActs[6] = new QAction(trUtf8("&6: Argument of \xce\xa6"), algoAG);
   algoActs[6]->setShortcut(QString("6"));
   algoActs[7] = new QAction(trUtf8("&7: Binary decomposition"), algoAG);
   algoActs[7]->setShortcut(QString("7"));
   algoActs[8] = new QAction(trUtf8("&8: Yoccoz puzzle 1/2 1/3 2/3"), algoAG);
   algoActs[8]->setShortcut(QString("8"));
   algoActs[9] = new QAction(trUtf8("&9: Zeros of q_n(c) [Q]"), algoAG);
   algoActs[9]->setShortcut(QString("9"));
   algoActs[10] = new QAction(trUtf8("&c: Cr. points of q_n(c) [Q]"), algoAG);
   algoActs[11] = new QAction(trUtf8("&n: Newton for q_n(c) [Q]"), algoAG);
   connect(algoAG, SIGNAL(triggered(QAction*)), this, SLOT(drawing(QAction*)));
// mndlbrot.cpp  by Wolf Jung (C) 2007-2019.
uint mndlbrot::pixcolor(mdouble x, mdouble y)
{ 

   uint j, cl = 0; mdouble a, b;
   
   if (sign > 0) { a = x; b = y; } //z = c, not z = 0 !
   else { a = A; b = B; }


   if (drawmode == 4) return marty(a, b, x, y);
   if (drawmode == 7) return binary(a, b, x, y);
   if (drawmode == 8) return yoccoz(a, b, x, y);
   if (sign < 0 ||
      ( (x*x + y*y - .0625)*(x*x + y*y - .0625) >
         ((x - .25)*(x - .25) + y*y)*.25  &&
      (x + 1)*(x + 1) + y*y > .0625 ) )  //not per 1, 2 in M
   {  mdouble x1 = x, y1 = y, u, v;
      for (j = 1; j <= maxiter; j++)
      {  u = x*x; v = y*y;
         if (u + v <= bailout) { y = 2*x*y + b; x = u - v + a; }
         else
         
         {  if (drawmode == 2)
            { cl = ((j % 24) >> 3); if (!cl) cl = 4; if (j & 1) cl |= 8; }
            else cl = j;
            break;
         }
      }
      x = x1; y = y1;
   }
   if (drawmode == 2 || drawmode == 3) return cl;
   if (cl) cl = dist(a, b, x, y);
   if (cl & 2) Period = 1;
   else { Period = 0; if (drawmode == 6) return 0; }
   if (drawmode == 6) return turn(a, b, x, y);
   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);
   if (drawmode == 11) return newton(a, b);
   return cl; //drawmode == 5
}

Level sets of escape time[edit | edit source]

Is it possibly, for algorithm 2 = ( level sets of ) escape time , to set parameters ( ? escape radius) such that boundaries of level sets go thru critical point ( and it's preimages ) ?

You can change the escape radius only by recompiling. But you can choose the parameter such that it is on an escape line and on the desired parameter ray, then the critical value will be on an escape line as well.


// mndlbrot.cpp  by Wolf Jung (C) 2007-2019
uint mndlbrot::esctime(mdouble x, mdouble y)
{  uint j = 0;
   mdouble a, b, x0, y0, u, v;
   if (sign > 0) { a = x; b = y; } //z = c, not z = 0 !
   else { a = A; b = B; }

   if (b == 0.0 && !drawmode && sign < 0
      && (a == 0.25 || a == -0.75)) return parabolic(x, y);
   if (drawmode > 255)
   {  if (sign > 0) { x = 0; y = 0; }
      return renormtime(a, b, x, y);
   }
   if (sign > 0 && (x*x + y*y - .0625)*(x*x + y*y - .0625) <
        ((x - .25)*(x - .25) + y*y)*.25 ) j = 1;
   else if (sign > 0 && (x + 1)*(x + 1) + y*y < .0625 ) j = 2;
   if (sign > 0 && -1.76875 < x && x < -1.74375 && y*y < .00015625)
   {  root(-4*a - 7, -4*b, x0, y0);
      u = a + 2 + a*x0 - b*y0; v = b + a*y0 + b*x0;
      if (u*u + v*v < .0625) j = 3;
   }
   if (sign > 0 && ( (x + .125)*(x + .125)+(y - .75)*(y - .75) < .015625 ||
         (x + .125)*(x + .125) + (y + .75)*(y + .75) < .015625 ) )
   {  root(-4*a - 7, -4*b, x0, y0);
      u = a + 2 - a*x0 + b*y0; v = b - a*y0 - b*x0;
      if (u*u + v*v < .0625) j = 3;
   }
   if (j)
   { if (drawmode) return 65280u | (drawmode & 15); else return 65280 | j; }
   /*if (drawmode && x*x + y*y < 5 && y < 0)
   {  long long int U, V, W = 1LL << 28, Bailout = 5LL << 56;
      long int X = (long int)(x*W), Y = (long int)(y*W),
         A = (long int)(a*W), B = (long int)(b*W);
      for (j = 1; j <= maxiter; j++)
      {  U = X; U *= X; V = Y; V *= Y; if (U + V >= Bailout) return j;
         W = X; W *= Y; W >>= 27; Y = (long int)(W) + B;
         U -= V; U >>= 28; X = (long int)(U) + A;
      }
      return 65284;//0 | (drawmode & 15);
   }//test integer speed up ... problem with Green*/
   for (j = 1; j <= maxiter; j++)
   {  u = x*x; v = y*y;
      if (u + v <= bailout)
      { y = 2*x*y + b; x = u - v + a; }
      else return j;
      if (sign < 0 && !drawmode
         && (x - temp[0])*(x - temp[0]) + (y - temp[1])*(y - temp[1]) < 1e-6)
      {  int cl = 1 + j % Period;
         if (Period <= 4 && j % (2*Period) >= Period ) cl |= 8;
         return 65280u | cl;
      }
   }
   if (drawmode) return 65280u | (drawmode & 15);
   if (sign > 0)
   {  mdouble x0 = x, y0 = y;
      for (j = 0; j < 60; j++)
      {  if (x*x + y*y <= bailout) f(a, b, x, y);
         else return maxiter + j;
         if ( (x - x0)*(x - x0) + (y - y0)*(y - y0) < 1e-6 )
           return 65281u + (j  % 12);
      }
   }
   return 65293u;
}

Distance estimation[edit | edit source]

// mndlbrot.cpp  by Wolf Jung (C) 2007-2023. ... are part of Mandel 5.19,
int mndlbrot::dist(mdouble a, mdouble b, mdouble x, mdouble y)
{  uint j; mdouble xp = 1, yp = 0, nz, nzp; //zp = 1
   for (j = 1; j <= maxiter; j++)
   {  nz = 2*(x*xp - y*yp); yp = 2*(x*yp + y*xp); xp = nz; //zp = 2*z*zp;
      if (sign > 0) xp++; //zp = 2*z*zp + 1
      nz = x*x - y*y + a; y = 2*x*y + b; x = nz; //z = z*z + c;
      nz = x*x + y*y; nzp = xp*xp + yp*yp;
      if (nzp > 1e40 || nz > bailout) break;
   }
   if (nz < bailout) return 1; //not escaping, rare
   if (nzp < nz) return 10; //includes escaping through 0
   x = log(nz); x = x*x*nz / (temp[1]*nzp); //4*square of dist/pixelwidth
   if (x < 0.04) return 1;
   if (x < 0.24) return 9;
   return 10;
} //dist

potential[edit | edit source]

Complex potential of Julia set

The green procedure is named after the British mathematician George Green

// mndlbrot.cpp  by Wolf Jung (C) 2007-2019
mdouble mndlbrot::green(int sg, mdouble x, mdouble y)
{  mdouble a = x, b = y; if (sg < 0) { a = A; b = B; }
   uint j; mdouble s = 1.0, dr = 0.5, u, v;
   for (j = 0; j <= maxiter; j++)
   {  s *= dr; u = x*x; v = y*y;
      if (u + v > 1e12) return s*log(u + v);
      y = 2*x*y + b; x = u - v + a;
   }
   return 0;
}


Simplified function

/* 
this function is based on function by W Jung 
it is used in : File:Cpmj.png
*/
double jlogphi(double zx0, double zy0, double cx, double cy)
{ 
  int j; 
  double 
    zx=zx0,
    zy=zy0,
    s = 0.5, 
    zx2=zx*zx,
    zy2=zy*zy,
    t;
  //
  for (j = 1; j < 400; j++)
    { s *= 0.5; 
      zy = 2 * zx * zy + cy;
      zx = zx2 - zy2 + cx; 
      zx2 = zx*zx; 
      zy2 = zy*zy;
      t = fabs(zx2 + zy2); // abs(z)
      if ( t > 1e24) break; 
    } 
  return s*log2(t);  // log(zn)* 2^(-n)
}//jlogphi

The equipotential line through the current point, or for a given potential, is drawn with the key g.

Equipotential lines can be defined everywhere, but they will be a figure-eight, when they go through the critcal point z = 0 or through one of its preimages.

External rays through 0 or its preimages branch as well.

if (act == greenAct)
   {  mdouble 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);
   }
// qmnplane.cpp by Wolf Jung (C) 2007-2023.
//the following code will be rewritten completely!!!
void QmnPlane::green(mndynamics *f, int st, mdouble g, int quality,
  QColor color) //5, white
{ //trace the boundary ...
  if (g <= 0 || type) return;
  uint m = 1; f->prepare(st, nmax, m, &g); //just to set nmax
  int i, k, i0, k0, i1, k1, i2, k2, j, vert, side, sg, start, T; stop();
  QPainter *p = new QPainter(buffer); p->setPen(color);
for (T = 20 - kmax; T < kmax - 20; T += imax*2/quality)
{ for (vert = 0; vert <= 1; vert++)
  { for (side = -1; side <= 1; side += 2)
    { start = 10000; //no startpoint
      if (vert) //vertical lines from above and below
      { if (f->green(st, hmid - ph*side*kmax + pw*T, vmid - pw*side*kmax + ph*T) <= g)
           continue;
        for (j = kmax-1; j >= -kmax/2; j--)
        { if (f->green(st, hmid + ph*side*j + pw*T, vmid - pw*side*j + ph*T) <= g)
          {start = side*j; break;} }
      }
      else//vert: horizontal lines from left and right
      { if (f->green(st, hmid + pw*side*imax - ph*T, vmid + ph*side*imax - pw*T) <= g)
          continue;
        for (j = imax-1; j >= -imax/2; j--)
        { if (f->green(st, hmid + pw*side*j - ph*T, vmid + ph*side*j - pw*T) <= g)
          {start = side*j; break;} }
      }//vert
      if (start == 10000) continue; //no startpoint
      for (sg = -1; sg <= 1; sg += 2) //go in both directions
      { if (vert)
        { k0 = start; k1 = k0 + side; k2 = k1; i0 = T; i1 = T; i2 = i1 + sg;}
        else
        { i0 = start; i1 = i0 + side; i2 = i1; k0 = T; k1 = T; k2 = k1 + sg;}
        for (j = 1; j < 8*imax; j++)
        { if (i0 < -imax || i0 >= imax || k0 < -kmax || k0 >= kmax) break;
          else p->drawLine(imax + i0, kmax + k0, imax + i0, kmax + k0);
          if (f->green(st, hmid + pw*i2 + ph*k2, vmid + ph*i2 - pw*k2) <= g)
          {i = i0; k = k0; i0 = i2; k0 = k2; i2 += (i1 - i); k2 += (k1 - k);}
          else
          {i = i1; k = k1; i1 = i2; k1 = k2; i2 += (i0 - i); k2 += (k0 - k);}
        }//for j
      }//for sg
    }//for side
  }//for vert
}//for T
  delete p; update();
} //green


Potential :

See also:


Binary decomposition[edit | edit source]

int mndlbrot::binary(mdouble a, mdouble b, mdouble x, mdouble y)
{  mdouble u, v; uint j;
   for (j = 1; j <= maxiter; j++)
   {  u = x*x; v = y*y;
      if (u + v <= 1e4) { y = 2*x*y + b; x = u - v + a; }
      else { if (j & 1) j = 2; else j = 10; if (y > 0) j += 2; return j; }
      //else return (x > 50.0 + y && x > 50.0 - y ? 12 : 10); //carrots
   }
   return 0;
} //binary

Renormalization[edit | edit source]

if (act == renormAct)
   {  uint n1, n2 = 0;
      QmnUIntDialog *dialog = new QmnUIntDialog(trUtf8(
         "Enter the renormalization period \xe2\x89\xa4 255.\n"
         "For embedded Julia sets, enter the preperiod\n"
         "and the period, separated with a comma:"),
         &n1, &n2, 65002u, 1, this);
      connect(dialog, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
      if (dialog->exec() && n1 <= 255 && n2 <= 255) n1 |= (n2 << 8);
      else return;
      if (signtype > 0) { drawLater = 0; dplane->stop(); }
      else pplane->stop();
      theplane->draw(f, signtype, &n1, 2);
   }

periodic and preperiodic points[edit | edit source]

  "Each hyperbolic component has a unique center, and when the current point is within a hyperbolic component, you can move it to the center with the keys x and Return."


Feigenbaum tuning[edit | edit source]

Steps:

  • Start with a center of period n
  • find an approximation for the center of period 2n and use Newton method to find it precisely The approximation can be done with Feigenbaum rescaling, but this did work only in the 1/2-limb, so I used a composition of the multiplier maps of periods 1 and 2 in the other sublimbs.

Keys: Ctrl+Alt+C

C++ code

// from the file qmnshell.cpp  by Wolf Jung (C) 2007-2018
const mdouble cFb = -1.40115518909205060052L, 
              dFb = 4.66920160910299067185L,
              aFb = 2.50290787509589282228L;  //Fibonacci -1.8705286321646448888906L
//
if (act == tuneAct && per && per <= 512)
   {  f->find(signtype, 0, per, x, y);
      if (x < -.5)
         { x = cFb + (x - cFb)/dFb; 
           y /= dFb; }
      else
      {  mndynamics::root(1/16.0 - x*.25, -y*.25, x, y);
         x = -.75 - x; 
          y = -y;
       }
      per *= 2; 
      f->find(signtype, 0, per, x, y);
      if (f->period(x, y) == per) theplane->setPoint(x, y);
   }


See also:

/*

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  by Wolf Jung (C) 
   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++ f.cpp -Wall -lm
   ./a.out
   
   ===========================================
   
Period =          1 	center =  0.000000000000000000
Period =          2 	center = -1.000000000000000000
Period =          4 	center = -1.310702641336832884
Period =          8 	center = -1.381547484432061470
Period =         16 	center = -1.396945359704560642
Period =         32 	center = -1.400253081214782798
Period =         64 	center = -1.400961962944841041
Period =        128 	center = -1.401113804939776124
Period =        256 	center = -1.401146325826946179
Period =        512 	center = -1.401153290849923882
Period =       1024 	center = -1.401154782546617839  // c = -1.401154782546620  +0.000000000000000 i    period = 1024
Period =       2048 	center = -1.401155102022463976
Period =       4096 	center = -1.401155170444411267
Period =       8192 	center = -1.401155185098297292
Period =      16384 	center = -1.401155188236710937
Period =      32768 	center = -1.401155188908863045
Period =      65536 	center = -1.401155189052817413
Period =     131072 	center = -1.401155189083648072
Period =     262144 	center = -1.401155189090251057
Period =     524288 	center = -1.401155189091665208
Period =    1048576 	center = -1.401155189091968106
Period =    2097152 	center = -1.401155189092033014
Period =    4194304 	center = -1.401155189092046745
Period =    8388608 	center = -1.401155189092049779
Period =   16777216 	center = -1.401155189092050532
Period =   33554432 	center = -1.401155189092051127
Period =   67108864 	center = -1.401155189092050572
Period =  134217728 	center = -1.401155189092050593
Period =  268435456 	center = -1.401155189092050599




*/


#include <iostream> // std::cout
#include <cmath>     // sqrt
#include <limits>
#include <cfloat>
  
typedef  unsigned int  uint;
typedef  long double  mdouble; // mdynamo.h


// from the file qmnshell.cpp  by Wolf Jung (C) 2007-2018
mdouble cFb = -1.40115518909205060052L; 
mdouble dFb = 4.66920160910299067185L;
    

mdouble bailout = 16.0L; // mdynamoi.h


// c = A+B*i
mdouble A= 0.0L;  
mdouble B = 0.0L;

/* 
   function from mndlbrot.cpp  by Wolf Jung (C) 2007-2017 ...
   part of Mandel 5.14, 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
   ----------------------------------------------
   
   it is used to find :
   * periodic or preperiodic points on dynamic plane 
   * on parameter plane 
   ** centers
   ** Misiurewicz points
   
   using Newton method
  
   
*/
int find(int sg, uint preper, uint per, mdouble &x, mdouble &y) 
{  mdouble a = A, b = B, fx, fy, px, py, w;
  uint i, j;
   
  for (i = 0; i < 30; i++)
    { if (sg > 0) // parameter plane
	{ a = x; b = y; } 
         
      if (!preper) // preperiod==0
	{  if (sg > 0) // parameter plane
	    { fx = 0; 
	      fy = 0; 
	      px = 0; 
	      py = 0; }
	       
	  else // dynamic plane
	    { fx = -x; 
	      fy = -y; 
	      px = -1; 
	      py = 0; }
	}
	
      else // preperiod > 0
	{ 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++; // parameter plane
	      w = fx*fx - fy*fy + a; 
	      fy = 2*fx*fy + b; 
	      fx = w;
	    }
	}
	
      mdouble 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++; // parameter plane
	  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()
{

  
	int plane = 1; // positive is parameter plane, negative is dynamic plane = signtype
  	uint preper = 0; // " 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 )
  	uint per ; // period
  	mdouble x ;
  	mdouble y = 0.0L;
  	int n;
  
  	// Starting with a center of period  n
  	per = 1;
	x = 0.0L;
	
	// find an approximation for the center of period  2n 
	for (n=1; n<30; n++){
     
     		 printf("Period = %10u \tcenter = %.18Lf\n", per, x);
        
        	// next center
      		per *= 2; // period doubling 
        	// approximate  of next value using Feigenbaum rescaling ( in the 1/2-limb )
        	x = cFb + (x - cFb)/dFb; 
        	// more precise value of x useing Newton method  
        	find(plane, preper, per, x, y); 
      
   }


  return 0;
}


C code

#include <stdio.h> // printf
#include <math.h>     // sqrt

  
typedef  unsigned int  uint;
typedef  long double  mdouble; // mdynamo.h


// from the file qmnshell.cpp  by Wolf Jung (C) 2007-2018
mdouble cFb = -1.40115518909205060052L; 
mdouble dFb = 4.66920160910299067185L;
    

mdouble bailout = 16.0L; // mdynamoi.h


// c = A+B*i
mdouble A= 0.0L;  
mdouble B = 0.0L;

/* 
   function from mndlbrot.cpp  by Wolf Jung (C) 2007-2017 ...
   part of Mandel 5.14, 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
   ----------------------------------------------
   
   it is used to find :
   * periodic or preperiodic points on dynamic plane 
   * on parameter plane 
   ** centers
   ** Misiurewicz points
   
   using Newton method
  
   
*/
int find(int sg, uint preper, uint per, mdouble *x, mdouble *y) 
{  
	mdouble a = A;
 	mdouble b = B;
 	mdouble tx;
 	mdouble ty;
 	mdouble fx;
 	mdouble fy;
 	mdouble px;
 	mdouble py;
 	mdouble w;
  	uint i, j;
   
  for (i = 0; i < 30; i++)
    { if (sg > 0) // parameter plane
	{ 	a = *x; 
		b = *y; } 
         
      if (!preper) // preperiod==0
	{  if (sg > 0) // parameter plane
	    { fx = 0; 
	      fy = 0; 
	      px = 0; 
	      py = 0; }
	       
	  else // dynamic plane
	    { fx = - *x; 
	      fy = - *y; 
	      px = -1; 
	      py = 0; }
	}
	
      else // preperiod > 0
	{ 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++; // parameter plane
	      w = fx*fx - fy*fy + a; 
	      fy = 2*fx*fy + b; 
	      fx = w;
	    }
	}
	
      mdouble 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++; // parameter plane
	  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; }
      	
      tx -= (fx*px + fy*py)/w; 
      ty += (fx*py - fy*px)/w;
    }
    x = &tx;
    y = &ty;
  return 0;
}







int main()
{

  
	int plane = 1; // positive is parameter plane, negative is dynamic plane = signtype
  	uint preper = 0; // " 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 )
  	uint per ; // period
  	mdouble x ;
  	mdouble y = 0.0L;
  	int n;
  
  	// Starting with a center of period  n
  	per = 1;
	x = 0.0L;
	
	// find an approximation for the center of period  2n 
	for (n=1; n<30; n++){
     
     		 printf("Period = %10u \tcenter = %.18Lf\n", per, x);
        
        	// next center
      		per *= 2; // period doubling 
        	// approximate  of next value using Feigenbaum rescaling ( in the 1/2-limb )
        	x = cFb + (x - cFb)/dFb; 
        	// more precise value of x useing Newton method  
        	find(plane, preper, per, &x, &y); 
      
   }


  return 0;
}

find[edit | edit source]

To find periodic or preperiodic point z use:

  • Menu/Points/Find point
  • key x

then enter:

  • period
  • preperiod,period

It is numerical procedure so:

  • it finds one and the nearest solution.
  • it finds nothing ( on parameter plane it goes to origin c=0) if current point is:
    • inside set ( Mandel of filled-in Julia)
    • not near good solution

For example on the parameter plane there are 3 period 3 centers. If:

  • current point is c=0 then result is c=0 ( nothing found)
  • current point is outside Mandelbrot set :
    • c = 0.350000000000000 +0.725000000000000 i then result is c = -0.122561166876654 +0.744861766619744 i period = 3 ( good result)
    • c = 0.580000000000000 +0.240000000000000 i period = 0 then result is c = 0.000000000000000 +0.000000000000000 i period = 1 ( bad result !!!)
    • c = -1.590000000000000 +0.170000000000000 i period = 0 then result is c = -1.754877666246693 +0.000000000000000 i period = 3 ( good result)


Use mode Del-Newton to show basins of attraction of Newton method, adjust period with q

It calls function find with

  • input:
    • sg ( "positive is parameter plane, negative is dynamic plane.")
    • preper = preperiod. ( "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 )
    • per = period
  • output:
    • x = creal(z)
    • y = cimag(z)
    • returning value ( type int)

Returning value:

  • 0 ( good result)
  • 2 when (px*px + py*py > 1e100)
  • 1 when (px*px + py*py > 1e100) = escaping in preperiodic loop
  • -1 when (w < 1e-100)
/* function from mndlbrot.cpp  by Wolf Jung (C) 2007-2017 ...
   part of Mandel 5.14, 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.
*/
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;
}

Example of use:

#include <iostream> // std::cout

typedef  unsigned int  uint;

// c = A+B*i
double A= -0.965000000000000;  
double B = 0.085000000000000;

/* 
   function from mndlbrot.cpp  by Wolf Jung (C) 2007-2017 ...
   part of Mandel 5.14, 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
   ----------------------------------------------
   
   it is used to find :
   * periodic or preperiodic points on dynamic plane 
   * on parameter plane 
   ** centers
   ** Misiurewicz points
   
   using Newton method
   -------------------------------------------
   to compile :
   g++ f.cpp -Wall 
   ./a.out
   
   
*/
int find(int sg, uint preper, uint per, double &x, double &y) 
{  double a = A, b = B, fx, fy, px, py, w;
  uint i, j;
   
  for (i = 0; i < 30; i++)
    { if (sg > 0) // parameter plane
	{ a = x; b = y; } 
         
      if (!preper) // preperiod==0
	{  if (sg > 0) // parameter plane
	    { fx = 0; 
	      fy = 0; 
	      px = 0; 
	      py = 0; }
	       
	  else // dynamic plane
	    { fx = -x; 
	      fy = -y; 
	      px = -1; 
	      py = 0; }
	}
	
      else // preperiod > 0
	{ 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++; // parameter plane
	      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++; // parameter plane
	  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()
{

  // input
  int plane = -1; // positive is parameter plane, negative is dynamic plane.
  uint preperiod =0; // " 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 )
  uint period = 3;

  // output
  double re,im;
  int r; // returnig value = error message

  r = find(plane, preperiod,period,re,im);
  
  
  
  if (r==0)
    {
      std::cout.precision( 15 ); 
      std::cout << "Preperiod = "<< preperiod << "\nPeriod = "<< period <<"\n";
      if (plane<0) {
	std::cout << " parameter c = "<< A << " ; " << B << ":\n";
	std::cout << " point z = "<< re << " ; " << im << ":\n";
      }
      if (plane>=0) std::cout << "point c = "<< re << " ; " << im << ":\n";
    }
   
  else { std::cout << "error \n";
    if (r==-1) std::cout << " w < 1e-100 \n";
    if (r== 2) std::cout << "(px*px + py*py > 1e100) = escaping in periodic    loop \n";
    if (r== 1) std::cout << "(px*px + py*py > 1e100) = escaping in preperiodic loop\n";
  }

  return 0;
}

Output:

Preperiod = 0
Period = 3
 parameter c = -0.965 ; 0.085:
 point z = -0.602943702541717 ; 0.0385332450804691:

algorithm 9 : zeros of qn(c)[edit | edit source]

Description:

  • "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 "allows locating zeroes of the iterated polynomial at a certain period where 4 colours meet." [6]


It can be used for:

  • on the dynamic plane
    • making parabolic checkerboard ( chessboard) when c is a parbolic parameter
    • showing center of components of Filled Julia sets, when c is a center of hyperbolic component and z=0 is a one of attracting periodic orbit

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

The four quadrants of a Cartesian coordinate system

4 colors are used because there are 4 quadrants :

  • re(z_n) > 0 and im(z_n) > 0 ( 1-st quadrant )
  • re(z_n) < 0 and im(z_n) > 0 ( 2-nd quadrant )
  • re(z_n) < 0 and im(z_n) < 0 ( 3-rd quadrant )
  • re(z_n) > 0 and im(z_n) < 0 (4-th quadrant ).

"... 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 [7]

// 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 :

  • there is no bailout test (for example for max value of abs(zn) ). It means that every point has the same number of iterations. For large n overflow is possible. Also one can't use test Iteration==IterationMax
  • Iteration limit is relatively small ( start for example IterationMax = 8 )

Compare :

  • code in Sage[8]

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 */

Newton for qn(c)[edit | edit source]

It can be used with:

  • Menu/Draw/Algorithm/Newton for Qn(c)
  • key n

Key q lets you change n ( period)


Let:

  • q_n(c) be the n-th iterate of z^2 + c with z = 0
  • N_n(c) the corresponding Newton method.

Then:

  • The zeros of q_n are the centers of period dividing n and at the same time the superattracting fixed points of N_n . There is another fixed point at infinity , which is repelling.
  • The critical points of q_n are poles of N_n . In many cases these are located in hyperbolic components of period less than n .
  • Each immediate attracting basin of N_n is connected to infinity with at least one channel. Therefore I do not see a relation to the atom domains.
  • It seems that each hyperbolic component of period n is contained completely in one attracting basin. I do not have a proof for this.

Misiurewicz point[edit | edit source]

"... 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

  • go to the parameter plane
  • mark the wake by drawing 2 external rays landing on the root point :
    • use Menu/Rays/Angles_of_the_wake
    • input fraction k/r in lowest turns ( k/r = rotation number = int_angle in turns )
  • click outside of the Mandelbrot set but inside a wake
  • find point ()
    • Menu/Points/Find point or key x
    • enter preperiod,period
  • if the wake is to small ( = you don't see a free space between rays ) then zoom in

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 )


Compare with:

critical orbit[edit | edit source]

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

  • on parameter plane choose internal angle ( key c or Manu/Poins/bifurcate)
  • go to dynamic plane
  • redraw image to remove critical orbit
  • press and keep on pressing keys Ctrl-F. Each key hit shall give 4000 iterations;

Functions :

  • in qmnshell.cpp in the function QmnShell::dFinished() see  : dplane->drawOrbit(f, x, y, 0, 4000);
  • in QmnShell::map(QAction *act) = Ctrl-F
// qmnplane.cpp by Wolf Jung (C) 2007-2017.
// part of Mandel 5.15,

void QmnPlane::drawOrbit(mndynamics *f, 
                         mdouble &x0, 
                         mdouble &y0,
                         int preiter, 
                         int plotiter)
{  if (type > 0) return;

   stop(); 
   int i, j, k;
   
   QPainter *p = new QPainter(buffer); 
   p->setPen(Qt::white);
   
   // iterate without polotting
   for (j = 0; j < preiter; j++)
        { f->iterate(x0, y0); 
          if (x0*x0 + y0*y0 > 1e6) break; }
   // iterate and plot
   for (j = 0; j < plotiter; j++)
   {  f->iterate(x0, y0); 
      if (x0*x0 + y0*y0 > 1e6) break;
      if (!pointToPos(x0, y0, i, k)) p->drawPoint(i, k);
   }

   delete p; update();
}

Tiling of parabolic Julia set interior[edit | edit source]

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

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

  • rotation number is 0 or 1/2 ( key C or MainMenu/Point/Bifurcation ), then c is the root of a satellite component.
  • Period of parent component is 1
  • map is complex quadratic polynomial ( MainMenu/ Map or key= Ctrl+0 )

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;
}

equipotential[edit | edit source]

  • key G
  • Main Menu/ Rays/Equipotewntial


"... to draw equipotential lines. ... Boundary scanning is slow and the Newton method has three problems: finding a starting point, drawing several lines around a disconnected Julia set, and choosing the number of points depending on the chosen subset of the plane. Therefore I am using boundary tracing with starting points on several lines through the image. See QmnPlane::green(). Still, sometimes not the whole line is drawn, or not all components are found."

" Hit the key g to draw the equipotential line through the current point z (outside of Kc) or c (outside of M), or to specify the potential level."

green(f, signtype, y, 10);
// qmnshell.cpp  by Wolf Jung (C) 2007-2023.
if (act == greenAct)
   {  mdouble 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);
   }



// qmnplane.cpp by Wolf Jung (C) 2007-2017.
// Mandel 5.15, 
//the following code will be rewritten completely!!!
void QmnPlane::green(mndynamics *f, int st, mdouble g, int quality,
  QColor color) //5, white
{ //trace the boundary ...
  if (g <= 0 || type) return;
  uint m = 1; f->prepare(st, nmax, m, &g); //just to set nmax
  int i, k, i0, k0, i1, k1, i2, k2, j, vert, side, sg, start, T; stop();
  QPainter *p = new QPainter(buffer); p->setPen(color);
for (T = 20 - kmax; T < kmax - 20; T += imax*2/quality)
{ for (vert = 0; vert <= 1; vert++)
  { for (side = -1; side <= 1; side += 2)
    { start = 10000; //no startpoint
      if (vert) //vertical lines from above and below
      { if (f->green(st, hmid - ph*side*kmax + pw*T, vmid - pw*side*kmax + ph*T) <= g)
           continue;
        for (j = kmax-1; j >= -kmax/2; j--)
        { if (f->green(st, hmid + ph*side*j + pw*T, vmid - pw*side*j + ph*T) <= g)
          {start = side*j; break;} }
      }
      else//vert: horizontal lines from left and right
      { if (f->green(st, hmid + pw*side*imax - ph*T, vmid + ph*side*imax - pw*T) <= g)
          continue;
        for (j = imax-1; j >= -imax/2; j--)
        { if (f->green(st, hmid + pw*side*j - ph*T, vmid + ph*side*j - pw*T) <= g)
          {start = side*j; break;} }
      }//vert
      if (start == 10000) continue; //no startpoint
      for (sg = -1; sg <= 1; sg += 2) //go in both directions
      { if (vert)
        { k0 = start; k1 = k0 + side; k2 = k1; i0 = T; i1 = T; i2 = i1 + sg;}
        else
        { i0 = start; i1 = i0 + side; i2 = i1; k0 = T; k1 = T; k2 = k1 + sg;}
        for (j = 1; j < 8*imax; j++)
        { if (i0 < -imax || i0 >= imax || k0 < -kmax || k0 >= kmax) break;
          else p->drawLine(imax + i0, kmax + k0, imax + i0, kmax + k0);
          if (f->green(st, hmid + pw*i2 + ph*k2, vmid + ph*i2 - pw*k2) <= g)
          {i = i0; k = k0; i0 = i2; k0 = k2; i2 += (i1 - i); k2 += (k1 - k);}
          else
          {i = i1; k = k1; i1 = i2; k1 = k2; i2 += (i0 - i); k2 += (k0 - k);}
        }//for j
      }//for sg
    }//for side
  }//for vert
}//for T
  delete p; update();
} //green

Argument of Phi or External angle[edit | edit source]

It can be applied with:

  • Menu/Draw/Algorithm/Argument of Phi
  • key 6


// mndlbrot.cpp  by Wolf Jung (C) 2007-2017 from Mandel 5.14
//uint mndlbrot::pixcolor(double x, double y)
if (drawmode == 6) return turn(a, b, x, y);

Here is the turn function for class mndlbrot

int mndlbrot::turn(double a, double b, double x, double y)
{  //Already checked that escaping. Requires z = c instead of z = 0
   int j; double s = 1.0, dr = 0.5, theta, u, v, r;
   if (x*x  + y*y < 1e-12) return 8; //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 8;
      u -= v; v = 2*x*y; x = u + a; y = v + b;
      //theta += s*u; Adjust 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; }
      theta += s*u; if (r > 1e18*s) break;
   }
   theta *= (.5/PI); theta -= floor(theta);
   theta *= (*temp); if (theta > 1e9) return 1;
   return 9 + (long int)(theta) % 4;
/* j = (long int)(theta); int cl = (j % 24) >> 3; if (!cl) cl = 4;
   if (j & 1) cl |= 8; return cl; //*/
} //turn

External rays[edit | edit source]

Methods :

  • Backwards iteration
  • Newton method
  • Argument tracing
// 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[edit | edit source]

Lamination[edit | edit source]

Hit o to draw an orbit portrait, which is shown as a lamination when you give the command from the parameter plane, and Ctrl+o for the precritical lamination.


If you are:

  • on the parameter plane then choosing lamination and angle 1/7 draws lamination of the unit disk in the right window
  • on the dynamic plane and choose lamination and angle 1/7 draws lamination of the dynamic plane in the right window


Lamination of:

// qmnshell.cpp  by Wolf Jung (C) 2007-2017. From Mandel 5.15
if (act == laminatAct)
   {  mdouble 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 = (mdouble)(N2);
      if (!k && p)
      { N = N1; mndAngle::conjugate(N, N2); v = ((mdouble)(N))/u; }
      u = ((mdouble)(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);
   }


//qmnplane.cpp by Wolf Jung (C) 2007-2017 from Mandel 5.15,
// It requests the conjugate angle and iterates the pair forward.
void QmnPlane::drawLamination(mdouble alpha, mdouble beta, uint n)
{  drawOrtho(alpha, beta);
   n--; if (!n) return;
   alpha *= 0.5; if (alpha < x) alpha += 0.5;
   beta *= 0.5; if (beta < x) beta += 0.5;
   drawLamination(alpha, beta, n);
   if (alpha < 0.5) alpha += 0.5; else alpha -= 0.5;
   if (beta < 0.5) beta += 0.5; else beta -= 0.5;
   drawLamination(alpha, beta, n);  
}
Orbit portrait[edit | edit source]

Steps

  • Main menu/Rays/Orbit portrait
  • qmnshell.cpp
  • portraitAct
// qmnshell.cpp  by Wolf Jung (C) 2007-2017. From Mandel 5.15
 if (act == portraitAct)
   {  mdouble 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 = mdouble(N2);
      for (k = 0; k < p; k++)
      {  dplane->drawOrtho(mdouble(N1)/x, mdouble(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);
   }
Spider algorithm[edit | edit source]

"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[edit | edit source]

In program mandel:

  • from Main Menu / Rays / Angles of the wake
  • keys: Ctrl+k

Original code from mandel 5.19

//qmnshell.cpp  by Wolf Jung (C) 2007-2023.  Defines class QmnShell

if (act == wakeAct)
   {  uint k = 1, r = 2; qulonglong n, d; mdouble x, y;
      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;
      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(((mdouble)(n))/((mdouble)(d)), 0.01*nn); return;
      }
      pplane->newtonRay(1, n, d, x, y, 5, Qt::darkMagenta, 1); n++;
      if (pplane->newtonRay(1, n, d, x, y, 5, Qt::darkMagenta, 2) <= 1
          && !f->find(1, 0, r, x, y)) pplane->setPoint(x, y);
   }
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[edit | edit source]

Backwards iteration[edit | edit source]
// 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[edit | edit source]
// 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[edit | edit source]
//  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[edit | edit source]

  • purple = RGB( 255,0,255)= HTML(ff00ff) = points which do not escape and do not fall into fixed attractor ( udecidable). For tiling of parabolic Julia set it can show repelling petals in case of low max number of iteration ( for example 100 )

In most drawmodes, Mandel uses a palette of colors 0 to 15.

In mode 3, a value from 0 to 255 is interpreted as a hue. See line 160 in qmnplane.cpp:

q->setPen(QColor::fromHsv(cl0 & 255, 255, 255));

To get grayscale, you may replace it with:

q->setPen(QColor::fromRgb(cl0, cl0, cl0));

4 colours:[9]

  • Cyan / Aqua = (0,255,255), green = ( 0,255,0), red = ( 255,0,0), blue = (0,0,255)


Change colors[edit | edit source]

This function can be started

  • by main menu File/Change colors
  • key F5
Replace a VGA-color by entering the numbers,
separated with a comma:
1   2   3   4   5   6   7   8
9 10 11 12 13 14 15 16

User may put 2 numbers separated by a comma

  • first number is the color which should change
  • second number is the color which should be

Code[edit | edit source]

classes[edit | edit source]

Classes types :

  • base
  • derived from base classes

Base classes :

  • mndynamics = abstract base class (escape to infinity)
  • mndScale  : scaling of the parameter plane at Misiurewicz points
  • mndAngle  : external angles
  • mndCombi  : kneading sequences or internal addresses
  • mndPath  : spider with legs or with path

mndynamics[edit | edit source]

mndynamics is an abstract base class (escape to infinity)

Classes derived from mndynamics

  • with critical relations and escape to infinity :
    • mndlbrot  : z^2 + c , Mandelbrot set
    • mndmulti  : z^q + c , Multibrot set
    • mndbfpoly  : cz(1 + z/q)^q (Branner-Fagella)
    • mndcubic  : various cubic polynomial families
    • mndquartic : various quartic polynomial families
    • mndmating  : various quadratic rational families (mating), some of these check only for periodicity, there is no escape.
    • mndmenage  : c(z + 0.5/z^2) , rotationally symmetric (Henriksen)
    • mndsingpert: z^2 + c/z^2 , singular perturbation of z^2
    • mndexpo  : various exponential families
    • mndtrigo  : various trigonometric families
    • mndsurge  : piecewise polynomial, modeling quasi-conformal surgery
    • mndtricorn : (z*)^q + c , anti-analytic, Tricorn or Multicorn
    • mndrealcubic: cubic poly., real-analytic para., 2 indep. cr.pt.
    • mndifs  : iterated function system for Liouville roots
  • mndsiegel = abstract derived class (persistent Siegel, two orbits, one may escape)
    • mndcubicsiegel : cubic polynomials (Zakeri-Buff-Henriksen & 2-per), derived from mndsiegel
    • mndquartsiegel : quartic polynomials (additional critical relation)
    • mndexposiegel  : exponential
    • mndtrigosiegel : trigonometric
    • mndmatesiegel  : quadratic rational
    • mndnewtonsiegel: quartic Newton mappings with Siegel cycle
    • mndnewtonpara  : cubic Newton maps with 2 roots and 1 parabolic basin
    • mndnewtonqara  : cubic Newton maps with 1 root and 2 parabolic basins
    • mndnewtonrara  : cubic Newton maps with 0 roots and 3 parabolic basins
  • mndnewton = Abstract derived class (cubic or quartic Newton mappings, with critical relations)
    • mndcubicnewton  : Newton for cubic polynomials, derived from mndnewton
    • mndquarticnewton : Newton for quartic polynomials, derived from mndnewton
  • mndherman = cubic rational with Herman ring. Derived from mndynamics (special case with escape to 0 or infinity)
  • mndhenon  : Henon mappingDerived from mndynamics (special case, escape to -infinity, no cr.pt.)
  • mndlambda  : Henon mapping. Derived from mndynamics (special case, escape to squares, no cr.pt.)

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

where :

  
 

mndAngle[edit | edit source]

The class mndAngle represents an angle (in turns) by a fraction of unsigned long long integers , which are <= 2^64 - 1 . They are normalized to the denominator 2^K * (2^P - 1) when possible, i.e., K + P <= 64 . Otherwise P = 0 . Many applications only need the static functions.

  • normalize() returns the period.
  • twice() doubles the angle modulo 1.
  • conjugate() computes the conjugate angle, returns period of the char.pt.
  • wake() computes the angles of a limb.
  • radians() gives the angle in radians.
  • lambda() computes e^{core entropy h} for dyadic angles.
  • realSpider() checks if the root is real and computes the center c .
  • truncatedTuning() truncates (u1, u2)*w to the next n-periodic angle.

mndCombi[edit | edit source]

The class mndCombi represents a *-periodic kneading sequence with period P <= 64 and the corresponding internal address. They are stored each in the bits of an unsigned long long integer . QString conversion is done by QmnCombiDialog . In the address, the first bit (2^63) is period 64 and the last bit (2^0) is period 1. Example:

    AABBAA*  is  0...01100111  

since the upper kneading sequence is AABBAAA

    1-3-4-7  is  0...01001101

The kneading sequence is always the upper one: setKneading() ignores the last entry, and the last entry from getKneading() shall be printed as a '*' . The set-functions compute the internal address from the kneading sequence or vice versa.

  • fromAngle() computes the combinatorics of a periodic angle.
  • renorm() gives the lowest period of simple renormalization.
  • count() returns the number of realizations, assuming admissibility.
  • failsBS() returns the first failure of the Bruin-Schleicher admissibility condition, see http://arXiv.org/abs/0801.4662 . 0 means that there is no evil orbit and the combinatorial data corresponds to a planar Hubbard tree with an orientation-preserving mapping, thus to a quadratic polynomial realizing these combinatorics. The corresponding angles are computed in toAngles() by backwards iteration for each period, always choosing the 1/q-limb when q > 2 . The simple algorithm fails for the first time at 1-2-6-7-13-14 , so the wake is subdivided by calling backAngles() recursively. See the explanation in mndcombi.cpp .

Global variables[edit | edit source]

  • drawmode uint. Defined in mndynamo.h
  • maxiter uint. Defined in mndynamo.h
  • a
  • b
  • A double. Defined in mndynamo.h
  • B double. Defined in mndynamo
  • sign int. Defined in mndynamo.h . "positive is parameter plane, negative is dynamic plane."
  • subtype int. Defined in mndynamo.h { subtype = subtype0; A = 0; B = 0; bailout = 16.0; temp = new double[5]; }
  • signtype = +- subtype : its modulus remembers the chosen subtype, and its sign denotes the current plane; pass its value for sign.
  • Period uint. Defined in mndynamo.h
  • *temp double. Defined in mndynamo.h
  • bailout double. Defined in mndynamo.h
  • k preperiod ( of critical value not critical point). It means that c= -2 has k=1 and p=1 so it is M_1,1 in this notation not M_2,1 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.

  • p period
  • plane parameters
    • mdouble hmid = horizontal mid = plane center x
    • vmid = vertical mid = plane center y
    • pw = plane width
    • ph = plane height
    • uint drawmode;
    • int imid ( horizontal pixel coordinate of center
    • int kmid ( vertical coorinate of the center
    • imax = horizontal maximal pixel coordinate
    • kmax = vertical maximal pixel coordinate
    • mode
    • alpha;

viewport[edit | edit source]

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

  • center of the rectangle : xmid + i ymid
  • the midpoint of the right side are given by : rewidth + i imwidth
/* 
   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 :

  • parameter plane
  • dynamic plane : -2 <= Re(z) <=2, -2 <= Im(z) <=2
// qmndemos.cpp  by Wolf Jung (C) 2007-2023
pplane = new QmnPlane(360, 360); // parameter plane
dplane = new QmnPlane(360, 360, 0); // dynamic plane

normalize[edit | edit source]

/*
qmnshell.cpp  by Wolf Jung (C) 2007-2019.  
mndynamo.h  by Wolf Jung (C) 2007-2019
mndcombi.cpp  by Wolf Jung (C) 2007-2019


   ... part of Mandel 5.17, 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.






g++ n.cpp -Wall -Wextra -lm


*/

#include <iostream> // std::cout
#include <cmath>     // sqrt
#include <limits>
#include <cfloat>
#include <bitset>
#include <string>
// #include <QString> // https://doc.qt.io/qt-5/qstring.html


using namespace std;


typedef  unsigned long long int  qulonglong;
typedef  long double  mdouble;






// mndAngle::twice
void mndAngle_twice(qulonglong &n, qulonglong &d)
{  if (n >= d) return;
   if (!(d & 1)) 
   	{ d >>= 1; if (n >= d) n -= d; return; }
   qulonglong large = 1LL; 
   large <<= 63; //avoid overflow:
   if (n < large) 
   	{ n <<= 1; if (n >= d) n -= d; return; }
   n -= large; 
   n <<= 1; 
   large -= (d - large); 
   n += large;
}





// "The function is computing the preperiod and period (of n/d under doubling map)
// and setting the denominator to  2^preperiod*(2^period - 1) if possible.
// So 1/5 becomes 3/15 and 2/10 becomes 3/15 as well.
// The period is returned as the value of the function, 
// n and d are changed ( Arguments passed to function by reference)
// and the preperiod is returned in k." (Wolf Jung)
// Question : if result is >=0 why do not use unsigneg char or unsigned int for type of result ???
// mndAngle::normalize
/*
input 
n = numerator
d = denominator

output
p = period
k = preperiod

*/
int mndAngle_normalize(qulonglong &n, qulonglong &d, int &k)
{  if (!d) return 0;
   n %= d; 
   while (!(n & 1) && !(d & 1)) 
   	{ n >>= 1; d >>= 1; }
   int p; //
   qulonglong n0;
   qulonglong n1 = n; 
   qulonglong d1 = d;
   qulonglong np;
   
   k = 0; 
   while (!(d1 & 1)) 
   	{ k++; d1 >>= 1; if (n1 >= d1) n1 -= d1; }
   n0 = n1;
   for (p = 1; p <= 65 - k; p++) 
   	{ 	mndAngle_twice(n1, d1); 
   		if (n1 == n0) break; }
   if (k + p > 64) return 0; // "The angle  N1/N2 has  preperiod + period > 64."
   
   np = 1LL << (p - 1); np--; 
   np <<= 1; 
   np++; //2^p - 1 for p <= 64
   n0 = np; 
   d >>= k; 
   n1 = d; 
   if (n1 > n0) { n1 = n0; n0 = d; }
   while (1) 
   	{ d1 = n0 % n1; 
   	if (!d1) break; 
   	n0 = n1; 
   	n1 = d1; } //gcd n1
   	
   n /= d/n1; 
   n *= np/n1; 
   d = np << k;
   return p;
}





// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// ++++++++++++++++++++ main +++++++++++++++++++++++++++++++++++++++++++++++++
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++



int main(void)
{  
	// input
	qulonglong n = 1;
	qulonglong d = 5;
	
	// output
	int k; // 
	int p; // 
	
	cout << "The input angle = n/d = " << n <<"/" << d << endl;
	p = mndAngle_normalize(n, d, k); 
	cout << "after normalisation angle = n/d = " << n <<"/" <<  d << endl;
	cout << "preperiod k = " << k <<" period p = " <<  p << endl;
	
	//if (!p) 
	//	{};
	//	else {}


	return 0;
}

compile

g++ n.cpp -Wall -Wextra -lm

run:

./a.out

result:

The input angle = n/d = 1/5
after normalisation angle = n/d = 3/15
preperiod = k = 0 period = p = 4

Images[edit | edit source]

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

opening images[edit | edit source]

Example image which was checked and modified with Mandel

One can open png image in Mandel and check it.

Image must be :

  • png type
  • square ( width = height )
  • size must be typical ( 2000 not 2001 pixels)
  • world coordinate should be :
    • parameter plane ( to do )
    • dynamic plane : -2 <= Re(z) <=2, -2 <= Im(z) <=2

Steps :

  • open Mandel
  • adjust image size ( Menu/File/Resize )
  • if you open image of dynamic plane change coordinate on parameter plane ( Menu/Points/Coordinates )
  • open image ( Menu/File/Load.png )

Then you can :

  • draw : orbit, ...
  • save

You can't:

  • zoom
  • redraw

Demos[edit | edit source]

See :

  • Main Menu/Help/
  • src code : qmndemos.cpp

Demo 2 : periodic points and bifurcations[edit | edit source]

page 7[edit | edit source]

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

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

  • start with internal angle t = Golden ratio conjugate[10]
  • compute c from t
  • go to the next t ( here y is increased by 0.01) with Go key

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; //  t = atan2()  is the argument in radians. It is increased by  0.01  so  1/628 turns
      x = 0.5*cos(y) - 0.25*cos(2*y); 
      y = 0.5*sin(y) - 0.25*sin(2*y);
      pplane->setPoint(x, y);}
}

page 10[edit | edit source]

elephants demo: there is just a rotation and rescaling on a small region instead of the non-linear transformation, like in unrolling


The method of parabolic renormalization[11] was used to show that the boundary of the Mandelbrot set has Hausdorff dimension 2,[12] and recently to construct a Julia set of positive area[13]. For another application of parabolic implosion, recall that the 1/r-limb is attached to the main cardioid at a parameter c, such that the fixed point is parabolic with combinatorial rotation number 1/r. These limbs converge to the root c = 0.25 of the main cardioid. Pierre Lavaurs has studied these limbs using the Fatou coordinates and the related horn map of Ecalle cylinders. It is conjectured[14] that the limbs converge to a limiting "elephant" under a suitable rescaling. The left image shows the 1/2-limb rotated by 90o. All limbs will be shown with increasing magnification and appropriate rotation, such that they point upwards from the boundary of the main cardioid. (Which is rotated such that it is horizontal at the root of the satellite component.) The right image shows the filled Julia set for the corresponding center of period r. Hit Return or push the Go-button to increase the period r of the limbs. This kind of scaling is not available in the main program. See page 9 of Chapter 5 for another scaling phenomenon at parabolic parameters.))


{  n++; mdouble t = 2*PI/n,
         x = 0.5*cos(t) - 0.25*cos(2*t), 
         y = 0.5*sin(t) - 0.25*sin(2*t),
         u = sin(t) - sin(2*t), 
         v = -cos(t) + cos(2*t);
      t = 2.5/(n*n); 
      u *= t; 
      v *= t; //move root below the center of the image:
      pplane->setPlane(x - 0.6*v, y + 0.6*u, u, v);
      //pplane->setPlane(x - 0.8*v - 0.2*u, y + 0.8*u - 0.2*v, 0.1*u, 0.1*v);
      pplane->draw(pf, 1, &mode);
      pf->find(1, 0, n, x, y); 
      pplane->setPoint(x, y);
   }

Here:

  • rotation number 1/r is t = 1/n
  • n is increasing by 1 at every go button click
  • setPlane(mdouble xmid, mdouble ymid, mdouble rewidth, mdouble imwidth, bool nowayBack)

Demo 3: external rays[edit | edit source]

page 1[edit | edit source]

Parametr plane : Mandelbrot set with external rays for angles where

void QmnDemoER::pFinished()
{  int j; mdouble x, y;
   if (page == 1)
   {  for (j = 0; j < 16; j++) pplane->traceRay(1, mdouble(j)/16.0, pf, x, y, 2);
      for (y = 1.35; y > 0.39; y *= 0.78)
         pplane->green(pf, 1, pf->green(1, -1.0, y), 3);
   }


Demo 5: renormalization[edit | edit source]

Videos[edit | edit source]

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

Mandelbrot[edit | edit source]

  • Bifurcation of Julia sets
  • Similarity and self-similarity of the Mandelbrot set
  • Visualization of the Thurston algorithm

Approaching the airplane root with Fatou-Lavaurs translation (1/σ ~ time)[edit | edit source]

The variable sigma is related to the multiplier of the approximately parabolic point, and 1/sigma is going to infinity when approaching the root. 1/sigma is also approximately the number of iterations from inside the approximate basin to the outside. In one video the time is proportional to 1/sigma, which means that the spiral is turning at constant speed. In the other video, the parameter c is moving at constant speed, which makes 1/sigma go faster and faster, such that it reaches infinity in a finite time.


Bifurcation of preperiodic points and rays at c = γM(5/12)[edit | edit source]

Video

The parameter c is moving to the left on the real axis. As it crosses the Misiurewicz parameter with angles 5/12 and 7/12 , the critical value has the same angles and the critical point has the angles 5/24 , 7/24 , 17/24 , 19/24 ; these rays are landing in two pairs in a different pattern for parameters right or left of this one.

Note that bifurcations happen all the time: parts of the Julia sets move down and up toward 0 and then to the left and right.


The angle  5/12  or  01p10 has  preperiod = 2  and  period = 2. The corresponding parameter ray lands at a Misiurewicz point of preperiod 2 and period dividing 2.
the landing point is c = -1.543689012692076  +0.000000000000000 i

Feigenbaum scaling[edit | edit source]

Video: Mandelbrot set and quadratic polynomials\ Slideshow of Feigenbaum scaling

Feigenbaum rescaling in the 1/2-limb

Description from the Madel demo 5 page 10:

The sequence of repeated period doubling leads to the infinitely renormalizable Feigenbaum parameter cF ≈ -1.401155189, the nested small Mandelbrot sets are scaled asymptotically by the first Feigenbaum constant  , and the small Julia sets around z = 0 (not z = c) by the second Feigenbaum constant . The decorations in the parameter plane are becoming more hairy, and converging to a dense subset of the plane. This conjecture by John Milnor was proved by Mikhail Lyubich, and Curt McMullen has shown convergence of blow-ups of the Feigenbaum Julia set.

Hit the Return key or push the Go-button a few times to rescale the parameter plane by δF around c = cF and the dynamic plane by αF around z = 0. In the main program, Feigenbaum scaling is available in the parameter plane with the key Ctrl+Alt+f. You will have to increase the number of iterations with n. The corresponding tuning of center parameters is obtained with the key Ctrl+Alt+c. Renormalization with the key r works best in the primitve case.


See also:

Rational period 2[edit | edit source]

Rational period 3[edit | edit source]

Bitransitive[edit | edit source]

Symmetric[edit | edit source]

Wishlist[edit | edit source]

  • showing window (viewport) parameters : center and radius
  • Main Menu/Rays. What about changing it to Main Menu/Curves ? It is uded to draw rays and equipotential curves
  • saving description of the image, maybe inside png image or in different txt file
  • making one image from various algorithms ( like adding boundary made by DEM or IIM to other images), layers of the image
  • scripting language form making complex images/operations
  • OpenMP , GPU
  • more description in the code and about algorithms
  • change the silent error checking to explicit information
  • options to do long operations without keeping the key preseed. For example try to draw critical orbit for bifurcation from period 1 using t = 34/89. It draws only some points.
  • dictionary
    • A Qt Linguist phrase book = qph file extension. These are human-readable XML files containing standard phrases and their translations. These files are created and updated by Qt Linguist and may be used by any number of projects and applications.
    • Qt Linguist Manual
  • not use of finite binary expansion but infinite version, because it describes better the dynamics. For example, now is: "The angle 1/4 or 01 has preperiod = 2 and period = 1.". Better version is :"The angle 1/4 or 00p1 has preperiod = 2 and period = 1."
  • description how to add a new fractal algorithm or coloring method
  • internal rays for parameter and dynamic plane


Bugs[edit | edit source]

  • "after loading certain png-images from another source, things like changing colors and drawing do not work any longer ..."

Translation[edit | edit source]

Example A Qt Linguist phrase book (qph file). It is a human-readable XML file containing standard phrases and their translations. This file is created and updated by Qt Linguist and may be used by any number of projects and applications.

<!DOCTYPE QPH>
<QPH sourcelanguage="en" language="pl_PL">
<phrase>
    <source>Arnold&apos;s Cat map </source>
    <target>Odwzorowanie kota Arnolda</target>
    <definition>In mathematics, Arnold&apos;s cat map is a chaotic map from the torus into itself, named after Vladimir Arnold, who demonstrated its effects in the 1960s using an image of a cat, hence the name.</definition>
</phrase>
<phrase>
    <source>lamination </source>
    <target>laminacja</target>
    <definition>Optional definition</definition>
</phrase>
<phrase>
    <source>map </source>
    <target>odwzorowanie </target>
    <definition>Optional definition</definition>
</phrase>
<phrase>
    <source>Boettcher map</source>
    <target>Przekształcenie Boettchera</target>
</phrase>
<phrase>
    <source>basin of attraction</source>
    <target>obszar przyciągania</target>
</phrase>
<phrase>
    <source>wake</source>
    <target>kilwater</target>
    <definition>podzbiór płaszczyny parametrów ograniczona przez 2 promienie zewnętrzne. Porównaj limb, subwake</definition>
</phrase>
<phrase>
    <source>dynamic plane </source>
    <target>płaszczyzna dynamiczna</target>
</phrase>
<phrase>
    <source>mating</source>
    <target>?parowanie?</target>
</phrase>
</QPH>

Papers with images made with Mandel[edit | edit source]

University[edit | edit source]

References[edit | edit source]

  1. Mandel by W Jung
  2. Mandelbrot set in C++
  3. Mandelbrot set in C++ and SDL
  4. ArgPhi
  5. symbolic dynamics
  6. periodicity_scan_revisited by Claude
  7. Multiplatform cpp program Mandel by Wolf Jung
  8. Formation of Escape-Time Fractals By Christopher Olah
  9. rapidtables: RGB Color Codes Chart
  10. wikipedia : Golden_ratio_conjugate
  11. Parabolic Renormalization by Hiroyuki Inou and Mitsuhiro Shishikura
  12. The Hausdorff dimension of the boundary of the Mandelbrot set and Julia sets by Mitsuhiro Shishikura
  13. Quadratic Julia Sets with Positive Area by Xavier Buff, Arnaud Cheritat
  14. Douady, Adrien : Does a Julia set depend continuously on the polynomial?
  15. Mandelbrot set and quadratic polynomials - videos by Wolf Jung