# Fractals/Iterations in the complex plane/periodic points

How to describe a cycle ( periodic orbit)

• period
• stability of a periodic Orbit ( cycle)
• (cyclic) permutation[1]
• basin of attraction

## period

Here period means period of iterated function (map).[2]

Point z has period p under f if :

${\displaystyle z:\ f^{p}(z)=z}$

or in other notation:

${\displaystyle z_{p}=z_{0}}$

The smallest positive integer p satisfying the above is called the prime period or least period of the point z.

## equation and polynomial F

So equation for periodic points:

${\displaystyle f^{p}(z)=z}$

or in other notation

${\displaystyle F^{p}(z)=f^{p}(z)-z=0}$

### degree

Degree d of polynomial is the highest of the degrees of the polynomial's monomials (individual terms) with non-zero coefficients.

Degree of rational function is the higher value from the degree of numerator and denominator

In case of iterated quadratic polynomial  :

${\displaystyle d=deg(f^{p})==2^{p}}$

# How to find periodic points for given period?

## Visual check

Visual check[3] of critical orbit helps to:

• understand dynamics
• find periodic points

## Symbolic methods

 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 Fp which roots are periodic z-points of period p and its divisors

• math notation : ${\displaystyle \ F_{p}(z,c)=f_{c}^{(p)}(z)-z}$  [4]
• ggplot[5]
• Maxima CAS function :
F(p, z, c) := fn(p, z, c) - z ;


### Function for computing reduced polynomial Gp which roots are periodic z-points of period p without its divisors

