Fractals/Mathematics/Period

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

Period of the angle under doubling map[edit]

How to find period of angle under doubling map

C++ version[edit]

/*
based on : 
   mndcombi.cpp  by Wolf Jung (C) 2010. 
   http://mndynamics.com/indexp.html 
   which is the part of Mandel 5.5 
   multiplatform C++ GUI program using QT 
   on the same licence as above 
 
 "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 ???
 
*/
int normalize(unsigned long long int &n, unsigned long long int &d, int &k)
{  if (!d) return 0; // d==0 error
   n %= d; 
   while (!(n & 1) && !(d & 1)) { n >>= 1; d >>= 1; }
   int p; 
   unsigned long long int n0, n1 = n, d1 = d, np;
   k = 0; 
   while (!(d1 & 1)) { k++; d1 >>= 1; if (n1 >= d1) n1 -= d1; }
   n0 = n1;
   for (p = 1; p <= 65 - k; p++) 
	{ twice(n1, d1); 
	  if (n1 == n0) break; }
   if (k + p > 64) return 0; // more then max unsigned long long int
   np = 1LL; 
   np <<= (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;
}

Lisp version[edit]

(defun give-period (ratio-angle)
	"gives period of angle in turns (ratio) under doubling map"
	(let* ((n (numerator ratio-angle))
	       (d (denominator ratio-angle))
	       (temp n)) ; temporary numerator
 
	  (loop for p from 1 to 100 do 
		(setq temp  (mod (* temp 2) d)) ; (2 x n) modulo d = doubling)
		when ( or (= temp n) (= temp 0)) return p )))

Maxima CAS version[edit]

DoublingMap(r):=
block([d,n],
 n:ratnumer(r),
 d:ratdenom(r),
 mod(2*n,d)/d)$

/*
Tests : 
GivePeriod (1/7)
3
GivePeriod (1/14) 
0
GivePeriod (1/32767)
15
GivePeriod (65533/65535)
16

Gives 0 if :
* not periodic ( preperiodic )
* period >pMax
*/


GivePeriod (r):=
block([rNew, rOld, period, pMax, p],
      pMax:100,
      period:0,
       
      p:1, 
      rNew:DoublingMap(r),
      while ((p<pMax) and notequal(rNew,r)) do
        (rOld:rNew,
         rNew:DoublingMap(rOld),
         p:p+1
        ),
      if equal(rNew,r) then period:p,
      period
);


Haskell version[edit]

Haskell version[1]

Conversion from an integer type (Int or Integer) to anything else is done by "fromIntegral". The target type is inferred automatically

-- by Claude Heiland-Allen
-- import Data.List (findIndex, groupBy)
-- type N = Integer
-- type Q = Rational
 period :: Q -> N
 period p =
  let Just i = (p ==) `findIndex` drop 1 (iterate double p)
  in  fromIntegral i + 1

Real quadratic map[edit]

Complex quadratic map[edit]

Finding periodic points of quadratic map[edit]

Complex quadratic map

 f(z,c):=z*z+c;
 fn(p, z, c) :=
  if p=0 then z
  elseif p=1 then f(z,c)
  else f(fn(p-1, z, c),c);

Standard polynomial F_p \, which roots are periodic z-points of period p and its divisors[edit]

  • math notation :  \ F_p(z,c) = f^{(p)} _c (z) - z [2]
  • Maxima CAS function :
F(p, z, c) := fn(p, z, c) - z ;

Function for computing reduced polynomial G_p\, which roots are periodic z-points of period p without its divisors[edit]

  • math definition : G_p(z,c) = \frac{F_p(z,c) }{ \prod_{m|p,m<p} G_m(z,c)}  \,
  • Maxima function:
