确实存在所谓的超操作序列。它的工作原理就像你构建乘法一样a*b=a+a+a+a...+a并添加了许多a重复b次。然后是求幂a^b = a*a*a*a*...*a与许多乘法a重复b次。然后，就出现了四联，表示为指数塔，就像a^^b == a^a^a^...^a， 重复b times.
我感兴趣如何编写这个函数，用于浮点数和复数？
我已经在 glsl 中编写了乘法和求幂函数：
// complex multiplication:
vec2 cmul(in vec2 a, in vec2 b) {
return vec2(a.x*b.x-a.y*b.y, a.x*b.y+a.y*b.x);
}
// complex exponent e^a
vec2 cexp(in vec2 a) {
float ea = exp(a.x);
float vl = a.y;
return ea * vec2( cos(vl), sin(vl) );
}
// complex natural logarithm ln(a)
vec2 cln(in vec2 a) {
float ql = length(a);
return vec2( log(ql), atan(a.y, a.x));
}
// complex power function a^b
vec2 cpow(in vec2 a, in vec2 b) {
return cexp(cmul(cln(a), b));
}
但后来我被困住了！我们怎样才能写ctet(in vec2 a, in vec2 b)四分函数，不仅适用于浮点数，还适用于整个复平面本身？
好吧，让我们开始吧真实域名 and integer b only:
a^^b = a^a^a^a^a...^a // a is there b times
这可以这样评估C++:
double tetration(double a,int b) // a^^b = a^a^a^a... b times
{
double c;
if (b<=0) return 0;
for (c=a;b>1;b--) c=pow(a,c);
return c;
}
因为你已经得到了pow对于复杂的领域，你也可以在那里做同样的事情......为了保持简单，我现在不会碰它......
这里有一些结果：
a\b| 1| 2| 3| 4
-------------------------------
1 | 1| 1| 1| 1
2 | 2| 4| 16|65536
3 | 3| 27|7625597484987|
4 | 4| 256| |
5 | 5|3125| |
btw.所有这些超级操作与阿克曼函数你可以在这里找到我的 C++ 迭代实现：
然而，由于极快的增长，即使double很快就会超出范围（因此缺少值）...
现在如何移动b to 真实域名？对代数方法一无所知，但几何方法是可能的。
简单来说就是“剧情”a^^b作为变量的函数b和常数a对于整数值b围绕你想要的真实b然后插值真实域名 b using 整数域 b作为控制点。它类似于获得函数的非整数阶导数。
So (X,Y)将是你的(a^^b,b)。现在使用任何插值来构造实域函数。
线性插值将如下所示：
y0 = a^^(int(b))
y1 = a^^(int(b)+1)
a^^b = y0 + (b-int(b))*(y1-y0)
然而，需要更高阶的插值，并且插值参数也应该缩放为非线性度量。欲了解更多信息，请参阅：
经过一些阐述立方体（t^3) and log^2事实证明比例足够了（使用我的 128 位浮点的 C++ 示例f128类只需将其重命名为double):
f128 tetration_fi(f128 a,int b) // a^^b = a^a^a^a... b times
{
f128 c;
if (b==-1) return 0.0; // first singularity
if (b== 0) return 1.0; // second singularity
if (b< -1) return 0.0; // uncomputed
for (c=a;b>1;b--) c=pow(a,c);
return c;
}
//---------------------------------------------------------------------------
f128 tetration_ff(f128 a,f128 b) // a^^b = a^a^a^a... b times
{
int bi;
f128 z0,z1,z2,z3,a0,a1,a2,a3,t,tt,ttt,o=2.0;
if (b==-1) return 0.0; // first singularity
if (b== 0) return 1.0; // second singularity
if (b< -1) return 0.0; // uncomputed
bi=b.ToInt(); b-=bi;
if (b.iszero()) return tetration_fi(a,bi);
z0=tetration_fi(a,bi-1); // known points around a^^b
z1=pow(a,z0);
z2=pow(a,z1);
z3=pow(a,z2);
z0=log2(log2(z0+o)+o); // log^2 scale
z1=log2(log2(z1+o)+o);
z2=log2(log2(z2+o)+o);
z3=log2(log2(z2+o)+o);
t =0.5*(z2-z0); // cubic interpolation coeff.
tt=0.5*(z3-z1);
a0=z1;
a1=t;
a2=(3.0*(z2-z1))-(2.0*t)-tt;
a3=t+tt+(2.0*(z1-z2));
t=b-bi; // cubic interpolation
tt=t*t;
ttt=tt*t;
z0=a0+(a1*t)+(a2*t*t)+(a3*t*t*t);
z0=exp2(exp2(z0)-o)-o; // linear scale
return z0;
}
//---------------------------------------------------------------------------
这是我将其与以下内容进行比较：
我选择相同的图表基础a from a^^b正如您所看到的，它是一个非常好的匹配，只有低于 1.0 的范围略有偏差。
让我们研究复杂域分形
现在当你想去复杂域你不能像在 Real 中那样做，因为结果对于插值来说太混乱了。所以我们只能坚持整数b或者使用Kneser算法来计算。
幸运的是，我们有更多的方法来显示分形......例如我们可以评估整数b来自a^^b其中只有a很复杂，并使用结果对输出进行着色。这里GLSL示例（基于我的曼德尔布罗特着色器和你复杂的数学）：
分段：
// Fragment
#version 450 core
uniform dvec2 p0=dvec2(0.0,0.0); // mouse position <-1,+1>
uniform double zoom=1.000; // zoom [-]
in smooth vec2 p32;
out vec4 col;
//---------------------------------------------------------------------------
// All components are in the range [0…1], including hue.
vec3 rgb2hsv(vec3 c)
{
vec4 K = vec4(0.0, -1.0 / 3.0, 2.0 / 3.0, -1.0);
vec4 p = mix(vec4(c.bg, K.wz), vec4(c.gb, K.xy), step(c.b, c.g));
vec4 q = mix(vec4(p.xyw, c.r), vec4(c.r, p.yzx), step(p.x, c.r));
float d = q.x - min(q.w, q.y);
float e = 1.0e-10;
return vec3(abs(q.z + (q.w - q.y) / (6.0 * d + e)), d / (q.x + e), q.x);
}
//---------------------------------------------------------------------------
// All components are in the range [0…1], including hue.
vec3 hsv2rgb(vec3 c)
{
vec4 K = vec4(1.0, 2.0 / 3.0, 1.0 / 3.0, 3.0);
vec3 p = abs(fract(c.xxx + K.xyz) * 6.0 - K.www);
return c.z * mix(K.xxx, clamp(p - K.xxx, 0.0, 1.0), c.y);
}
//---------------------------------------------------------------------------
vec3 spectral_color(float l) // RGB <0,1> <- lambda l <400,700> [nm]
{
float t; vec3 c=vec3(0.0,0.0,0.0);
if ((l>=400.0)&&(l<410.0)) { t=(l-400.0)/(410.0-400.0); c.r= +(0.33*t)-(0.20*t*t); }
else if ((l>=410.0)&&(l<475.0)) { t=(l-410.0)/(475.0-410.0); c.r=0.14 -(0.13*t*t); }
else if ((l>=545.0)&&(l<595.0)) { t=(l-545.0)/(595.0-545.0); c.r= +(1.98*t)-( t*t); }
else if ((l>=595.0)&&(l<650.0)) { t=(l-595.0)/(650.0-595.0); c.r=0.98+(0.06*t)-(0.40*t*t); }
else if ((l>=650.0)&&(l<700.0)) { t=(l-650.0)/(700.0-650.0); c.r=0.65-(0.84*t)+(0.20*t*t); }
if ((l>=415.0)&&(l<475.0)) { t=(l-415.0)/(475.0-415.0); c.g= +(0.80*t*t); }
else if ((l>=475.0)&&(l<590.0)) { t=(l-475.0)/(590.0-475.0); c.g=0.8 +(0.76*t)-(0.80*t*t); }
else if ((l>=585.0)&&(l<639.0)) { t=(l-585.0)/(639.0-585.0); c.g=0.84-(0.84*t) ; }
if ((l>=400.0)&&(l<475.0)) { t=(l-400.0)/(475.0-400.0); c.b= +(2.20*t)-(1.50*t*t); }
else if ((l>=475.0)&&(l<560.0)) { t=(l-475.0)/(560.0-475.0); c.b=0.7 -( t)+(0.30*t*t); }
return c;
}
//---------------------------------------------------------------------------
// complex domain math
vec3 color_wheel(vec2 a) // complex -> polar -> HSV -> RGB
{
float an=(atan(-a.y,-a.x)*0.15915494309189533576888376337251)+0.5;
float r=length(a); r-=floor(r); r*=0.75; r+=0.25;
return hsv2rgb(vec3(an,1.0,r));
}
vec3 color_spectral(vec2 a) // complex -> wavelength -> RGB
{
float r=length(a); r-=floor(r);
return spectral_color(400.0+(300.0*r));
}
vec2 cadd(vec2 a,vec2 b) // a+b
{
return a+b;
}
vec2 csub(vec2 a,vec2 b) // a-b
{
return a-b;
}
vec2 cmul(vec2 a,vec2 b) // a*b
{
return vec2((a.x*b.x)-(a.y*b.y),(a.x*b.y)+(a.y*b.x));
}
vec2 cdiv(vec2 a,vec2 b) // a/b
{
float an=atan(-a.y,-a.x)-atan(-b.y,-b.x);
float r=length(a)/length(b);
return r*vec2(cos(an),sin(an));
}
vec2 csqr(vec2 a) // a^2
{
return cmul(a,a);
}
vec2 cexp(vec2 a) // e^a
{
// e^(x+y*i)= e^x * e^(y*i) = e^x * ( cos(y) + i*sin(y) )
return exp(a.x)*vec2(cos(a.y),sin(a.y));
}
vec2 cln(vec2 a) // ln(a)
{
return vec2(log(length(a)),atan(a.y,a.x));
}
vec2 cpow(vec2 a,vec2 b) // a^b
{
return cexp(cmul(cln(a),b));
}
vec2 ctet(vec2 a,int b) // a^^b
{
vec2 c=vec2(1.0,0.0);
for (;b>0;b--) c=cpow(a,c);
return c;
}
//---------------------------------------------------------------------------
void main()
{
// poistion (double)
dvec2 p=dvec2(p32);
p=(p/zoom)-p0; // x,y (-1.0, 1.0)
// position (float)
vec2 pp=vec2(p);
// [chose function]
// complex domain test function 1 (color wheel)
// vec2 a=cdiv(cmul(csub(cmul(pp,pp),vec2(1.0,0.0)),csqr(csub(pp,vec2(2.0,1.0)))),cadd(csqr(pp),vec2(2.0,2.0)));
// complex domain test function 2 (color wheel)
// vec2 a=pp; a=cln(a);
// exponentiation escape fractal 1 (color wheel)
// vec2 a=cpow(pp,vec2(100,0));
// exponentiation escape fractal 2 (color wheel)
// vec2 a=vec2(1.0,1.0); for (int i=0;i<100;i++) a=cpow(a,pp);
// exponentiation escape fractal 3 (color wheel)
// vec2 a=vec2(0.0,0.0),b=vec2(1.0,0.0); float r=0.5,rr=1.0,wt=0.1; for (int i=0;i<20;i++){ a+=rr*cexp(vec2(-b.y,b.x)*wt); b=cmul(b,pp); rr*=r; } a*=(1.0-r);
// tetration escape fractal (grayscale)
// vec2 a=ctet(pp,100);
// pentation escape fractal (grayscale)
vec2 a=pp; for (int i=0;i<20;i++) a=ctet(a,20); a*=100.0;
// [chose coloring method]
// grayscale based on escape
float r=0.2*length(a); r-=floor(r); r=0.25+0.75*r; col=vec4(r,r,r,1.0);
// RGB based on result
// col=vec4(a,a.x+a.y,1.0);
// result -> wavelength+intensity
// col=vec4(color_wheel(a),1.0);
// result -> spectral color
// col=vec4(color_spectral(a),1.0);
}
和四步预览：
这是我对比的：
它符合我的结果只是反映在两个x,y
所以我所做的是计算a^^100 where a是屏幕上片段的复杂域位置<-1,+1>和一些panning and zooming并渲染根据结果构造的颜色......
我还在那里留下了一个测试函数（不是分形），我用来测试着色方法和复杂的数学here第一个来自 Wiki，第二个是着色器结果（色轮）：
您可以像曼德尔布罗那样进行逃逸测试或其他任何算法来显示分形。
这里的着色选项是 tetration 的屏幕截图（我喜欢灰度）zoom=500.0 pos=-0.188418+0.234466i
最后是五旬节：
