分形/数学/周期
< 分形
所有落在相同周期点的射线都具有相同的周期：射线的公共周期是其着陆点周期的（可能适当的）倍数；因此，区分了：射线周期与轨道周期。[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 [email protected]: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://wikibooks.cn/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://wikibooks.cn/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;
}
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 绘制的曼德布罗特集