G[p,z,c]:=
block(
[f:divisors(p),
t:1], /* t is temporary variable = product of Gn for (divisors of p) other than p */
f:delete(p,f), /* delete p from list of divisors */
if p=1
then return(F(p,z,c)),
for i in f do 
 t:t*G[i,z,c],
g: F(p,z,c)/t,
return(ratsimp(g))
)$
Finding equations for periodic points[edit]
Standard equations[edit]

Equation for periodic point for period p \, is :

 \  F_p(z,c) = 0 \,

so in Maxima :

for i:1 thru 4 step 1 do disp(i,expand(Fn(i,z,c)=0));
1
z^2-z+c=0
2
z^4+2*c*z^2-z+c^2+c=0
3
z^8+4*c*z^6+6*c^2*z^4+2*c*z^4+4*c^3*z^2+4*c^2*z^2-z+c^4+2*c^3+c^2+c=0
4
z^16+8*c*z^14+28*c^2*z^12+4*c*z^12+56*c^3*z^10+24*c^2*z^10+70*c^4*z^8+60*c^3*z^8+6*c^2*z^8+2*c*z^8+56*c^5*z^6+80*c^4*z^6+
24*c^3*z^6+8*c^2*z^6+28*c^6*z^4+60*c^5*z^4+36*c^4*z^4+16*c^3*z^4+4*c^2*z^4+8*c^7*z^2+24*c^6*z^2+24*c^5*z^2+16*c^4*z^2+
8*c^3*z^2- z+c^8+4*c^7+6*c^6+6*c^5+5*c^4+2*c^3+c^2+c=0

Degree of standard equations is 2^p \,

These equations give periodic points for period p and its divisors. It is because these polynomials are product of lower period polynomials G_p\,:

for i:1 thru 4 step 1 do disp(i,factor(Fz(i,z,c)));
1
z^2-z+c
2
(z^2-z+c)*(z^2+z+c+1)
3
(z^2-z+c)*(z^6+z^5+3*c*z^4+z^4+2*c*z^3+z^3+3*c^2*z^2+3*c*z^2+z^2+c^2*z+2*c*z+z+c^3+2*c^2+c+1)
4
(z^2-z+c)*(z^2+z+c+1)(z^12+6*c*z^10+z^9+15*c^2*z^8+3*c*z^8+4*c*z^7+20*c^3*z^6+12*c^2*z^6+z^6+6*c^2*z^5+2*c*z^5+15*c^4*z^4+18*c^3*z^4+3*
c^2*z^4+4*c*z^4+4*c^3*z^3+4*c^2*z^3+z^3+6*c^5*z^2+12*c^4*z^2+6*c^3*z^2+5*c^2*z^2+c*z^2+c^4*z+2*c^3*z+c^2*z+2*c*z+c^6+
3*c^5+3*c^4+3*c^3+2*c^2+1)

so

 \  F_1(z,c) = z^2-z+c = G_1 \,

 \  F_2(z,c) = (z^2-z+c)*(z^2+z+c+1) = G_1 * G_2\,

 \  F_3(z,c) =  (z^2-z+c)*(z^6+z^5+3*c*z^4+z^4+2*c*z^3+z^3+3*c^2*z^2+3*c*z^2+z^2+c^2*z+2*c*z+z+c^3+2*c^2+c+1) = G_1 * G_3\,

 \  F_4(z,c) =  G_1 * G_2 * G_4 \,

 \  F_4(z,c) =  G_1 * G_5 \,

Standard equations can be reduced to equations G_p\, giving periodic points for period p without its divisors.

Reduced equations[edit]

Reduced equation is

G_p(z,c) = 0 \,

so in Maxima gives :

