朱莉娅用实数集计算

Question

我正在尝试制作一个 julia 集图像并用 C 进行编程。我曾尝试研究一种算法来创建 Julia 集的图像，但对大多数在线示例感到困惑（因为大多数示例似乎是复制的，但对逐步发生的事情几乎没有解释）。

我已经多次重写我的代码以适应我在网上找到的算法，但收效甚微。目前，我正在使用迭代函数来定义以下每个像素的颜色：

// Loop through each pixel to work out new image set
int x, y;
for(x = 0; x < SIZE; x++){
    for(y = 0; y < SIZE; y++){

        int i = iterate(x, y, maxIterations, c);
        image[x][y] = i % 256 + 255 + 255 * (i < maxIterations);

    }
}

/* Iterate function */
int iterate(int x, int y, int maxI, double c[]){
    double z, zx, zy, oldRe, oldIm;   //real and imaginary parts of new and old
    double xmin = -1.0, xmax = 1.0, ymin = -1.0, ymax = 1.0;
    int k; // number of times iterated

    //calculate the initial real and imaginary part of z
    // z0 = (x + yi)^2 + c = (0 + 0i) + c = c
    zx = 1.5*(x - SIZE/2)/(0.5*SIZE);
    zy = 1.0*(y - SIZE/2)/(0.5*SIZE);

    //start the iteration process
    for(k = 1; k < maxI; k++){
      //remember value of previous iteration
      oldRe = zx;
      oldIm = zy;
      z = zx*zx - zy*zy + c[0];

      //the actual iteration, the real and imaginary part are calculated
      zx = oldRe * oldRe - oldIm * c[1] + c[0];
      zy = 2 * oldRe * oldIm + c[1];
      zy = 2.0*zx*zy + c[1];
      zx = z;

      //if the point is outside the circle with radius 2: stop
      if((zx * zx + zy * zy) > 4) break;
}

/*
while(((zx * zx + zy * zy) > 4) && (k < maxI)) {
    //remember value of previous iteration
    z = zx*zx - zy*zy + c[0];
    zy = 2*x*y-c[1];
    zx = z;
    //if the point is outside the circle with radius 2: stop
}
return k;

} */ 注释掉的 while 循环是我在网上找到的第二种算法，但也是不正确的。我不确定我在等式中遗漏了什么。有什么猜测吗？

（我不是在寻找解决方案，只是一些可以将实数、虚数复数转换为我的边界框中的位置的伪代码，在本例中为：-1,1）。

Answer 1

这是我刚刚得到验证的一个学校项目的可读代码的重新转录

YOUR_CR 和 YOUR_CI 是 Julia 的两个可能变量， LOWER_LIMIT 和 UPPER_LIMIT（X 和 Y）对应于您要在其中绘制的窗口。代码应正确地将输出缩放到您要在其中绘制分形的窗口和缩放级别（限制比例是分形的数学限制，不是像素）。

# define X 0
# define Y 1

# define ZR 0
# define ZI 1
# define CR 2
# define CI 3

static int          iterations(long double coords[])
{
    long double     pos[4];
    int             i;
    double          tmp;

    pos[CR] = YOUR_CR;
    pos[CI] = YOUR_CI;
    pos[ZR] = coords[X];
    pos[ZI] = coords[Y];
    i = 0;
    while (i < MAX_ITERATIONS && (PYTHAGORE(pos[ZR], pos[ZI]) < 4))
    {
        tmp = pos[ZR];
        pos[ZR] = pos[ZR] * pos[ZR] - pos[ZI] * pos[ZI] + pos[CR];
        pos[ZI] = 2 * pos[ZI] * tmp + pos[CI];
        i++;
    }
    return (i);
}

# define TRANSLATE_X_PIXEL(a, b) ((int)(((a - LOWER_LIMIT_X) / (UPPER_LIMIT_X - LOWER_LIMIT_X)) * WIDTH))
# define TRANSLATE_Y_PIXEL(a, b) ((int)(((a - LOWER_LIMIT_Y) / (UPPER_LIMIT_Y - LOWER_LIMIT_Y)) * HEIGHT))

void                julia()
{
    long double     coords[2];
    int             itera;

    coords[X] = LOWER_LIMIT_X;
    while ((coords[X] += ((UPPER_LIMIT_X - LOWER_LIMIT_X) /
                    (double)(WIDTH + 1))) <= UPPER_LIMIT_X)
    {
        coords[Y] = LOWER_LIMIT_Y;
        while ((coords[Y] += ((UPPER_LIMIT_Y - LOWER_LIMIT_Y) /
                        (double)(HEIGHT + 1))) <= UPPER_LIMIT_Y)
        {
                itera = iterations(coords);
                set_pixel(TRANSLATE_X_PIXEL(coords[Y]), TRANSLATE_Y_PIXEL(coords[Y]), itera * DESIRED_COLOR);
        }
    }
}

Answer 2

注释掉的while 代码是错误的，因为它在虚部计算中使用x,y 而不是zx,zy。您的代码看起来好像您采用了一种方法并在其中粘贴了另一种方法。

你想要的经典二次Julia迭代，它是实部的实现

znext = z*z + c = (zx+*zy)*(zx+i*zy)+(cx+*cy)
      = (zx*zx-zy*zy+cx) + i*(2*zx*zy+cy) 

您不能进行就地替换，因为zx 和zy 都在两个公式中使用，并且需要保留它们的旧值。正确的代码应该是

for(k = 1; k < maxI; k++){
    // compute the new real part, store temporarily in helper variable
    double z = zx*zx - zy*zy + c[0];
    // compute the imaginary part, since no further operations, do it in place
    zy = 2.0*zx*zy + c[1];
    //store the new real part in the real part variable
    zx = z;
    //if the point is outside the circle with radius 2: stop
    if((zx * zx + zy * zy) > 4) break;
}

如果你存储正方形，你将需要更少的乘法，

for(k = 1; k < maxI; k++){
    double sqzx = zx*zx, sqzy = zy*zy;
    // compute the imaginary part first, since the real part only uses the squares, not zy
    zy = 2.0*zx*zy + c[1];
    // compute the new real part using the pre-computed squares
    zx = sqzx - sqzy + c[0];
    //if the point is outside the circle with radius 2: stop
    if( ( sqzx + sqzy ) > 4 ) break;
}