Fractals/Mathematics/Period

external angle edit

All rays landing at the same periodic point have the same period: the common period of the rays is a (possibly proper) multiple of the period of their landing point; therefore one distinguishes: the ray period from the orbit period.[1]


The doubling map edit

How to find the period of angle under doubling map

  • visual
  • numerical
    • read period from denominator of decimal fraction ( reduced rational fraction m/n )
    • find period/preperiod in the binary expansion ( binary sequence)
    • read it from the itinerary of angle under doubling map

Period of binary expansion of reduced rational fraction m/n is equal to the multiplicative order of 2 modulo n:


 


C version edit

double precision: forward and inverse doubling map edit

/*


doubling map 
en.wikipedia.org/wiki/Dyadic_transformation

2*t mod 1 



how to invert doubling map  



Inverse of doubling map is multivalued function: 2 preimages
t/2 and  (t+1)/2
to choose proper preimage one needs extra information from forward iteration = itinerary





itinerary : list of symbols 
for  coding the orbits of a given dynamical system 
by partitioning the space X and forming an itinerary 



http://www.maths.qmul.ac.uk/~sb/cf_chapter4.pdf


see also how to convert proper decimal fraction to binary fraction

commons.wikimedia.org/wiki/File:Binary_decomposition_of_dynamic_plane_for_f0(z)_%3D_z%5E2.png



---------- git -------------------- 
cd existing_folder
git init
git remote add origin git@gitlab.com:adammajewski/doubling_map.git
git add .
git commit -m "Initial commit"
git push -u origin master



*/
#include <stdio.h> // printf
#include <math.h> // fabs


#define iMax  8 //

int main(){

	double t0 ;
	double t ;
	double ti; // final t after iMax iterations
	double tr; //
	double dt;
	
	int itinerary[iMax]= {0};
	
	
	
	
	
	int i;
	
	
	t0 = (double) 1/7;
	t = t0;
	
	// check the input : it should be   0.0 <= t < 1.0
	if (t>1.0) {printf("t is > 1.0\n"); return 1;}
	if (t<0.0) {printf("t is < 0.0\n"); return 1;}
	
	
	printf("forward iteration of doubling map\n");
	for(i=0; i<iMax; i++){
	        
	        printf("t%d = %f", i, t);
	        // https://en.wikipedia.org/wiki/Dyadic_transformation
		t = t*2.0; // doubling 
		if (t>1.0) {	
			
			itinerary[i]= 1;
			t = t - 1.0; 
			printf(" wrap\n");} // modulo 1 
			else printf("\n");
		}
	printf("t%d = %f\n", i, t);	
		
	//		
	ti = t;	
	
	printf("\nbackward iteration of doubling map = halving map \n");
	
	//
	for(i=iMax; i>0; i--){ // reverse counting
	        
	        printf("t%d = %f", i, t);
	        
	        if (itinerary[i-1]==1) { // i-1 !!! 
	        	 
	        	t = t + 1.0; 
	        	printf(" unwrap\n");} // modulo 1 
			else printf("\n");
	        t = t/2.0; // halving 
		
		}
	printf("t%d = %f\n", i, t);	
		
			
         tr = t;		
		
		
	//
	printf("\n\nresults \n");
	printf("t0 = %f\n", t0);	
	printf("t%d = %f\n", iMax, ti);
	
	dt = fabs(t0- tr);
	printf("tr = %f\n", tr);
	printf("dt = fabs(t0- tr) = %f\n", dt );
	printf("\nitinerary:\n");
	for(i=0; i<iMax; i++) printf("itinerary[%d] = %d \n", i, itinerary[i]);
	
	
	printf("\ndecimal %f has binary expansion = 0.", t0);
	for(i=0; i<iMax; i++) printf("%d", itinerary[i]);
	printf("\n");
	
	if (dt < 0.0000000001) printf("program works good !\n");
		else printf("program fails !\n");
	
		

	return 0;}

arbitrary precision edit

// gcc d.c -lgmp -Wall

#include <stdio.h>
#include <gmp.h>

