Curtis T McMullen 的 C 程序^[1]

如何运行

julia.tar

描述^[2]

步骤

运行 tcsh
make
将 ps2gif 脚本复制到每个子目录
将 ./ 添加到每个运行脚本
运行

名称

"大象的眼睛" 停留在主心形的内部角 p/q = 主 Misiurewicz 点

图片

julia.tar
在斐波那契点 c = -1.8705286321646448888906 处的曼德勃罗集
飞机 Julia 集。C 是实轴上周期 3 分量的中心
在周期 16 分量附近（尾流 1/16）的断开连接的 Julia 集
来自尾流 1/19 的断开连接的（康托尔）Julia 集

kleinian.tar
球体上的克莱因群极限集
平面上肖特基（克莱因）群极限集

代码

复数

// from cx.c  programs by Curtis T McMullen
// http://www.math.harvard.edu/~ctm/programs/home/prog/julia/src/julia.tar

/* Compute sqrt(z) in the half-plane perpendicular to w. */
complex contsqrt(z,w)
complex z,w;
{
	complex cx_sqrt(), t;

	t = cx_sqrt(z);
	if (0 > (t.x*w.x + t.y*w.y)) 
		{t.x = -t.x; t.y = -t.y;}
	return(t);
}


/* Values in [-pi,pi]. */
double arg(z)  
	complex z;
{
	return(atan2(z.y,z.x));
}


complex cx_exp(z)
	complex z;
{
	complex w;
	double m;

	m = exp(z.x);
	w.x = m * cos(z.y);
	w.y = m * sin(z.y);
	return(w);
}

complex cx_log(z)
	complex z;
{
	complex w;

	w.x = log(cx_abs(z));
	w.y = arg(z);
	return(w);
}

complex cx_sin(z)
	complex z;
{
	complex w;

	w.x = sin(z.x) * cosh(z.y);
	w.y = cos(z.x) * sinh(z.y);
	return(w);
}

complex cx_cos(z)
	complex z;
{
	complex w;

	w.x = cos(z.x) * cosh(z.y);
	w.y = -sin(z.x) * sinh(z.y);
	return(w);
}

complex cx_sinh(z)
	complex z;
{
	complex w;

	w.x = sinh(z.x) * cos(z.y);
	w.y = cosh(z.x) * sin(z.y);
	return(w);
}

complex cx_cosh(z)
	complex z;
{
	complex w;

	w.x = cosh(z.x) * cos(z.y);
	w.y = sinh(z.x) * sin(z.y);
	return(w);
}

complex power(z,t)   /* Raise z to a real power t */
	complex z;
	double t;
{	double arg(), cx_abs(), pow();
	complex polar();

	return(polar(pow(cx_abs(z),t), t*arg(z)));
}

/*  Map points in the unit disk onto the lower hemisphere of the
    Riemann sphere by inverse stereographic projection.
    Projecting, r -> s = 2r/(r^2 + 1);  inverting this,
    s -> r = (1 - sqrt(1-s^2))/s.   */
 
complex disk_to_sphere(z)
	complex z;
{       complex w;
        double cx_abs(), r, s, sqrt();
 
        s = cx_abs(z);
        if (s == 0) return(z);
                else r = (1 - sqrt(1-s*s))/s;
        w.x = (r/s)*z.x;
        w.y = (r/s)*z.y;
        return(w);
}

complex mobius(a,b,c,d,z)
        complex a,b,c,d,z;

{       return(divide(add(mult(a,z),b),
                   add(mult(c,z),d)));
}

逃逸时间和势

// from quad.c  programs by Curtis T McMullen
// http://www.math.harvard.edu/~ctm/programs/home/prog/julia/src/julia.tar

/* F : \cx-K(f) -> {Re z > 0} is the multivalued
   function given by lim 2^{-n} log f^n(z).
   Re F(z) = rate, 
   (Re F(z))/|F'(z)| = dist.
   Returns 1 if z escaped, 0 otherwise.
*/

escape(z,c,iterlim,rate,dist)
	complex z,c;
	int iterlim;
	double *rate, *dist;
{
	int i;
	double fp, r, s;

	s = 1.0;
	fp = 1.0;
	for(i=0; i<iterlim; i++)
	{	r = cx_abs(z);
		if(r > LARGE)
		{	*rate = log(r)/s;
			*dist = r*log (r)/fp;
			return(1);
		}
		fp = fp * 2.0 * r;
		z = add(mult(z,z),c);
		s = s * 2.0;
	}
	return(0);
}

外部射线

// from quad.c  programs by Curtis T McMullen
// http://www.math.harvard.edu/~ctm/programs/home/prog/julia/src/julia.tar