• math definition : ${\displaystyle 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 points ###### Standard equations Equation for periodic point for period ${\displaystyle p\,}$ is : ${\displaystyle \ 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 ${\displaystyle 2^{p}\,}$ These equations give periodic points for period p and its divisors. It is because these polynomials are product of lower period polynomials ${\displaystyle 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 ${\displaystyle \ F_{1}(z,c)=z^{2}-z+c=G_{1}\,}$ ${\displaystyle \ F_{2}(z,c)=(z^{2}-z+c)*(z^{2}+z+c+1)=G_{1}*G_{2}\,}$ ${\displaystyle \ 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}\,}$ ${\displaystyle \ F_{4}(z,c)=G_{1}*G_{2}*G_{4}\,}$ ${\displaystyle \ F_{5}(z,c)=G_{1}*G_{5}\,}$ Standard equations can be reduced to equations ${\displaystyle G_{p}\,}$ giving periodic points for period p without its divisors. See also : prime factors in wikipedia ###### Reduced equations ${\displaystyle 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  period degree Fp degree Gp 1 2^1 = 2 2 2 2^2=4 2 3 2^3 = 8 6 4 2^4 = 16 12 p d = 2^p ##### Computing periodic points Definition in wikipedia [6] ###### Fixed points ( period = 1 ) Equation defining fixed points :  z^2-z+c = 0  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 :  complex double GiveFixed(complex double c){ /* Equation defining fixed points : z^2-z+c = 0 z*2+c = z z^2-z+c = 0 coefficients of standard form ax^2+ bx + c a = 1 , b = -1 , c = c The discriminant d is d=b^2- 4ac d = 1 - 4c alfa = (1-sqrt(d))/2 */ complex double d = 1-4*c; complex double z = (1-csqrt(d))/2.0; return z; }  double complex GiveAlfaFixedPoint(double complex c) { // Equation defining fixed points : z^2-z+c = 0 // a = 1, b = -1 c = 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 Using non-reduced equation in Maxima CAS : (%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  Reduced equation :  ${\displaystyle z^{2}+z+c+1=0}$  so standard coefficient names of single-variable quadratic polynomial[7] ${\displaystyle ax^{2}+bx+c}$ are :  a = 1 b = 1 c = c+1  One can use the quadratic formula[8] to find the exact solution ( roots of quadratic equation):[9] ${\displaystyle z_{1,2}={\frac {-b\pm {\sqrt {d}}}{2a}}={\frac {-1\pm {\sqrt {d}}}{2}}}$  where : ${\displaystyle d=b^{2}-4ac=1-4(c+1)=-4c-3}$  Another formula[10] can be used to solve square root of complex number " ${\displaystyle s={\sqrt {d}}={\sqrt {r}}{\frac {d+r}{|d+r|}}}$ where : ${\displaystyle r=|d|}$ so ###### Period 3 points 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;  ## Numerical methods ### Newton method #### find one root ##### algorithm For given ( and fixed = constant values): • period p • parameter c • initial aproximation ${\displaystyle z_{0}}$ = initial value of Newton iteration find one periodic point ${\displaystyle z}$ of complex quadratic polynomial ${\displaystyle f}$ :  ${\displaystyle z:f^{p}(z)=z}$  where  ${\displaystyle f(z)=z^{2}+c}$  Note that periodic point ${\displaystyle z}$ is only one of p points of periodic cycle. To find periodic point one have to solve equation ( = find a single zero of a function F): ${\displaystyle f^{p}(z)-z=0}$  Lets call left side of this equation (arbitrary name) function ${\displaystyle F}$ : ${\displaystyle F_{p}(z)=f^{p}(z)-z}$  Derivative of function ${\displaystyle F}$ with respect to variable z :  ${\displaystyle F'_{p}(z)=f^{p'}(z)-1}$  It can be checked using Maxima CAS: (%i1) f:z*z+c; (%o1) z^2 + c (%i2) F:f-z; (%o2) z^2 - z + c (%i3) diff(F,z,1); (%o3) 2 z - 1  Newton iteration is : ${\displaystyle z_{n+1}=z_{n}-{\frac {F_{p}(z_{n})}{F'_{p}(z_{n})}}=z_{n}-{\frac {f_{p}(z_{n})-z_{n}}{f'_{p}(z_{n})-1}}=z_{n}-N_{p}(z_{n})}$  gives sequence of ${\displaystyle z_{n}}$ points: ${\displaystyle z_{0}}$ ${\displaystyle z_{1}=z_{0}-N_{p}(z_{0})}$ ${\displaystyle z_{2}=z_{1}-N_{p}(z_{1})}$ ${\displaystyle ...}$ ${\displaystyle z_{n+1}=z_{n}-N_{p}(z_{n})}$ where: • n is a number of Newton iteration • p is a period ( fixed positive integer) • ${\displaystyle z_{0}}$ is an initial aproximation of periodic point • ${\displaystyle z_{n}}$ is a final aproximation of periodic point • ${\displaystyle f_{p}(z_{n})}$ is computed using quadratic iteration with ${\displaystyle z_{n}}$ as a starting point • ${\displaystyle f'_{p}(z_{n})}$ is a value of first derivative . It computed using iteration • ${\displaystyle N_{p}}$ is called Newton function. Sometimes different and shorter notation is used : ${\displaystyle N_{p}(z_{n})=z_{n}-{\frac {F_{p}(z_{n})}{F'_{p}(z_{n})}}}$ ##### stoping criteria • ${\displaystyle \epsilon _{succes}}$ : predefined accuracy threshold ( succes = converged )[11] • ${\displaystyle n_{max}}$ : maximal number of iterations ( fail to converge) ${\displaystyle \epsilon =|z_{n+1}-z_{n}|}$ Usually one starts with :  ${\displaystyle \epsilon _{succes}=10^{-16}}$  for long double floating point number format  ${\displaystyle n_{max}=10^{d}}$  ##### precision Precision is proportional to the degree of polynomial F, which is proportional to the period p  ${\displaystyle d=d(p)=2^{p-1}}$  so  ${\displaystyle period\to degree\to accuracy\to numberformat}$  IEEE 754 - 2008 Common name C++ data type Base ${\displaystyle b}$ Precision ${\displaystyle r}$ Machine epsilon = ${\displaystyle b^{-(r-1)}}$ binary16 half precision short 2 11 2−10 ≈ 9.77e-04 binary32 single precision float 2 24 2−23 ≈ 1.19e-07 binary64 double precision double 2 53 2−52 ≈ 2.22e-16 extended precision _float80[12] 2 64 2−63 ≈ 1.08e-19 binary128 quad(ruple) precision _float128[12] 2 113 2−112 ≈ 1.93e-34 where : • precision = number of fraction's binary digits • accuracy < machine epsilon  precision = -log2(accuracy)  Relation between accuracy and degree d:[13] • accuracy = 1/d • accuracy = 1/(2 d log d) Compare: ##### number of iterations ##### code ###### c /* gcc -std=c99 -Wall -Wextra -pedantic -O3 m.c -lm ./a.out */ #include <complex.h> #include <math.h> #include <stdio.h> #include <stdlib.h> const double pi = 3.141592653589793; static inline double cabs2(complex double z) { return (creal(z) * creal(z) + cimag(z) * cimag(z)); } /* newton function N(z) = z - (fp(z)-z)/fp'(z)) used for computing periodic point of complex quadratic polynomial f(z) = z*z +c */ complex double N( complex double c, complex double zn , int pMax, double er2){ complex double z = zn; complex double d = 1.0; /* d = first derivative with respect to z */ int p; for (p=0; p < pMax; p++){ d = 2*z*d; /* first derivative with respect to z */ z = z*z +c ; /* complex quadratic polynomial */ if (cabs(z) >er2) break; } z = zn - (z - zn)/(d - 1) ; return z; } /* compute periodic point of complex quadratic polynomial using Newton iteration = numerical method */ complex double GivePeriodic(complex double c, complex double z0, int period, double eps2, double er2){ complex double z = z0; complex double zPrev = z0; // prevoiuus value of z int n ; // iteration const int nMax = 64; for (n=0; n<nMax; n++) { z = N( c, z, period, er2); if (cabs2(z - zPrev)< eps2) break; zPrev = z; } return z; } int main() { const double ER2 = 2.0 * 2.0; // ER*ER const double EPS2 = 1e-100 * 1e-100; // EPS*EPS int period = 3; complex double c = -0.12565651+0.65720*I ; // https://commons.wikimedia.org/wiki/File:Fatou_componenets_3.png complex double z0 = 0.0 ; // initial aproximation for Newton method; usuallly crotical point complex double zp; // periodic point zp = GivePeriodic( c , z0, period, EPS2, ER2); // printf (" c = %f ; %f \n", creal(c), cimag(c)); printf (" period = %d \n", period); printf (" z0 = %f ; %f \n", creal(z0), cimag(z0)); printf (" zp = %.16f ; %.16f \n", creal(zp), cimag(zp)); return 0; }  ###### Maxima CAS /* batch file for Maxima CAS find periodic point of iterated function solve : c0: -0.965 + 0.085*%i$ //     period of critical orbit = 2
2						1 ( alfa ) repelling 						2						1 ( beta) repelling
[0.08997933216560844 %i - 0.972330689471866, 	0.0385332450804691 %i - 0.6029437025417168, 	(- 0.0899793321656081 %i) - 0.02766931052813359, 	1.602943702541716 - 0.03853324508046944 %i] // periodic z points
[1.952970301777068, 				1.208347498712429, 				0.1882750218384916, attr					3.206813573935709] // stability ( 1)
[0.3676955262170045, 				1.460103677644584, 				0.3676955262170032, attr					10.28365329797831] // stability (2)

Newton :
z0: 0 				 zn : - 0.08997933216560859 %i - 0.02766931052813337
z0: -0.86875-0.03125*%i 	 zn: 0.08997933216560858*%i-0.9723306894718666

(%o21) (-0.08997933216560859*%i)-0.02766931052813337
(%i22) zn2:Periodic(period,c0,(-0.86875)-0.03125*%i,ER,eps)
(%o22) 0.08997933216560858*%i-0.9723306894718666

*/

kill(all); /* if you want to delete all variables and functions. Careful: this does not completely reset the environment to its initial state. */
remvalue(z);
reset(); /*  Resets many (global) system variables */
ratprint : false; /* remove "rat :replaced  " */
display2d:false$/* display in one line */ mode_declare(ER, real)$
ER: 1e100$mode_declare(eps,real)$
eps:1e-100$mode_declare( c0, complex)$
c0: -0.965 + 0.085*%i$/* period of critical orbit = 2*/ mode_declare(z0, complex)$
z0:0$mode_declare(period, integer)$
period:2$/* ----------- 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);

Fn(p,z,c):=  fn(p, z, c) - z$/* find periodic point using solve = symbolic method */ S(p,c):= ( s: solve(Fn(p,z,c)=0,z), s: map( rhs, s), s: float(rectform(s)))$

/* newton function : N(z) = z - (fp(z)-z)/(f'(z)-1) */

N(zn, c, iMax, _ER):=block(

[z, d], /* d = first derivative with respect to z */
/* mode_declare([z,zn, c,d], complex), */

d  : 1+0*%i,
z  : zn,

catch( for i:1 thru iMax step 1 do (

/*  print("i =",i,", z=",z, ", d = ",d),   */

d : float(rectform(2*z*d)), /* first derivative with respect to z */
z : float(rectform(z^2+c)), /* complex quadratic polynomial */

if ( cabs(z)> _ER)  /* z is escaping */
then  throw(z)

)), /* catch(for i */

if (d!=0) /* if not : divide by 0 error */
then  z:float(rectform( zn - (z - zn)/(d-1))),

return(z) /*  */

)$/* compute periodic point of complex quadratic polynomial using Newton iteration = numerical method */ Periodic(p, c, z0, _ER, _eps) :=block ( [z, n, temp], /* local variables */ if (p<1) then print("to small p ") else ( n:0, z: z0, /* */ temp:0, /* print("n = ", n, "zn = ", z ), */ catch( while (1) do ( z : N(z, c, p, _ER), n:n+1, /* print("n = ", n, "zn = ", z ), */ if(cabs(temp-z)< _eps) then throw(z), temp:z, if(n>64) then throw(z) ))), /* catch(while) */ return(z) )$

compile(all)$/* ------------------ run ---------------------------------------- */ zn1: Periodic(period, c0, z0, ER, eps); zn2 : Periodic(period, c0,-0.868750000000000-0.031250000000000*%i, ER, eps);  output: (%i21) zn1:Periodic(period,c0,z0,ER,eps) (%o21) (-0.08997933216560859*%i)-0.02766931052813337 (%i22) zn2:Periodic(period,c0,(-0.86875)-0.03125*%i,ER,eps) (%o22) 0.08997933216560858*%i-0.9723306894718666  #### find all roots Methods: • The iterated refinement Newton method[13] • the Ehrlich-Aberth iteration ( Mpsolve[14] ) Number: • number of all roots = degree of polynomial ( Fundamental Theorem of Algebra ) • degree = ${\displaystyle 2^{p}}$ where p is a period period degree number of all roots 1 2^1 = 2 2 2 2^2=4 4 3 2^3 = 8 8 p d = 2^p d ##### test "Viete tests[15] on the roots found: • if p is any polynomial of degree d and normalized (i.e. the leading coefficient is 1), then the sum of all roots of p must be equal to the negative of the coefficient of degree d − 1. • Moreover, the product of the roots must be equal to (−1)^d times the constant coefficient. If one root vanishes, then the product of the remaining d − 1 roots is equal to (−1)^(d−1) times the linear term. " ${\displaystyle z_{1}+z_{2}+\dots +z_{d-1}+z_{d}=-{\dfrac {a_{d-1}}{a_{d}}}}$ ${\displaystyle z_{1}z_{2}\dots z_{d}=(-1)^{d}{\dfrac {a_{0}}{a_{d}}}}$ "The sum of all roots should be 0 = the coefficient a of the degree d − 1 term of our polynomial. "[16] ${\displaystyle a_{d-1}=0}$ ${\displaystyle a_{d}=1}$ The constant term is a0. Its value depends on the period (degree). ###### code /* batch file for Maxima CAS find all periodic point of iterated function using solve */ kill(all); /* if you want to delete all variables and functions. Careful: this does not completely reset the environment to its initial state. */ remvalue(all); reset(); /* Resets many (global) system variables */ ratprint : false; /* remove "rat :replaced " */ display2d:false$   /* display in one line */

mode_declare(c0, complex)$c0: -0.965 + 0.085*%i$ /*     period of critical orbit = 2*/
mode_declare(z, complex)$mode_declare(period, integer)$
period:2  $/* ----------- 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);

Fn(p,z,c):=  fn(p, z, c) - z$/* find periodic point using solve = symbolic method */ S(p,c):= ( s: solve(Fn(p,z,c)=0,z), s: map( rhs, s), s: float(rectform(s)))$

/* gives a list of coefficients */
GiveCoefficients(P,x):=block
( [a, degree],
P:expand(P), /* if not expanded */
degree:hipow(P,x),
a:makelist(coeff(P,x,degree-i),i,0,degree));

GiveConstCoefficient(period, z, c):=(
[l:[]],
l :  GiveCoefficients( Fn(period,z,c), z),
l[length(l)]
)$/* https://en.wikipedia.org/wiki/Vieta%27s_formulas sum of all roots , here it should be zero */ VietaSum(l):=( abs(sum(l[i],i, 1, length(l))) )$

/*

it should be = (-1)^d* a0 = constant coefficient of Fn

*/

VietaProduct(l):=(
rectform(product(l[i],i, 1, length(l)))
)$compile(all); s: S(period, c0)$

vs: VietaSum(s);

vp : VietaProduct(s);

a0: GiveConstCoefficient( period, z,c0);


###### example results

The equation for period p=2 cycle is :

 z^4+2*c*z^2-z+c^2+c = 0


so d-1 = 3 coefficient is 0

${\displaystyle a_{d-1}=a_{3}=0}$

${\displaystyle a_{0}=c^{2}+c}$

Degree d :

${\displaystyle d=2^{p-1}=2^{2-1}=2^{2}=4}$


The four roots for c= -0.965 + 0.085*i are :

[0.08997933216560844*%i-0.972330689471866, 0.0385332450804691*%i-0.6029437025417168, (-0.0899793321656081*%i)-0.02766931052813359, 1.602943702541716-0.03853324508046944*%i]


their sum is :

 ${\displaystyle sum=6.938893903907228e-18\approx 0}$


their product :

  ${\displaystyle product=(-0.07904999999999943*\%i)-0.04100000000000011\approx a_{0}=c^{2}+c\approx (-0.07905*\%i)-0.04100000000000004}$


So results have passed Viete's tests

# How to compute stability of periodic point ?

## stability of periodic points (orbit) - multiplier

The multiplier (or eigenvalue, derivative) ${\displaystyle m(f^{p},z_{0})=\lambda \,}$  of a rational map ${\displaystyle f\,}$  iterated ${\displaystyle p}$  times at cyclic point ${\displaystyle z_{0}\,}$ is defined as:

${\displaystyle m(f^{p},z_{0})=\lambda ={\begin{cases}f^{p\prime }(z_{0}),&{\mbox{if }}z_{0}\neq \infty \\{\frac {1}{f^{p\prime }(z_{0})}},&{\mbox{if }}z_{0}=\infty \end{cases}}}$

where ${\displaystyle f^{p\prime }(z_{0})}$  is the first derivative of${\displaystyle \ f^{p}}$  with respect to ${\displaystyle z\,}$  at ${\displaystyle z_{0}}$ .

Because the multiplier is the same at all periodic points on a given orbit, it is called a multiplier of the periodic orbit.

The multiplier is:

• a complex number;
• invariant under conjugation of any rational map at its fixed point;[17]
• used to check stability of periodic (also fixed) points with stability index ${\displaystyle abs(\lambda ).\,}$

A periodic point is[18]

• attracting when ${\displaystyle abs(\lambda )<1;\,}$
• super-attracting when ${\displaystyle abs(\lambda )=0;\,}$
• attracting but not super-attracting when ${\displaystyle 0
• indifferent ( neutral) when ${\displaystyle abs(\lambda )=1;\,}$
• rationally indifferent or parabolic or weakly attractive[19] if ${\displaystyle \lambda \,}$  is a root of unity;
• irrationally indifferent if ${\displaystyle abs(\lambda )=1\,}$  but multiplier is not a root of unity;
• repelling when ${\displaystyle abs(\lambda )>1.\,}$

Periodic points

• that are attracting are always in the Fatou set
• that are repelling are in the Julia set;
• that are indifferent fixed points may be in one or the other.[20] A parabolic periodic point is in the Julia set.

## code

### c

/*
gcc -std=c99 -Wall -Wextra -pedantic -O3 m.c -lm

c = -0.965000 ; 0.085000
period = 2
z0 = 0.000000 ; 0.000000
zp = -0.027669 ; -0.089979
m = 0.140000 ; 0.340000
Internal Angle = 0.187833
internal radius = 0.367696

-----------------------------------------
for (int p = 1; p < pMax; p++){

if( cabs(c)<=2){
w =  MultiplierMap(c,p);
if ( iRadius <= 1.0) { // attractring
iAngle = GiveTurn(w);
period=p;
break;}
}   }

*/

#include <complex.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

const double pi = 3.141592653589793;

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

/* newton function : N(z) = z - (fp(z)-z)/f'(z)) */

complex double N( complex double c, complex double zn , int pMax, double er2){

complex double z = -zn;
complex double d = -1.0; /* d = first derivative with respect to z */
int p;

for (p=0; p < pMax; p++){

d = 2*z*d; /* first derivative with respect to z */
z = z*z +c ; /* complex quadratic polynomial */
if (cabs(z) >er2) break;

}

if ( cabs2(d) > 2)
z = zn - (z - zn)/d ;

return z;
}

/*
compute periodic point
of complex quadratic polynomial
using Newton iteration = numerical method

*/

complex double GivePeriodic(complex double c, complex double z0, int period, double eps2, double er2){

complex double z = z0;
complex double zPrev = z0; // prevoiuus value of z
int n ; // iteration
const int nMax = 64;

for (n=0; n<nMax; n++) {

z = N( c, z, period, er2);
if (cabs2(z - zPrev)< eps2) break;

zPrev = z; }

return z;
}

complex double AproximateMultiplierMap(complex double c, int period, double eps2, double er2){

complex double z;  // variable z
complex double zp ; // periodic point
complex double zcr = 0.0; // critical point
complex double d = 1;

int p;

// first find periodic point
zp =  GivePeriodic( c, zcr, period,  eps2, er2); // Find periodic point z0 such that Fp(z0,c)=z0 using Newton's method in one complex variable

// Find w by evaluating first derivative with respect to z of Fp at z0
if ( cabs2(zp)<er2) {

z = zp;
for (p=0; p < period; p++){
d = 2*z*d; /* first derivative with respect to z */
z = z*z +c ; /* complex quadratic polynomial */

}}
else d= 10000; //

return d;
}

complex double MultiplierMap(complex double c, int period, double eps2, double er2){

complex double w;
switch(period){
case 1: w = 1.0 - csqrt(1.0-4.0*c); 	break; // explicit
case 2: w = 4.0*c + 4; 			break; //explicit
default:w = AproximateMultiplierMap(c, period, eps2, er2); 	break; //

}

return w;

}

/*
gives argument of complex number in turns
https://en.wikipedia.org/wiki/Turn_(geometry)

*/

double GiveTurn( double complex z){
double t;

t =  carg(z);
t /= 2*pi; // now in turns
if (t<0.0) t += 1.0; // map from (-1/2,1/2] to [0, 1)
return (t);
}

int main() {

const double ER2  = 2.0 * 2.0; // ER*ER
const double EPS2 = 1e-100 * 1e-100; // EPS*EPS

int period = 3;
complex double c = -0.12565651+0.65720*I ; // https://commons.wikimedia.org/wiki/File:Fatou_componenets_3.png
complex double z0 = 0.0 ; // initial aproximation for Newton method; usuallly crotical point

complex double zp; // periodic point
complex double m; // multiplier

zp = GivePeriodic( c , z0, period, EPS2, ER2); //

m = MultiplierMap ( c, period,  EPS2, ER2 );

printf (" c = %f ; %f \n",  creal(c), cimag(c));
printf (" period = %d  \n",  period);
printf (" z0 = %f ; %f \n",  creal(z0), cimag(z0));
printf (" zp = %.16f ; %.16f \n",  creal(zp), cimag(zp));
printf (" m = %.16f ; %.16f \n",  creal(m), cimag(m));
printf (" Internal Angle = %.16f \n internal radius = %.16f \n", GiveTurn( m), cabs(m));

return 0;
}


### Maxima CAS

/*
Computes periodic points
z: z=f^n(z)
and its stability index ( abs(multiplier(z))

for complex quadratic polynomial :
f(z)=z*z+c

batch file for Maxima cas
fraktal.republika.pl
*/

/* basic function  = monic and centered complex quadratic polynomial

*/
f(z,c):=z*z+c $/* iterated function */ fn(n, z, c) := if n=1 then f(z,c) else f(fn(n-1, z, c),c)$
/* roots of Fn are periodic point of  fn function */
Fn(n,z,c):=fn(n, z, c)-z $/* gives irreducible divisors of polynomial Fn[p,z,c] which roots are periodic points of period p */ G[p,z,c]:= block( [f:divisors(p),t:1], g, f:delete(p,f), if p=1 then return(Fn(p,z,c)), for i in f do t:t*G[i,z,c], g: Fn(p,z,c)/t, return(ratsimp(g)) )$

/* use :
roots:GiveRoots_bf(G[3,z,c]);

*/
GiveRoots_bf(g):=
block(

[cc:bfallroots(expand(g)=0)],
cc:map(rhs,cc),/* remove string "c=" */
return(cc)
)$/* -------------- main ----------------- */ load(cpoly); /* c should be inside Mandelbrot set to make attracting periodic points */ /* c:0.74486176661974*%i-0.12256116687665; center of period 3 component */ coeff:0.37496784+%i*0.21687214; /* -0.11+0.65569999*%i;*/ disp(concat("c=",string(coeff))); for p:1 step 1 thru 5 do /* period of z-cycle */ ( disp(concat("period=",string(p))), /* multiplier */ define(m(z),diff(fn(p,z,c),z,1)), /* G should be found when c is variable without value */ eq:G[p,z,c], /* here should be c as a symbol not a value */ c:coeff, /* put a value to a symbol here */ /* periodic points */ roots[p]:GiveRoots_bf(ev(eq)), /* ev puts value instead of c symbol */ /* multiplier of periodic points */ for z in roots[p] do disp(concat( "z=",string(float(z)),"; abs(multiplier(z))=",string(cabs(float(m(z)))) ) ), kill(c) )$

/*



results:

c:0.74486176661974*%i-0.12256116687665
periodic z-points:
supperattracting cycle (multiplier=0):
3.0879115871966895*10^-15*%i-1.0785250533182843*10^-14=0
0.74486176661974*%i-0.12256116687665=c
0.5622795120623*%i-0.66235897862236=c*c+c

=== repelling cycle
0.038593416246119-1.061217891258986*%i,
0.99358185708018-0.90887310643243*%i,

0.66294971900937*%i-1.247255127827274]
--------------------------------------------

for c:-0.11+0.65569999*%i;
abs(multiplier) = 0.94840556944351

z1a=0.17434140717278*%i-0.23080799291989,
z2a=0.57522120945524*%i-0.087122596659277,
z3a=0.55547045915754*%i-0.4332890929585,

abs(multiplier) =15.06988063154376
z1r=0.0085747730232722-1.05057855305735*%i,
z2r=0.95628667892607-0.89213756745704*%i,
z3r=0.63768304472883*%i-1.213641769411674


# How to compute period of given periodic point ?

How is it possible to calculate the true period?

• Robert Munafo's polygon iteration method described in Mu-Ency
• Knighty's Taylor ball iteration method

Both find the lowest period nucleus in a region. If there is no nucleus in the region (either 100% exterior, or 100% interior not surrounding a nucleus) then I don't know what happens.

The polygon method can be adapted to pre-periodic Misiurewicz points. Probably the Taylor ball method can be too. A version of the ball method can be adapted to Burning Ship etc using Jacobians, where the polygon method can sometimes fail due to early folds.

int GivePeriod(complex double c){

if (cabs(c)>2.0) {return 0;} // exterior
if (cabs(1.0 - csqrt(1.0-4.0*c))<=1.0 ) {return 1;} // main cardioid
if (cabs(4.0*c + 4)<=1.0){return 2;} // period 2 component

int period =  GivePeriodByIteration(c);

if ( period < 0) // last chance
{
iUnknownPeriod +=1;
//fprintf(stderr,"unknown period = %d < 0 c = %.16f %+.16f \n", period, creal(c), cimag(c) ); // print for analysis
//period = m_d_box_period_do(c, 0.5, iterMax_LastIteration);
//fprintf(stderr,"box period = %d  \n", period);

}

// period > 0 means is periodic
// period = 0 means is not periodic = exterior = escaping to infinity
// period < 0 means period not found, maybe increase global variable iterMax_Period ( see local_setup)

return period;

}


## marcm200

A heuristics to get to a point where the period detection might work correctly could also be: Compute the rolling derivative from the start on, i.e 2^n*product over all iterates. If the orbit approaches an attracting cycle of  length p, cyclic portions of that overall product will become smaller than one in magnitude. And repeatedly multiplying those values will eventually let the |rolling derivative| fall below any ever so small threshold. So if one  were to define beforehand such a small value like 10^-40, if the rolling derivative is first below, then doing a periodicity check by number equaltiy up to an epsilon should give the correct period.  I only use it with unicritical polynomials so chances are that during iteration one does not fall onto a critical point and then the rolling derivative gets stuck at 0 forever (so for z^2+c, starting at c is usually feasible, as +-sqrt(-c) is then almost always irrationally complex and not representable in software. However, with bad luck, its approximation might still lead to computer-0 within machine-epsilon).


( marcm200)

## Knighty's Taylor ball iteration method

Finding the roots of the SA polynomial as suggested by knightly is a totally different approach. Which of the many polynomial root finding algorithms do you use?

ball interval arithmetic" is also called "midpoint-radius interval".

Divide and conquer, using ball interval arithmetic to isolate the roots and Newton-Raphson to refine them.
My approach to Ball interval arithmetic is  very simple:
To a given number z we associate a disc representing an estimation of the errors. A ball interval may be represented like this: [z] = z + r [E]

Here z is the center of the ball, r is the radius and [E] is the error symbol representing the unit disc centred at 0.

Being not a mathematician, be aware that these definitions are very very sloppy. In particular, I don't care about rounding stuff.

One can also think about [E] as a function in some "mute" variables. for example, for complex numbers we can write:
[E] = f(a,b) = a * exp(i *b), where a is in [0,1] and b any real number.
it can also be written like this:
[E]= [0,1] * exp(i * [-inf, +inf])

A nice property of [E] is that raising it to a (positive) power will yield [E].
Another one is: [E] = -[E]
... Another one: z[E] = |z| [E]

Some basic operations on ball interval are:
[v] + [w] = v + rv [E] + w + rw [E] = (v+w) + (rv+rw) [E]
-Subtraction:
[v] - [w] = (v-w) + (rv+rw) [E]; because [E] will absorb the '-' sign.
-Multiplication:
[v] [w] = (v + rv [E]) (w + rw [E]) = v*w + v*rw [E] + w*rv [E] + rv*rw [E]? = v*w + (|v|*rw + |w|*rv + rv*rw) [E]; because [E]? = [E]
The division is more complicated: It is not always defined (when 0 is inside the ball) and needs special care.
-Raising to a power:
use successive multiplications... but that may give too large balls... there are other methods.

With these operations one can evaluate polynomials using Horner method for example. It can also be used to iterate zn+1 = zn? + c for example and be used to find nucleus and its period instead of the box method.


You can take a look at the codes attached to below programs

See this note