for i:1 thru 5 step 1 do disp(i,expand(G[i,z,c]=0));
1
z^2-z+c=0
2
z^2+z+c+1=0
3
z^6+z^5+3*c*z^4+z^4+2*c*z^3+z^3+3*c^2*z^2+3*c*z^2+z^2+c^2*z+2*c*z+z+c^3+2*c^2+c+1=0
4
z^12+6*c*z^10+z^9+15*c^2*z^8+3*c*z^8+4*c*z^7+20*c^3*z^6+12*c^2*z^6+z^6+6*c^2*z^5+2*c*z^5+15*c^4*z^4+18*c^3*z^4+
3*c^2*z^4+4*c*z^4+4*c^3*z^3+4*c^2*z^3+z^3+6*c^5*z^2+12*c^4*z^2+
6*c^3*z^2+5*c^2*z^2+c*z^2+c^4*z+2*c^3*z+c^2*z+2*c*z+c^6+3*c^5+3*c^4+3*c^3+2*c^2+1=0
Computing periodic points[edit]

Definition in wikipedia [3]


Fixed points ( period = 1 )[edit]

in Maxima :

Define c value

(%i1) c:%i;
(%o1) %i

Compute fixed points

(%i2) p:float(rectform(solve([z*z+c=z],[z])));
(%o2) [z=0.62481053384383*%i-0.30024259022012,z=1.30024259022012-0.62481053384383*%i]

Find which is repelling :

abs(2*rhs(p[1]));


Show that sum of fixed points is 1

(%i10) p:solve([z*z+c=z], [z]);
(%o10) [z=-(sqrt(1-4*%i)-1)/2,z=(sqrt(1-4*%i)+1)/2]
(%i14) s:radcan(rhs(p[1]+p[2]));
(%o14) 1

Draw points :

(%i15) xx:makelist (realpart(rhs(p[i])), i, 1, length(p));
(%o15) [-0.30024259022012,1.30024259022012]
(%i16) yy:makelist (imagpart(rhs(p[i])), i, 1, length(p));
(%o16) [0.62481053384383,-0.62481053384383]
(%i18) plot2d([discrete,xx,yy], [style, points]);

One can add point z=1/2 to the lists:

(%i31) xx:cons(1/2,xx);
(%o31) [1/2,-0.30024259022012,1.30024259022012]
(%i34) yy:cons(0,yy);
(%o34) [0,0.62481053384383,-0.62481053384383]
(%i35) plot2d([discrete,xx,yy], [style, points]);


In C :

double complex GiveAlfaFixedPoint(double complex c)
{
  double dx, dy; //The discriminant is the  d=b^2- 4ac = dx+dy*i
  double r; // r(d)=sqrt(dx^2 + dy^2)
  double sx, sy; // s = sqrt(d) = sqrt((r+dx)/2)+-sqrt((r-dx)/2)*i = sx + sy*i
  double ax, ay; // alfa = ax+ay*I
 
 // d=1-4c = dx+dy*i
 dx = 1 - 4*creal(c);
 dy = -4 * cimag(c);
 // r(d)=sqrt(dx^2 + dy^2)
 r = sqrt(dx*dx + dy*dy);
 //sqrt(d) = s =sx +sy*i
 sx = sqrt((r+dx)/2);
 sy = sqrt((r-dx)/2);
 // alfa = ax +ay*i = (1-sqrt(d))/2 = (1-sx + sy*i)/2
 ax = 0.5 - sx/2.0;
 ay =  sy/2.0;
 
 
return ax+ay*I;
}
period 2 points[edit]
(%i2) solve([(z*z+c)^2+c=z], [z]);
(%o2) [z=-(sqrt(-4*c-3)+1)/2,z=(sqrt(-4*c-3)-1)/2,z=-(sqrt(1-4*c)-1)/2,z=(sqrt(1-4*c)+1)/2]

Show that z1+z2 = -1

(%i4) s:radcan(rhs(p[1]+p[2]));
(%o4) -1

Show that z2+z3=1

(%i6) s:radcan(rhs(p[3]+p[4]));
(%o6) 1
Period 3 points[edit]

In Maxima CAS

kill(all);
remvalue(z);

f(z,c):=z*z+c; /* Complex quadratic map */
 finverseplus(z,c):=sqrt(z-c);
 finverseminus(z,c):=-sqrt(z-c); 