/* CAUTION is 2^{1/n} */
#define LOG2 .69314718055994530941
#define SQRT2 1.414213562373095
#define LARGE 200.0
#define LENGTH 200
#define CAUTION 0.9330329915368074
#define CIRCLEPTS 1000
#define EPS 0.000001
#define RAYEPS 1e-20
#define JFAC 1.0

/* Defaults */
#define C {0, 1}
#define CENTER {0.0,0.0}
#define RADIUS 1.5
#define PCLIM 0
#define SUBDIVIDE 8
#define ITERLIM 50  


/* Analytically continue Riemann mapping */
/* F : unit disk -> \cx-K(f) */
/* Assuming F(r,ang) is approx. z, find exact value*/

complex riemann(r,ang,z)
	double r;
	rational ang;
	complex z;
{
	complex orb[LENGTH], neworb[LENGTH];
	int i, n;

	orb[0] = z;
	n = 0;
	for (i=0; i<LENGTH-1; i++)
	{	if(cx_abs(orb[i]) > LARGE) break;
		orb[i+1] =
		add(mult(orb[i],orb[i]),c);
		n++;
		r = 2*r;
		ang.num = (2*ang.num)%ang.den;
	}
	neworb[n] = polar(pow(2.0,r),2*PIVAL*ang.num/ang.den);

	for(i=n; i>0; i--)
	{	neworb[i-1] =
		contsqrt(sub(neworb[i],c),orb[i-1]);
	}
	return(neworb[0]);
}


/* Set the color on a line joining z to w */
scan_line(z,w,col)
	complex z, w;
	int col;
{	
	int i, n;
	complex u;

	n = 4*cx_abs(sub(z,w))/pixrad;
	if(n == 0)
	{	scan_set(z,col);
		return;
	}
	for(i=0; i<=n; i++)
	{	u.x = (i*z.x + (n-i)*w.x)/n; 	
		u.y = (i*z.y + (n-i)*w.y)/n; 	
		scan_set(u,col);
	}
}


/* Draw an external ray */
draw_ray(level,ang)
	double level;
	rational ang;
{
	complex lastz, z;
	double r;

	r = level;
	while(pow(2.0,r) < LARGE) r *= 2;
	z = polar(pow(2.0,r),2*PIVAL*ang.num/ang.den);
	do
	{	r = r*CAUTION;
		lastz = z;
		z = riemann(r,ang,lastz);
		if(r < level*0.99) 
			scan_line(lastz,z,RAY);
	} while(r > RAYEPS);
}

使用 Curtis McMullen 的逆 Boettcher 映射绘制动态外部射线

该方法基于 Curtis McMullen^[3] 的 C 程序及其由 Matjaz Erat^[4] 编写的 Pascal 版本

这里有两个平面^[5]

w 平面（或 f0 平面）
z 平面（fc 平面的动态平面）

该方法包含 3 个主要步骤

计算圆的外部射线在 $f_{0}\,$ 平面上的某些 w 点，这些点对应于角度 $\theta \,$ 和不同的半径（光栅化）

w_{0}=R*e^{2\pi it}\,

其中

1<R<\infty \,

使用逆 Boettcher 映射 $\Psi _{c}=\Phi _{c}^{-1}\,$ 将 w 点映射到 z 点

z_{0}=\Psi _{c}(w_{0})=f_{c}^{-n}(w_{0}^{2^{n}})\,

绘制 z 点（并使用线段连接它们（线段是直线上由两个不同的端点界定的部分^[6]））

第一步和最后一步很简单，但第二步并非如此，需要更多解释。

光栅化

对于 $f_{0}\,$ 平面中给定的外部射线，射线上的每个点都有

常数 $t\,$ （以转数表示的外部角度）
可变半径 $R\,$

所以 $w\,$ 射线的点由半径 $R\,$ 参数化，可以使用复数的指数形式

w=F_{t}(R)=R*e^{2*\Pi *i*t}\,

可以使用线性尺度沿着射线前进

t:1/3; /* example value */
R_Max:4;
R_Min:1.1;
for R:R_Max step -0.5 thru R_Min do w:R*exp(2*%pi*%i*t);
/* Maxima allows non-integer values in for statement */

它会给出一些 w 点，它们之间具有相等的距离。

另一种方法是使用非线性尺度。

为了做到这一点，我们引入浮点指数 $r\,$ ，使得

R=2^{r}\,

和

w=G_{t}(r)=2^{r}*e^{2*\Pi *i*t}\,

为了在 $f_{0}\,$ 平面上计算角度为 $t\,$ 的外部射线的 w 点，使用以下 Maxima 代码