//  a multiple precision integer, as defined by the GMP library. The C data type for such integers is mpz_t

int print_z(mpz_t  z, int base, char *s){
  printf("%s= ", s);
  mpz_out_str (stdout, 10, z);
  printf (" for base = %d\n", base);
  return 0;
}

// rop = (2*op) mod 1
// wikipedia : dyadic_transformation or doubling map
void mpq_doubling(mpq_t rop, const mpq_t op)
{
  mpz_t n; // numerator
  mpz_t d; // denominator
  mpz_inits(n, d, NULL);

 
  //  
  mpq_get_num (n, op); // 
  mpq_get_den (d, op); 
 
  // n = (n * 2 ) % d
  mpz_mul_ui(n, n, 2); 
  mpz_mod( n, n, d);
  
      
  // output
  mpq_set_num(rop, n);
  mpq_set_den(rop, d);
    
  mpz_clears(n, d, NULL);

}

int main ()
{

        int i;
        //
        unsigned long int e = 89; // exponent is also a period of doubling map 
        unsigned long int b = 2;
       
        // arbitrary precision variables from GMP library
        mpz_t  n ; // numerator of q
        mpz_t  d ; // denominator of q
        mpq_t q;   // rational number q = n/d

        // init and set variables 
        mpz_init_set_ui(n, 1);

        // d = (2^e) -1 
        // http://fraktal.republika.pl/mset_external_ray.html
        mpz_init(d);
        mpz_ui_pow_ui(d, b, e) ;  // d = b^e
        mpz_sub_ui(d, d, 1);   // d = d-1

        //   q = n/d
        mpq_init (q); //
        mpq_set_num(q,n);
        mpq_set_den(q,d);
        mpq_canonicalize (q); // It is the responsibility of the user to canonicalize the assigned variable before any arithmetic operations are performed on that variable.

        // print 
        //print_z(d, 10, "d ");
        //print_z(n, 10, "n ");
        gmp_printf ("q = %Qd\n",q); //  
        
        // 
        for (i=0; i<(1+2*e) ; i++){  
          mpq_doubling(q, q);
          gmp_printf ("q = %Qd\n",q); // 
        }   
         
       
        
        
        // clear memory
        mpq_clear (q);
        mpz_clears(n, d, NULL);
        
        
        return 0;
}

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[3]

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


What c parameters are hard cases when computing period ? edit

  • boundary and near boundary points ( parabolic, Siegel, Cremer)
  • "for example, c=1/4−10^{−10} takes over 800000 iterations to reach a fixed point in double precision floating point, but 20 iterations of Newton's method suffice to reach a fixed point. The two fixed points are slightly different, but both are fixed (each to themselves) in double precision." - Claude Heiland-Allen
  • " You may need very small epsilon and very large n, otherwise for example c=−3/4+10{−10} will probably give an incorrect period of 2 instead of the correct period of 1, which error will compound to an incorrect interior distance estimate (for example distance 3.8e-8 with your method (epsilon 1e-12, n 79,573,343) instead of 2e-10 with my method). " Claude Heiland-Allen
  • "one can iterate a "whole lot" of times, then see if the last iteration equals any of the previous iterations. However, such methods are never foolproof. Typically, they fail when used on points that are comparatively close to the boundary of their mu-atom. In such cases the algorithm will usually "find" an integer multiple of the actual period. This happens because as you get close to a child mu-atom, the points approaching the limit cycle converge in a spiral or pinheel-like manner. For example, imagine the period is actually 5, and our point is close to a child mu-atom with period 35. If you look at each 5th iterate, you'll typically find a pattern that spirals in on a limit point, taking 7 steps to go each time around the spiral. Because of the spiral pattern, the 7th step (35th iteration) is a lot closer than the 1st step (5th iteration), and simple period-finding algorithms will catch the closer one." Robert Munafo From the Mandelbrot Set Glossary and Encyclopedia


Solution:

finding periodic points for a given period edit

checking the period using position of parameter = periodicity checking edit


Jordan curve method edit


There is a greatly enhanced period detection algorithm (using the jordan curve method) that locates islands in the view even if they are far too small to see at the current magnification level. Robert Munafo


