[Math] Random binary array shows patterns around prime numbers

algorithmsprime numbers

First post, so please let me know if I'm doing something wrong or if this question does not belong here.

I have been toying with java to visualize an interesting 2D binary array I thought of today in my math class. My thought was to start with one origin cell, and set all cells in an infinite grid equal to 1. Then, starting with the 8 cells surrounding the origin and moving outward, each cell would be evaluated by the following rule:

-If the cell at (x,y) == 1, then any cell at (nx,ny) = 0, for all positive integers n > 1

Basically I was thinking that if (1,0) was on, then all of the other cells behind it – (2,0), (3,0), etc… would be "in it's shadow" and therefore equal to zero. I figured the result would have some significance with prime numbers, and it did.

I wrote a program to plot all the cells in the 4th quadrant, up to about x=2000, y=1000 and took screenshots at various zoom levels. Here are the results:

enter image description here

The top left cell (0,0) is arbitrarily set to 1. Here's the algorithm:

    boolean [][] grid = new boolean[height][width];

    for(int x=0; x<width; x++)
    {
      for(int y=0; y<height; y++)
      {
          grid[x][y] = true;
      }
    }

    for(int x=0; x<width; x++)
    {
      for(int y=0; y<height; y++)
      {
          if(grid[x][y])
          {
              int xt = 2*x;
              int yt = 2*y;

              while((xt < width && yt < height) && (xt != 0 || yt != 0))
              {
                  grid[xt][yt] = false;

                  xt += x;
                  yt += y;
              }
          }
      }
    }

