# Period of the angle under doubling mapEdit

How to find period of angle under doubling map

## C++ versionEdit

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

(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 versionEdit

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 versionEdit 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 mapEdit # Complex quadratic mapEdit ## Finding periodic points of quadratic mapEdit 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 divisorsEdit • 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 divisorsEdit • math definition : $G_p(z,c) = \frac{F_p(z,c) }{ \prod_{m|p,m • 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 pointsEdit
###### Standard equationsEdit

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 equationsEdit

$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 pointsEdit

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]
(%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 pointsEdit
(%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 pointsEdit

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 orbitEdit

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 orbitEdit

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

#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

*/

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