分形/數學/週期
< 分形
所有落在相同週期點的射線都具有相同的週期：射線的公共週期是其著陸點週期的（可能適當的）倍數；因此，區分了：射線週期與軌道週期。[1]
- 視覺
- 數值
簡化有理分數 m/n 的二進位制展開的週期等於2 模n的乘法階
-
角度 5/31 在倍增對映下的軌道
/*
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;}
// 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;
}
/*
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;
}
(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 )))
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 版本[3]
從整數型別（Int 或 Integer）到任何其他型別的轉換是透過“fromIntegral”完成的。目標型別是自動推斷的
-- 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
- 邊界和近邊界點（拋物線、西格爾、克雷默）
- "例如，c=1/4−10^{−10} 需要超過 800000 次迭代才能在雙精度浮點數中達到固定點，但牛頓法的 20 次迭代就足以達到固定點。這兩個固定點略有不同，但在雙精度中都是固定的（每個都固定到自身）。” - 克勞德·海蘭德-艾倫
- "你可能需要非常小的 epsilon 和非常大的 n，否則例如 c=−3/4+10{−10} 可能給出錯誤的週期 2 而不是正確的週期 1，該錯誤將累加到錯誤的內部距離估計（例如，你的方法的距離 3.8e-8（epsilon 1e-12，n 79,573,343）而不是我的方法的 2e-10）。" 克勞德·海蘭德-艾倫
- "人可以進行很多次迭代，然後檢視最後一次迭代是否等於之前任何一次迭代。但是，這種方法從來都不是萬無一失的。通常，它們在用於靠近其微原子的邊界點的點時會失敗。在這種情況下，演算法通常會“找到”實際週期的整數倍。這是因為當您靠近子微原子時，逼近極限環的點以螺旋或針腳狀方式收斂。例如，假設週期實際上是 5，並且我們的點靠近週期為 35 的子微原子。如果您檢視每第 5 次迭代，您通常會發現一種以螺旋形收斂到極限點的模式，每次圍繞螺旋形轉動需要 7 步。由於螺旋模式，第 7 步（第 35 次迭代）比第 1 步（第 5 次迭代）更接近，簡單的週期查詢演算法會捕捉到更接近的那一個。"Robert Munafo 來自曼德布羅特集合詞彙表和百科全書
解決方案
- 不要使用樸素方法透過迭代計算週期
- 使用 Robert P. Munafo 的 Jordan 曲線方法
- 來自 Robert Munafo 的曼德布羅特集合詞彙表和百科全書，版權所有 1987-2017 年。
- 使用 Claude Heiland-Allen 的盒子週期 Jordan 曲線方法進行週期性掃描
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
來自 Claude Heiland-Allen 的 mandelbrot-numerics 庫 的 C 程式碼
// 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;
}
/* 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;
}
-
具有 3 週期迴圈的臨界軌道 -
M 集的週期性成分 -
繪製 Julia 集和臨界軌道。計算週期 -
計算用於繪製曼德布羅特集的週期
方法
查詢週期用於
使用臨界點的正向迭代查詢臨界軌道的週期
/*
b batch file for maxima
*/
kill(all);
remvalue(all);
/* =================== functions ============ */
/*
https://wikibook.tw/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://wikibook.tw/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;
比較兩個查詢週期的函式
/*
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;
}
非線性比例顯示更大的週期（沿著曼德布羅特集的實數切片）
/*
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;
}
使用者：Simpsons_contributor：週期性檢查
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;
}
}
- ↑ H. Bruin 和 D. Schleicher，二次多項式的符號動力學，Mittag-Leffler 研究所，瑞典皇家科學院，7。
- ↑ 數學堆疊交換問題：有限二進位制序列的週期
- ↑ Claude Heiland-Allen 使用 Haskell 和 SVG 輸出的 lavaurs 演算法
- ↑ 使用區間方法對電子電路中的週期軌道進行嚴格研究，作者：Zbigniew Galias
- ↑ Milan 繪製的曼德布羅特集