The following C++ code "knows" the cubic Julia set z_n+1=z_n³+(0.5625+0.625i).

It asks the user to enter a complex number and outputs with mathematical certainty:

- interior - the point is within the interior of the filled-in set
- exterior - the point escapes to infinity
- unknown - that's the price to pay for the certainty in the other cases.

The routine can be used as a filter if one e.g. wants to point-sample this set but is only interested in a specific region of the complex plane.

The "knowledge" is based upon the brilliant true shape algorithm by Figueiredo et. al. (see literature). I implemented their algorithm, computed the Julia set to a trustworthy image and derived a number of polygons describing the filled-in set. (I omit the details to avoid double-posting, a description can be found at: https://fractalforums.org/fractal-mathematics-and-new-theories/28/julia-sets-true-shape-and-escape-time/2725/msg14770)

Technical notes:

• I did my best to ensure the correctness of the implementation to the question at hand: mathematical certainty, and performed every computational and manual check I could think of.
• All calculations were done in a conservative manner to ensure no falsely classified number (rounding, neighbours etc).
  
• I use double precision in the code, however the test will work correctly for any number of any precision and still retain the mathematical guarantee.
  
• The image was computed in a spatial resolution of 6*2^-14 per pixel.
• The program reads all files named intpolyNNNN (151 interior polygons) and extpolyNNNN (1 exterior), attached as a zip file.
  
• The code was intended to be concise and does not perform checks for file correctness, pointer validity, array index etc.
  
• The main routine is called "jsoracle(x,y)" which performs the point-in-polygon tests for the complex number x+i*y.
  
• There are no special functions used. I believe it should be compilable (maybe after changing header names) on any platform.
  

Literature:

• Figueiredo, L., Nehab, D. "Images of Julia sets that you can trust", 2013.
• Figueiredo, L., Nehab, D. "Rigorous bounds for polynomial Julia sets", 2016.
• Galetzka, G., Glauner, P. "A Simple and Correct Even-Odd Algorithm for the Problem for Complex Polygons", 2017.

#include "stdio.h"
#include "math.h"

enum { PIP_INTERIOR=1, PIP_EXTERIOR, PIP_UNKNOWN, PIP_BOUNDARY=4 };
const int MAXPOLYGON=1024;

class Polygon {
public:
int pointcount,*X,*Y;

Polygon();
virtual ~Polygon();
void setsize(const int);
int read(const char*);
};

double globalcx0,globalcx1,globalcy0,globalcy1;
Polygon intp[MAXPOLYGON];
Polygon extp[MAXPOLYGON];
int globaldenominator,intpcount=0,extpcount=0;

int Polygon::read(const char* afn) {
FILE *f=fopen(afn,"rt");
if (!f) return 0;
fscanf(f,"%i",&globaldenominator);
fscanf(f,"%lf,%lf,%lf,%lf",&globalcx0,&globalcx1,&globalcy0,&globalcy1);
fscanf(f,"%i",&pointcount);
setsize(pointcount);
for(int i=0;i<pointcount;i++) fscanf(f,"%i,%i",&X[i],&Y[i]);
fclose(f);

return 1;
}

Polygon::Polygon() {
pointcount=0;
X=Y=NULL;
}

Polygon::~Polygon() {
if (X) delete[] X;
if (Y) delete[] Y;
}

void Polygon::setsize(const int a) {
if (X) delete[] X;
if (Y) delete[] Y;
X=new int[a];
Y=new int[a];
pointcount=a;
}

int minimum(const int a,const int b) {
if (a<b) return a;
return b;
}

int maximum(const int a,const int b) {
if (a>b) return a;
return b;
}

int point_in_polygon(Polygon& polygon,const int ax,const int ay) {
int even=1;
for(int i=1;i<polygon.pointcount;i++) {
if (polygon.X[i] == polygon.X[i-1]) {
int miy=minimum(polygon.Y[i],polygon.Y[i-1]);
int may=maximum(polygon.Y[i],polygon.Y[i-1]);

if (
(polygon.X[i] == ax) &&
(miy <= ay) && (ay <= may)
) return PIP_BOUNDARY;

if (
(ax < polygon.X[i]) &&
(miy <= ay) && (ay <= may)
) {
if (ay == polygon.Y[i]) {
int y0,y1,y2;
y1=polygon.Y[i];
if (i>0) y0=polygon.Y[i-1]; else y0=polygon.Y[polygon.pointcount-2];
if (i<(polygon.pointcount-1)) y2=polygon.Y[i+1]; else y2=polygon.Y[1];

if (
( (y0 < y1) && (y1 < y2) ) ||
( (y0 > y1) && (y1 > y2) )
) even=1-even;
} else
if ((miy < ay) && (ay < may)) even=1-even;
}
} else
if (polygon.Y[i] == polygon.Y[i-1]) {
int minx=minimum(polygon.X[i],polygon.X[i-1]);
int maxx=maximum(polygon.X[i],polygon.X[i-1]);

if ((polygon.Y[i] == ay) && (minx <= ax) && (ax <= maxx)) return PIP_BOUNDARY;

if ((ay == polygon.Y[i]) && (minx > ax)) {
int y0,y1,y2;
if (i > 1) {
y0=polygon.Y[i-2];
y1=polygon.Y[i];
} else {
if (i==1) {
y0=polygon.Y[polygon.pointcount-2];
y1=polygon.Y[i];
} else {
y0=polygon.Y[polygon.pointcount-2];
y1=polygon.Y[0];
}
}

if (i < (polygon.pointcount-1)) y2=polygon.Y[i+1]; else y2=polygon.Y[1];

if (( (y0 < y1) && (y1 < y2) ) || ( (y0 > y1) && (y1 > y2) )) even=1-even;
}
}
}

if (even > 0) return PIP_EXTERIOR;
return PIP_INTERIOR;
}

int jsorakel(const double ax,const double ay) {
int px=(int)floor(ax*globaldenominator);
int py=(int)floor(ay*globaldenominator);
int mx=2+(int)floor(globalcx1-globalcx0);
int my=2+(int)floor(globalcx1-globalcx0);

int ergint=PIP_UNKNOWN;
for(int i=0;i<intpcount;i++) {
int ic=0;
for(int dy=-my;((ic>=0)&&(dy<=my));dy++) {
for(int dx=-mx;((ic>=0)&&(dx<=mx));dx++) {
if (point_in_polygon(intp[i],px+dx,py+dy) == PIP_INTERIOR) ic++;
else ic=-1;
}
}

if (ic==((mx+mx+1)*(my+my+1)) ) return PIP_INTERIOR;
}

int ergext=PIP_UNKNOWN;
for(int i=0;i<extpcount;i++) if (extp[i].pointcount>0) {
int ic=0;
for(int dy=-my;((ic>=0)&&(dy<=my));dy++) {
for(int dx=-mx;((ic>=0)&&(dx<=mx));dx++) {
if (point_in_polygon(extp[i],px+dx,py+dy) == PIP_EXTERIOR) ic++;
else ic=-1;
}
}

if (ic!=((mx+mx+1)*(my+my+1)) ) return PIP_UNKNOWN; else ergext=PIP_EXTERIOR;
}

if ((ergext == PIP_EXTERIOR) && (extpcount>0)) return PIP_EXTERIOR;

return PIP_UNKNOWN;
}

int main() {
char tmp[1024];
intpcount=0;
for(int i=0;i<MAXPOLYGON;i++) {
sprintf(tmp,"intpoly%04i",intpcount);
if (intp[intpcount].read(tmp) <= 0) break;
intpcount++;
}

extpcount=0;
for(int i=0;i<MAXPOLYGON;i++) {
sprintf(tmp,"extpoly%04i",extpcount);
if (extp[extpcount].read(tmp) <= 0) break;
extpcount++;
}

double re,im;
printf("Enter real part> "); scanf("%lf",&re);
printf("Enter imag part> "); scanf("%lf",&im);

switch (jsorakel(re,im)) {
case PIP_INTERIOR: printf("definitely interior.\n"); break;
case PIP_EXTERIOR: printf("definitely exterior.\n"); break;
default: printf("unknown\n"); break;
};

return 0;
}


Linkback: https://fractalforums.org/index.php?topic=2852.0

So... you create "outline" polygons of the shape to a certain resolution and then using a bounds check on the polygon instead of the actual formula. Quite ingenious even if a bit error-prone. Congrats on having some success with it!

Thanks!
I'm currently working on speeding the construction up. It takes about one whole day of computation to derive the polygons, and the outside one quite often makes some loops which I have to remove manually. I hope I can get rid of that to completely automate this process.

...even if a bit error-prone.

What do you refer to exactly by "error-prone"?

The algorithm's answers "inside" and "outside" come with a mathematical guarantee because of the underlying calculation process: cell mapping and interval arithmetics, so a point in the image is actually a whole square and its color is applicable to the fate of every number therein (provided my implementation is correct, of course). So there's no error in the sense that a user-entered number is falsely classified as interior or exterior.

Or do you mean that there is always a (not so tiny) set of complex numbers that cannot be judged and the algorithm's output is "unknown"? That's unfortunately the price one has to pay to be mathematically certain in the other two cases.

Relevant to the algorithm's output is "outside the outside polygon" (currently I'm working only with connected shapes so I only need one outside polygon) and "inside one inside polygon", all other point check results lead to the answer "unknown".

