分形/曼德尔
< 分形
Mandel:[1] 软件是 Wolf Jung 开发的用于实数和复数动力学的软件。这里可以找到关于它的非官方文档。
- 下载
- 解压
- qmake
- make
- ./mandel
Mandel 程序使用自己的二进制分数表示法
Enter the angle with Notation: decimal 1/7 or binary p001 for (binary periodic angles) decimal 9/56 or binary 001p010 for (binary preperiodic angles) decimal 1/4 or binary 01 for 1/4 = 0.01 (binary dyadic angles with both finite and equal infinite preperiodic binary expansion)
函数族（= 映射）
- 曼德尔布罗特集合 z^2+c（键 Ctrl+0）。这是默认映射。这里的映射度数 q = 2
- 多重布罗特集合 z^q+c（键 Ctrl+1）。使用键 q 更改映射度数
- 布兰纳-法格拉 cz(1+z/q)^q（键 Ctrl+2）
- 三次多项式
- 四次多项式
- 有理映射
- 牛顿映射
- 超越函数
- 非解析映射
这些族由 类 定义
复二次多项式 的单项式和中心形式
"曼德尔布罗特集合基于二次多项式 fc(z) = z2 + c 的单参数族。"（来自帮助）
// from mndynamo.cpp by Wolf Jung (C) 2007-2014.
void mndynamics::root(double x, double y, double &u, double &v)
{ v = sqrt(x*x + y*y);
if (x > 0.0) { u = sqrt(0.5*(v + x)); v = 0.5*y/u; return; }
if (x < 0.0)
{ v = sqrt(0.5*(v - x)); if (y < 0.0) v = -v; u = 0.5*y/v; return; }
if (y >= 0.0) { u = sqrt(0.5*y); v = u; return; }
u = sqrt(-0.5*y); v = -u;
}
// from mndlbrot.cpp by Wolf Jung (C) 2007-2014. part of Mandel 5.12
// input : c = a+b*i and z = x+y*i
// output : z = x+y*i
void mndlbrot::f(double a, double b, double &x, double &y) const
{ double u = x*x - y*y + a; // u is temporary variablke
y = 2*x*y + b;
x = u; }
它可用于
- 前向迭代 ，键 f = 模式 0
- 2 个逆（或反向）迭代（多值，幅角调整）
- 第一个逆函数 = 键 a =
- 第二个逆函数 = 键 b =
"Use the keys a and b for the inverse mapping. (The two branches are approximately mapping to the parts A and B of the itinerary.)
- 与 kneading 序列相关的动态平面的划分
-
1/4 -
1/6 预周期 = 1，周期 = 2 -
9/56 预周期 = 3，周期 = 3 -
129/16256
另请参见 示例 c 代码
// from mndlbrot.cpp by Wolf Jung (C) 2007-2014. part of Mandel 5.12
// input : c, z , mode
// c = A+B*i where A and B are global variables defined in mndynamo.h
// z = x+y*i
//
// output : z = x+y*i
void mndlbrot::iterate(double &x, double &y, int mode) const // = 0
{ // Mapping f = forward itearation = f(z) = key f = mode 0
if (!mode) f(A, B, x, y);
if (mode > 0) { x = 0; y = 0; }
if (mode >= 0 || mode < -2) return;
// f^{-1}(z) = inverse with principal value
if (A*A + B*B < 1e-20)
{ root(x - A, y - B, x, y); // 2-nd inverse function = key b
if (mode & 1) { x = -x; y = -y; } // 1-st inverse function = key a
return;
}
/*f^{-1}(z) = inverse with argument adjusted
"When you write the real and imaginary parts in the formulas as complex numbers again,
you see that it is sqrt( -c / |c|^2 ) * sqrt( |c|^2 - conj(c)*z ) ,
so this is just sqrt( z - c ) except for the overall sign:
the standard square-root takes values in the right halfplane, but this is rotated by the squareroot of -c .
The new line between the two planes has half the argument of -c .
(It is not orthogonal to c ... )"
...
"the argument adjusting in the inverse branch has nothing to do with computing external arguments. It is related to itineraries and kneading sequences, ...
Kneading sequences are explained in demo 4 or 5, in my slides on the stripping algorithm, and in several papers by Bruin and Schleicher.
W Jung " */
double u, v, w = A*A + B*B;
root(-A/w, -B/w, u, v);
root(w - A*x - B*y, B*x - A*y, x, y);
w = u*x - v*y;
y = u*y + v*x;
x = w;
// 2-nd inverse function = key b
if (mode & 1) // mode = -1
{ x = -x; y = -y; } // 1-st inverse function = key a
}
它被调用为
// forward iteration
fAct = new QAction(trUtf8("Mapping &f"), mapAG);
fAct->setShortcut(QString("f"));
// inverse iteration
inv1Act = new QAction(trUtf8("&1st inverse"), mapAG);
inv1Act->setShortcut(QString("a"));
inv2Act = new QAction(trUtf8("&2nd inverse"), mapAG);
inv2Act->setShortcut(QString("b"));
fAct->setEnabled(ftype != 48);
inv1Act->setEnabled(!ftype || (ftype == 28 && signtype*signtype == 4) || ftype == 38 || ftype == 48 || ftype == 58);
inv2Act->setEnabled(!ftype || (ftype == 28 && signtype*signtype == 4) || ftype == 48 || ftype == 58);
void QmnShell::map(QAction *act)
{ double x, y;
dplane->getPoint(x, y);
if (act == fAct) f->iterate(x, y); // mode = 0 ??
if (act == inv1Act) { f->iterate(x, y, -1); /*if (!ftype)*/ dplane->Move(12, f); }
if (act == inv2Act) { f->iterate(x, y, -2); /*if (!ftype)*/ dplane->Move(13, f); }
dplane->setPoint(x, y);
}
内部地址
[edit | edit source]- 这些数字（周期）必须按定义递增
算法
[edit | edit source]- 绘制
- 平面（绘制模式或算法操作）
- LSM 逃逸时间（算法 2）
- DEM（算法 5）
- 平面（绘制模式或算法操作）
- 曲线
- 点
- 轨道（迭代）
- 临界轨道
- 正向和逆向轨道
- 轨道（迭代）
- 计算
使用 Mandel 中代码的示例程序
// qmnshell.cpp by Wolf Jung (C) 2007-2019.
QActionGroup *algoAG = new QActionGroup(this);
algoActs[0] = new QAction(QString("&0"), algoAG);
algoActs[0]->setShortcut(QString("0"));
algoActs[1] = new QAction(trUtf8("&1: Comparing neighbors"), algoAG);
algoActs[1]->setShortcut(QString("1"));
algoActs[2] = new QAction(trUtf8("&2: Escape time"), algoAG);
algoActs[2]->setShortcut(QString("2"));
algoActs[3] = new QAction(trUtf8("&3: Escape time, rainbow colors"),algoAG);
algoActs[3]->setShortcut(QString("3"));
algoActs[4] = new QAction(trUtf8("&4: Marty normality"), algoAG);
algoActs[4]->setShortcut(QString("4"));
algoActs[5] = new QAction(trUtf8("&5: Distance estimate"), algoAG);
algoActs[5]->setShortcut(QString("5"));
algoActs[6] = new QAction(trUtf8("&6: Argument of \xce\xa6"), algoAG);
algoActs[6]->setShortcut(QString("6"));
algoActs[7] = new QAction(trUtf8("&7: Binary decomposition"), algoAG);
algoActs[7]->setShortcut(QString("7"));
algoActs[8] = new QAction(trUtf8("&8: Yoccoz puzzle 1/2 1/3 2/3"), algoAG);
algoActs[8]->setShortcut(QString("8"));
algoActs[9] = new QAction(trUtf8("&9: Zeros of q_n(c) [Q]"), algoAG);
algoActs[9]->setShortcut(QString("9"));
algoActs[10] = new QAction(trUtf8("&c: Cr. points of q_n(c) [Q]"), algoAG);
algoActs[11] = new QAction(trUtf8("&n: Newton for q_n(c) [Q]"), algoAG);
connect(algoAG, SIGNAL(triggered(QAction*)), this, SLOT(drawing(QAction*)));
// mndlbrot.cpp by Wolf Jung (C) 2007-2019.
uint mndlbrot::pixcolor(mdouble x, mdouble y)
{
uint j, cl = 0; mdouble a, b;
if (sign > 0) { a = x; b = y; } //z = c, not z = 0 !
else { a = A; b = B; }
if (drawmode == 4) return marty(a, b, x, y);
if (drawmode == 7) return binary(a, b, x, y);
if (drawmode == 8) return yoccoz(a, b, x, y);
if (sign < 0 ||
( (x*x + y*y - .0625)*(x*x + y*y - .0625) >
((x - .25)*(x - .25) + y*y)*.25 &&
(x + 1)*(x + 1) + y*y > .0625 ) ) //not per 1, 2 in M
{ mdouble x1 = x, y1 = y, u, v;
for (j = 1; j <= maxiter; j++)
{ u = x*x; v = y*y;
if (u + v <= bailout) { y = 2*x*y + b; x = u - v + a; }
else
{ if (drawmode == 2)
{ cl = ((j % 24) >> 3); if (!cl) cl = 4; if (j & 1) cl |= 8; }
else cl = j;
break;
}
}
x = x1; y = y1;
}
if (drawmode == 2 || drawmode == 3) return cl;
if (cl) cl = dist(a, b, x, y);
if (cl & 2) Period = 1;
else { Period = 0; if (drawmode == 6) return 0; }
if (drawmode == 6) return turn(a, b, x, y);
if (drawmode == 9 && sign < 0) return quadrant(a, b, x, y);
if (drawmode == 9) return quadrantqn(a, b);
if (drawmode == 10) return quadrantqnp(a, b);
if (drawmode == 11) return newton(a, b);
return cl; //drawmode == 5
}
逃逸时间的等值线
[edit | edit source]对于算法 2 =（等值线）逃逸时间，是否可以设置参数（？逃逸半径），使得等值线的边界穿过临界点（及其前像）？
您只能通过重新编译来更改逃逸半径。但是您可以选择参数，使其位于逃逸线上和所需的参数射线上，那么临界值也将位于逃逸线上。
// mndlbrot.cpp by Wolf Jung (C) 2007-2019
uint mndlbrot::esctime(mdouble x, mdouble y)
{ uint j = 0;
mdouble a, b, x0, y0, u, v;
if (sign > 0) { a = x; b = y; } //z = c, not z = 0 !
else { a = A; b = B; }
if (b == 0.0 && !drawmode && sign < 0
&& (a == 0.25 || a == -0.75)) return parabolic(x, y);
if (drawmode > 255)
{ if (sign > 0) { x = 0; y = 0; }
return renormtime(a, b, x, y);
}
if (sign > 0 && (x*x + y*y - .0625)*(x*x + y*y - .0625) <
((x - .25)*(x - .25) + y*y)*.25 ) j = 1;
else if (sign > 0 && (x + 1)*(x + 1) + y*y < .0625 ) j = 2;
if (sign > 0 && -1.76875 < x && x < -1.74375 && y*y < .00015625)
{ root(-4*a - 7, -4*b, x0, y0);
u = a + 2 + a*x0 - b*y0; v = b + a*y0 + b*x0;
if (u*u + v*v < .0625) j = 3;
}
if (sign > 0 && ( (x + .125)*(x + .125)+(y - .75)*(y - .75) < .015625 ||
(x + .125)*(x + .125) + (y + .75)*(y + .75) < .015625 ) )
{ root(-4*a - 7, -4*b, x0, y0);
u = a + 2 - a*x0 + b*y0; v = b - a*y0 - b*x0;
if (u*u + v*v < .0625) j = 3;
}
if (j)
{ if (drawmode) return 65280u | (drawmode & 15); else return 65280 | j; }
/*if (drawmode && x*x + y*y < 5 && y < 0)
{ long long int U, V, W = 1LL << 28, Bailout = 5LL << 56;
long int X = (long int)(x*W), Y = (long int)(y*W),
A = (long int)(a*W), B = (long int)(b*W);
for (j = 1; j <= maxiter; j++)
{ U = X; U *= X; V = Y; V *= Y; if (U + V >= Bailout) return j;
W = X; W *= Y; W >>= 27; Y = (long int)(W) + B;
U -= V; U >>= 28; X = (long int)(U) + A;
}
return 65284;//0 | (drawmode & 15);
}//test integer speed up ... problem with Green*/
for (j = 1; j <= maxiter; j++)
{ u = x*x; v = y*y;
if (u + v <= bailout)
{ y = 2*x*y + b; x = u - v + a; }
else return j;
if (sign < 0 && !drawmode
&& (x - temp[0])*(x - temp[0]) + (y - temp[1])*(y - temp[1]) < 1e-6)
{ int cl = 1 + j % Period;
if (Period <= 4 && j % (2*Period) >= Period ) cl |= 8;
return 65280u | cl;
}
}
if (drawmode) return 65280u | (drawmode & 15);
if (sign > 0)
{ mdouble x0 = x, y0 = y;
for (j = 0; j < 60; j++)
{ if (x*x + y*y <= bailout) f(a, b, x, y);
else return maxiter + j;
if ( (x - x0)*(x - x0) + (y - y0)*(y - y0) < 1e-6 )
return 65281u + (j % 12);
}
}
return 65293u;
}
距离估计
[edit | edit source]// mndlbrot.cpp by Wolf Jung (C) 2007-2023. ... are part of Mandel 5.19,
int mndlbrot::dist(mdouble a, mdouble b, mdouble x, mdouble y)
{ uint j; mdouble xp = 1, yp = 0, nz, nzp; //zp = 1
for (j = 1; j <= maxiter; j++)
{ nz = 2*(x*xp - y*yp); yp = 2*(x*yp + y*xp); xp = nz; //zp = 2*z*zp;
if (sign > 0) xp++; //zp = 2*z*zp + 1
nz = x*x - y*y + a; y = 2*x*y + b; x = nz; //z = z*z + c;
nz = x*x + y*y; nzp = xp*xp + yp*yp;
if (nzp > 1e40 || nz > bailout) break;
}
if (nz < bailout) return 1; //not escaping, rare
if (nzp < nz) return 10; //includes escaping through 0
x = log(nz); x = x*x*nz / (temp[1]*nzp); //4*square of dist/pixelwidth
if (x < 0.04) return 1;
if (x < 0.24) return 9;
return 10;
} //dist
势
[edit | edit source]Julia 集的复势
绿色过程以英国数学家乔治·格林命名
// mndlbrot.cpp by Wolf Jung (C) 2007-2019
mdouble mndlbrot::green(int sg, mdouble x, mdouble y)
{ mdouble a = x, b = y; if (sg < 0) { a = A; b = B; }
uint j; mdouble s = 1.0, dr = 0.5, u, v;
for (j = 0; j <= maxiter; j++)
{ s *= dr; u = x*x; v = y*y;
if (u + v > 1e12) return s*log(u + v);
y = 2*x*y + b; x = u - v + a;
}
return 0;
}
简化函数
/*
this function is based on function by W Jung
it is used in : File:Cpmj.png
*/
double jlogphi(double zx0, double zy0, double cx, double cy)
{
int j;
double
zx=zx0,
zy=zy0,
s = 0.5,
zx2=zx*zx,
zy2=zy*zy,
t;
//
for (j = 1; j < 400; j++)
{ s *= 0.5;
zy = 2 * zx * zy + cy;
zx = zx2 - zy2 + cx;
zx2 = zx*zx;
zy2 = zy*zy;
t = fabs(zx2 + zy2); // abs(z)
if ( t > 1e24) break;
}
return s*log2(t); // log(zn)* 2^(-n)
}//jlogphi
等势线通过当前点，或对于给定势，用键 g 绘制。
等势线可以在任何地方定义，但当它们穿过临界点 z = 0 或其前像时，它们将是8 字形。
外部射线通过 0 或其前像分支。
if (act == greenAct)
{ mdouble x, y; theplane->getPoint(x, y); y = f->green(signtype, x, y);
QmnDoubleDialog *dialog = new QmnDoubleDialog(trUtf8(
"To draw an equipotential line, enter the\n"
"potential log|\xce\xa6(%1)| > 0 :\n"
"The suggested value is the potential of the\n"
"current point. Just hit Return to accept it."
).arg((signtype > 0 ? QChar('c') : QChar('z'))), &y, 0, 3, this);
connect(dialog, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
if (!dialog->exec() || y <= 0 || y > 10.0) return;
if (signtype > 0) { drawLater = 0; dplane->stop(); }
else pplane->stop();
theplane->green(f, signtype, y, 10);
}
// qmnplane.cpp by Wolf Jung (C) 2007-2023.
//the following code will be rewritten completely!!!
void QmnPlane::green(mndynamics *f, int st, mdouble g, int quality,
QColor color) //5, white
{ //trace the boundary ...
if (g <= 0 || type) return;
uint m = 1; f->prepare(st, nmax, m, &g); //just to set nmax
int i, k, i0, k0, i1, k1, i2, k2, j, vert, side, sg, start, T; stop();
QPainter *p = new QPainter(buffer); p->setPen(color);
for (T = 20 - kmax; T < kmax - 20; T += imax*2/quality)
{ for (vert = 0; vert <= 1; vert++)
{ for (side = -1; side <= 1; side += 2)
{ start = 10000; //no startpoint
if (vert) //vertical lines from above and below
{ if (f->green(st, hmid - ph*side*kmax + pw*T, vmid - pw*side*kmax + ph*T) <= g)
continue;
for (j = kmax-1; j >= -kmax/2; j--)
{ if (f->green(st, hmid + ph*side*j + pw*T, vmid - pw*side*j + ph*T) <= g)
{start = side*j; break;} }
}
else//vert: horizontal lines from left and right
{ if (f->green(st, hmid + pw*side*imax - ph*T, vmid + ph*side*imax - pw*T) <= g)
continue;
for (j = imax-1; j >= -imax/2; j--)
{ if (f->green(st, hmid + pw*side*j - ph*T, vmid + ph*side*j - pw*T) <= g)
{start = side*j; break;} }
}//vert
if (start == 10000) continue; //no startpoint
for (sg = -1; sg <= 1; sg += 2) //go in both directions
{ if (vert)
{ k0 = start; k1 = k0 + side; k2 = k1; i0 = T; i1 = T; i2 = i1 + sg;}
else
{ i0 = start; i1 = i0 + side; i2 = i1; k0 = T; k1 = T; k2 = k1 + sg;}
for (j = 1; j < 8*imax; j++)
{ if (i0 < -imax || i0 >= imax || k0 < -kmax || k0 >= kmax) break;
else p->drawLine(imax + i0, kmax + k0, imax + i0, kmax + k0);
if (f->green(st, hmid + pw*i2 + ph*k2, vmid + ph*i2 - pw*k2) <= g)
{i = i0; k = k0; i0 = i2; k0 = k2; i2 += (i1 - i); k2 += (k1 - k);}
else
{i = i1; k = k1; i1 = i2; k1 = k2; i2 += (i0 - i); k2 += (k0 - k);}
}//for j
}//for sg
}//for side
}//for vert
}//for T
delete p; update();
} //green
势
另请参阅
二进制分解
[edit | edit source]int mndlbrot::binary(mdouble a, mdouble b, mdouble x, mdouble y)
{ mdouble u, v; uint j;
for (j = 1; j <= maxiter; j++)
{ u = x*x; v = y*y;
if (u + v <= 1e4) { y = 2*x*y + b; x = u - v + a; }
else { if (j & 1) j = 2; else j = 10; if (y > 0) j += 2; return j; }
//else return (x > 50.0 + y && x > 50.0 - y ? 12 : 10); //carrots
}
return 0;
} //binary
重整化
[edit | edit source]if (act == renormAct)
{ uint n1, n2 = 0;
QmnUIntDialog *dialog = new QmnUIntDialog(trUtf8(
"Enter the renormalization period \xe2\x89\xa4 255.\n"
"For embedded Julia sets, enter the preperiod\n"
"and the period, separated with a comma:"),
&n1, &n2, 65002u, 1, this);
connect(dialog, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
if (dialog->exec() && n1 <= 255 && n2 <= 255) n1 |= (n2 << 8);
else return;
if (signtype > 0) { drawLater = 0; dplane->stop(); }
else pplane->stop();
theplane->draw(f, signtype, &n1, 2);
}
周期点和前周期点
[edit | edit source]"Each hyperbolic component has a unique center, and when the current point is within a hyperbolic component, you can move it to the center with the keys x and Return."
费根鲍姆调谐
[edit | edit source]步骤
- 从周期为 n 的中心开始
- 找到周期为 2n 的中心的近似值，并使用牛顿法精确地找到它。可以使用费根鲍姆重新缩放来进行近似，但这只在 1/2 分支中有效，因此我在其他子分支中使用了周期 1 和 2 的乘子映射的组合。
键：Ctrl+Alt+C
C++ 代码
// from the file qmnshell.cpp by Wolf Jung (C) 2007-2018
const mdouble cFb = -1.40115518909205060052L,
dFb = 4.66920160910299067185L,
aFb = 2.50290787509589282228L; //Fibonacci -1.8705286321646448888906L
//
if (act == tuneAct && per && per <= 512)
{ f->find(signtype, 0, per, x, y);
if (x < -.5)
{ x = cFb + (x - cFb)/dFb;
y /= dFb; }
else
{ mndynamics::root(1/16.0 - x*.25, -y*.25, x, y);
x = -.75 - x;
y = -y;
}
per *= 2;
f->find(signtype, 0, per, x, y);
if (f->period(x, y) == per) theplane->setPoint(x, y);
}
另请参阅
- 视频：Feigenbaum_scaling
- Mandel 演示（重整化）第 10 页
/*
This is not official program by W Jung,
but it usess his code ( I hope in a good way)
These functions are part of Mandel by Wolf Jung (C)
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.
http://www.mndynamics.com/indexp.html
to compile :
g++ f.cpp -Wall -lm
./a.out
===========================================
Period = 1 center = 0.000000000000000000
Period = 2 center = -1.000000000000000000
Period = 4 center = -1.310702641336832884
Period = 8 center = -1.381547484432061470
Period = 16 center = -1.396945359704560642
Period = 32 center = -1.400253081214782798
Period = 64 center = -1.400961962944841041
Period = 128 center = -1.401113804939776124
Period = 256 center = -1.401146325826946179
Period = 512 center = -1.401153290849923882
Period = 1024 center = -1.401154782546617839 // c = -1.401154782546620 +0.000000000000000 i period = 1024
Period = 2048 center = -1.401155102022463976
Period = 4096 center = -1.401155170444411267
Period = 8192 center = -1.401155185098297292
Period = 16384 center = -1.401155188236710937
Period = 32768 center = -1.401155188908863045
Period = 65536 center = -1.401155189052817413
Period = 131072 center = -1.401155189083648072
Period = 262144 center = -1.401155189090251057
Period = 524288 center = -1.401155189091665208
Period = 1048576 center = -1.401155189091968106
Period = 2097152 center = -1.401155189092033014
Period = 4194304 center = -1.401155189092046745
Period = 8388608 center = -1.401155189092049779
Period = 16777216 center = -1.401155189092050532
Period = 33554432 center = -1.401155189092051127
Period = 67108864 center = -1.401155189092050572
Period = 134217728 center = -1.401155189092050593
Period = 268435456 center = -1.401155189092050599
*/
#include <iostream> // std::cout
#include <cmath> // sqrt
#include <limits>
#include <cfloat>
typedef unsigned int uint;
typedef long double mdouble; // mdynamo.h
// from the file qmnshell.cpp by Wolf Jung (C) 2007-2018
mdouble cFb = -1.40115518909205060052L;
mdouble dFb = 4.66920160910299067185L;
mdouble bailout = 16.0L; // mdynamoi.h
// c = A+B*i
mdouble A= 0.0L;
mdouble B = 0.0L;
/*
function from mndlbrot.cpp by Wolf Jung (C) 2007-2017 ...
part of Mandel 5.14, 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.
http://www.mndynamics.com/indexp.html
----------------------------------------------
it is used to find :
* periodic or preperiodic points on dynamic plane
* on parameter plane
** centers
** Misiurewicz points
using Newton method
*/
int find(int sg, uint preper, uint per, mdouble &x, mdouble &y)
{ mdouble a = A, b = B, fx, fy, px, py, w;
uint i, j;
for (i = 0; i < 30; i++)
{ if (sg > 0) // parameter plane
{ a = x; b = y; }
if (!preper) // preperiod==0
{ if (sg > 0) // parameter plane
{ fx = 0;
fy = 0;
px = 0;
py = 0; }
else // dynamic plane
{ fx = -x;
fy = -y;
px = -1;
py = 0; }
}
else // preperiod > 0
{ fx = x;
fy = y;
px = 1.0;
py = 0;
for (j = 1; j < preper; j++)
{ if (px*px + py*py > 1e100) return 1;
w = 2*(fx*px - fy*py);
py = 2*(fx*py + fy*px);
px = w;
if (sg > 0) px++; // parameter plane
w = fx*fx - fy*fy + a;
fy = 2*fx*fy + b;
fx = w;
}
}
mdouble Fx = fx, Fy = fy, Px = px, Py = py;
for (j = 0; j < per; j++)
{ if (px*px + py*py > 1e100) return 2;
w = 2*(fx*px - fy*py);
py = 2*(fx*py + fy*px);
px = w;
if (sg > 0) px++; // parameter plane
w = fx*fx - fy*fy + a;
fy = 2*fx*fy + b;
fx = w;
}
fx += Fx;
fy += Fy;
px += Px;
py += Py;
w = px*px + py*py;
if (w < 1e-100) return -1;
x -= (fx*px + fy*py)/w;
y += (fx*py - fy*px)/w;
}
return 0;
}
int main()
{
int plane = 1; // positive is parameter plane, negative is dynamic plane = signtype
uint preper = 0; // " the usual convention is to use the preperiod of the critical value. This has the advantage, that the angles of the critical value have the same preperiod under doubling as the point, and the same angles are found in the parameter plane." ( Wolf Jung )
uint per ; // period
mdouble x ;
mdouble y = 0.0L;
int n;
// Starting with a center of period n
per = 1;
x = 0.0L;
// find an approximation for the center of period 2n
for (n=1; n<30; n++){
printf("Period = %10u \tcenter = %.18Lf\n", per, x);
// next center
per *= 2; // period doubling
// approximate of next value using Feigenbaum rescaling ( in the 1/2-limb )
x = cFb + (x - cFb)/dFb;
// more precise value of x useing Newton method
find(plane, preper, per, x, y);
}
return 0;
}
C 代码
#include <stdio.h> // printf
#include <math.h> // sqrt
typedef unsigned int uint;
typedef long double mdouble; // mdynamo.h
// from the file qmnshell.cpp by Wolf Jung (C) 2007-2018
mdouble cFb = -1.40115518909205060052L;
mdouble dFb = 4.66920160910299067185L;
mdouble bailout = 16.0L; // mdynamoi.h
// c = A+B*i
mdouble A= 0.0L;
mdouble B = 0.0L;
/*
function from mndlbrot.cpp by Wolf Jung (C) 2007-2017 ...
part of Mandel 5.14, 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.
http://www.mndynamics.com/indexp.html
----------------------------------------------
it is used to find :
* periodic or preperiodic points on dynamic plane
* on parameter plane
** centers
** Misiurewicz points
using Newton method
*/
int find(int sg, uint preper, uint per, mdouble *x, mdouble *y)
{
mdouble a = A;
mdouble b = B;
mdouble tx;
mdouble ty;
mdouble fx;
mdouble fy;
mdouble px;
mdouble py;
mdouble w;
uint i, j;
for (i = 0; i < 30; i++)
{ if (sg > 0) // parameter plane
{ a = *x;
b = *y; }
if (!preper) // preperiod==0
{ if (sg > 0) // parameter plane
{ fx = 0;
fy = 0;
px = 0;
py = 0; }
else // dynamic plane
{ fx = - *x;
fy = - *y;
px = -1;
py = 0; }
}
else // preperiod > 0
{ fx = *x;
fy = *y;
px = 1.0;
py = 0;
for (j = 1; j < preper; j++)
{ if (px*px + py*py > 1e100) return 1;
w = 2*(fx*px - fy*py);
py = 2*(fx*py + fy*px);
px = w;
if (sg > 0) px++; // parameter plane
w = fx*fx - fy*fy + a;
fy = 2*fx*fy + b;
fx = w;
}
}
mdouble Fx = fx, Fy = fy, Px = px, Py = py;
for (j = 0; j < per; j++)
{ if (px*px + py*py > 1e100) return 2;
w = 2*(fx*px - fy*py);
py = 2*(fx*py + fy*px);
px = w;
if (sg > 0) px++; // parameter plane
w = fx*fx - fy*fy + a;
fy = 2*fx*fy + b;
fx = w;
}
fx += Fx;
fy += Fy;
px += Px;
py += Py;
w = px*px + py*py;
if (w < 1e-100)
{return -1; }
tx -= (fx*px + fy*py)/w;
ty += (fx*py - fy*px)/w;
}
x = &tx;
y = &ty;
return 0;
}
int main()
{
int plane = 1; // positive is parameter plane, negative is dynamic plane = signtype
uint preper = 0; // " the usual convention is to use the preperiod of the critical value. This has the advantage, that the angles of the critical value have the same preperiod under doubling as the point, and the same angles are found in the parameter plane." ( Wolf Jung )
uint per ; // period
mdouble x ;
mdouble y = 0.0L;
int n;
// Starting with a center of period n
per = 1;
x = 0.0L;
// find an approximation for the center of period 2n
for (n=1; n<30; n++){
printf("Period = %10u \tcenter = %.18Lf\n", per, x);
// next center
per *= 2; // period doubling
// approximate of next value using Feigenbaum rescaling ( in the 1/2-limb )
x = cFb + (x - cFb)/dFb;
// more precise value of x useing Newton method
find(plane, preper, per, &x, &y);
}
return 0;
}
找到
[edit | edit source]要找到周期点或 前周期点 z，请使用
- 菜单/点/查找点
- 键 x
然后输入
- 周期
- 前周期，周期
这是一个数值过程，因此
- 它找到一个最接近的解。
- 如果当前点位于
- 集合内（Mandel 或填充的 Julia）
- 不在好的解附近，则它什么也找不到（在参数平面上它会到原点 c=0）
例如，在参数平面上有 3 个周期为 3 的中心。如果
- 当前点是 c=0，则结果是 c=0（未找到）
- 当前点在 Mandelbrot 集之外
- c = 0.350000000000000 +0.725000000000000 i，则结果是 c = -0.122561166876654 +0.744861766619744 i，周期 = 3（好的结果）
- c = 0.580000000000000 +0.240000000000000 i，周期 = 0，则结果是 c = 0.000000000000000 +0.000000000000000 i，周期 = 1（不好的结果 ！！！）
- c = -1.590000000000000 +0.170000000000000 i，周期 = 0，则结果是 c = -1.754877666246693 +0.000000000000000 i，周期 = 3（好的结果）
使用 Del-Newton 模式显示牛顿法的吸引盆，用 q 调整周期
它调用 find 函数，
- 输入
- sg（“正数是参数平面，负数是动态平面。”）
- preper = 前周期。（“通常的做法是使用临界值的预周期。这样做的优势是临界值的角在加倍时与点的预周期相同，并且在参数平面上会找到相同的角。”沃尔夫·荣格）
- per = 周期
- 输出
- x = creal(z)
- y = cimag(z)
- 返回值（类型为 int）
返回值
- 0（好的结果）
- 2 当 (px*px + py*py > 1e100)
- 1 当 (px*px + py*py > 1e100) = 在前周期循环中逃逸
- -1 当 (w < 1e-100)
/* function from mndlbrot.cpp by Wolf Jung (C) 2007-2017 ...
part of Mandel 5.14, 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.
*/
int mndlbrot::find(int sg, uint preper, uint per, double &x, double &y) const
{ double a = A, b = B, fx, fy, px, py, w; uint i, j;
for (i = 0; i < 30; i++)
{ if (sg > 0) { a = x; b = y; }
if (!preper)
{ if (sg > 0) { fx = 0; fy = 0; px = 0; py = 0; }
else { fx = -x; fy = -y; px = -1; py = 0; }
}
else
{ fx = x; fy = y; px = 1.0; py = 0;
for (j = 1; j < preper; j++)
{ if (px*px + py*py > 1e100) return 1;
w = 2*(fx*px - fy*py); py = 2*(fx*py + fy*px);
px = w; if (sg > 0) px++;
w = fx*fx - fy*fy + a; fy = 2*fx*fy + b; fx = w;
}
}
double Fx = fx, Fy = fy, Px = px, Py = py;
for (j = 0; j < per; j++)
{ if (px*px + py*py > 1e100) return 2;
w = 2*(fx*px - fy*py); py = 2*(fx*py + fy*px);
px = w; if (sg > 0) px++;
w = fx*fx - fy*fy + a; fy = 2*fx*fy + b; fx = w;
}
fx += Fx; fy += Fy; px += Px; py += Py;
w = px*px + py*py; if (w < 1e-100) return -1;
x -= (fx*px + fy*py)/w; y += (fx*py - fy*px)/w;
}
return 0;
}
使用示例
#include <iostream> // std::cout
typedef unsigned int uint;
// c = A+B*i
double A= -0.965000000000000;
double B = 0.085000000000000;
/*
function from mndlbrot.cpp by Wolf Jung (C) 2007-2017 ...
part of Mandel 5.14, 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.
http://www.mndynamics.com/indexp.html
----------------------------------------------
it is used to find :
* periodic or preperiodic points on dynamic plane
* on parameter plane
** centers
** Misiurewicz points
using Newton method
-------------------------------------------
to compile :
g++ f.cpp -Wall
./a.out
*/
int find(int sg, uint preper, uint per, double &x, double &y)
{ double a = A, b = B, fx, fy, px, py, w;
uint i, j;
for (i = 0; i < 30; i++)
{ if (sg > 0) // parameter plane
{ a = x; b = y; }
if (!preper) // preperiod==0
{ if (sg > 0) // parameter plane
{ fx = 0;
fy = 0;
px = 0;
py = 0; }
else // dynamic plane
{ fx = -x;
fy = -y;
px = -1;
py = 0; }
}
else // preperiod > 0
{ fx = x;
fy = y;
px = 1.0;
py = 0;
for (j = 1; j < preper; j++)
{ if (px*px + py*py > 1e100) return 1;
w = 2*(fx*px - fy*py);
py = 2*(fx*py + fy*px);
px = w;
if (sg > 0) px++; // parameter plane
w = fx*fx - fy*fy + a;
fy = 2*fx*fy + b;
fx = w;
}
}
double Fx = fx, Fy = fy, Px = px, Py = py;
for (j = 0; j < per; j++)
{ if (px*px + py*py > 1e100) return 2;
w = 2*(fx*px - fy*py);
py = 2*(fx*py + fy*px);
px = w;
if (sg > 0) px++; // parameter plane
w = fx*fx - fy*fy + a;
fy = 2*fx*fy + b;
fx = w;
}
fx += Fx;
fy += Fy;
px += Px;
py += Py;
w = px*px + py*py;
if (w < 1e-100) return -1;
x -= (fx*px + fy*py)/w;
y += (fx*py - fy*px)/w;
}
return 0;
}
int main()
{
// input
int plane = -1; // positive is parameter plane, negative is dynamic plane.
uint preperiod =0; // " the usual convention is to use the preperiod of the critical value. This has the advantage, that the angles of the critical value have the same preperiod under doubling as the point, and the same angles are found in the parameter plane." ( Wolf Jung )
uint period = 3;
// output
double re,im;
int r; // returnig value = error message
r = find(plane, preperiod,period,re,im);
if (r==0)
{
std::cout.precision( 15 );
std::cout << "Preperiod = "<< preperiod << "\nPeriod = "<< period <<"\n";
if (plane<0) {
std::cout << " parameter c = "<< A << " ; " << B << ":\n";
std::cout << " point z = "<< re << " ; " << im << ":\n";
}
if (plane>=0) std::cout << "point c = "<< re << " ; " << im << ":\n";
}
else { std::cout << "error \n";
if (r==-1) std::cout << " w < 1e-100 \n";
if (r== 2) std::cout << "(px*px + py*py > 1e100) = escaping in periodic loop \n";
if (r== 1) std::cout << "(px*px + py*py > 1e100) = escaping in preperiodic loop\n";
}
return 0;
}
输出
Preperiod = 0 Period = 3 parameter c = -0.965 ; 0.085: point z = -0.602943702541717 ; 0.0385332450804691:
算法 9 : qn(c) 的零点
[edit | edit source]-
参数平面根据第 5 次迭代的符号进行颜色编码。它使用算法 9，阶数 n=5。 -
动态平面根据第 8 次迭代的符号进行颜色编码。它使用算法 9，阶数 q=8 -
动态平面 q = 5
描述
- "它旨在显示预临界点的位置，即 0 的原像。因此它用四种不同的颜色显示，f^{n-1}(z) 在哪个象限。这意味着实轴被拉回了 n 次，因此当参数为实数时它给出了棋盘格的近似值。" (沃尔夫·荣格)
- 它“允许在 4 种颜色相遇的特定周期内定位迭代多项式的零点”。[6]
它可以用于
- 在动力学平面上
- 制作 抛物线棋盘格（棋盘） 当 c 是抛物线参数时
- 显示填充的朱利亚集的组件中心，当 c 是双曲组件的中心并且 z=0 是吸引周期轨道之一时
这是平面的径向 4 次分解（将其与 LSM/M 的第 n 次分解进行比较）
使用 4 种颜色是因为有 4 个象限
- re(z_n) > 0 且 im(z_n) > 0 (第一象限)
- re(z_n) < 0 且 im(z_n) > 0 (第二象限)
- re(z_n) < 0 且 im(z_n) < 0 (第三象限)
- re(z_n) > 0 且 im(z_n) < 0 (第四象限)。
"... 当四种颜色在某个地方相遇时，您就有了 q_n(c) 的零点，即周期除以 n 的中心。此外，浅色或深色表示 c 是否在 M 中。" (沃尔夫·荣格) 参见 mndlbrot.cpp 文件中 mndlbrot 类中的 quadrantqn 函数 [7]
// mndlbrot.cpp by Wolf Jung (C) 2007-2014. part of Mandel 5.12
int mndlbrot::quadrant(double a, double b, double x, double y)
{ //exterior checked in Period (dist)
int cl = 1, j;
for (j = 1; j < subtype; j++)
{ f(a, b, x, y); if (x*x + y*y > 1e100) return 0; }
if (x > 0) cl = 3; if (y < 0) cl++;
if (Period) cl |= 8; return cl;
}//quadrant
int mndlbrot::quadrantqn(double a, double b)
{ //exterior checked in Period (dist)
double x = a, y = b; int cl = 1, j;
for (j = 1; j < subtype; j++)
{ f(a, b, x, y); if (x*x + y*y > 1e100) return 0; }
if (x > 0) cl = 3; if (y < 0) cl++;
if (Period) cl |= 8; return cl;
}//quadrantqn
int mndlbrot::quadrantqnp(double a, double b)
{ //exterior checked in Period (dist)
int cl = 1, j;
mndplex C(a, b), Q = C, P = 1.0;
for (j = 1; j < subtype; j++)
{ P = 2.0*Q*P + 1.0; Q = Q*Q + C;
if (norm(Q) > 1e100 || norm (P) > 1e100) return 0;
}
if (real(P) > 0) cl = 3; if (imag(P) < 0) cl++;
if (Period) cl |= 8; return cl;
} //quadrantqnp
// mndlbrot.cpp by Wolf Jung (C) 2007-2014. part of Mandel 5.12,
uint mndlbrot::pixcolor(double x, double y)
{ uint j, cl = 0; double a, b;
if (sign > 0) { a = x; b = y; } //z = c, not z = 0 !
else { a = A; b = B; }
if (drawmode == 9 && sign < 0) return quadrant(a, b, x, y);
if (drawmode == 9) return quadrantqn(a, b);
if (drawmode == 10) return quadrantqnp(a, b);
与标准循环的区别是
- 没有逃逸测试（例如 abs(zn) 的最大值）。这意味着每个点都具有相同数量的迭代。对于较大的 n，可能会发生溢出。同样，也不能使用测试 Iteration==IterationMax
- 迭代限制相对较小（例如，从 IterationMax = 8 开始）
比较
- Sage 中的代码[8]
以下是用于计算点 c = cx + cy * i 的 8 位颜色的 c 代码片段
/* initial value of orbit = critical point Z= 0 */
Zx=0.0;
Zy=0.0;
Zx2=Zx*Zx;
Zy2=Zy*Zy;
/* the same number of iterations for all points !!! */
for (Iteration=0;Iteration<IterationMax; Iteration++)
{
Zy=2*Zx*Zy + Cy;
Zx=Zx2-Zy2 +Cx;
Zx2=Zx*Zx;
Zy2=Zy*Zy;
};
/* compute pixel color (8 bit = 1 byte) */
if ((Zx>0) && (Zy>0)) color[0]=0; /* 1-st quadrant */
if ((Zx<0) && (Zy>0)) color[0]=10; /* 2-nd quadrant */
if ((Zx<0) && (Zy<0)) color[0]=20; /* 3-rd quadrant */
if ((Zx>0) && (Zy<0)) color[0]=30; /* 4-th quadrant */
用于 qn(c) 的牛顿法
[edit | edit source]它可以与以下内容一起使用
- 菜单/绘制/算法/用于 Qn(c) 的牛顿法
- 键 n
键 q 允许您更改 n（周期）
令
- q_n(c) 是 z^2 + c 的第 n 次迭代，其中 z = 0
- N_n(c) 是相应的牛顿方法。
然后
- q_n 的零点是周期除以 n 的中心，同时也是 N_n 的超吸引不动点。在无穷远处还有另一个不动点，它是排斥的。
- q_n 的临界点是 N_n 的极点。在许多情况下，这些位于周期小于 n 的双曲组件中。
- N_n 的每个直接吸引盆地都至少通过一个通道连接到无穷远。因此，我认为它与原子域无关。
- 似乎每个周期为 n 的双曲组件都完全包含在一个吸引盆地中。对此我还没有证明。
米修列维奇点
[edit | edit source]"... 通常的做法是使用临界值的预周期。这样做的好处是，临界值的角在加倍下具有与该点相同的预周期，并且在参数平面上发现了相同的角。" (沃尔夫·荣格)
查找主心形唤醒 k/r 的主要米修列维奇点
- 转到参数平面
- 通过绘制落在根点上的两条外部射线来标记唤醒
- 使用菜单/射线/唤醒的角度
- 以最低的圈数输入分数 k/r（k/r = 旋转数 = 以圈数表示的积分角）
- 点击曼德布罗特集之外，但位于唤醒内
- 查找点（）
- 菜单/点/查找点或键 x
- 输入预周期、周期
- 如果唤醒太小（= 您在射线之间看不到空隙），则放大
int_angle prep,period c e_angles_of_the_wake e_engles_of_Misiurewicz
1/2 2,1 c = -1.543689012692076 +0.000000000000000 i ( 1/3 or p01 and 2/3 or p10 )
1/3 3,1 c = -0.101096363845622 +0.956286510809142 i ( 1/7 or p001 and 2/7 or p010 ) 9/56, 11/56 and 15/56
1/4 4,1 c = 0.366362983422764 +0.591533773261445 i
1/5 5,1 c = 0.437924241359463 +0.341892084338116 i
1/6 6,1 c = 0.424512719050040 +0.207530228166745 i
1/7 7,1 c = 0.397391822296541 +0.133511204871878 i
1/8 8,1 c = 0.372137705449577 +0.090398233158173 i
1/9 9,1 c = 0.351423759052522 +0.063866559813293 i (1/511= p000000001 ; 2/511 = p000000010 )
1/10 10,1 c = 0.334957506651529 +0.046732666062027 i (1/1023 or p0000000001 and 2/1023 or p0000000010 )
--------------------------------------------------------------------------------------------------------------------------------------------
1/11 11,1 c = 0.321911396847221 +0.035204463294452 i (1/2047 or p00000000001 and 2/2047 or p00000000010 )
1/12 12,1 c = 0.311507660281508 +0.027173719501342 i (1/4095 or p000000000001 and 2/4095 or p000000000010 )
1/13 13,1 c = 0.303127979909719 +0.021411628038965 i (1/8191 or p0000000000001 and 2/8191 or p0000000000010 )
1/14 14,1 c = 0.296304475349759 +0.017171379707062 i (1/16383 or p00000000000001 and 2/16383 or p00000000000010 )
1/15 15,1 c = 0.290687752431041 +0.013982147106557 i (1/32767 or p000000000000001 and 2/32767 or p000000000000010 )
1/16 16,1 c = 0.286016666695566 +0.011537401448441 i (1/65535 or p0000000000000001 and 2/65535 or p0000000000000010 )
1/17 17,1 c = 0.282094678248954 +0.009631861589570 i (1/131071 or p00000000000000001 and 2/131071 or p00000000000000010 )
1/18 18,1 c = 0.278772459129383 +0.008124579648410 i ( 1/262143 or p000000000000000001 and 2/262143 or p000000000000000010 )
1/19 19,1 c = 0.275935362441682 +0.006916613801737 i (1/524287 or p0000000000000000001 and 2/524287 or p0000000000000000010 )
1/20 20,1 c = 0.273494431676492 +0.005937124623590 i ( 1/1048575 or p00000000000000000001 and 2/1048575 or p00000000000000000010 )
-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------
1/21 21,1 c = 0.271379927384804 +0.005134487480794 i ( 1/2097151 or p000000000000000000001 and 2/2097151 or p000000000000000000010 )
1/22 22,1 c = 0.269536632623270 +0.004470475898698 i ( 1/4194303 or p0000000000000000000001 and 2/4194303 or p0000000000000000000010 )
1/23 23,1 c = 0.267920417711850 +0.003916372840385 i ( 1/8388607 or p00000000000000000000001 and 2/8388607 or p00000000000000000000010 )
1/24 24,1 c = 0.266495701963622 +0.003450320181976 i ( 1/16777215 or p000000000000000000000001 and 2/16777215 or p000000000000000000000010 )
1/25 25,1 c = 0.265233559524147 +0.003055480072447 i ( 1/33554431 or p0000000000000000000000001 and 2/33554431 or p0000000000000000000000010 )
1/26 26,1 c = 0.264110291981537 +0.002718738820085 i ( 1/67108863 or p00000000000000000000000001 and 2/67108863 or p00000000000000000000000010 )
1/27 27,1 c = 0.263106342463072 +0.002429779655182 i ( 1/134217727 or p000000000000000000000000001 and 2/134217727 or p000000000000000000000000010 )
1/28 28,1 c = 0.262205461953178 +0.002180410330127 i ( 1/268435455 or p0000000000000000000000000001 and 2/268435455 or p0000000000000000000000000010 )
1/29 29,1 c = 0.261394063652992 +0.001964069379345 i ( 1/536870911 or p00000000000000000000000000001 and 2/536870911 or p00000000000000000000000000010 )
1/30 30,1 c = 0.260660718810682 +0.001775459345952 i ( 1/1073741823 or p000000000000000000000000000001 and 2/1073741823 or p000000000000000000000000000010 )
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
1/35 35,1 c = 0.257872123091544 +0.001121742345611 i ( 1/34359738367 or p00000000000000000000000000000000001 and 2/34359738367 or p00000000000000000000000000000000010 )
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
1/62 62,1 c = 0.252537923584907 +0.000203841219482 i ( 1/4611686018427387903 = p00000000000000000000000000000000000000000000000000000000000001 ;
2/4611686018427387903 = p00000000000000000000000000000000000000000000000000000000000010 )
1/63 63,1 c = 0.252458520819920 +0.000194332543868 i ( 1/9223372036854775807 = p000000000000000000000000000000000000000000000000000000000000001 ;
2/9223372036854775807 = p000000000000000000000000000000000000000000000000000000000000010 )
1/64 64,1 c = 0.252382784694904 +0.000185406363701 i ( 1/18446744073709551615 = p0000000000000000000000000000000000000000000000000000000000000001 ;
2/18446744073709551615 = p0000000000000000000000000000000000000000000000000000000000000010 )
与以下内容进行比较
临界轨道
[edit | edit source]如果您想绘制比 Mandel 作为标准绘制更多的迭代的临界轨道
- 在参数平面上选择内部角（键 c 或 Manu/Poins/bifurcate）
- 转到动力学平面
- 重新绘制图像以删除临界轨道
- 按住并持续按下键 Ctrl-F。每次击键将进行 4000 次迭代；
功能
- 在 qmnshell.cpp 中的 QmnShell::dFinished() 函数中，请参见 : dplane->drawOrbit(f, x, y, 0, 4000);
- 在 QmnShell::map(QAction *act) = Ctrl-F 中
// qmnplane.cpp by Wolf Jung (C) 2007-2017.
// part of Mandel 5.15,
void QmnPlane::drawOrbit(mndynamics *f,
mdouble &x0,
mdouble &y0,
int preiter,
int plotiter)
{ if (type > 0) return;
stop();
int i, j, k;
QPainter *p = new QPainter(buffer);
p->setPen(Qt::white);
// iterate without polotting
for (j = 0; j < preiter; j++)
{ f->iterate(x0, y0);
if (x0*x0 + y0*y0 > 1e6) break; }
// iterate and plot
for (j = 0; j < plotiter; j++)
{ f->iterate(x0, y0);
if (x0*x0 + y0*y0 > 1e6) break;
if (!pointToPos(x0, y0, i, k)) p->drawPoint(i, k);
}
delete p; update();
}
抛物线朱利亚集内部的平铺
[edit | edit source]
它适用于抛物线朱利亚集的内部，仅适用于以下情况
- 旋转数为 0 或 1/2（键 C 或主菜单/点/分叉），则 c 是卫星组件的根。
- 父组件的周期为 1
- 映射是复二次多项式（主菜单/映射或键= Ctrl+0）
它制作了 抛物线棋盘格
// from mndlbrot.cpp by Wolf Jung (C) 2007-2015.
uint mndlbrot::parabolic(double x, double y)
{ uint j; double u, v;
for (j = 1; j <= maxiter; j++)
{ u = x*x; v = y*y;
if (u + v <= bailout)
{ y = 2*x*y; x = u - v + A; }
else return j;
if (A > 0 && x >= 0 && x <= 0.5 && (y > 0 ? y : -y) <= 0.5 - x)
return (y > 0 ? 65281u : 65289u);
if (A < 0 && x >= -0.5 && x <= 0 && (y > 0 ? y : -y) <= 0.3 + 0.6*x)
{ if (j & 1) return (y > 0 ? 65282u : 65290u);
else return (y > 0 ? 65281u : 65289u);
}
// c para 1/3 ; not working now
// if (x > -0.25 && y > 0 && x + y < 0.183012701892219)
// { j %= 3; if (!j) return 65290;
// else if (j & 1) return 65291; else return 65289;}
//
}
return 65293u;
}
等势线
[edit | edit source]- 键 G
- 主菜单/射线/等势线
"... 绘制等势线。... 边界扫描速度很慢，牛顿方法存在三个问题：找到起点，在断开的朱利亚集周围绘制多条线，以及根据所选平面的子集选择点数。因此，我使用边界跟踪，在穿过图像的多条线上找到起点。参见 QmnPlane::green()。 尽管如此，有时不会绘制整条线，或者不会找到所有组件。"
"按下键 g 以绘制穿过当前点 z（Kc 外部）或 c（M 外部）的等势线，或指定电势水平。"
green(f, signtype, y, 10);
// qmnshell.cpp by Wolf Jung (C) 2007-2023.
if (act == greenAct)
{ mdouble x, y; theplane->getPoint(x, y); y = f->green(signtype, x, y);
QmnDoubleDialog *dialog = new QmnDoubleDialog(trUtf8(
"To draw an equipotential line, enter the\n"
"potential log|\xce\xa6(%1)| > 0 :\n"
"The suggested value is the potential of the\n"
"current point. Just hit Return to accept it."
).arg((signtype > 0 ? QChar('c') : QChar('z'))), &y, 0, 3, this);
connect(dialog, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
if (!dialog->exec() || y <= 0 || y > 10.0) return;
if (signtype > 0) { drawLater = 0; dplane->stop(); }
else pplane->stop();
theplane->green(f, signtype, y, 10);
}
// qmnplane.cpp by Wolf Jung (C) 2007-2017.
// Mandel 5.15,
//the following code will be rewritten completely!!!
void QmnPlane::green(mndynamics *f, int st, mdouble g, int quality,
QColor color) //5, white
{ //trace the boundary ...
if (g <= 0 || type) return;
uint m = 1; f->prepare(st, nmax, m, &g); //just to set nmax
int i, k, i0, k0, i1, k1, i2, k2, j, vert, side, sg, start, T; stop();
QPainter *p = new QPainter(buffer); p->setPen(color);
for (T = 20 - kmax; T < kmax - 20; T += imax*2/quality)
{ for (vert = 0; vert <= 1; vert++)
{ for (side = -1; side <= 1; side += 2)
{ start = 10000; //no startpoint
if (vert) //vertical lines from above and below
{ if (f->green(st, hmid - ph*side*kmax + pw*T, vmid - pw*side*kmax + ph*T) <= g)
continue;
for (j = kmax-1; j >= -kmax/2; j--)
{ if (f->green(st, hmid + ph*side*j + pw*T, vmid - pw*side*j + ph*T) <= g)
{start = side*j; break;} }
}
else//vert: horizontal lines from left and right
{ if (f->green(st, hmid + pw*side*imax - ph*T, vmid + ph*side*imax - pw*T) <= g)
continue;
for (j = imax-1; j >= -imax/2; j--)
{ if (f->green(st, hmid + pw*side*j - ph*T, vmid + ph*side*j - pw*T) <= g)
{start = side*j; break;} }
}//vert
if (start == 10000) continue; //no startpoint
for (sg = -1; sg <= 1; sg += 2) //go in both directions
{ if (vert)
{ k0 = start; k1 = k0 + side; k2 = k1; i0 = T; i1 = T; i2 = i1 + sg;}
else
{ i0 = start; i1 = i0 + side; i2 = i1; k0 = T; k1 = T; k2 = k1 + sg;}
for (j = 1; j < 8*imax; j++)
{ if (i0 < -imax || i0 >= imax || k0 < -kmax || k0 >= kmax) break;
else p->drawLine(imax + i0, kmax + k0, imax + i0, kmax + k0);
if (f->green(st, hmid + pw*i2 + ph*k2, vmid + ph*i2 - pw*k2) <= g)
{i = i0; k = k0; i0 = i2; k0 = k2; i2 += (i1 - i); k2 += (k1 - k);}
else
{i = i1; k = k1; i1 = i2; k1 = k2; i2 += (i0 - i); k2 += (k0 - k);}
}//for j
}//for sg
}//for side
}//for vert
}//for T
delete p; update();
} //green
Phi 的幅角或外部角
[edit | edit source]它可以与以下内容一起使用
- 菜单/绘制/算法/Phi 的幅角
- 键 6
// mndlbrot.cpp by Wolf Jung (C) 2007-2017 from Mandel 5.14
//uint mndlbrot::pixcolor(double x, double y)
if (drawmode == 6) return turn(a, b, x, y);
以下是 mndlbrot 类的转动函数
int mndlbrot::turn(double a, double b, double x, double y)
{ //Already checked that escaping. Requires z = c instead of z = 0
int j; double s = 1.0, dr = 0.5, theta, u, v, r;
if (x*x + y*y < 1e-12) return 8; //prevent atan2(0, 0) if disconnected
theta = atan2(y, x); root(.25 - a, -b, u, v); double X = .5 - u, Y = -v;
for (j = 1; j <= 63; j++)
{ s *= dr; u = x*x; v = y*y; r = u + v; if (r < 1e-12) return 8;
u -= v; v = 2*x*y; x = u + a; y = v + b;
//theta += s*u; Adjust in triangle:
u = atan2(u*y - v*x, u*x + v*y);
if ( (y*a - x*b)*(Y*a - X*b) > 0
&& (y*X - x*Y)*(b*X - a*Y) > 0
&& ((b-y)*(a-X) - (a-x)*(b-Y))*(a*Y - b*X) > 0)
{ if (u < 0) u += 2*PI; else u -= 2*PI; }
theta += s*u; if (r > 1e18*s) break;
}
theta *= (.5/PI); theta -= floor(theta);
theta *= (*temp); if (theta > 1e9) return 1;
return 9 + (long int)(theta) % 4;
/* j = (long int)(theta); int cl = (j % 24) >> 3; if (!cl) cl = 4;
if (j & 1) cl |= 8; return cl; //*/
} //turn
外部射线
[edit | edit source]方法
- 反向迭代
- 牛顿法
- 幅角跟踪
// qmnshell.cpp by Wolf Jung (C) 2007-2014. Mandel 5.12
void QmnShell::setRay(QAction *act)
{ int error = 0;
QString enterAngleString = trUtf8(
"<b></b>Enter the angle θ "
"with 0 ≤ θ < 1. Notation :<br>"
"\"1/7\" or \"p001\" for 1/7 = "
".<nobr style=\"text-decoration:overline\">001</nobr> "
"(periodic angles),<br>"
"\"9/56\" or \"001p010\" for 9/56 = "
".001<nobr style=\"text-decoration:overline\">010</nobr> "
"(preperiodic angles),<br>"
"\"1/4\" or \"01\" for 1/4 = .01 (dyadic angles).");
QString pullBackString = trUtf8(
"<br>Then hit Return repeatedly to perform the iteration<br>"
"of the Thurston Algorithm. A connecting path between<br>"
"point configurations is used instead of spider legs.<br>"
"Use Home or enter 0 to quit.");
if (act == greenAct)
{ double x, y; theplane->getPoint(x, y); y = f->green(signtype, x, y);
QmnDoubleDialog *dialog = new QmnDoubleDialog(trUtf8(
"To draw an equipotential line, enter the\n"
"potential log|\xce\xa6(%1)| > 0 :\n"
"The suggested value is the potential of the\n"
"current point. Just hit Return to accept it."
).arg((signtype > 0 ? QChar('c') : QChar('z'))), &y, 0, 3, this);
connect(dialog, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
if (!dialog->exec() || y <= 0 || y > 10.0) return;
if (signtype > 0) { drawLater = 0; dplane->stop(); }
else pplane->stop();
theplane->green(f, signtype, y, 10);
}
if (act == rayAct)
{ uint method, q = 0; double x, y; f->getParameter(x, y);
qulonglong N1 = 1LL, N2 = 2LL;
if (signtype < 0) { N1 = 0LL; if (dplane->getType() < 0) N2 = 1LL; }
if (ftype == 1) { N1 = 2LL; N2 = 2LL; }
if (ftype == 28) { N1 = 1LL; N2 = 1LL; }
QmnRayDialog *dialog = new QmnRayDialog(enterAngleString, trUtf8(
"Adjust the quality, 2...10. It is the number of intermediate\n"
"points, or it concerns the distance to the starting point."),
&N1, &N2, &method, &q, 3, this);
connect(dialog, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
if (!dialog->exec()) return; if (ftype == 28) method = 5;
if (!method && dplane->backRay(N1, N2, x, y, q, Qt::white, 1))
method = 1;
if (method & 1)
theplane->newtonRay(signtype, N1, N2, x, y, q, Qt::white, method);
if (method & 2)
{ if (signtype > 0) { drawLater = 0; dplane->stop(); }
else pplane->stop();
theplane->traceRay(signtype, double(N1)/double(N2), f, x, y, q);
}
}
if (act == rayPointAct)
{ double x, y; f->getParameter(x, y);
qulonglong N1 = 0LL, N2;
QmnCombiDialog *dialog = new QmnCombiDialog(enterAngleString,
&N1, &N2, 3, this);
connect(dialog, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
if (!dialog->exec()) return;
int k, p = mndAngle::normalize(N1, N2, k), q = 0; qulonglong n1 = N1;
if (!p)
{ QmnUIntDialog *dialog1 = new QmnUIntDialog(trUtf8(
"The angle %1/%2\n"
"has preperiod + period > 64.").arg(N1).arg(N2),
0, 0, 0u, 3, this, 0);
connect(dialog1,SIGNAL(needHelp(int)), helptabs,SLOT(showPage(int)));
dialog1->exec(); return;
}
QString bin; QmnCombiDialog::numbersToBinary(N1, N2, bin); double l;
if (p == 1) l = mndAngle::lambda(N1, N2, 1.0e-10, 1000);
else
{ qulonglong L = 1ULL; L <<= 60;
l = ((double)(L))*((double)(N1))/((double)(N2));
l = -mndAngle::lambda(((qulonglong)(l)), L, 1.0e-10, 1000);
}
QString text = trUtf8(
"The angle %1/%2 or %3\n"
"has preperiod = %4 and period = %5.\n"
).arg(N1).arg(N2).arg(bin).arg(k).arg(p);
if (l > 0.0 && signtype > 0) text += trUtf8(
"Entropy: e^h = 2^B = \xce\xbb = %1\n").arg(
QString::number(l, 'f', 8));
if (k && signtype < 0) text += trUtf8(
"The corresponding dynamic ray is landing\n"
"at a preperiodic point of preperiod %1 and\n"
"period dividing %2.\n"
"Do you want to draw the ray and to shift z\n"
"to the landing point?").arg(k).arg(p);
if (k && signtype > 0) text += trUtf8(
"The corresponding parameter ray is landing\n"
"at a Misiurewicz point of preperiod %1 and\n"
"period dividing %2.\n"
"Do you want to draw the ray and to shift c\n"
"to the landing point?").arg(k).arg(p);
if (!k && signtype < 0) text += trUtf8(
"The dynamic ray is landing at a repelling\n"
"or parabolic point of period dividing %1.\n"
"Do you want to draw the ray and to shift\n"
"z to the landing point?").arg(p);
if (!k && signtype > 0)
{ q = mndAngle::conjugate(n1, N2);
QmnCombiDialog::numbersToBinary(n1, N2, bin); text += QString(trUtf8(
"The conjugate angle is %1/%2 or %3 .\n"
)).arg(n1).arg(N2).arg(bin);
if (q < p) bin = trUtf8("satellite"); else bin = trUtf8("primitive");
QString t1, t2; mndCombi c; c.fromAngle(N1, N2); qulonglong b;
c.getKneading(b); QmnCombiDialog::codeToKneading(b, t1);
c.getAddress(b); QmnCombiDialog::codeToAddress(b, t2);
text += trUtf8(
"The kneading sequence is %1 and\n"
"the internal address is %2 .\n").arg(t1).arg(t2);
text += trUtf8(
"The corresponding parameter rays are landing\n"
"at the root of a %1 component of period %2.\n").arg(bin).arg(p);
if (q && q < p) text += trUtf8(
"It is bifurcating from period %1.\n").arg(q);
text += trUtf8(
"Do you want to draw the rays and to shift c\n"
"to the corresponding center?");
}
QmnUIntDialog *dialog1 = new QmnUIntDialog(text, 0, 0, 0u, 3, this);
connect(dialog1, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
if (!dialog1->exec()) return;
if (ftype == 58)
{ x = ((double)(N1))/((double)(N2)); if (l < 0.0) l = -l;
uint n; pplane->getNmax(n);
pplane->setPoint(x, l*0.01*n); return;
}
if (q) pplane->newtonRay(1, n1, N2, x, y, 5, Qt::white, 1);
if (signtype < 0) q = dplane->backRay(N1, N2, x, y, 5, Qt::white, 2);
else q = pplane->newtonRay(1, N1, N2, x, y, 5, Qt::white, 2);
if (q <= 1 && !f->find(signtype, k, p, x, y)) theplane->setPoint(x, y);
}
if (act == portraitAct)
{ double x, y; f->getParameter(x, y);
qulonglong n1, N1 = 0ULL, N2;
QmnCombiDialog *dialog = new QmnCombiDialog(trUtf8(
"Enter the characteristic angle \xce\xb8 with 0 < \xce\xb8 < 1,\n"
"2 \xe2\x89\xa4 period \xe2\x89\xa4 64. "
"Notation : \"2/7\" or \"p010\" ."),
&N1, &N2, 3, this);
connect(dialog, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
if (!dialog->exec()) return;
int k, p = mndAngle::normalize(N1, N2, k); if (k || p <= 1) return;
n1 = N1; int q = mndAngle::conjugate(n1, N2);
if (signtype < 0)
{ for (k = 0; k < p; k++)
{ dplane->backRay(N1, N2, x, y); mndAngle::twice(N1, N2); }
if (q == p) for (k = 0; k < p; k++)
{ dplane->backRay(n1, N2, x, y); mndAngle::twice(n1, N2); }
return;
}
bool sphere = (dplane->getType() < 0); if (sphere) dplane->setType(0);
dplane->setPlane(0, 0, 1.2, 0); dplane->draw(f, 0, themode);
dplane->drawEllipse(0, 0, 1.0, 1.0, Qt::green); x = double(N2);
for (k = 0; k < p; k++)
{ dplane->drawOrtho(double(N1)/x, double(n1)/x);
mndAngle::twice(N1, N2); mndAngle::twice(n1, N2);
}
updateRegion = true; drawLater = 0; dplane->setPoint(0, 0);
if (sphere) dplane->setType(-1); else dplane->setPlane(0, 0, 2.0, 0);
}
if (act == laminatAct)
{ double x, y, u, v; dplane->getPoint(x, y);
qulonglong N = 0ULL, N1 = 0ULL, N2;
QmnCombiDialog *dialog = new QmnCombiDialog(enterAngleString,
&N1, &N2, 3, this);
connect(dialog, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
if (!dialog->exec() || !N1) return;
uint n; int k, p = mndAngle::normalize(N1, N2, k); u = (double)(N2);
if (!k && p)
{ N = N1; mndAngle::conjugate(N, N2); v = ((double)(N))/u; }
u = ((double)(N1))/u; dplane->setPoint(0.5*u, 0.0);
if (N && (N1 & 1)) dplane->setPoint(0.5*v, 0.0);
//problem 1/7 avoided, why? new problem 41/127, need rational angles?
pplane->getNmax(n); if (n > 20u) n = 12u;
if (signtype < 0) //no crash when n + k + p ~ 64
{ dplane->setPoint(x, y); f->getParameter(x, y); qulonglong K, D = 1ULL;
for (k = 0; k <= n; k++)
{ for (K = 0ULL; K < D; K++)
{ dplane->backRay(N1 + K*N2, D*N2, x, y);
if (N) dplane->backRay(N + K*N2, D*N2, x, y);
}
D <<= 1;
}
return;
}
bool sphere = (dplane->getType() < 0); if (sphere) dplane->setType(0);
dplane->setPlane(0, 0, 1.2, 0); dplane->draw(f, 0, themode);
dplane->drawEllipse(0, 0, 1.0, 1.0, Qt::green);
if (N)
{ dplane->drawOrtho(u, v);
dplane->drawLamination(0.5*u, 0.5*v + 0.5, n);
dplane->drawLamination(0.5*v, 0.5*u + 0.5, n);
}
else if (N2 == 56ULL && (N1 == 9ULL || N1 == 11ULL || N1 == 15ULL))
{ dplane->drawOrtho(9.0/56.0, 11.0/56.0);
dplane->drawOrtho(11.0/56.0, 15.0/56.0);
dplane->drawOrtho(9.0/56.0, 15.0/56.0);
dplane->drawLamination(9.0/112.0, 11.0/112.0, n);
dplane->drawLamination(11.0/112.0, 15.0/112.0, n);
dplane->drawLamination(15.0/112.0, 65.0/112.0, n);
dplane->drawLamination(65.0/112.0, 67.0/112.0, n);
dplane->drawLamination(67.0/112.0, 71.0/112.0, n);
dplane->drawLamination(71.0/112.0, 9.0/112.0, n);
}
else dplane->drawLamination(0.5*u, 0.5*u + 0.5, n);
updateRegion = true; drawLater = 0; dplane->setPoint(0, 0);
if (sphere) dplane->setType(-1); else dplane->setPlane(0, 0, 2.0, 0);
}
if (act == wakeAct)
{ uint k = 1, r = 2;
QmnUIntDialog *dialog = new QmnUIntDialog(trUtf8(
"Determine the wake of a limb at the main cardioid.\n"
"Enter a fraction k/r for the rotation number,\n"
"in lowest terms, with 1 \xe2\x89\xa4 k < r \xe2\x89\xa4 64 :"),
&k, &r, 65001u, 3, this);
connect(dialog, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
if (!dialog->exec()) return;
qulonglong n, d = mndAngle::wake(((int) (k)), ((int) (r)), n);
if (!d) return;
QString b1, b2; QmnCombiDialog::numbersToBinary(n, d, b1);
QmnCombiDialog::numbersToBinary(n + 1LL, d, b2);
QmnUIntDialog *dialog1 = new QmnUIntDialog(trUtf8(
"The %1/%2-wake of the main cardioid is\n"
"bounded by the parameter rays with the angles\n"
"%3/%4 or %5 and\n"
"%6/%4 or %7 .\n"
"Do you want to draw the rays and to shift c\n"
"to the center of the satellite component?").arg(k).arg(r).arg(
n).arg(d).arg(b1).arg(n + 1LL).arg(b2), 0, 0, 0u, 3, this);
connect(dialog1, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
if (!dialog1->exec()) return;
if (ftype == 58)
{ uint nn; pplane->getNmax(nn);
pplane->setPoint(((double)(n))/((double)(d)), 0.01*nn); return;
}
double x, y; pplane->newtonRay(1, n, d, x, y, 5, Qt::white, 1); n++;
if (pplane->newtonRay(1, n, d, x, y, 5, Qt::white, 2) <= 1
&& !f->find(1, 0, r, x, y)) pplane->setPoint(x, y);
}
if (act == kneadAct)
{ qulonglong N1 = 1LL, N2, n1, n2, d;
QmnCombiDialog *dialog = new QmnCombiDialog(trUtf8(
"Enter a *-periodic kneading sequence,\n"
"e.g, \"ABAAB*\" or \"10110*\".\n"
"Or enter an internal address,\n"
"e.g., \"1-2-4-5-6\".\n"
"(The periods \xe2\x89\xa4 64 are increasing.)"), &N1, &N2, 3, this);
connect(dialog, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
if (!dialog->exec()) return;
int p, q; QString t1, t2, text; mndCombi c;
if (!N2) { p = c.setKneading(N1); c.getAddress(N2); }
else { p = c.setAddress(N2); c.getKneading(N1); }
QmnCombiDialog::codeToKneading(N1, t1);
QmnCombiDialog::codeToAddress(N2, t2);
text = trUtf8(
"The kneading sequence %1 corresponds\n"
"to the internal address %2 .\n"
"The period is %3.\n").arg(t1).arg(t2).arg(p);
q = c.renorm();
if (q <= 1) text += trUtf8(
"It is not simply renormalizable.\n");
else text += trUtf8(
"It is simply renormalizable with lowest period %1.\n").arg(q);
q = c.failsBS();
if (q) text += trUtf8(
"This combinatorics is not realized by quadratic polynomials,\n"
"since it fails the Bruin-Schleicher admissibility condition:\n"
"the abstract Hubbard tree has an evil branch point of period %1."
).arg(q);
else
{ q = c.count();
if (q == 1)
{ text += trUtf8(
"This combinatorics is realized once on the real axis.\n");
t1 = trUtf8("external");
}
else
{ text += trUtf8(
"This combinatorics is realized for %1 complex parameters.\n"
).arg(q);
t1 = trUtf8("smallest");
}
q = c.toAngles(n1, n2, d);
if (q) text += QString("Angles not computed: Error %1").arg(q);
else text += trUtf8(
"The %4 angles are %1/%3 and %2/%3 .\n"
"Do you want to draw the rays and to shift c\n"
"to the corresponding center?").arg(n1).arg(n2).arg(d).arg(t1);
}
QmnUIntDialog *dialog1
= new QmnUIntDialog(text, 0, 0, 0u, 3, this, (q ? 0 : 1) );
connect(dialog1, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
if (!dialog1->exec() || q) return;
if (ftype == 58)
{ double l; qulonglong L = 1ULL; L <<= 60;
l = ((double)(L))*((double)(n1))/((double)(d));
l = mndAngle::lambda(((qulonglong)(l)), L, 1.0e-12, 500);
uint n; pplane->getNmax(n);
pplane->setPoint(((double)(n1))/((double)(d)), l*0.01*n); return;
}
q = 0; double x, y; pplane->newtonRay(1, n1, d, x, y, 5, Qt::white, 1);
if (pplane->newtonRay(1, n2, d, x, y, 5, Qt::white, 2) <= 1
&& !f->find(1, 0, p, x, y)) pplane->setPoint(x, y);
while (1)
{ q++; N2 -= 1LL << (p - 1); p = c.setAddress(N2); if (p <= 1) break;
c.toAngles(n1, n2, d);
pplane->newtonRay(1, n1, d, x, y, 5,
(q & 1 ? Qt::yellow : Qt::white), 1);
pplane->newtonRay(1, n2, d, x, y, 5,
(q & 1 ? Qt::yellow : Qt::white), 1);
}
pplane->newtonRay(1, 0LL, 1LL, x, y, 5,
(q & 1 ? Qt::yellow : Qt::white), 1);
}
if (act == spiderAct || act == slowspiderAct)
{ if (act == spiderAct) pullBackString = trUtf8(
"<br>Then hit Return repeatedly to perform the iteration of<br>"
"the Spider Algorithm. This tentative implementation is not<br>"
"refining the discretization when spider legs get twisted.<br>"
"Use Home or enter 0 to quit.");
qulonglong N1 = 0LL, N2;
QmnCombiDialog *dialog = new QmnCombiDialog(
enterAngleString + pullBackString, &N1, &N2, 3, this);
connect(dialog, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
if (!dialog->exec()) return;
int k, p = mndAngle::normalize(N1, N2, k);
if (!k && p == 1)
{ if (gamma < 0.0) setRay(0); else pplane->setPoint(0.0, 0.0); return; }
if (!p)
{ QmnUIntDialog *dialg = new QmnUIntDialog(trUtf8(
"The angle %1/%2\n"
"has preperiod + period > 64.").arg(N1).arg(N2),
0, 0, 0u, 3, this, 0);
connect(dialg, SIGNAL(needHelp(int)), helptabs,SLOT(showPage(int)));
dialg->exec(); return;
}
if (act == spiderAct)
{ path->spiderInit(N1, N2); double t = mndAngle::radians(N1, N2);
dplane->setPoint(cos(t), sin(t));
}
else path->angleInit(N1, N2);
error = 1;
}
if (act == pathAct)
{ updateRegion = true; lsppp = 0; uint n, m = 0;
if (dplane->getType()) resizeWin(sphereAct);
double x, y; pplane->getPoint(x, y); n = f->period(x, y);
if (n < 3 || n > 64)
{ n = 3; if (gamma < -1.0) { n = (uint)(-gamma); gamma = 0.0; } }
pplane->stop(); if (gamma < 0.0) setRay(0); else pMoved();
gamma = 0.0;
QmnUIntDialog *dialog = new QmnUIntDialog(trUtf8(
"<b></b>A polynomial may be deformed such that the resulting<br>"
"postcritically finite branched covering is equivalent<br>"
"to a polynomial again:"
"<ul><li>"
"You may shift the n-periodic critical value along a<br>"
"closed path back to itself, moving around and between<br>"
"the other points of the critcal orbit. The resulting<br>"
"branched covering is equivalent to a polynomial with<br>"
"critical period n again."
"</li><li>"
"If the path is enclosing a single postcritcal point,<br>"
"this gives a Dehn twist about that point and the<br>"
"critical value. Note that a right-handed path gives a<br>"
"left-handed twist. You may turn around several times<br>"
"to define a Dehn twist with higher winding number."
"</li><li>"
"Or shift the critcal value c to another point z<sub>1</sub> ,<br>"
"which is mapped to c in n iterations. (You can find<br>"
"it with the keys a and b, or approximately from your<br>"
"intuition of the dynamics.) E.g., if the internal<br>"
"address 1-...-k-n is realized in the 1/2-sublimb of<br>"
"period k, you may take the center with 1-...-k as the<br>"
"initial parameter c and choose z<sub>1</sub> by iterating<br>"
"0 backwards according to the kneading sequence.<br>"
"For a direct path, the branched covering will be<br>"
"equivalent to the polynomial realizing 1-...-k-n.<br>"
"But recapture surgery is defined in more general<br>"
"cases: the initial c may be any parameter except 0."
"</ul>"
"Examples of period n = 3 are available with the key Ctrl+d.<br>"
"Now enter the period 3 ≤ n ≤ 64 and draw the path by<br>"
"dragging the mouse. (You may cancel it by releasing the<br>"
"left button outside of the image. To zoom in, you may<br>"
"turn the path off temporarily with the context menu.)")
+ pullBackString, &n, &m, 64u, 3, this, 1);
connect(dialog, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
if (!dialog->exec()) return;
if (n < 3) { setRay(0); return; }
gamma = -((double)(n)); if (signtype > 0) paraDyn();
dplane->setPoint(x, y); dplane->drawOrbit(f, x, y, 0, 4000);
dplane->move(9); return;
}
if (act == redrawAct) //from dMoved(), from user path
{ if (dplane->hasPathLength() <= 0 || gamma > -2.0) return;
ldouble *X = new ldouble[101], *Y = new ldouble[101];
dplane->getUserPath(100, X, Y); int M = (int)(X[0]);
double x, y; f->getParameter(x, y);
X[0] = (ldouble)(x); Y[0] = (ldouble)(y); dplane->move(10); paraDyn();
path->shiftInit((int)(-gamma), M, X, Y); delete[] X; delete[] Y;
error = 1; act = pathAct;
}
if (act == twistAct)
{ const int M = 100; int m, p = 0;
ldouble a, b, A, B, u, v; double w, z;
a = -0.122561166876654L; b = 0.744861766619744L; //rabbit
//a = -1.266423235901739L; b = 0.357398971748218L; //6347/16383
pplane->getPoint(w, z); A = (ldouble)(w); B = (ldouble)(z); w = 1.0;
if (dplane->getType()) dplane->setType(0); if (gamma < 0.0) setRay(0);
updateRegion = true; lsppp = 0; pplane->stop(); pplane->setPoint(a, b);
QmnDoubleDialog *dialog = new QmnDoubleDialog(trUtf8(
"<b></b>To perform a Dehn twist around the Rabbit's ears,<br>"
"enter the winding number: 1 ... 10 is left handed<br>"
"and -1 ... -10 is right-handed.<br>"
"Or enter 400 or 500 to see examples of recapture<br>"
"surgery: the critical value is shifted to a<br>"
"preimage along an arc.") + pullBackString, &w, 0, 3, this);
connect(dialog, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
if (!dialog->exec() || w < -10.0 || w > 500.0) p = 1;
else { m = (int)(w); if(!m || (m > 10 && m != 400 && m != 500)) p = 1; }
if (p) { pplane->setPoint(A, B); return; }
dplane->setPoint(a, b);
ldouble *X = new ldouble[M+1], *Y = new ldouble[M+1];
if (m >= 400) //shift to preimage along straight line segment
{ if (m == 400)
{ p = 4; A = -0.069858045513714L; B = 0.978233344895305L; }
else
{ p = 5; A = 0.027871676743202L; B = 0.902736237344649L; }
A -= a; B -= b; A /= M; B /= M;
for (m = 0; m <= M; m++) { X[m] = a + m*A; Y[m] = b + m*B; }
}
else //shift to itself along circle around z = c^2 + c
{ w = (-2.0*m)*PI/M; //left-handed twist: exterior left, interior right
p = 3; A = -0.662358978622373L; B = 0.562279512062301L;//rabbit
//p = 14; A = -1.253762663411539L; B = 0.368547951368429L;//6347/16383
A = 0.7*A + 0.3*a; B = 0.7*B + 0.3*b; a -= A; b -= B;
for (m = 0; m <= M; m++)
{ v = w*m; u= cos(v); v = 0.3L*sin(v);
X[m] = A + u*a - v*b; Y[m] = B + u*b + v*a;
}
}
path->shiftInit(p, M, X, Y); delete[] X; delete[] Y; error = 1;
}
if (error)
{ gamma = -1.0; drawLater = 0; lsppp = 0;
pplane->setCursorPix(spiderPix); dplane->setCursorPix(spiderPix);
updateRegion = true; if (dplane->getType()) dplane->setType(0);
dplane->setNmax(0); updateActions();
dplane->setPlane(0.0, 0.0, 2.0, 0.0); dplane->draw(f, 0, themode);
if (act == slowspiderAct) act = stepAct;
else dplane->drawPathSegment(path);
}
error = 0;
if (act == stepAct && gamma == -1.0)
{ ldouble x, y; error = path->pathStep(x, y);
dplane->setPoint(x, y); pplane->setPoint(x, y);
f->setParameter(x, y); dplane->drawPathSegment(path);
}
/*if (error > 0)
{ QmnUIntDialog *dialog = new QmnUIntDialog(trUtf8(
"The homotopy type has changed because the\n"
"discretization is too coarse.\n"
"The algorithm may converge to a wrong value."),
0, 0, 0u, 3, this, 0);
connect(dialog, SIGNAL(needHelp(int)), helptabs,SLOT(showPage(int)));
dialog->exec(); return;
} //optional quit?*/
if (error < 0)
{ QmnUIntDialog *dialog = new QmnUIntDialog(trUtf8(
"Further iteration is pointless due\n"
"to floating-point cancellations."),
0, 0, 0u, 3, this, 0);
connect(dialog, SIGNAL(needHelp(int)), helptabs,SLOT(showPage(int)));
dialog->exec(); act = 0;
}
if (act == 0 && gamma < 0.0)
{ gamma = 0.0; path->deactivate(); dplane->move(10);
pplane->setCursorPix(0); dplane->setCursorPix(0);
updateActions(); pMoved(); //twice when called from homeAct
}
} //setRay
算法
[edit | edit source]- 射线到点
- 外部射线
层压
[edit | edit source]按下 o 绘制轨道肖像，当您从参数平面发出命令时，该肖像将显示为层压，而 Ctrl+o 用于预临界层压。
如果您是
- 在参数平面上，则选择层压和角 1/7 会在右侧窗口中绘制单位圆盘的层压
- 在动力学平面上，则选择层压和角 1/7 会在右侧窗口中绘制动力学平面的层压
层压
// qmnshell.cpp by Wolf Jung (C) 2007-2017. From Mandel 5.15
if (act == laminatAct)
{ mdouble x, y, u, v; dplane->getPoint(x, y);
qulonglong N = 0ULL, N1 = 0ULL, N2;
QmnCombiDialog *dialog = new QmnCombiDialog(enterAngleString,
&N1, &N2, 3, this);
connect(dialog, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
if (!dialog->exec() || !N1) return;
uint n; int k, p = mndAngle::normalize(N1, N2, k); u = (mdouble)(N2);
if (!k && p)
{ N = N1; mndAngle::conjugate(N, N2); v = ((mdouble)(N))/u; }
u = ((mdouble)(N1))/u; dplane->setPoint(0.5*u, 0.0);
if (N && (N1 & 1)) dplane->setPoint(0.5*v, 0.0);
//problem 1/7 avoided, why? new problem 41/127, need rational angles?
pplane->getNmax(n); if (n > 20u) n = 12u;
if (signtype < 0) //no crash when n + k + p ~ 64
{ dplane->setPoint(x, y); f->getParameter(x, y); qulonglong K, D = 1ULL;
for (k = 0; k <= n; k++)
{ for (K = 0ULL; K < D; K++)
{ dplane->backRay(N1 + K*N2, D*N2, x, y);
if (N) dplane->backRay(N + K*N2, D*N2, x, y);
}
D <<= 1;
}
return;
}
bool sphere = (dplane->getType() < 0); if (sphere) dplane->setType(0);
dplane->setPlane(0, 0, 1.2, 0); dplane->draw(f, 0, themode);
dplane->drawEllipse(0, 0, 1.0, 1.0, Qt::green);
if (N)
{ dplane->drawOrtho(u, v);
dplane->drawLamination(0.5*u, 0.5*v + 0.5, n);
dplane->drawLamination(0.5*v, 0.5*u + 0.5, n);
}
else if (N2 == 56ULL && (N1 == 9ULL || N1 == 11ULL || N1 == 15ULL))
{ dplane->drawOrtho(9.0/56.0, 11.0/56.0);
dplane->drawOrtho(11.0/56.0, 15.0/56.0);
dplane->drawOrtho(9.0/56.0, 15.0/56.0);
dplane->drawLamination(9.0/112.0, 11.0/112.0, n);
dplane->drawLamination(11.0/112.0, 15.0/112.0, n);
dplane->drawLamination(15.0/112.0, 65.0/112.0, n);
dplane->drawLamination(65.0/112.0, 67.0/112.0, n);
dplane->drawLamination(67.0/112.0, 71.0/112.0, n);
dplane->drawLamination(71.0/112.0, 9.0/112.0, n);
}
else dplane->drawLamination(0.5*u, 0.5*u + 0.5, n);
updateRegion = true; drawLater = 0; dplane->setPoint(0, 0);
if (sphere) dplane->setType(-1); else dplane->setPlane(0, 0, 2.0, 0);
}
//qmnplane.cpp by Wolf Jung (C) 2007-2017 from Mandel 5.15,
// It requests the conjugate angle and iterates the pair forward.
void QmnPlane::drawLamination(mdouble alpha, mdouble beta, uint n)
{ drawOrtho(alpha, beta);
n--; if (!n) return;
alpha *= 0.5; if (alpha < x) alpha += 0.5;
beta *= 0.5; if (beta < x) beta += 0.5;
drawLamination(alpha, beta, n);
if (alpha < 0.5) alpha += 0.5; else alpha -= 0.5;
if (beta < 0.5) beta += 0.5; else beta -= 0.5;
drawLamination(alpha, beta, n);
}
轨道肖像
[edit | edit source]步骤
- 主菜单/射线/轨道肖像
- qmnshell.cpp
- portraitAct
// qmnshell.cpp by Wolf Jung (C) 2007-2017. From Mandel 5.15
if (act == portraitAct)
{ mdouble x, y; f->getParameter(x, y);
qulonglong n1, N1 = 0ULL, N2;
QmnCombiDialog *dialog = new QmnCombiDialog(trUtf8(
"Enter the characteristic angle \xce\xb8 with 0 < \xce\xb8 < 1,\n"
"2 \xe2\x89\xa4 period \xe2\x89\xa4 64. "
"Notation : \"2/7\" or \"p010\" ."),
&N1, &N2, 3, this);
connect(dialog, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
if (!dialog->exec()) return;
int k, p = mndAngle::normalize(N1, N2, k); if (k || p <= 1) return;
n1 = N1; int q = mndAngle::conjugate(n1, N2);
if (signtype < 0)
{ for (k = 0; k < p; k++)
{ dplane->backRay(N1, N2, x, y); mndAngle::twice(N1, N2); }
if (q == p) for (k = 0; k < p; k++)
{ dplane->backRay(n1, N2, x, y); mndAngle::twice(n1, N2); }
return;
}
bool sphere = (dplane->getType() < 0); if (sphere) dplane->setType(0);
dplane->setPlane(0, 0, 1.2, 0); dplane->draw(f, 0, themode);
dplane->drawEllipse(0, 0, 1.0, 1.0, Qt::green); x = mdouble(N2);
for (k = 0; k < p; k++)
{ dplane->drawOrtho(mdouble(N1)/x, mdouble(n1)/x);
mndAngle::twice(N1, N2); mndAngle::twice(n1, N2);
}
updateRegion = true; drawLater = 0; dplane->setPoint(0, 0);
if (sphere) dplane->setType(-1); else dplane->setPlane(0, 0, 2.0, 0);
}
蜘蛛算法
[edit | edit source]"瑟斯顿算法正在反向迭代临界轨道以计算对应于给定外部角的参数 c。为了选择正确的原像，Chéritat 在慢配对中使用点配置之间的连接路径。 蜘蛛算法 使用射线到无穷远。使用 s 或 Ctrl+s 来说明这些算法。使用 d，您可以沿某些路径移动临界值以定义德恩扭转或重新捕获。使用 Ctrl+d 可以获得特殊示例。重复按下回车键以观察迭代。使用 Home 或输入 0 退出。" (来自帮助)
唤醒
[edit | edit source]在程序 mandel 中
- 从主菜单/射线/唤醒的角度
- 键：Ctrl+k
来自 mandel 5.19 的原始代码
//qmnshell.cpp by Wolf Jung (C) 2007-2023. Defines class QmnShell
if (act == wakeAct)
{ uint k = 1, r = 2; qulonglong n, d; mdouble x, y;
QmnUIntDialog *dialog = new QmnUIntDialog(trUtf8(
"Determine the wake of a limb at the main cardioid.\n"
"Enter a fraction k/r for the rotation number,\n"
"in lowest terms, with 1 \xe2\x89\xa4 k < r \xe2\x89\xa4 64 :"),
&k, &r, 65001u, 3, this);
connect(dialog, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
if (!dialog->exec()) return;
d = mndAngle::wake(((int) (k)), ((int) (r)), n);
if (!d) return;
QString b1, b2; QmnCombiDialog::numbersToBinary(n, d, b1);
QmnCombiDialog::numbersToBinary(n + 1LL, d, b2);
QmnUIntDialog *dialog1 = new QmnUIntDialog(trUtf8(
"The %1/%2-wake of the main cardioid is\n"
"bounded by the parameter rays with the angles\n"
"%3/%4 or %5 and\n"
"%6/%4 or %7 .\n"
"Do you want to draw the rays and to shift c\n"
"to the center of the satellite component?").arg(k).arg(r).arg(
n).arg(d).arg(b1).arg(n + 1LL).arg(b2), 0, 0, 0u, 3, this);
connect(dialog1, SIGNAL(needHelp(int)), helptabs, SLOT(showPage(int)));
if (!dialog1->exec()) return;
if (ftype == 58)
{ uint nn; pplane->getNmax(nn);
pplane->setPoint(((mdouble)(n))/((mdouble)(d)), 0.01*nn); return;
}
pplane->newtonRay(1, n, d, x, y, 5, Qt::darkMagenta, 1); n++;
if (pplane->newtonRay(1, n, d, x, y, 5, Qt::darkMagenta, 2) <= 1
&& !f->find(1, 0, r, x, y)) pplane->setPoint(x, y);
}
示例
/*
This is not official program by W Jung,
but it usess his code ( I hope in a good way)
These functions are part of Mandel 5.9 by Wolf Jung (C) 2007-2013,
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.
http://www.mndynamics.com/indexp.html
to compile :
g++ w.cpp -Wall -lm
./a.out
*/
#include <iostream> // std::cout
#include <complex> // std::complex, std::norm
#define PI 3.1415926535897932385L //from mndynamo.h
// type qulonglong = unsigned long long int
// n is a numerator of external angle that land on root point of the wake k/r
// d is a denominator
// funcion mndAngle::wake from mndcombi.cpp by Wolf Jung (C) 2007-2013
unsigned long long int wake(int k, int r, unsigned long long int &n)
{
if (k <= 0 || k >= r || r > 64) return 0LL; //
unsigned long long int d = 1LL; int j, s = 0; n = 1LL;
for (j = 1; j < r; j++)
{ s -= k;
if (s < 0) s += r;
if (!s) return 0LL;
if (s > r - k) n += d << j;
}
d <<= (r - 1);
d--;
d <<= 1;
d++; //2^r - 1 for r <= 64
return d;
}
// from mndynamo.cpp by Wolf Jung (C) 2007-2013
void root(double x, double y, double &u, double &v)
{ v = sqrt(x*x + y*y);
if (x > 0.0) { u = sqrt(0.5*(v + x)); v = 0.5*y/u; return; }
if (x < 0.0)
{ v = sqrt(0.5*(v - x)); if (y < 0.0) v = -v; u = 0.5*y/v; return; }
if (y >= 0.0) { u = sqrt(0.5*y); v = u; return; }
u = sqrt(-0.5*y); v = -u;
}
// int mndlbrot::bifurcate from mndlbrot.cpp by Wolf Jung (C) 2007-2013
// type mndplex = complex
int bifurcate(double t, double &a, double &b)
{ int per = 1; if (a < -0.75) per = 2;
if (a < -1.5 || b > sqrt(27/64.0L) || b < -sqrt(27/64.0L) ) per = 3;
if (t <= -1.0) return per;
t *= (2*PI);
if (per == 1)
{ a = 0.5*cos(t) - 0.25*cos(2*t); b = 0.5*sin(t) - 0.25*sin(2*t); }
if (per == 2) { a = 0.25*cos(t) - 1.0; b = 0.25*sin(t); }
if (per <= 2) return per;
std::complex<double> u, c, c1, l = std::complex<double>(cos(t), sin(t));
if (a < -1.54) c1 = -1.754877666;
else
{ c1 = std::complex<double>(-.122561167, .744861767); if (b < 0) c1 = conj(c1); }
c = 81.0*l*l-528.0*l+4416.0; root(real(c), imag(c), a, b);
u = pow(-13.5*l*l+144.0*l-800.0 + (-1.5*l+12.0)*std::complex<double>(a, b), 1/3.0);
c = (1/3.0)*(.25*u + (1.5*l+4.0)/u - 2.0);
if (norm(c - c1) > .25)
{ u *= std::complex<double>(-0.5, sqrt(0.75)); c = (1/3.0)*(.25*u + (1.5*l+4.0)/u - 2.0); }
if (norm(c - c1) > .25)
{ u *= std::complex<double>(-0.5, sqrt(0.75)); c = (1/3.0)*(.25*u + (1.5*l+4.0)/u - 2.0); }
a = real(c); b = imag(c); return 3;
} //bifurcate
// mndlbrot::find from mndlbrot.cpp by Wolf Jung (C) 2007-2013
// sign int. Defined in mndynamo.h positive is parameter plane, negative is dynamic plane."
int find(int sg, unsigned int preper, unsigned int per, double cx, double cy, double &x, double &y)
{ double a = cx, b = cy, fx, fy, px, py, w;
unsigned int i, j;
for (i = 0; i < 30; i++)
{ if (sg > 0) { a = x; b = y; }
if (!preper)
{ if (sg > 0) { fx = 0; fy = 0; px = 0; py = 0; }
else { fx = -x; fy = -y; px = -1; py = 0; }
}
else
{ fx = x; fy = y; px = 1.0; py = 0;
for (j = 1; j < preper; j++)
{ if (px*px + py*py > 1e100) return 1;
w = 2*(fx*px - fy*py); py = 2*(fx*py + fy*px);
px = w; if (sg > 0) px++;
w = fx*fx - fy*fy + a; fy = 2*fx*fy + b; fx = w;
}
}
double Fx = fx, Fy = fy, Px = px, Py = py;
for (j = 0; j < per; j++)
{ if (px*px + py*py > 1e100) return 2;
w = 2*(fx*px - fy*py); py = 2*(fx*py + fy*px);
px = w; if (sg > 0) px++;
w = fx*fx - fy*fy + a; fy = 2*fx*fy + b; fx = w;
}
fx += Fx; fy += Fy; px += Px; py += Py;
w = px*px + py*py; if (w < 1e-100) return -1;
x -= (fx*px + fy*py)/w; y += (fx*py - fy*px)/w;
}
return 0;
}
// =====================================================================================================================
// ====================================================================================================================
// ========================================================================================================================
int main()
{
unsigned long long int p, q;
unsigned long long int num, denom;
double cx,cy;
double zx,zy;
double t;
// --------------------------------------------------------------------------------------------------------------------
// -------------------- initial value ------------------------------------------------------------------------------
// ------------------------------------------------------------------------------------------------------------------
std::cout << "Determine the wake of a limb at the main cardioid of Mandelbrot set. " << "\n";
std::cout << "Enter a fraction k/r for the rotation number, in lowest terms, with 1 <= k < r < 64 " << "\n";
while (1)
{
std::cout << " Enter numerator of the rotation number, it is integer 1 <= k < 64 : " << "\n";
std::cin >> p ;
if (std::cin.fail()) // no extraction took place
{
std::cin.clear(); // reset the state bits back to goodbit so we can use ignore()
std::cin.ignore(1000, '\n'); // clear out the bad input from the stream
continue; // try again
}
std::cin.ignore(1000, '\n'); // clear out any additional input from the stream
if (p >= 0) break; // if input value is in good range then exit
}
while (1)
{
std::cout << "Enter the denominator of the rotation number, it is integer 1 <= r < 64 : " << "\n";
std::cin >> q ;
if (std::cin.fail()) // no extraction took place
{
std::cin.clear(); // reset the state bits back to goodbit so we can use ignore()
std::cin.ignore(1000, '\n'); // clear out the bad input from the stream
continue; // try again
}
std::cin.ignore(1000, '\n'); // clear out any additional input from the stream
if (q > 0) break; // if input value is in good range then exit
}
std::cout.precision( 15 );
std::cout << "The rotation number is " << p << "/" << q << "\n";
denom = wake(p,q,num);
if (denom!=0LL)
{
std::cout << "The "<< p << "/" <<q <<" wake of the main cardioid is bounded by the parameter rays with the angles :\n";
std::cout << num << "/" << denom << " and "<< num+1 << "/" << denom << "\n";
}
else std::cout << "(k <= 0 || k >= r || r > 64) \n";
t=(double)p/((double)q);
bifurcate(t, cx,cy);
std::cout << "The root point of wake is c = "<< cx << " ; " << cy << ":\n";
find(-1,0,1,cx,cy,zx,zy);
std::cout << "The fixed point alfa is z = "<< zx << " ; " << zy << ":\n";
return 0;
}
方法
[edit | edit source]// qmnplane.cpp by Wolf Jung (C) 2007-2014. Mandel 5.12,
//rewrite double array allocation as in mndPath
int QmnPlane::backRay(qulonglong num, qulonglong denom, double &a, double &b,
const int quality, QColor color, int mode) // = 5, white, 1
{ //Draw a dynamic ray with angle num/denom by backwards iteration with
//quality points per step. At present only for quadratic polynomials.
//mode : 0 all rays, 1 one ray, 2 return endpoint in a, b.
if (type > 0) return 2; //works on sphere and on plane
if (quality <= 1) return 3;
int i, j, k, preper = 0, pppp = 0; //pppp = preper + per
mndAngle t; pppp = t.setAngle(num, denom, preper); if (!pppp) return 4;
//pppp += preper; double c, s, X[quality][], Y[quality][];
//X = new double[quality][pppp + 1]; Y = new double[quality][pppp + 1];
pppp += preper; double c, s,
X[quality][pppp + 1], Y[quality][pppp + 1];
//The second index corresponds to the i-th iterate of the angle. Initial
//radii between 2^12 and 2^24 : 2^(24*2^(-k/quality))
s = log(2); c = 24*s; s = exp(-s/double(quality));
for (k = 0; k < quality; k++)
{ c *= s; X[k][pppp] = exp(c); Y[k][pppp] = 0.5/X[k][pppp]; }
//Using the approximation Phi_c^{-1}(w) ~ w - 0.5*c/w :
for (i = 0; i < pppp;i++)
{ s = t.radians(); c = cos(s); s = sin(s);
for (k = 0; k < quality; k++)
{ X[k][i] = c*X[k][pppp] - (b*s + a*c)*Y[k][pppp];
Y[k][i] = s*X[k][pppp] + (a*s - b*c)*Y[k][pppp];
}
t.twice();
}
for (k = 0; k < quality; k++)
{ X[k][pppp] = X[k][preper]; Y[k][pppp] = Y[k][preper]; }
stop(); QPainter *p = new QPainter(buffer); p->setPen(color);
//2 more for bailout |z| = 4 :
for (j = nmax + 2; j; j--) for (k = 0; k < quality; k++)
{ for (i = 0; i < pppp; i++)
{ c = X[(k ? k - 1 : quality - 1)][i];
s = Y[(k ? k - 1 : quality - 1)][i];
mndynamics::root(X[k][i + 1] - a, Y[k][i + 1] - b, X[k][i], Y[k][i]);
if (c*X[k][i] + s*Y[k][i] < 0) //choose closest preimage
{ X[k][i] = -X[k][i]; Y[k][i] = -Y[k][i]; }
//if (k & 1) p->setPen(Qt::red); else p->setPen(Qt::blue);
int i1, k1, i2, k2;
if ( (!i || !mode) && pointToPos(c, s, i1, k1) <= 1
&& pointToPos(X[k][i], Y[k][i], i2, k2) <= 1)
p->drawLine(i1, k1, i2, k2);
}
X[k][pppp] = X[k][preper]; Y[k][pppp] = Y[k][preper];
}
if (mode == 2) { a = X[quality - 1][0]; b = Y[quality - 1][0]; }
delete p; /*delete[] Y; delete[] X;*/ update(); return 0;
} //backRay
// qmnplane.cpp by Wolf Jung (C) 2007-2014. Mandel 5.12,
//Time ~ nmax^2 , therefore limited nmax .
int QmnPlane::newtonRay(int signtype, qulonglong N, qulonglong D,
double &a, double &b, int q, QColor color, int mode) //5, white, 1
{ uint n; int j, i, k, i0, k0; mndAngle t; t.setAngle(N, D, j);
double logR = 14.0, x, y, u;
u = exp(0.5*logR); y = t.radians();
x = u*cos(y); y = u*sin(y);
if (pointToPos(x, y, i0, k0) > 1) i0 = -10;
stop(); QPainter *p = new QPainter(buffer); p->setPen(color);
for (n = 1; n <= (nmax > 5000u ? 5000u : nmax + 2); n++)
{ t.twice();
for (j = 1; j <= q; j++)
{ if (mode & 4 ? tricornNewton(signtype, n, a, b, x, y,
exp(-j*0.69315/q)*logR, t.radians())
: rayNewton(signtype, n, a, b, x, y,
exp(-j*0.69315/q)*logR, t.radians()) )
{ n = (n <= 20 ? 65020u : 65010u); break; }
if (pointToPos(x, y, i, k) > 1) i = -10;
if (i0 > -10)
{ if (i > -10) { if (!(mode & 8)) p->drawLine(i0, k0, i, k); }
else { n = 65020u; break; }
}
i0 = i; k0 = k;
}
}
//if rayNewton fails after > 20 iterations, endpoint may be accepted
delete p; update(); if (n >= 65020u) return 2;
if (mode & 2) { a = x; b = y; }
if (n >= 65010u) return 1; else return 0;
} //newtonRay
int QmnPlane::rayNewton(int signtype, uint n, double a, double b,
double &x, double &y, double rlog, double ilog)
{ uint k, l; double fx, fy, px, py, u, v = 0.0, d = 1.0 + x*x + y*y, t0, t1;
for (k = 1; k <= 60; k++)
{ if (signtype > 0) { a = x; b = y; }
fx = cos(ilog); fy = sin(ilog);
t0 = exp(rlog)*fx - 0.5*exp(-rlog)*(a*fx + b*fy);
t1 = exp(rlog)*fy + 0.5*exp(-rlog)*(a*fy - b*fx);
fx = x; fy = y; px = 1.0; py = 0.0;
for (l = 1; l <= n; l++)
{ u = 2.0*(fx*px - fy*py); py = 2.0*(fx*py + fy*px);
px = u; if (signtype > 0) px++;
u = fx*fx; v = fy*fy; fy = 2.0*fx*fy + b; fx = u - v + a;
u += v; v = px*px + py*py; if (u + v > 1.0e100) return 1;
}
fx -= t0; fy -= t1; if (v < 1.0e-50) return 2;
u = (fx*px + fy*py)/v; v = (fx*py - fy*px)/v;
px = u*u + v*v; if (px > 9.0*d) return -1;
x -= u; y += v; d = px; if (px < 1.0e-28 && k >= 5) break;
}
return 0;
} //rayNewton
// mndlbrot.cpp by Wolf Jung (C) 2007-2014. Mandel 5.12,
int mndlbrot::turnsign(double x, double y)
{/*Calculate the turn, i.e., the argument of Phi. Return +-1 by comparing
temp[0] and the turn, 0 for failure or far from the ray. Using two
tricks to reduce the ambiguity from the multi-valued argument:
First, the argument should jump on the Julia set instead of the
line [0, c]. Approximate the Julia set by the lines [0, alpha] and
[alpha, c] and change the argument accordingly within the triangle.
Second, keep track of the arguments in temp[j] to detect jumps.
Before searching a starting point for drawing the external ray,
calling prepare(201) sets temp[0] = theta and temp[j] = 0,
and it disables comparison by setting drawmode = 0. Later on, before
tracing the ray, prepare(202) enables comparison by drawmode = 1.
*/
double a = x, b = y; if (sign < 0) { a = A; b = B; }
int j; double s = 1.0, dr = 0.5, theta, u, v, r, eps = 0.004;
if (x*x + y*y < 1e-12) return 0; //prevent atan2(0, 0) if disconnected
theta = atan2(y, x); root(.25 - a, -b, u, v); double X = .5 - u, Y = -v;
for (j = 1; j <= 63; j++)
{ s *= dr; u = x*x; v = y*y; r = u + v; if (r < 1e-12) return 0;
u -= v; v = 2*x*y; x = u + a; y = v + b;
//theta += s*u; First adjust argument in triangle:
u = atan2(u*y - v*x, u*x + v*y);
if ( (y*a - x*b)*(Y*a - X*b) > 0
&& (y*X - x*Y)*(b*X - a*Y) > 0
&& ((b-y)*(a-X) - (a-x)*(b-Y))*(a*Y - b*X) > 0)
{ if (u < 0) u += 2*PI; else u -= 2*PI; }
//Second compare and shift. 3.6 is ok for initial value 0:
if (drawmode)
{ if (u > temp[j] + 3.6) u -= 2*PI;
if (u < temp[j] - 3.6) u += 2*PI;
temp[j] = u;
}
theta += s*u; if (r > 1e18*s) break;
}
//Problem: j larger is inaccuarte, but thus ray ends at esctime 64:
if (r < 100) return 0; //prevent strong inaccuracy. Or r < 1e10*s ?
theta *= (.5/PI);
theta -= temp[0]; theta -= floor(theta); // 0 <= theta < 1
if (theta < eps) return 1; if (1 - eps < theta) return -1;
return 0;
} //turnsign
// qmnplane.cpp by Wolf Jung (C) 2007-2014. Mandel 5.12,
//rewrite with posToPoint?
int QmnPlane::traceRay(int signtype, double t, mndynamics *f,
double &x0, double &y0, int quality, QColor color) //5, white
{ //Draw the external ray for the turn t with the argument tracing method.
//Return 0 if the endpoint was found, 1 no startpoint, 2 or 3 no endpoint.
//f->turnsign() checks if the turn is in a small intervall around t,
//returns +-1 for the side and 0 otherwise, may adjust jumps.
if (type) return 4;
int i, i0, i1, i2, k, k0, k1, k2, j, sign, draw = 0, iold, kold, u, v;
//Set signtype and t in mndynamics, disable adjusting of jumps:
uint m = 201; f->prepare(signtype, 0, m, &t);
//First, search the startpoint on the boundary enlarged by quality >= 1:
i = quality*imax; k = -quality*kmax - 1; i2 = -2; k2 = 0; //go left on top
j = f->turnsign(hmid + pw*i + ph*k, vmid + ph*i - pw*k);
while (1)
{ i += i2; k += k2;
sign = f->turnsign(hmid + pw*i + ph*k, vmid + ph*i - pw*k);
if (j < 0 && sign > 0)
{ iold = i2/2; kold = k2/2;
i1 = i; k1 = k;
i0 = i - iold; k0 = k - kold;
if (f->turnsign(hmid + pw*i0 + ph*k0, vmid + ph*i0 - pw*k0) > 0)
{ i1 = i0; k1 = k0; i0 -= iold; k0 -= kold; }
i2 = i0 + kold; k2 = k0 - iold; u = i1; v = k1;
break; //found startpoint
}
j = sign;
if (i < -quality*imax && i2 < 0) { i2 = 0; k2 = 2; } //down on left line
if (k >= quality*kmax && k2 > 0) { i2 = 2; k2 = 0; } //right on bottom
if (i >= quality*imax && i2 > 0) { i2 = 0; k2 = -2; } //up on right line
if (k < -quality*kmax && k2 < 0) return 1;
}
//Second, trace the ray by triangles with sign(z0) < 0, sign(z1) > 0:
stop(); QPainter *p = new QPainter(buffer); p->setPen(color);
m = 202; f->prepare(signtype, 0, m, &t); //adjusting jumps enabled
for (j = 1; j < quality*3*(imax + kmax); j++)
{ sign = f->turnsign(hmid + pw*i2 + ph*k2, vmid + ph*i2 - pw*k2);
if (sign > 0) { i = i1; k = k1; i1 = i2; k1 = k2; }
else { i = i0; k = k0; i0 = i2; k0 = k2; }
//New reflected point z2:
i2 = i0 + i1 - i; k2 = k0 + k1 - k;
//If not reflected at diagonal edge, take maximal distance to (u, v) :
if (i0 == i1) { if (mabs(k1-v) > mabs(k0-v)) k2 = k1; else k2 = k0; }
if (k0 == k1) { if (mabs(i1-u) > mabs(i0-u)) i2 = i1; else i2 = i0; }
u = i; v = k;
//Check viewport and draw at z1 (t = 0 | 1/2 too low with z0; z2 rough):
if (-imax <= i1 && i1 < imax && -kmax <= k1 && k1 < kmax)
{ if (draw)
{ // >= 8 is less smooth but not going in circles at 399/992:
if ((i1-iold)*(i1-iold)+(k1-kold)*(k1-kold) >= 5)
{ p->drawLine(imax+iold, kmax+kold, imax+i1, kmax+k1);
iold = i1; kold = k1;
}
if (!sign)
{ x0 = hmid + pw*i1 + ph*k1; y0 = vmid + ph*i1 - pw*k1;
delete p; update(); return 0;
}
}
else
{ draw = 1; iold = i1; kold = k1; }
}
else if (draw) //has left the viewport, no endpoint
{ delete p; update(); return 2; }
}
delete p; update(); return 3;
} //traceRay
- 紫色 = RGB( 255,0,255)= HTML(ff00ff) = 不逃逸也不落入固定吸引子的点（无法判定）。对于抛物型朱利亚集的平铺，它可以在迭代次数较少的情况下（例如 100 次）显示排斥花瓣。
在大多数绘制模式中，Mandel 使用 0 到 15 的调色板。
在模式 3 中，0 到 255 之间的值被解释为色调。请参见 qmnplane.cpp 中的第 160 行。
q->setPen(QColor::fromHsv(cl0 & 255, 255, 255));
为了获得灰度，您可以将其替换为：
q->setPen(QColor::fromRgb(cl0, cl0, cl0));
4 种颜色：[9]
- 青色/水色 = (0,255,255)，绿色 = ( 0,255,0)，红色 = ( 255,0,0)，蓝色 = (0,0,255)
此功能可以通过以下方式启动：
- 主菜单文件/更改颜色
- F5 键
Replace a VGA-color by entering the numbers, separated with a comma: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
用户可以输入两个数字，用逗号隔开
- 第一个数字是要更改的颜色
- 第二个数字是要更改为的颜色
类类型
- 基本类
- 从基本类派生
基本类
- mndynamics = 抽象基本类（逃逸到无穷大）
- mndScale ：Misiurewicz 点处参数平面的缩放
- mndAngle ：外部角
- mndCombi ：捏合序列或内部地址
- mndPath ：带有腿或路径的蜘蛛
mndynamics 是一个抽象基本类（逃逸到无穷大）
从 mndynamics 派生的类
- 具有临界关系和逃逸到无穷大
- mndlbrot ：z^2 + c，曼德博集合
- mndmulti ：z^q + c，多重曼德博集合
- mndbfpoly ：cz(1 + z/q)^q（布兰纳-法格拉）
- mndcubic ：各种三次多项式族
- mndquartic ：各种四次多项式族
- mndmating ：各种二次有理函数族（交配），其中一些仅检查周期性，没有逃逸。
- mndmenage ：c(z + 0.5/z^2)，旋转对称（亨里克森）
- mndsingpert：z^2 + c/z^2，z^2 的奇异扰动
- mndexpo ：各种指数族
- mndtrigo ：各种三角函数族
- mndsurge ：分段多项式，模拟拟共形手术
- mndtricorn ：(z*)^q + c，反解析，三叉星或多叉星
- mndrealcubic：三次多项式，实解析参数，2 个独立的临界点。
- mndifs ：用于刘维尔根的迭代函数系统
- mndsiegel = 抽象派生类（持久西格尔，两个轨道，一个可能逃逸）
- mndcubicsiegel ：三次多项式（扎凯里-巴夫-亨里克森 & 2 周期），从 mndsiegel 派生
- mndquartsiegel ：四次多项式（附加临界关系）
- mndexposiegel ：指数
- mndtrigosiegel ：三角函数
- mndmatesiegel ：二次有理函数
- mndnewtonsiegel：具有西格尔循环的四次牛顿映射
- mndnewtonpara ：具有 2 个根和 1 个抛物线盆的立方牛顿映射
- mndnewtonqara ：具有 1 个根和 2 个抛物线盆的立方牛顿映射
- mndnewtonrara ：具有 0 个根和 3 个抛物线盆的立方牛顿映射
- mndnewton = 抽象派生类（立方或四次牛顿映射，具有临界关系）
- mndcubicnewton ：三次多项式的牛顿法，从 mndnewton 派生
- mndquarticnewton ：四次多项式的牛顿法，从 mndnewton 派生
- mndherman = 具有赫曼环的立方有理函数。从 mndynamics 派生（逃逸到 0 或无穷大的特殊情况）
- mndhenon ：Henon 映射从 mndynamics 派生（特殊情况，逃逸到负无穷大，没有临界点）
- mndlambda ：Henon 映射。从 mndynamics 派生（特殊情况，逃逸到方块，没有临界点）
大多数类描述了复函数的一个参数族
其中
mndAngle 类使用无符号长长整数表示一个角度（以转表示），这些整数 <= 2^64 - 1。它们在可能的情况下被归一化为 2^K * (2^P - 1) 的分母，即 K + P <= 64。否则 P = 0。许多应用程序只需要静态函数。
- normalize() 返回周期。
- twice() 将角度加倍，模 1。
- conjugate() 计算共轭角，返回特征点的周期。
- wake() 计算肢体的角度。
- radians() 以弧度给出角度。
- lambda() 为二元角度计算 e^{核心熵 h}。
- realSpider() 检查根是否为实数，并计算中心 c。
- truncatedTuning() 将 (u1, u2)*w 截断到下一个 n 周期角度。
mndCombi 类表示一个 *-周期捏合序列，周期为 P <= 64，以及相应的内部地址。它们都存储在一个无符号长长整数的位中。QString 转换由 QmnCombiDialog 完成。在地址中，第一个位 (2^63) 是周期 64，最后一个位 (2^0) 是周期 1。示例
AABBAA* is 0...01100111
由于上捏合序列为 AABBAAA
1-3-4-7 is 0...01001101
捏合序列始终是上捏合序列：setKneading() 忽略最后一个条目，而 getKneading() 中的最后一个条目应打印为 '*'。set 函数根据捏合序列计算内部地址，反之亦然。
- fromAngle() 计算周期角度的组合。
- renorm() 给出简单重整化的最低周期。
- count() 返回实现的数量，假设可接受性。
- failsBS() 返回布鲁因-施莱舍可接受性条件的第一次失败，请参见 http://arXiv.org/abs/0801.4662。0 表示没有邪恶轨道，组合数据对应于具有方向保持映射的平面哈伯德树，因此对应于实现这些组合的二次多项式。相应的角度通过对每个周期进行反向迭代在 toAngles() 中计算，始终在 q > 2 时选择 1/q 肢。简单的算法第一次在 1-2-6-7-13-14 处失败，因此通过递归调用 backAngles() 对醒觉进行细分。请参见 mndcombi.cpp 中的解释。
- drawmode uint。在 mndynamo.h 中定义。
- maxiter uint。在 mndynamo.h 中定义。
- a
- b
- A double。在 mndynamo.h 中定义。
- B double。在 mndynamo 中定义。
- sign int。在 mndynamo.h 中定义。"正数是参数平面，负数是动力学平面。"。
- subtype int。在 mndynamo.h 中定义 { subtype = subtype0; A = 0; B = 0; bailout = 16.0; temp = new double[5]; }
- signtype = +- subtype ：它的模量记住所选的子类型，它的符号表示当前的平面；将它的值传递给 sign。
- Period uint。在 mndynamo.h 中定义。
- *temp double。在 mndynamo.h 中定义。
- bailout double。在 mndynamo.h 中定义。
- k 前周期（临界值的不是临界点）。这意味着 c= -2 具有 k=1 和 p=1，所以它在这个符号中是 M_1,1 而不是 M_2,1。这有一个优点，即临界值的角在加倍下具有与点相同的预周期，
并且在参数平面中找到了相同的角。
- p 周期
- 平面参数
- mdouble hmid = 水平中间 = 平面中心 x
- vmid = 垂直中间 = 平面中心 y
- pw = 平面宽度
- ph = 平面高度
- uint drawmode;
- int imid （水平像素坐标为中心
- int kmid （垂直坐标为中心
- imax = 水平最大像素坐标
- kmax = 垂直最大像素坐标
- 模式
- alpha;
沃尔夫·荣格使用中心和宽度来描述复平面的矩形视窗
- 矩形中心：xmid + i ymid
- 右边的中点由以下给出：rewidth + i imwidth
/*
from mndlbrot.cpp by Wolf Jung (C) 201
These classes are part of Mandel 5.7, 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::startPlane(int sg, double &xmid, double &rewidth) const
{ if (sg > 0)
{ xmid = -0.5; rewidth = 1.6; } // parameter plane
else { xmid = 0; rewidth = 2.0; } // dynamic plane
所以角点的初始值是
- 参数平面
- 动态平面：-2 <= Re(z) <=2, -2 <= Im(z) <=2
// qmndemos.cpp by Wolf Jung (C) 2007-2023
pplane = new QmnPlane(360, 360); // parameter plane
dplane = new QmnPlane(360, 360, 0); // dynamic plane
/*
qmnshell.cpp by Wolf Jung (C) 2007-2019.
mndynamo.h by Wolf Jung (C) 2007-2019
mndcombi.cpp by Wolf Jung (C) 2007-2019
... part of Mandel 5.17, 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.
g++ n.cpp -Wall -Wextra -lm
*/
#include <iostream> // std::cout
#include <cmath> // sqrt
#include <limits>
#include <cfloat>
#include <bitset>
#include <string>
// #include <QString> // https://doc.qt.ac.cn/qt-5/qstring.html
using namespace std;
typedef unsigned long long int qulonglong;
typedef long double mdouble;
// mndAngle::twice
void mndAngle_twice(qulonglong &n, qulonglong &d)
{ if (n >= d) return;
if (!(d & 1))
{ d >>= 1; if (n >= d) n -= d; return; }
qulonglong large = 1LL;
large <<= 63; //avoid overflow:
if (n < large)
{ n <<= 1; if (n >= d) n -= d; return; }
n -= large;
n <<= 1;
large -= (d - large);
n += large;
}
// "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 ???
// mndAngle::normalize
/*
input
n = numerator
d = denominator
output
p = period
k = preperiod
*/
int mndAngle_normalize(qulonglong &n, qulonglong &d, int &k)
{ if (!d) return 0;
n %= d;
while (!(n & 1) && !(d & 1))
{ n >>= 1; d >>= 1; }
int p; //
qulonglong n0;
qulonglong n1 = n;
qulonglong d1 = d;
qulonglong np;
k = 0;
while (!(d1 & 1))
{ k++; d1 >>= 1; if (n1 >= d1) n1 -= d1; }
n0 = n1;
for (p = 1; p <= 65 - k; p++)
{ mndAngle_twice(n1, d1);
if (n1 == n0) break; }
if (k + p > 64) return 0; // "The angle N1/N2 has preperiod + period > 64."
np = 1LL << (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;
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// ++++++++++++++++++++ main +++++++++++++++++++++++++++++++++++++++++++++++++
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
int main(void)
{
// input
qulonglong n = 1;
qulonglong d = 5;
// output
int k; //
int p; //
cout << "The input angle = n/d = " << n <<"/" << d << endl;
p = mndAngle_normalize(n, d, k);
cout << "after normalisation angle = n/d = " << n <<"/" << d << endl;
cout << "preperiod k = " << k <<" period p = " << p << endl;
//if (!p)
// {};
// else {}
return 0;
}
编译
g++ n.cpp -Wall -Wextra -lm
运行
./a.out
结果
The input angle = n/d = 1/5 after normalisation angle = n/d = 3/15 preperiod = k = 0 period = p = 4
在commons上有一个类别：Category:Fractals created with Mandel
可以在Mandel中打开png图像并检查它。
图像必须是
- png类型
- 正方形（宽度 = 高度）
- 大小必须是典型的（2000而不是2001像素）
- 世界坐标应该是
- 参数平面（待办事项）
- 动态平面：-2 <= Re(z) <=2, -2 <= Im(z) <=2
步骤
- 打开Mandel
- 调整图像大小（菜单/文件/调整大小）
- 如果您打开动态平面的图像，请更改参数平面的坐标（菜单/点/坐标）
- 打开图像（菜单/文件/加载.png）
然后你可以
- 绘制：轨道，...
- 保存
你不能
- 缩放
- 重绘
参见
- 主菜单/帮助/
- 源代码：qmndemos.cpp
查找具有西格尔盘的 c 参数（在主心形边界上）
- 从内部角度 t = 黄金分割共轭开始[10]
- 从 t 计算 c
- 使用 Go 键转到下一个 t（此处 y 增加 0.01）
这是原始代码
// qmndemos.cpp by Wolf Jung (C) 2007-2013. part of Mandel 5.12,
//
pplane->draw(pf, 1, &mode); double x = 0, y = 0;
pf->bifurcate(.61803398874989484820L, x, y); pplane->setPoint(x, y);
void QmnDemoPB::go()
{
// .....
if (page == 7)
{ double x, y; pplane->getPoint(x, y); // 2alpha :
mndynamics::root(1 - 4*x, -4*y, x, y); //
x = 1 - x;
y = -y;
y = atan2(y, x) + 0.01; // t = atan2() is the argument in radians. It is increased by 0.01 so 1/628 turns
x = 0.5*cos(y) - 0.25*cos(2*y);
y = 0.5*sin(y) - 0.25*sin(2*y);
pplane->setPoint(x, y);}
}
大象演示：在一个小区域上只进行旋转和缩放，而不是像展开中那样的非线性变换。
抛物线重整化方法[11]被用来证明曼德勃罗集的边界具有 2 的豪斯多夫维数，[12]并且最近用来构建一个具有正面积的朱莉娅集。 [13]对于抛物线内爆的另一个应用，请回想一下 1/r-肢体连接到主心形在参数 c 处，使得固定点 是抛物线的，具有组合旋转数 1/r。这些肢体收敛到主心形的根 c = 0.25。皮埃尔·拉瓦尔使用法图坐标和相关的埃卡勒圆柱体的角映射研究了这些肢体。据推测[14]，在适当的缩放下，这些肢体收敛于一个极限“大象”。左图显示了旋转 90o 的 1/2-肢体。所有肢体将以越来越大的放大倍数和适当的旋转显示，使得它们从主心形的边界向上指向。（主心形被旋转，使得它在卫星分量的根部是水平的。）右图显示了对应于周期 r 的中心的填充朱莉娅集。按回车键或按下 Go 按钮以增加肢体的周期 r。这种缩放在主程序中不可用。有关抛物线参数的另一种缩放现象，请参阅第 5 章第 9 页。）
{ n++; mdouble t = 2*PI/n,
x = 0.5*cos(t) - 0.25*cos(2*t),
y = 0.5*sin(t) - 0.25*sin(2*t),
u = sin(t) - sin(2*t),
v = -cos(t) + cos(2*t);
t = 2.5/(n*n);
u *= t;
v *= t; //move root below the center of the image:
pplane->setPlane(x - 0.6*v, y + 0.6*u, u, v);
//pplane->setPlane(x - 0.8*v - 0.2*u, y + 0.8*u - 0.2*v, 0.1*u, 0.1*v);
pplane->draw(pf, 1, &mode);
pf->find(1, 0, n, x, y);
pplane->setPoint(x, y);
}
这里
- 旋转数 1/r 是 t = 1/n
- 每次按下 Go 按钮时，n 都会增加 1
- setPlane(mdouble xmid, mdouble ymid, mdouble rewidth, mdouble imwidth, bool nowayBack)
参数平面：曼德勃罗集，具有角度为 的外部射线，其中
void QmnDemoER::pFinished()
{ int j; mdouble x, y;
if (page == 1)
{ for (j = 0; j < 16; j++) pplane->traceRay(1, mdouble(j)/16.0, pf, x, y, 2);
for (y = 1.35; y > 0.39; y *= 0.78)
pplane->green(pf, 1, pf->green(1, -1.0, y), 3);
}
Wolf Jung 使用 FFmpeg 制作的视频，将使用 Mandel 制作的图像组合在一起[15]
- 朱莉娅集的分岔
- 曼德勃罗集的相似性和自相似性
- 瑟斯顿算法的可视化
变量 σ 与近似抛物线点的乘数有关，当接近根时，1/σ 将趋于无穷大。1/σ 也近似等于从近似盆地内部到外部的迭代次数。在一个视频中，时间与 1/σ 成正比，这意味着螺旋以恒定速度旋转。在另一个视频中，参数 c 以恒定速度移动，这使得 1/σ 越来越快，以至于在有限时间内达到无穷大。
- Julia 集的分岔：c = 时前周期点和射线的分岔（c 由落在临界值的外部射线的角度 θ 参数化）
参数 c 正在实轴上向左移动。当它穿过角度为 5/12 和 7/12 的 Misiurewicz 参数时，临界值具有相同的角度，临界点具有角度 5/24、7/24、17/24、19/24；这些射线以不同的模式成对地落在该参数的左侧或右侧。
请注意，分岔一直在发生：Julia 集的部分上下移动，朝向 0，然后向左和向右移动。
The angle 5/12 or 01p10 has preperiod = 2 and period = 2. The corresponding parameter ray lands at a Misiurewicz point of preperiod 2 and period dividing 2. the landing point is c = -1.543689012692076 +0.000000000000000 i
视频：Mandelbrot 集和二次多项式\ 费根鲍姆标度的幻灯片
1/2 分支中的费根鲍姆重标度
Madel 演示 5 第 10 页的描述
重复的周期加倍序列会导致无限可重整化的费根鲍姆参数 cF ≈ -1.401155189，嵌套的小 Mandelbrot 集按第一个费根鲍姆常数 渐近地缩放，而围绕 z = 0（而不是 z = c）的小 Julia 集按第二个费根鲍姆常数 缩放。参数平面中的装饰变得越来越毛茸茸，并收敛到平面的密集子集。John Milnor 的这个猜想被 Mikhail Lyubich 证明，Curt McMullen 证明了费根鲍姆 Julia 集的放大收敛。
按 Return 键或多次按下 Go 按钮，将参数平面按 δF 围绕 c = cF 重标度，将动态平面按 αF 围绕 z = 0 重标度。在主程序中，费根鲍姆标度可以使用 Ctrl+Alt+f 在参数平面上进行。您需要使用 n 增加迭代次数。相应的中心参数调整可以使用 Ctrl+Alt+c 获得。使用 r 键进行的重整化在原始情况下效果最好。
另请参阅
- 显示窗口（视口）参数：中心和半径
- 主菜单/射线。是否将其更改为主菜单/曲线？它用于绘制射线和等势线
- 保存图像描述，可能在 png 图像内或在不同的 txt 文件中
- 从各种算法制作一个图像（例如，将 DEM 或 IIM 生成的边界添加到其他图像），图像的层
- 用于制作复杂图像/操作的脚本语言表单
- OpenMP、GPU
- 代码中和有关算法的更多描述
- 将静默错误检查更改为显式信息
- 在不保持键预播的情况下执行长时间操作的选项。例如，尝试绘制从周期 1 使用 t = 34/89 进行的分岔的临界轨道。它只绘制一些点。
- 词典
- Qt Linguist 词典 = qph 文件扩展名。这些是包含标准短语及其翻译的人类可读的 XML 文件。这些文件由 Qt Linguist 创建和更新，可供任何数量的项目和应用程序使用。
- Qt Linguist 手册
- 不使用有限二进制展开，而是使用无限版本，因为它更好地描述了动力学。例如，现在是：“角度 1/4 或 01 的前周期 = 2，周期 = 1。”。更好的版本是：“角度 1/4 或 00p1 的前周期 = 2，周期 = 1。”
- 描述如何添加新的分形算法或着色方法
- 参数平面和动态平面的内部射线
- “从另一个来源加载某些 png 图像后，更改颜色和绘图等操作将不再起作用……”
示例 Qt Linguist 词典（qph 文件）。它是一个包含标准短语及其翻译的人类可读的 XML 文件。此文件由 Qt Linguist 创建和更新，可供任何数量的项目和应用程序使用。
<!DOCTYPE QPH>
<QPH sourcelanguage="en" language="pl_PL">
<phrase>
<source>Arnold's Cat map </source>
<target>Odwzorowanie kota Arnolda</target>
<definition>In mathematics, Arnold's cat map is a chaotic map from the torus into itself, named after Vladimir Arnold, who demonstrated its effects in the 1960s using an image of a cat, hence the name.</definition>
</phrase>
<phrase>
<source>lamination </source>
<target>laminacja</target>
<definition>Optional definition</definition>
</phrase>
<phrase>
<source>map </source>
<target>odwzorowanie </target>
<definition>Optional definition</definition>
</phrase>
<phrase>
<source>Boettcher map</source>
<target>Przekształcenie Boettchera</target>
</phrase>
<phrase>
<source>basin of attraction</source>
<target>obszar przyciągania</target>
</phrase>
<phrase>
<source>wake</source>
<target>kilwater</target>
<definition>podzbiór płaszczyny parametrów ograniczona przez 2 promienie zewnętrzne. Porównaj limb, subwake</definition>
</phrase>
<phrase>
<source>dynamic plane </source>
<target>płaszczyzna dynamiczna</target>
</phrase>
<phrase>
<source>mating</source>
<target>?parowanie?</target>
</phrase>
</QPH>
- Mary E. Wilkerson 的《二次函数交配的有限细分规则：存在和构造》
- D. GHIOCA、H. KRIEGER、K. NGUYEN 和 H. YE 的《动力学 Andre-Oort 猜想：单临界多项式》
- ↑ W Jung 的 Mandel
- ↑ C++ 中的曼德尔布罗集合
- ↑ C++ 和 SDL 中的曼德尔布罗集合
- ↑ ArgPhi
- ↑ 符号动力学
- ↑ 克劳德的 periodicity_scan_revisited
- ↑ 沃尔夫·容格的多平台 C++ 程序 Mandel
- ↑ 克里斯托弗·奥拉赫关于逃逸时间分形的形成
- ↑ rapidtables：RGB 颜色代码表
- ↑ wikipedia : 黄金分割共轭
- ↑ 井上裕之和宍仓光弘的抛物线重整化
- ↑ 宍仓光弘关于曼德尔布罗集合边界和朱利亚集的豪斯多夫维数
- ↑ Xavier Buff 和 Arnaud Cheritat 关于具有正面积的二次朱利亚集
- ↑ 阿德里安·杜阿迪 : 朱利亚集是否连续依赖于多项式？
- ↑ 沃尔夫·容格关于曼德尔布罗集合和二次多项式的视频