t:1/3; /* external angle in turns */
/* range for computing  R ; as r tends to 0 R tends to 1 */
rMax:2; /* so Rmax=2^2=4 /
rMin:0.1; /* rMin > 0   */
caution:0.93; /* positive number < 1 ; r:r*caution  gives smaller r */
r:rMax;
unless r<rMin do
( 
 r:r*caution, /* new smaller r */ 
 R:2^r, /* new smaller R */
 w:R*exp(2*%pi*%i*t) /* new point w in f0 plane */
);

在这种方法中，点之间的距离并不相等，而是与填充的朱利亚集边界到点的距离成反比。

这很好，因为这里射线具有更大的曲率，因此曲线将更加平滑。

映射

从 $f_{0}\,$ -平面到 $f_{c}\,$ -平面的映射包括 4 个小的步骤

在 $f_{0}\,$ 平面上进行前向迭代

w_{n}=w_{0}^{2^{n}}=R^{2^{n}}*e^{2\pi i2^{n}t}\,

直到 $w_{n}\,$ 接近无穷大

abs(w_{n})>R_{limit}\,

切换平面（从 $f_{0}\,$ 到 $f_{c}\,$ )

z_{n}=w_{n}\,

（因为在这里，靠近无穷大： $w=z\,$ ）

在 $f_{c}\,$ 平面上进行反向迭代，迭代次数相同（ $n\,$ ）
最后一个点 $z_{0}\,$ 位于我们的外部射线上

步骤 1、2 和 4 很容易。第三步不容易。

z_{n-1}=f_{c}^{-1}(z_{n})={\sqrt {z_{n}-c}}\,

反向迭代使用复数的平方根。它是一个双值函数，因此反向迭代会产生二叉树。

在没有额外信息的情况下，无法在这样的树中选择好的路径。为了解决这个问题，我们将使用以下两点

无穷大的吸引盆的等度连续性
$f_{0}\,$ 和 $f_{c}\,$ 平面之间的共轭关系

无穷大的吸引盆的等度连续性

无穷大的吸引盆（填充的 Julia 集的补集）包含在正向迭代下趋于无穷大的所有点。

$A_{f_{c}}(\infty )\ {\overset {\underset {\mathrm {def} }{}}{=}}\ \{z\in \mathbb {C} :f_{c}^{(k)}(z)\to \infty \ as\ k\to \infty \}.$

无穷大是一个超吸引不动点，所有点的轨道都具有相似的行为。换句话说，如果两个点的初始位置很近，则假设它们的轨道会保持靠近。

这就是等度连续性（与正规性相比较）。

在 $f_{c}\,$ 平面上，可以使用射线上先前点的正向轨道来计算下一个点的反向轨道。

算法的详细版本

计算射线上的第一个点（从无穷大附近开始，朝向 Julia 集移动）

w=2^{r}*e^{2\Pi it}\,

其中

r=r_{max}\,

这里可以轻松地切换平面

z=w\,

这是我们射线的第一个 z 点。

计算射线的下一个 z 点
- 为 $r=r*caution\,$ 计算射线的下一个 w 点
- 计算 2 个点的向前迭代：之前的 z 点和实际的 w 点。保存 z 轨道和最后一个 w 点
- 切换平面并将最后一个 w 点用作起点： $z_{n}=w_{last}\,$
- 使用之前 z 点的向前轨道，对新的 $z_{n}\,$ 进行向后的迭代，以获得新的 $z_{0}\,$
- $z_{0}\,$ 是我们射线的下一个 z 点
依此类推（下一个点），直到 $r<r_{Min}\,$

Maxima CAS 源代码

/* gives a list of z-points of external ray for angle t in turns and coefficient c */
GiveRay(t,c):=
block(
 [r],
 /* range for drawing  R=2^r ; as r tends to 0 R tends to 1 */
 rMin:1E-20, /* 1E-4;  rMin > 0  ; if rMin=0 then program has infinity loop !!!!! */
 rMax:2,  
 caution:0.9330329915368074, /* r:r*caution ; it gives smaller r */
 /* upper limit for iteration */
 R_max:300,
 /* */
 zz:[], /* array for z points of ray in fc plane */
 /*  some w-points of external ray in f0 plane  */
 r:rMax,
 while 2^r<R_max do r:2*r, /* find point w on ray near infinity (R>=R_max) in f0 plane */
 R:2^r,
 w:rectform(ev(R*exp(2*%pi*%i*t))),
 z:w, /* near infinity z=w */
 zz:cons(z,zz), 
 unless r<rMin do
 ( /* new smaller R */
  r:r*caution,  
  R:2^r,
  /* */
  w:rectform(ev(R*exp(2*%pi*%i*t))),
  /* */
  last_z:z,
  z:Psi_n(r,t,last_z,R_max), /* z=Psi_n(w) */
  zz:cons(z,zz)
 ),
 return(zz)
)$

