sandbox/Antoonvh/primes.c

A 4^2 Morton curve. Basilisk employs a diagnoally mirrored version, i.e. (capital) N-order curve to be precise. Image courtesy of Asger Hoedt

A 42 Morton curve. Basilisk employs a diagnoally mirrored version, i.e. (capital) N-order curve to be precise. Image courtesy of Asger Hoedt

Prime numbers along a Z-order space-filling curve

Stanislaw Ulam once dicided that is was a good idea to order the prime numbers (primes) along a spiralling space-filling curve. Doing so, he discovered what is now known as the Ulam’s spiral. It is attractive to use the Basilisk toolbox to study the behaviour of primes along a Z-order space-filling curve. If you are curious what role this curve plays within Baslisk, have a look here.

#include "grid/quadtree.h" //<- For it's 2D Z-order indexing iterator
#include "utils.h" //<- For visualization purposes
#include "tag.h" //<- For finding connected regions. 
int i=9;

Find primes

We need to find all primes upto n and store them in an array b of length n. This is done by using Eratosthenes’ sieve, an ancient algorithm that was not developed with computational efficientcy nor paralellization in mind.

void getprimes(int b[],int n){
  b[0]=0; //zero is not a prime
  b[1]=0; //one is not a prime
  for (int j=2;j<n;j++){
    b[j]=j;
  }
  int j=2;
  while(j<=ceil(sqrt(n))){
    int num = b[j];
    if (b[j]!=0){
      int ind = 2*num;
      while (ind<=n){
        b[ind]=0;
        ind+=num;
      }
    }
    j++;
  }
}

The loop

The Z-order (or (capital)N-order) space-filling curve seems to be most suitable to study N×N grids where N is a power of 2 (e.g. 2,4,128 etc.). In order to study the behaviour of the prime-number locations in the Z-order indexed grid we perform our analysis on an increasingly larger grid.

int main(){
  char name[100];
  FILE * fp2 = fopen("connectedregions","w");
  static FILE * fp1 = popen ("ppm2gif --delay 200 > g.gif ", "w");

Loop over increasingly larger grids

  for (int maxlevel=1;maxlevel<=i;maxlevel++){  
    init_grid(1<<maxlevel);
    int d[1<<(maxlevel*dimension)];
    getprimes(d,1<<(maxlevel*dimension));
    int m = 1;

Mark cells at prime locations along the space-filling curve:

    scalar field[];
    foreach(){
      field[]=d[m++];
    }

We can view the result by using this line of output.

    output_ppm (field, fp1,n=pow(2,i),min=0,max=1,map = gray);

Here it is, for all the iterations:

White indicates locations of the prime numbers

White indicates locations of the prime numbers

Connected regions

We see that there exist connected regions, these could potentially have interesting properties. Hence, we tag the connected regions with a unique tag.

    int am = tag(field);
    int regsize[am];
    for (m=0;m<am;m++)
      regsize[m]=0;

And we store the length (i.e. size) of each region.

    foreach(){
      if (field[]>0)
        regsize[(int)field[]]++;
    }
    sprintf(name,"prime%d.dat",maxlevel);
    FILE * fp = fopen (name, "w");
    for (m=1;m<am;m++)
      fprintf(fp,"%d\n",regsize[m]);

For the largest grid (i.e. 512×512 points), we plot some statistics on the lengths of the connected regions.

A truly remarkable feature is that the frequency (f(i)) of the region’s lengths (i) appears to scale with the region’s length, according to;

f(i)e2i.

We also log the number of connected regions for each grid-refinement iteration.

    fprintf(fp2,"%d\t%d\n",maxlevel,am+1);
  }

The result of this procedure are plotted below:

Again we obtain a result that I did not expect, but I am not an expert.

}

The next step

The next step is to increase the dimensionality of our Z-order-indexing curve and do the same analysis in 3 dimensions. The results are presented here.