C code from mandelbrot-numerics library by Claude Heiland-Allen

// mandelbrot-numerics -- numerical algorithms related to the Mandelbrot set
// Copyright (C) 2015-2018 Claude Heiland-Allen
// License GPL3+ http://www.gnu.org/licenses/gpl.html

/*
gcc b.c -std=c99 -Wall -Wextra -pedantic -lm
*/
#include <stdio.h> // fprintf
#include <stdbool.h>
#include <stdint.h>
#include <stdlib.h>
#include <complex.h>
#include <math.h>


static inline int sgn(double z) {
  if (z > 0) { return  1; }
  if (z < 0) { return -1; }
  return 0;
}

static inline bool odd(int a) {
  return a & 1;
}

static inline double cabs2(double _Complex z) {
  return creal(z) * creal(z) + cimag(z) * cimag(z);
}

static inline bool cisfinite(double _Complex z) {
  return isfinite(creal(z)) && isfinite(cimag(z));
}

//****************************************************************
//**************** Box period *************************************
//*****************************************************************



static double cross(double _Complex a, double _Complex b) {
  return cimag(a) * creal(b) - creal(a) * cimag(b);
}

static bool crosses_positive_real_axis(double _Complex a, double _Complex b) {
  if (sgn(cimag(a)) != sgn(cimag(b))) {
    double _Complex d = b - a;
    int s = sgn(cimag(d));
    int t = sgn(cross(d, a));
    return s == t;
  }
  return false;
}

static bool surrounds_origin(double _Complex a, double _Complex b, double _Complex c, double _Complex d) {
  return odd
    ( crosses_positive_real_axis(a, b)
    + crosses_positive_real_axis(b, c)
    + crosses_positive_real_axis(c, d)
    + crosses_positive_real_axis(d, a)
    );
}

typedef struct  {
  double _Complex c[4];
  double _Complex z[4];
  int p;
} m_d_box_period ;

m_d_box_period *m_d_box_period_new(double _Complex center, double radius) {
  m_d_box_period *box = (m_d_box_period *) malloc(sizeof(*box));
  if (! box) {
    return 0;
  }
  box->z[0] = box->c[0] = center + ((-radius) + I * (-radius));
  box->z[1] = box->c[1] = center + (( radius) + I * (-radius));
  box->z[2] = box->c[2] = center + (( radius) + I * ( radius));
  box->z[3] = box->c[3] = center + ((-radius) + I * ( radius));
  box->p = 1;
  return box;
}

void m_d_box_period_delete(m_d_box_period *box) {
  if (box) {
    free(box);
  }
}

bool m_d_box_period_step(m_d_box_period *box) {
  if (! box) {
    return false;
  }
  bool ok = true;
  for (int i = 0; i < 4; ++i) {
    box->z[i] = box->z[i] * box->z[i] + box->c[i];
    ok = ok && cisfinite(box->z[i]);
  }
  box->p = box->p + 1;
  return ok;
}

bool m_d_box_period_have_period(const m_d_box_period *box) {
  if (! box) {
    return true;
  }
  return surrounds_origin(box->z[0], box->z[1], box->z[2], box->z[3]);
}

int m_d_box_period_get_period(const m_d_box_period *box) {
  if (! box) {
    return 0;
  }
  return box->p;
}

int m_d_box_period_do(double _Complex center, double radius, int maxperiod) {
  m_d_box_period *box = m_d_box_period_new(center, radius);
  if (! box) {
    return 0;
  }
  int period = 0;
  for (int i = 0; i < maxperiod; ++i) {
    if (m_d_box_period_have_period(box)) {
      period = m_d_box_period_get_period(box);
      break;
    }
    if (! m_d_box_period_step(box)) {
      break;
    }
  }
  m_d_box_period_delete(box);
  return period;
}



//****************************************************************
//**************** Main  *************************************
//*****************************************************************


