分形/曼德爾
< 分形
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.io/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 關於具有正面積的二次朱利亞集
- ↑ 阿德里安·杜阿迪 : 朱利亞集是否連續依賴於多項式？
- ↑ 沃爾夫·容格關於曼德爾布羅集合和二次多項式的影片