Unexpected error in Julia setup

I am playing with Mandelbrot and Julia sets and I am facing an interesting problem. The Mandelbrot set can be rendered in double precision as long as the scale is in the region of 2 ^ 56 anywhere. However, the Julia kit sometimes produces artifacts much earlier than around 2 ^ 20 magnification (see image below).

Julia set error

Interestingly, this only happens in some areas, but the rest of the image is fine and you can even keep zooming. This usually occurs in the center of clusters and near the source [0,0].

  • Is this really caused by double arithmetic errors?
  • Can these artifacts be avoided without using arbitrary precision arithmetic?
  • Why does this only happen in some areas? Are these areas somewhat special?

Coords

I don't have the coordinates of the image above, but another one can be found here (pictured below):

  • Center = [0, 0]
  • Scale 2 ^ 26
  • C = [-8.01030596311150589e-01, 1.56495138793530941e-01]

Arbitrary precision

Arbitrary precision code seems to help, but only slightly. 82 bits of arbitrary precision is better than just double, but then 232 bits are exactly the same as 82. Maybe my implementation is flawed? The only number that has a lower precision is the C argument, taken from the Mandelbrot set, and it has the precision needed to capture it at the depth where it was installed. I don't think it matters as adding extra zeros to create a high precision won't change the result. Immediate variables that compute the result are highly accurate.

Double precision:

Double precision

Arbitrary precision 82 bits:

Arbitrary precision 82 bits

Arbitrary precision of 232 bits:

Arbitrary precision 232 bits

code

This is the core of my code (unnecessary parts omitted):

public void UberCompute(UberComplex origin, UberComplex size,
        UberComplex initialCoord, int maxIters, double[,] result,
        ref bool stop) {

    int wid = result.GetLength(1);
    int hei = result.GetLength(0);
    int area = wid * hei;

    int precision = Math.Max(origin.Precision,
        initialCoord == null ? 0 : initialCoord.Precision);
    UberComplex coord = new UberComplex(precision);
    UberComplex step = new UberComplex(precision);
    UberFloat.Div(step.Re, size.Re, wid);
    UberFloat.Div(step.Im, size.Im, hei);

    UberFloat imed1 = new UberFloat(precision);
    UberFloat imed2 = new UberFloat(precision);
    UberFloat imed3 = new UberFloat(precision);

    for (int y = 0; y < hei; ++y) {
        // double yt = (double)y / hei;
        // double im = origin.Im + yt * size.Im;
        UberFloat.Mul(imed1, step.Im, y);
        UberFloat.Add(coord.Im, imed1, origin.Im);

        if (stop) {
            break;
        }

        for (int x = 0; x < wid; ++x) {
            // double xt = (double)x / wid;
            // Complex coord = new Complex(origin.Re + xt * size.Re, im);
            UberFloat.Mul(imed1, step.Re, x);
            UberFloat.Add(coord.Re, imed1, origin.Re);

            result[y, x] = UberIterate(initialCoord ?? coord,
                coord, maxIters, smooth, imed1, imed2, imed3);
        }
    }
}

public double UberIterate(UberComplex coord, UberComplex initialCoord, int maxIters,
        UberFloat imed1, UberFloat imed2, UberFloat imed3) {

    Contract.Requires(imed1.Precision == initialCoord.Precision);
    Contract.Requires(imed2.Precision == initialCoord.Precision);
    Contract.Requires(imed3.Precision == initialCoord.Precision);

    int precision = coord.Precision;
    UberFloat re = new UberFloat(precision, initialCoord.Re);
    UberFloat im = new UberFloat(precision, initialCoord.Im);

    int i = 0;
    do {
        // re * re + im * im > maxRadiusSq
        UberFloat.Mul(imed1, re, re);
        UberFloat.Mul(imed2, im, im);
        UberFloat.Add(imed3, imed1, imed2);
        if (imed3 > MAX_RADIUS_SQ) {
            break;
        }

        // newRe = re * re - im * im + coord.Re;
        UberFloat.Sub(imed3, imed1, imed2);
        UberFloat.Add(imed1, imed3, coord.Re);

        // im = 2.0 * re * im + coord.Im;
        UberFloat.Mul(imed2, re, im);
        UberFloat.Mul(imed3, imed2, 2);
        UberFloat.Add(im, imed3, coord.Im);

        // re = newRe;
        UberFloat.Swap(re, imed1);
    } while (++i < maxIters);

    if (i == maxIters) {
        return Double.NaN;  // Did not diverged.
    }

    return i;
}

      

Other software

I have tested Ultra Fractal and it has this problem:

Ultra fractal

+3


source to share





All Articles