int main(void){
	
	complex double c = -0.120972062945854+0.643951893407125*I; //
	complex double dc = 1.0;
	int maxiters = 1000;
	
	/* find the period of a nucleus within a large box uses Robert P. Munafo's Jordan curve method  */
	int p = m_d_box_period_do(c, 4.0 * cabs(dc), maxiters);
	fprintf(stdout, "c = %.16f %+.16f \t period = %d\n", creal(c), cimag(c),p);
      

	return 0;

}

Finding period of the orbit edit

/* mndynamo.cpp  by Wolf Jung (C) 2007-2015.  Defines classes:
   mndynamics, mndsiegel, mndcubesiegel, mndquartsiegel, mndexposiegel,
   mndtrigosiegel, mndexpo, mndtrigo, mndmatesiegel, mndmating, mndsingpert,
   mndherman, mndnewtonsiegel, mndnewton, mndcubicnewton, mndquarticnewton

   These classes are part of Mandel 5.13, 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.
*/
uint mndynamics::period(double &a, double &b, int cycle) // = 0
{  //determine the period, if cycle then set a, b to periodic point.
   uint j; double x, y, x0, y0; critical(a, b, x, y);
   for (j = 1; j <= 1000; j++)
   { if (x*x + y*y <= bailout) f(a, b, x, y); else return 0; }
   x0 = x; y0 = y;
   for (j = 1; j <= 1024; j++)
   {  if (x*x + y*y <= bailout) f(a, b, x, y); else return 0;
      if ( (x - x0)*(x - x0) + (y - y0)*(y - y0) < 1e-16)
      {  if (cycle) { a = x; b = y; }
         return j;
      }
   }
   return 10000;
}

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 :

Maxima CAS edit


/*
 b batch file for maxima
*/

kill(all);
remvalue(all);

/* =================== functions ============ */

/*
 https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/qpolynomials
complex quadratic polynomial
*/

f(z,c):=z*z+c;

/* iterated map */

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

/* 
https://en.wikibooks.org/wiki/Fractals/Mathematics/Period#Complex_quadratic_map

period of c under complex quadratic polynomial f
*/
GivePeriod(c):=block(

[z: 0.0,
 k2Max:200, /* to big values couse bind stack overflow */
 k1Max:100,
 ER:2.0,
 dMax:0.0003, /* if too low then gives smaller period then */
 period:0 /* no period found = (period > k2Max) or ..... ????  */

],

/* to remove non periodic points , iterate and do not use it */
for k1:1 thru k1Max  do 
  (z: f(z,c),
   if  (cabs(z)>ER) then  (period : -1, /* escaping */ 
                           go(exit))
  ),

/* after k1Max iterations z SHOULD BE inside periodic orbit   */   
zOld:z,

for k2:1 thru k2Max  do  
 ( z: f(z,c), 
   if  (cabs(z)>ER) then  (period : -1,  go(exit)), /* escaping */
   if  (cabs(zOld-z)<dMax) then  (period : k2,  go(exit)) /* periodic */
  ),

exit,

return(period)

)$

/*

Tests :

good
G(0)
G(-1.75)
G(-1.77)
G(-1.778)
G(-0.155+0.75*%i)    period = 3
G(-1.7577+0.0138*%i)     period = 9
G(-0.615341000000000  +0.423900000000000*%i);    period = 7
G(-1.121550281113895  +0.265176187855967*%i);    period = 18)

Tuning : 
0 period ( when true period > k2Max
G(-1.119816337988403  +0.264371090395906*%i); 
gives 0 when k2Max =100
gives 108 when dMax = 0.003
but  true period = 162  ( set k2Max = 200 and dMax= 0.0003

-------------------
G(0.37496784+%i*0.21687214);
http://fraktal.republika.pl/period.html
gives 0

*/

G(c):=GivePeriod(c);

compile(all);
/* -------- input value ------ */

c : 0.25  +0.5 * %i$

/* ============== compute ===============  */
p:GivePeriod(c)$
p;

c edit

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

c# edit