/* */
fn(p, z, c) :=
  if p=0 then z
  elseif p=1 then f(z,c)
  else f(fn(p-1, z, c),c);

/*Standard polynomial F_p \, which roots are periodic z-points of period p and its divisors */
F(p, z, c) := fn(p, z, c) - z ;

/* Function for computing reduced polynomial G_p\, which roots are periodic z-points of period p without its divisors*/
G[p,z,c]:=
block(
[f:divisors(p),
t:1], /* t is temporary variable = product of Gn for (divisors of p) other than p */
f:delete(p,f), /* delete p from list of divisors */
if p=1
then return(F(p,z,c)),
for i in f do 
 t:t*G[i,z,c],
g: F(p,z,c)/t,
return(ratsimp(g))
)$


 GiveRoots(g):=
 block(
 [cc:bfallroots(expand(%i*g)=0)],
 cc:map(rhs,cc),/* remove string "c=" */
 cc:map('float,cc),
 return(cc)
  )$ 


compile(all);



k:G[3,z,c]$
c:1;
s:GiveRoots(ev(k))$
s;

Finding period of the orbit[edit]

Methods :

  • direct period detection from iterations
  • the spider algorithm
  • "methods based on interval arithmetic when implemented properly are capable of finding all period-n cycles for considerable large n." (ZBIGNIEW GALIAS )[4]
  • Floyd's cycle-finding algorithm [5]

Finding period is used to :

Period of critical orbit[edit]

Finding period of critical orbit using forward iteration of critical point :


z: 0.0,

for k1:1 thru k1Max  while (abs(z)<2) do z: z*z+c; /* to remove non periodic points from image ; compute and do not draw it */

k2:1,
while (abs(zOld-z)>dc) and k<=k2Max do 
 ( z:  z * z+c, 
   k2:k2+1
  );

if (abs(z)>2) /* bailout test */ 
  then period=0 /* if escaping */
  else (if (abs(zOld-z)<dc )
             then period : k2
             else period : -1), /* no period found = (period > k2Max) or ( escaping but after >k1Max+k2Max iterations ) */

Comparisen of 2 functions for finding a period :

/*
gcc p.c -Wall -lm 
time ./a.out
numerical approximation of period of limit cycle 
Adam Majewski
*/
 
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
 
 
 
long double ER2 = 4.0L;
unsigned int jMax = 1000; // iteration max = Max period 
 
// mndynamics::period(double &a, double &b, int cycle)
// mndynamo.cpp  by Wolf Jung (C) 2007-2014
// part of Mandel 5.10 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::f(double a, double b, double &x, double &y) const
{ double u = x*x - y*y + a; y = 2*x*y + b; x = u; }
*/
unsigned int GivePeriodJung(long double cx, long double cy, long double ER2, unsigned int jMax, long double precision2, long double Zp[2]) 
{  //determine the period, then set Zp to periodic point.
   // bailout = ER2 = (EscapeRadius)^2
   unsigned int j;
  // unsigned int jMax = 500000; 
   long double x=0.0L;
   long double y=0.0L; // z 
   long double x0, y0; // z0 inside periodic orbit
   long double t; // temp
   //long double precision = 1e-16;
 
   // iterate until z fall into periodic cycle ( = limit cycle) 
   for (j = 1; j <= jMax; j++)
   { 
     if (x*x + y*y <= ER2) 
       {t = x*x - y*y + cx; 
        y = 2*x*y + cy; 
        x = t;}
       else return 0; //escaping = definitely not periodic 
   } 
   // after jMax iterations z SHOULD BE inside periodic orbit 
   x0 = x; y0 = y; // z = z0
 
   // find a period 
   for (j = 1; j <= jMax; j++)
   {  
      if (x*x + y*y <= ER2) 
        {t = x*x - y*y + cx; 
        y = 2*x*y + cy; 
        x = t;}
        else return 0; // escaping = definitely not periodic
 
     if ( (x - x0)*(x - x0) + (y - y0)*(y - y0) < precision2) // periodic 
      {   Zp[0] = x; 
          Zp[1] = y; 
         return j;  // period = j 
      }
   }
   return (2*jMax+3); // (not escaping after 2*jMax = maybe periodic but period > jMax) or  
    // (maybe escaping but slow dynamics, so need more iterations then 2*jMax) 
}
 
int SameComplexValue(long double Z1x,long double Z1y,long double Z2x,long double Z2y, long double precision)
{
    if (fabsl(Z1x-Z2x)<precision && fabs(Z1y-Z2y)<precision) 
       return 1; /* true */
       else return 0; /* false */
    }
 
/*-------------------------------*/
// this function is based on program:
// Program MANCHAOS.BAS  
// http://sprott.physics.wisc.edu/chaos/manchaos.bas
// (c) 1997 by J. C. Sprott 
//
unsigned int GivePeriodS(long double Cx,long double Cy, unsigned int iMax, long double precision, long double Zp[2])
{  
 
 
  long double Zx2, Zy2, /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
         ZPrevieousX,ZPrevieousY,
         ZNextX,ZNextY;
 
     unsigned int i; 
     unsigned int  period = iMax+3; // not periodic or period > iMax
 
     /* dynamic 1D arrays for  x, y of z points   */
    long double *OrbitX; // zx
    long double *OrbitY;  // zy 
     int iLength = iMax; // length of arrays ;  array elements are numbered from 0 to iMax-1 
  //  creates dynamic arrays and checks if it was done properly
  OrbitX = malloc( iLength * sizeof(long double) );
  OrbitY = malloc( iLength * sizeof(long double) );
  if (OrbitX == NULL || OrbitY ==NULL)
    {
      printf("Could not allocate memory \n");
      return 1; // error
    }
 
 
  Zp[0] = 0.0;
  Zp[1] = 0.0; 
 
  /* starting point is critical point  */
   ZPrevieousX=0.0;
   ZPrevieousY=0.0;
   OrbitX[0] =0.0;
   OrbitY[0] =0.0;  
   Zx2=ZPrevieousX*ZPrevieousX;
   Zy2=ZPrevieousY*ZPrevieousY;
 
   /* iterate and save points to the array */
   for (i=0;i<iMax ;i++)
        {
            ZNextY=2*ZPrevieousX*ZPrevieousY + Cy;
            ZNextX=Zx2-Zy2 +Cx;
            Zx2=ZNextX*ZNextX;
            Zy2=ZNextY*ZNextY;
            if ((Zx2+Zy2)>ER2) return 0; /* basin of atraction to infinity */
            //if (SameComplexValue(ZPrevieousX,ZPrevieousY,ZNextX,ZNextY,precision))
            //   return 1; /* fixed point , period =1 */
            ZPrevieousX=ZNextX;
            ZPrevieousY=ZNextY;
            /* */
            OrbitX[i] = ZNextX;
            OrbitY[i] = ZNextY;   
 
        };
 
    /* find   */    
     for(i=iMax-2;i>0;i--) 
      if (SameComplexValue(OrbitX[iMax-1],OrbitY[iMax-1],OrbitX[i],OrbitY[i],precision))
        { 
          Zp[0] = OrbitX[i];
          Zp[1] = OrbitY[i]; 
          period = iMax-i-1; // compute period 
          break; // the loop 
        }
 
  // free memmory
  free(OrbitX);
  free(OrbitY);
 
  return period ; 
}
 
 
 
unsigned int GivePeriodReal(long double Cx,long double Cy)
{
 // check
 
  if ( -0.75L<Cx && Cx<0.25L ) return 1;
  if ( -1.25L<Cx && Cx<-0.75L ) return 2;
  if ( -1.368089448988708L<Cx && Cx<-1.25L ) return 4; // numerical approximation = maybe wrong
  if ( -1.394040000725660L<Cx && Cx<-1.368089448988708L ) return 8; // numerical approximation = maybe wrong
 return 0; // -1.36809742955000002314
 
}
 
int main()
 
 
{
// THE REAL SLICE OF THE MANDELBROT SET 
long double CxMin = -1.4011551890L; // The Feigenbaum Point = the limit of the period doubling cascade of bifurcations 
long double CxMax = -0.74L;  
long double Cx;
long double Cy = 0.0L;
long double PixelWidth = (CxMax-CxMin)/10000.0L;
long double precisionS = PixelWidth / 100.0L;
long double precisionJ = 1e-16; 
unsigned int periodS, periodJ, periodR;
long double Zp[2]; // periodic z points on dynamic plane
unsigned int iMax = 1000000; // iteration max = Max period 
 
// text file 
FILE * fp;  // result is saved to text file 
fp = fopen("data2p10pz.txt","w"); // create new file,give it a name and open it in binary mode  
fprintf(fp," periods of attracting orbits ( c points ) on real axis of parameter plane = real slice of the Mandelbrot set  \n");
fprintf(fp," from Cmin = %.20Lf to Cmax = %.20Lf \n", CxMin, CxMax);
fprintf(fp," dC = CxMax-CxMin = %.20Lf \n", CxMax- CxMin);
fprintf(fp," PixelWidth       = %.20Lf \n", PixelWidth);
fprintf(fp," precisionS        = %.20Lf ; precisionJ =  %.20Lf\n", precisionS, sqrtl(precisionJ));
fprintf(fp," iMaxS = %u ; iMaxJ = %u\n", iMax, 2*jMax);
fprintf(fp," \n\n\n");
 
// go along real axis from CxMin to CxMax using linear scale 
Cx = CxMin;
while (Cx<CxMax)
{ 
  // compute 
  periodR = GivePeriodReal(Cx,Cy);
  periodS = GivePeriodS(Cx, Cy, iMax, precisionS, Zp);
  periodJ = GivePeriodJung(Cx, Cy, ER2, jMax, precisionJ, Zp);
  // check and save 
  if (periodR>0)
    {
      if (periodJ==periodS && periodS==periodR ) // all periods are the same and real period is known 
         fprintf(fp," c = %.20Lf ; period = %u ; \n", Cx, periodS );
         else fprintf(fp," c = %.20Lf ; period = %u ; periodS = %u ; periodJ = %u ; difference !!! \n", Cx, periodR, periodS, periodJ );
    }
    else // PeriodR==00
     {
       if (periodJ==0 && periodS==0 ) 
         fprintf(fp," c = %.20Lf ; period = %u ; \n", Cx, periodS );// all periods are the same and real period is known 
         else { if (periodS==periodJ)
                fprintf(fp," c = %.20Lf ; periodJ = periodS = %u ; \n", Cx, periodS );
                else fprintf(fp," c = %.20Lf ; periodS = %u ; periodJ = %u ; difference !!! \n", Cx, periodS, periodJ );
              }
      }  
  // info message
  printf("c = %.20Lf \n",Cx);
  // next c point 
  Cx += PixelWidth;
}
 
 fclose(fp);
 printf(" result is saved to text file \n");
return 0;
}


Non-linear scale shows bigger periods ( along real slice of Mandelbrot set ) :

/*
 
gcc p.c -Wall -lm 
time ./a.out
 
 
 
numerical approximation  of limit cycle's period  
along real slice of Mandelbrot set  
 
Adam Majewski
 
 
*/
 
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
 
// part of THE REAL SLICE OF THE MANDELBROT SET where period doubling cascade is  
long double CxMin = -1.4011552; // 1890L; // > The Feigenbaum Point = the limit of the period doubling cascade of bifurcations 
long double CxMax = 0.26L;  
long double Cx;
long double Cy = 0.0L; // constant value 
long double PixelWidth ; // = (CxMax-CxMin)/10000.0L;
//long double precisionS ; //precisionS = PixelWidth / 100.0L;//= PixelWidth / 100.0L;
long double f= 4.669201609102990671853203820466L; // The Feigenbaum delta constant 
long double precisionJ = 1e-20; 
unsigned int periodJ, periodR;
long double Zp[2]; // periodic z points on dynamic plane
 
 
long double ER2 = 4.0L;
unsigned int jMax = 5000000; // iteration max = Max period 
unsigned int iNoPeriod;
//unsigned int iMax ; //= 2*jMax; // 1000000; // iteration max = Max period 
 
// mndynamics::period(double &a, double &b, int cycle)
// mndynamo.cpp  by Wolf Jung (C) 2007-2014
// part of Mandel 5.10 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::f(double a, double b, double &x, double &y) const
{ double u = x*x - y*y + a; y = 2*x*y + b; x = u; }
 
code with small changes 
 
*/
unsigned int GivePeriodJung(long double cx, long double cy, long double ER2, unsigned int jMax, long double precision2, long double Zp[2]) 
{  //determine the period, then set Zp to periodic point.
   // bailout = ER2 = (EscapeRadius)^2
   unsigned int j;
  // unsigned int jMax = 500000; 
   long double x=0.0L;
   long double y=0.0L; // z 
   long double x0, y0; // z0 inside periodic orbit
   long double t; // temp
   //long double precision = 1e-16;
 
   // iterate until z fall into periodic cycle ( = limit cycle) 
   for (j = 1; j <= jMax; j++)
   { 
     if (x*x + y*y <= ER2) 
       {t = x*x - y*y + cx; 
        y = 2*x*y + cy; 
        x = t;}
       else return 0; //escaping = definitely not periodic 
   } 
   // after jMax iterations z SHOULD BE inside periodic orbit 
   x0 = x; y0 = y; // z = z0
 
   // find a period 
   for (j = 1; j <= jMax; j++)
   {  
      if (x*x + y*y <= ER2) 
        {t = x*x - y*y + cx; 
        y = 2*x*y + cy; 
        x = t;}
        else return 0; // escaping = definitely not periodic
 
     if ( (x - x0)*(x - x0) + (y - y0)*(y - y0) < precision2) // periodic 
      {   Zp[0] = x; 
          Zp[1] = y; 
         return j;  // period = j 
      }
   }
   return (iNoPeriod); // (not escaping after 2*jMax = maybe periodic but period > jMax) or  
    // (maybe escaping but slow dynamics, so need more iterations then 2*jMax) 
}
 
 
 
// http://classes.yale.edu/Fractals/MandelSet/MandelScalings/CompDiam/CompDiam.html
unsigned int GivePeriodReal(long double Cx,long double Cy)
{
 long double Cx0= 0.25L; 
 long double Cx1= -0.75L; 
 long double Cx2= -1.25L; 
 long double Cx3= -1.368089448988708L; // numerical approximation = maybe wrong
 long double Cx4= -1.394040000725660L; // numerical approximation = maybe wrong
 
  if ( Cx1<Cx && Cx<Cx0 ) return 1;
  if ( Cx2<Cx && Cx<Cx1 ) return 2;
  if ( Cx3<Cx && Cx<Cx2 ) return 4; // numerical approximation = maybe wrong
  if ( Cx4<Cx && Cx<Cx3 ) return 8; // numerical approximation = maybe wrong
 return 0; // -1.36809742955000002314
 
}
 
 
// try to have the same number of the pixels = n
// inside each hyperbolic component of Mandelbrot set along real axis
// width of components 
 
long double GivePixelWidth(unsigned int period, unsigned int n)
{
 
  long double w ;
  unsigned int k;
 
 switch ( period )
 {  // A SCALING CONSTANT EQUAL TO UNITY IN 1D QUADRATIC MAPS M. ROMERA, G. PASTOR and F. MONTOYA
   case      0 : w=(CxMax-CxMin)/n;      break;
   case      1 : w=1.000000000000L/n;    break; // exact value
   case      2 : w=0.310700264133L/n;    break; // numerical approximation , maybe wrong 
   case      4 : w=0.070844843095L/n;    break; // w(2*p) = w(p)/f  ; f =  Feigenbaum constant
   case      8 : w=0.015397875272L/n;    break;
   case     16 : w=0.003307721510L/n;    break;
   case     32 : w=0.000708881730L/n;    break;
   case     64 : w=0.000151841994935L/n; break;
   case    128 : w=0.000032520887170L/n; break;
   case    256 : w=0.00000696502297L/n;  break;
   case    512 : w=0.000001491696694L/n; break;
   case   1024 : w=0.000000319475846L/n; break;
   case   2048 : w=0.000000068421948L/n; break;
   case   4096 : w=0.000000015L/n;       break;
   case   8192 : w=0.000000004L/n;       break;
   case  16384 : w=0.000000001L/n;       break;
   default : if (period == 2*jMax+3)  w=(CxMax-CxMin)/10.0L; // period not found or period > jMax
                else { k=period/16384; w = 0.000000001L; while (k>2) { w /=f; k /=2;};  w /=n;} // feigenbaum scaling  
 }
 
 return w;
}
 
int main()
 
 
{
 
PixelWidth = (CxMax-CxMin)/1000.0L;
precisionJ = PixelWidth/10000000.0L;
iNoPeriod = 2*jMax+3;
 
// text file 
FILE * fp;  // result is saved to text file 
fp = fopen("data64_50ff.txt","w"); // create new file,give it a name and open it in binary mode  
fprintf(fp," periods of attracting orbits ( c points ) on real axis of parameter plane = real slice of the Mandelbrot set  \n");
fprintf(fp," from  Cmax = %.20Lf to Cmin = %.20Lf \n", CxMax, CxMin);
fprintf(fp," dC = CxMax-CxMin = %.20Lf \n", CxMax- CxMin);
fprintf(fp," non-inear scale with varied step = PixelWidth       \n");
fprintf(fp," precisionJ =  %.20Lf\n", sqrtl(precisionJ));
fprintf(fp,"  jMax = %u\n",  2*jMax);
fprintf(fp," \n\n\n");
 
// go along real axis from CxMin to CxMax using linear scale 
Cx = CxMax;
while (Cx>CxMin)
{ 
  // compute 
  //periodR = GivePeriodReal(Cx,Cy);
  periodJ = GivePeriodJung(Cx, Cy, ER2, jMax, PixelWidth/10000000.0L, Zp);
  // check and save 
   if (periodJ == iNoPeriod) 
      fprintf(fp," c = %.20Lf ; periodJ = %u ; PixelWidth = %.20LF Period not found : error !!! \n", Cx, periodJ, PixelWidth );
      else fprintf(fp," c = %.20Lf ; periodJ = %u ; PixelWidth = %.20LF \n", Cx, periodJ, PixelWidth );
  printf("c = %.20Lf ; period = %u \n",Cx, periodJ);  // info message
  // next c point 
  PixelWidth =GivePixelWidth( periodJ, 50);
  Cx -= PixelWidth;
}
 
 fclose(fp);
 printf(" result is saved to text file \n");
 
return 0;
}

References[edit]

  1. lavaurs' algorithm in Haskell with SVG output by Claude Heiland-Allen
  2. Numerical Methods for Finding Periodic Orbits at Scholarpedia
  3. wikipedia : Periodic_points_of_complex_quadratic_mappings#Finite_fixed_points
  4. Rigorous Investigations Of Periodic Orbits In An Electronic Circuit By Means Of Interval Methods by Zbigniew Galias
  5. Mandelbrot set drawing by Milan

See also[edit]