On a large scale, it looks somewhat uniform, but on a small scale it looks random (or at least I can't see the pattern).

Here's my questions:

1 – Has this "function" been researched before? It seems interesting to me, but I can't find anything online about it.

2 – Where is this randomness coming from? The code is so simple, but the result looks very complex.

3 – Is there any mathematical way for me to search this grid for patterns more efficiently than the naked eye? I can see that the vertical and horizontal lines are most prominent along the prime numbers, which is interesting.

Thank you, and sorry if this is a dumb question. I couldn't think of any other place to ask it.

Best Answer

This pattern is usually described in the following way: Suppose you stand at the origin, and at each lattice point $(x,y)\in\mathbb{Z}^2$ there is a thin pole stuck into the ground. Which lattice points are visible, i.e. which poles are not hidden behind other poles? The answer is that $(x,y)$ is visible if and only if $x$ and $y$ have no common factor $>1$, since if $d$ is a common factor, then $(x,y)$ is hidden by $(x/d, y/d)$.

If you search for visible lattice points, you will find articles ranging from introduction to research papers.

Related Solutions

[Math] Prime numbers in binary.

In case the exact number helps, Mathematica can compute PrimePi[2^16] - PrimePi[2^15 - 1] to be $3030$. Choosing one-hundred odd integers uniformly at random from $[2^{15},2^{16}]$, the expected number of primes among them is $20200/10923\approx 18.4931$.

Calculating the approximation with the prime number theorem, I get approximately $16.8315$. The logarithms are supposed to be base-$e$: if I redo the calculation with base-$10$ logarithms, I get approximately $38.7559$.

Optimize table layout

I tried to implement Rahul's suggestion to view it as a convex optimization problem. The results are mixed. I can easily do small tables like 30 by 30, but 300 by 300 can be done with only about 1% precision if you are willing to wait 1 minute and getting down from there will take eternity. That is primarily due to the inefficiency of the solution finder I'm using (which is more or less just cycling over variables and optimizing over certain subsets of them; I wish I could find a better way or, at least, accelerate convergence somewhat). Nevertheless it is a good exercise in convex programming, so I'll post the details here. The algorithm can be modified to take into account "natural" restrictions of the kind $w_j\ge W_j$ or $h_i\ge H_i$ (width/height should not be too small) and the modification has pretty much the same rate of performance as far as I can tell from simulations, but I'll restrict myself to the original question here.

Let $w_j$ be the unknown widths and $a_{ij}$ be the known areas. We want to minimize $\sum_i\max_j \frac{a_{ij}}{w_j}$. It is useful to consider the dual problem. I will spare you from the general theory of duality and will just note that $$ \max_j \frac{a_{ij}}{w_j}=\max\left\{\sum_j b_{ij}\frac{a_{ij}}{w_j}:b_{ij}\ge 0, \sum_j b_{ij}=1\right\} $$ so if we consider all admissible vectors $w=(w_1,\dots,w_n)$ (non-negative entries, total sum $1$) and all admissible matrices $b=(b_{ij})$ (non-negative entries, all row sums equal to $1$), we can write our problem as that of finding $$ \min_w\max_b \sum_{i,j} b_{ij}\frac{a_{ij}}{w_j}\,. $$ The dual problem to that is finding $$ \max_b \min_w\sum_{i,j} b_{ij}\frac{a_{ij}}{w_j}\,. $$ The inner $\min_w$ is here easy to find: if we denote $S_j=\sum_i b_{ij}a_{ij}$, then it is just $(\sum_j \sqrt{S_j})^2$ with unique optimal $w_j$ proportional to $\sqrt{S_j}$.

There are two things one should understand about duality. The first one is that every admissible matrix $b$ (computed or just taken from the ceiling) can serve as the certificate of the impossibility to do better than a certain number in the original problem, i.e., the minimax is never less than the maximin. This is pretty trivial: just use the given $b$ to estimate the minimax from below. The second one is that the true value of minimax is actually the same as the true value of maximin (under some mild assumptions that certainly hold in our case). This is a somewhat non-trivial statement.

Together these two observations allow one to use the following strategy. We shall try to solve the dual problem. For every approximation $b$ to the solution, we will look at the easily computable lower bound $(\sum_j\sqrt{S_j})^2$ it produces and at the corresponding minimizer $w$. For that $w$ we can easily compute the original function $\sum_j\max_i\frac{a_{i,j}}{w_j}$. If its value is reasonably close to the lower bound, we know that we should look no further.

Now, of course, the question is how to maximize $\sum_j\sqrt S_j$ under our constraints on $b$. It doesn't look like an attractive problem because the number of unknowns increased from $n$ to $mn$. Still, one can notice that if we fix all rows of $b$ except, say, the $i$'th one, then the optimization of the $i$'th row is rather straightforward. Indeed, the corresponding problem is of the following kind:

**Find $\max\sum_j\sqrt{a_jb_j+c_j}$ where $a_j,c_j\ge 0$ are given and $b_j\ge 0$ are the unknowns subject to the constraint $\sum_j b_j=1$. Using the standard Lagrange multiplier mumbo-jumbo, we conclude that the optimal $b_j$ must satisfy the equations $\frac{a_{j}}{\sqrt{a_jb_j+c_j}}=\lambda$ whenever $b_j>0$ and the inequalities $\frac{a_{j}}{\sqrt{a_jb_j+c_j}}\le \lambda$ whenever $b_j=0$. Thus, the optimizer is just a vector $b_j=\max(\Lambda a_{j}-\frac{c_j}{a_j},0)$ with an unknown $\Lambda=\frac 1{\lambda^2}>0$ that should be found from the constraint $\sum_j b_j=1$. This is a one-variable equation for the root of a monotone function, so it can be easily solved in various ways.

Thus, we can optimize each row of $b$ with other rows fixed rather quickly. The natural idea is then just to cycle over rows optimizing each one in turn. It takes about 20 full cycles to get the lower bound and the value of the function within 1% range from each other on a random matrix (structured matrices seem to be even better) up to the size of 300 by 300 at least.

This is the description. The code (in Asymptote) is below.

srand(seconds());

int m=50, n=55;

real[][] a, b;
for(int i=0;i<m;++i)
{
    a[i]=new real[]; b[i]=new real[];
    for(int j=0; j<n; ++j)
    {
        a[i][j]=unitrand();
        b[i][j]=1/n;
    }
    //a[i][rand()%n]=2;
    a[i]*=unitrand();
}

real[] p, S;

for(int k=0;k<101;++k)
{
    for(int j=0;j<n;++j)
    {
        real s=0;
        for(int i=0;i<m;++i)
            s+=a[i][j]*b[i][j];
        S[j]=s;
        p[j]=sqrt(S[j]);
    }
    if(k%10==0)
    {
        write("*** Iteration "+string(k)+" ***");
        write(sum(map(sqrt,S))^2);
    }

    p/=sum(p);

    real s=0; 
    for(int i=0;i<m;++i)
    {
        real M=0; 
        for(int j=0;j<n;++j)
        {
            real h=a[i][j]/p[j];
            if(h>M)
                M=h;
        }
        s+=M;
    }
    if(k%10==0)
        write(s);
    //pause();

    for(int i=0;i<m;++i)
    {
        real[] A,V,C,B;
        for(int j=0;j<n;++j)
        {
            A[j]=a[i][j];
            V[j]=S[j]-a[i][j]*b[i][j];
            C[j]=V[j]/a[i][j];
        }
        real aa=(sum(C)+1)/sum(A);
        real da=1;
        while(da>1/10^10)
        {
            for(int j=0;j<n;++j)
            {
                B[j]=aa*A[j]-C[j];
                if(B[j]<0)
                {
                    A[j]=0;
                    B[j]=0;
                }
            }
            da=sum(B)-1; aa-=da/sum(A); 
        }
        for(int j=0;j<n;++j)
        {
            b[i][j]=B[j];
            S[j]=V[j]+a[i][j]*B[j];
        }
    }
}

write("************");

pause();

Best Answer

Related Solutions

[Math] Prime numbers in binary.

Optimize table layout

Related Question