逃逸时间的等高线

// from quad.c  programs by Curtis T McMullen
// http://www.math.harvard.edu/~ctm/programs/home/prog/julia/src/julia.tar

/* Color a level curve of escape */
draw_level(level)
	double level;
{
	double r;
	complex lastz, w, z;
	static rational zero={0,1};
	rational ang;

	r = level;
	while(pow(2.0,r) < LARGE) r *= 2;
	z = polar(pow(2.0,r),0.0);
	do
	{	r = r * CAUTION;
		lastz = z;
		z = riemann(r,zero,lastz);
	} while(r > level);

	lastz = z;
	ang.den = CIRCLEPTS;
	for(ang.num=0; ang.num<=CIRCLEPTS; ang.num++)
	{	z = riemann(level,ang,lastz);
		if(ang.num>0) scan_line(z,lastz,LEVEL);
		lastz = z;
	}
}

颜色

// from quad.c  programs by Curtis T McMullen
// http://www.math.harvard.edu/~ctm/programs/home/prog/julia/src/julia.tar

/* Colors */
#define FILLED_JULIA 0
#define JULIA 1
#define PC_SET 2 
#define RAY 3
#define LEVEL 4
#define RESERVED 5
#define FREE (256-RESERVED)

set_gray(g)
	double g;
{
	int i;

	scan_gray(FILLED_JULIA,g);
	scan_gray(JULIA,0.0);
	scan_gray(PC_SET,0.0);
	scan_gray(RAY,0.0);
	scan_gray(LEVEL,0.0);
	for(i=0; i<FREE; i++)
		scan_gray(RESERVED+i,1.0);
}

set_colors()
{
	int i;

	scan_hls(FILLED_JULIA,0.0,70.0,0.0);
	scan_rgb(JULIA,0,0,0);
	scan_rgb(PC_SET,0,255,0);
	scan_rgb(RAY,0,0,0);
	scan_rgb(LEVEL,0,0,0);
	for(i=0; i<FREE; i++)
		scan_hls(RESERVED+i,
	((FREE-i)*360.0/FREE),50.0,100.0);
}

color(z)
	complex z;
{
	complex grad;
	double rate, dist;
	int i;

	if(escape(z,c,iterlim,&rate,&dist))
	{	if(dist < jfac*SQRT2*pixrad) 
			return(JULIA);
		i = subdivide * log(rate)/LOG2;
		while (i<0) i += FREE;
		return(RESERVED + i%FREE);
	}
	return(FILLED_JULIA);
}

// from scan.c  programs by Curtis T McMullen
// http://www.math.harvard.edu/~ctm/programs/home/prog/julia/src/julia.tar

scan_hls(i,h,l,s)
	int i;
	double h, l, s;
{
	int r, g, b;

	hlsrgb(h,l,s,&r,&g,&b);
	red[i] = r; green[i] = g; blue[i] = b;
	cmap = COLOR;
}



/* Convert hue/lightness/saturation to rgb */
/* Domain is [0,360] x [0,100] x [0,100]. */
/* Range is [0,255] x [0,255] x [0,255]. */

hlsrgb(h, l, s, r, g, b)
	double h, l, s;
	int *r, *g, *b;
{
	double M, m;
	double floor(), bound(), gun();

	h = h - floor(h/360) * 360;
	l = bound(0.0, l/100, 1.0);
	s = bound(0.0, s/100, 1.0);
	M = l + s*(l > 0.5 ? 1.0 - l : l);
	m = 2*l - M;
	*r = 255*gun(m, M, h);
	*g = 255*gun(m, M, h-120);
	*b = 255*gun(m, M, h-240);
}

double gun(m, M, h)
double m, M, h;
{
	if(h < 0)
		h += 360;
	if(h < 60)
		return(m + h*(M-m)/60);
	if(h < 180)
		return(M);
	if(h < 240)
		return(m + (240-h)*(M-m)/60);
	return(m);
}

double bound(low, x, high)
double low, x, high;
{
	if(x < low)
		return(low);
	if(x > high)
		return(high);
	return(x);
}

参考文献

[1] Curtis T McMullen 编写的程序

[2] 如何运行程序

[3] Curtis T McMullen 编写的 c 程序（Julia.tar.gz 中的 quad.c）

[4] Matjaz Erat 编写的二次多项式

[5] 维基百科：复杂二次多项式 / 平面 / 动态平面

[6] 维基百科：线段

[1]

[2]

[3]

[4]

[5]

[6]