Have you checked the speed of this method vs just performing the math calculations outright?

I guess you refer to distance estimation? Since this method uses some sort of arbitrary parameter (large enough radius) and is mathematically only true in the limit at infinity its calculation would not be certain in the Figueiredo sense, falsely classified numbers will come up - but could not be detected (and I haven't found a reference as to how big the error of distance estimation is).

The point test itself is quite slow if checking a large set of numbers (it took 15 minutes to check about a quarter million points). But that could be sped up using a better data structure - I simply use an array and walk through it again and again. A tree sorted by x-coordinate or a hash-table would greatly reduce that time. And, of course, parallelization is possible as checking different numbers is independent from each other.

From a practical point of view one could always use the trustworthily computed Julia set image itself and simply calculate in which pixel the complex number of interest lies and then check its color (so direct array access - nothing beats that!). For the images I currently can compute polygons for (2^14 x 2^14 pixels) that would only need about 16 MBytes (assuming 4 pixels per byte), so not a problem at all.

But from a theoretical point of view I like the polygons better - and why take the easy road if there's an intersting detour?

@claude: I've read that article a couple of times but no luck so far putting those theoretical results into code for me. In fact, I'm a little confused, especially by the last statement about "For every epsilon...": This sounds to me like a tautology: One has to know that the true distance is between two numbers (how?) and then can deduce that it lies between two (tighter?) values. But the part about the error estimate for the escape iteration count is very interesting.

I uploaded the code to construct the polygons. Binaries (win32, win64, wine/ubuntu 12) - both with and without compiler optimization flags (sse4.1, fastmath), C++ source code and example images can be found here:

https://github.com/marcm200/julia-polygon

Features:

• command-line tool as always
• construct interior and exterior polygons
• quality checking the polygons
• oracle function: use the polygons to determine whether a complex number is in the set or not (or unknown)
• construction is very fast, quality checking needs time: 2k images are done in seconds, 4k need some minutes, 16k half an hour, up to 64k is possible

So far, with the new construction routine and the standard parameters, I haven't encountered a polygon that did not pass the tests. I hope that stays this way.