User:Simpsons_contributor: periodicity_checking

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Mandelbrot2
{
    public struct MandelbrotData
    {
        private double px;
        private double py;

        private double zx;
        private double zy;

        private long iteration;
        private bool inSet;
        private HowFound found;

        public MandelbrotData(double px, double py)
        {
            this.px = px;
            this.py = py;
            this.zx = 0.0;
            this.zy = 0.0;
            this.iteration = 0L;
            this.inSet = false;
            this.found = HowFound.Not;
        }

        public MandelbrotData(double px, double py,
                              double zx, double zy,
                              long iteration,
                              bool inSet,
                              HowFound found)
        {
            this.px = px;
            this.py = py;
            this.zx = zx;
            this.zy = zy;
            this.iteration = iteration;
            this.inSet = inSet;
            this.found = found;
        }

        public double PX
        {
            get { return this.px; }
        }

        public double PY
        {
            get { return this.py; }
        }

        public double ZX
        {
            get { return this.zx; }
        }

        public double ZY
        {
            get { return this.zy; }
        }

        public long Iteration
        {
            get { return this.iteration; }
        }

        public bool InSet
        {
            get { return this.inSet; }
        }

        public HowFound Found
        {
            get { return this.found; }
        }
    }

    public enum HowFound { Max, Period, Cardioid, Bulb, Not }

    class MandelbrotProcess
    {
        private long maxIteration;
        private double bailout;

        public MandelbrotProcess(long maxIteration, double bailout)
        {
            this.maxIteration = maxIteration;
            this.bailout = bailout;
        }

        public MandelbrotData Process(MandelbrotData data)
        {
            double zx;
            double zy;
            double xx;
            double yy;

            double px = data.PX;
            double py = data.PY;
            yy = py * py;

            #region Cardioid check

            //Cardioid
            double temp = px - 0.25;
            double q = temp * temp + yy;
            double a = q * (q + temp);
            double b = 0.25 * yy;
            if (a < b)
                return new MandelbrotData(px, py, px, py, this.maxIteration, true, HowFound.Cardioid);

            #endregion

            #region Period-2 bulb check

            //Period-2 bulb
            temp = px + 1.0;
            temp = temp * temp + yy;
            if (temp < 0.0625)
                return new MandelbrotData(px, py, px, py, this.maxIteration, true, HowFound.Bulb);

            #endregion

            zx = px;
            zy = py;

            int check = 3;
            int checkCounter = 0;

            int update = 10;
            int updateCounter = 0;

            double hx = 0.0;
            double hy = 0.0;

            for (long i = 1; i <= this.maxIteration; i++)
            {
                //Calculate squares
                xx = zx * zx;
                yy = zy * zy;

                #region Bailout check

                //Check bailout
                if (xx + yy > this.bailout)
                    return new MandelbrotData(px, py, zx, zy, i, false, HowFound.Not);

                #endregion

                //Iterate
                zy = 2.0 * zx * zy + py;
                zx = xx - yy + px;

                #region Periodicity check

                //Check for period
                double xDiff = Math.Abs(zx - hx);
                if (xDiff < this.ZERO)
                {
                    double yDiff = Math.Abs(zy - hy);
                    if (yDiff < this.ZERO)
                        return new MandelbrotData(px, py, zx, zy, i, true, HowFound.Period);
                } //End of on zero tests

                //Update history
                if (check == checkCounter)
                {
                    checkCounter = 0;

                    //Double the value of check
                    if (update == updateCounter)
                    {
                        updateCounter = 0;
                        check *= 2;
                    }
                    updateCounter++;

                    hx = zx;
                    hy = zy;
                } //End of update history
                checkCounter++;

                #endregion
            } //End of iterate for

            #region Max iteration

            return new MandelbrotData(px, py, zx, zy, this.maxIteration, true, HowFound.Max);

            #endregion
        }

        private double ZERO = 1e-17;
    }
}

References edit

  1. H. Bruin and D. Schleicher, Symbolic dynamics of quadratic polynomials, Institut Mittag-Leffler, The Royal Swedish Academy of Sciences, 7.
  2. math stackexchange question: period-of-a-finite-binary-sequence
  3. lavaurs' algorithm in Haskell with SVG output by Claude Heiland-Allen
  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