Thursday, February 2, 2017

How to build a projection (Schlegel diagram) of a tesseract, and show a four-dimensional point inside it

This is an explanation about how to build a projection (Schlegel diagram) of a tesseract, and to show a point inside it (this is a mirror of my question at MSE). Regarding the methodology applied to visualize the 4D point, basically, if we want to show a point inside the tesseract, we need to project the tesseract first, and then project the desired point as well, following the same projection rules.

1. The definition of the tesseract is as follows (please see my post at MSE for the credits):

The tesseract is a four dimensional cube. It has 16 edge points $v=(a,b,c,d)$, with $a,b,c,d$ either equal to $+1$ or $-1$. Two points are connected, if their distance is $2$. Given a projection $P(x,y,z,w)=(x,y,z)$ from four dimensional space to three dimensional space, we can visualize the cube as an object in familiar space. The effect of a linear transformation like a rotation
$$
R(t)=\pmatrix{1&0&0&0\\0&1&0&0&\\0&0&\cos(t)&\sin(t)\\0&0&-\sin(t)&\cos(t)}
$$
in $4D$ space can be visualized in $3D$ by viewing the points $v(t) = P R(t) v$ in $\mathbb R^3$.

2. The definition of the projection (please see my post at MSE for the credits):

$$
P(x, y, z, w) = \frac{h}{h - w}(x, y, z).
$$

3. And finally, the definition of the distance between two four-dimensional points is calculated as follows (this is used to show the edges of the tesseract properly, making lines between the correct projected vertices):

$$d=\sqrt{(x_0-x'_0)^2+(x_1-x'_1)^2+(x_2-x'_2)^2+(x_3-x'_3)^2}$$

4. I have prepared a Python code snippet that creates the frames (jpg) of an animation of a tesseract including an internal point. In this case, the length of the edge is $1000$ so the distance between the vertices is not $2$, but $1000$. For the projection I have used a light source located at three times the length of the edge, this is $h=3000$. Finally, I have applied a rotation as defined above. The star is marking the location of the point $(\frac{3}{4} \cdot \frac{1000}{2}, \frac{3}{4} \cdot \frac{1000}{2}, \frac{3}{4} \cdot \frac{1000}{2}, \frac{3}{4} \cdot \frac{1000}{2})$ while we rotate the tesseract. Be aware that the position of the camera in the animation is lateral. The typical location of the camera is from above, which is the usual view of a square inside a square. But for visualization purposes (we want to see clearly the movement of the projection of the point due to the rotation of the tesseract) the camera in this case was located in a lateral position. Please use and modify it freely (for instance instead of one point it is possible to show a set of points and verify if there is symmetry, etc.):






from math import pi, sin , cos, sqrt
import matplotlib.pyplot as plt
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D

edge_length=1000
edge_half_length= int(edge_length/2)

lotuples=[]
list_of_loxt_lists=[]
list_of_loyt_lists=[]
list_of_lozt_lists=[]

rotation_accuracy=100
filled_once=False
for ratio in range(0,rotation_accuracy):
    angle= ((2*pi)*ratio)/rotation_accuracy
    loxt=[]
    loyt=[]
    lozt=[]
  
    #t=edge_half_length (positive)
    a=-edge_half_length
    b=edge_half_length
    ret0=-edge_half_length
    ret1=edge_half_length
    finala=a
    finalb=b
    finalret0=(ret0*cos(angle))+(ret1*sin(angle))
    finalret1=(ret0*(-sin(angle)))+(ret1*cos(angle))
  
    light_projection_factor = ((edge_length*3)/((edge_length*3)-(finalret1)))
    loxt.append(light_projection_factor*finala)
    loyt.append(light_projection_factor*finalb)
    lozt.append(light_projection_factor*finalret0)
    if filled_once==False:
        lotuples.append([a,b,ret0,ret1])
      
    a=-edge_half_length
    b=-edge_half_length
    ret0=-edge_half_length
    ret1=edge_half_length
    finala=a
    finalb=b
    finalret0=(ret0*cos(angle))+(ret1*sin(angle))
    finalret1=(ret0*(-sin(angle)))+(ret1*cos(angle))  

    light_projection_factor = ((edge_length*3)/((edge_length*3)-(finalret1)))
    loxt.append(light_projection_factor*finala)
    loyt.append(light_projection_factor*finalb)
    lozt.append(light_projection_factor*finalret0)
    if filled_once==False:
        lotuples.append([a,b,ret0,ret1])

    a=edge_half_length
    b=edge_half_length
    ret0=-edge_half_length
    ret1=edge_half_length
    finala=a
    finalb=b
    finalret0=(ret0*cos(angle))+(ret1*sin(angle))
    finalret1=(ret0*(-sin(angle)))+(ret1*cos(angle))  

    light_projection_factor = ((edge_length*3)/((edge_length*3)-(finalret1)))
    loxt.append(light_projection_factor*finala)
    loyt.append(light_projection_factor*finalb)
    lozt.append(light_projection_factor*finalret0)
    if filled_once==False:
        lotuples.append([a,b,ret0,ret1])

    a=edge_half_length
    b=-edge_half_length
    ret0=-edge_half_length
    ret1=edge_half_length
    finala=a
    finalb=b
    finalret0=(ret0*cos(angle))+(ret1*sin(angle))
    finalret1=(ret0*(-sin(angle)))+(ret1*cos(angle))  

    light_projection_factor = ((edge_length*3)/((edge_length*3)-(finalret1)))
    loxt.append(light_projection_factor*finala)
    loyt.append(light_projection_factor*finalb)
    lozt.append(light_projection_factor*finalret0)
    if filled_once==False:
        lotuples.append([a,b,ret0,ret1])

    a=edge_half_length
    b=edge_half_length
    ret0=edge_half_length
    ret1=edge_half_length
    finala=a
    finalb=b
    finalret0=(ret0*cos(angle))+(ret1*sin(angle))
    finalret1=(ret0*(-sin(angle)))+(ret1*cos(angle))  

    light_projection_factor = ((edge_length*3)/((edge_length*3)-(finalret1)))
    loxt.append(light_projection_factor*finala)
    loyt.append(light_projection_factor*finalb)
    lozt.append(light_projection_factor*finalret0)
    if filled_once==False:
        lotuples.append([a,b,ret0,ret1])

    a=edge_half_length
    b=-edge_half_length
    ret0=edge_half_length
    ret1=edge_half_length
    finala=a
    finalb=b
    finalret0=(ret0*cos(angle))+(ret1*sin(angle))
    finalret1=(ret0*(-sin(angle)))+(ret1*cos(angle))  

    light_projection_factor = ((edge_length*3)/((edge_length*3)-(finalret1)))
    loxt.append(light_projection_factor*finala)
    loyt.append(light_projection_factor*finalb)
    lozt.append(light_projection_factor*finalret0)
    if filled_once==False:
        lotuples.append([a,b,ret0,ret1])

    a=-edge_half_length
    b=edge_half_length
    ret0=edge_half_length
    ret1=edge_half_length
    finala=a
    finalb=b
    finalret0=(ret0*cos(angle))+(ret1*sin(angle))
    finalret1=(ret0*(-sin(angle)))+(ret1*cos(angle))  

    light_projection_factor = ((edge_length*3)/((edge_length*3)-(finalret1)))
    loxt.append(light_projection_factor*finala)
    loyt.append(light_projection_factor*finalb)
    lozt.append(light_projection_factor*finalret0)
    if filled_once==False:
        lotuples.append([a,b,ret0,ret1])

    a=-edge_half_length
    b=-edge_half_length
    ret0=edge_half_length
    ret1=edge_half_length
    finala=a
    finalb=b
    finalret0=(ret0*cos(angle))+(ret1*sin(angle))
    finalret1=(ret0*(-sin(angle)))+(ret1*cos(angle))  

    light_projection_factor = ((edge_length*3)/((edge_length*3)-(finalret1)))
    loxt.append(light_projection_factor*finala)
    loyt.append(light_projection_factor*finalb)
    lozt.append(light_projection_factor*finalret0)
    if filled_once==False:
        lotuples.append([a,b,ret0,ret1])

    #t=-edge_half_length (negative)
    a=-edge_half_length
    b=edge_half_length
    ret0=-edge_half_length
    ret1=-edge_half_length
    finala=a
    finalb=b
    finalret0=(ret0*cos(angle))+(ret1*sin(angle))
    finalret1=(ret0*(-sin(angle)))+(ret1*cos(angle))  

    light_projection_factor = ((edge_length*3)/((edge_length*3)-(finalret1)))
    loxt.append(light_projection_factor*finala)
    loyt.append(light_projection_factor*finalb)
    lozt.append(light_projection_factor*finalret0)
    if filled_once==False:
        lotuples.append([a,b,ret0,ret1])

    a=-edge_half_length
    b=-edge_half_length
    ret0=-edge_half_length
    ret1=-edge_half_length
    finala=a
    finalb=b
    finalret0=(ret0*cos(angle))+(ret1*sin(angle))
    finalret1=(ret0*(-sin(angle)))+(ret1*cos(angle))  

    light_projection_factor = ((edge_length*3)/((edge_length*3)-(finalret1)))
    loxt.append(light_projection_factor*finala)
    loyt.append(light_projection_factor*finalb)
    lozt.append(light_projection_factor*finalret0)
    if filled_once==False:
        lotuples.append([a,b,ret0,ret1])

    a=edge_half_length
    b=edge_half_length
    ret0=-edge_half_length
    ret1=-edge_half_length
    finala=a
    finalb=b
    finalret0=(ret0*cos(angle))+(ret1*sin(angle))
    finalret1=(ret0*(-sin(angle)))+(ret1*cos(angle))  

    light_projection_factor = ((edge_length*3)/((edge_length*3)-(finalret1)))
    loxt.append(light_projection_factor*finala)
    loyt.append(light_projection_factor*finalb)
    lozt.append(light_projection_factor*finalret0)
    if filled_once==False:
        lotuples.append([a,b,ret0,ret1])

    a=edge_half_length
    b=-edge_half_length
    ret0=-edge_half_length
    ret1=-edge_half_length
    finala=a
    finalb=b
    finalret0=(ret0*cos(angle))+(ret1*sin(angle))
    finalret1=(ret0*(-sin(angle)))+(ret1*cos(angle))  

    light_projection_factor = ((edge_length*3)/((edge_length*3)-(finalret1)))
    loxt.append(light_projection_factor*finala)
    loyt.append(light_projection_factor*finalb)
    lozt.append(light_projection_factor*finalret0)
    if filled_once==False:
        lotuples.append([a,b,ret0,ret1])

    a=edge_half_length
    b=edge_half_length
    ret0=edge_half_length
    ret1=-edge_half_length
    finala=a
    finalb=b
    finalret0=(ret0*cos(angle))+(ret1*sin(angle))
    finalret1=(ret0*(-sin(angle)))+(ret1*cos(angle))  

    light_projection_factor = ((edge_length*3)/((edge_length*3)-(finalret1)))
    loxt.append(light_projection_factor*finala)
    loyt.append(light_projection_factor*finalb)
    lozt.append(light_projection_factor*finalret0)
    if filled_once==False:
        lotuples.append([a,b,ret0,ret1])

    a=edge_half_length
    b=-edge_half_length
    ret0=edge_half_length
    ret1=-edge_half_length
    finala=a
    finalb=b
    finalret0=(ret0*cos(angle))+(ret1*sin(angle))
    finalret1=(ret0*(-sin(angle)))+(ret1*cos(angle))  

    light_projection_factor = ((edge_length*3)/((edge_length*3)-(finalret1)))
    loxt.append(light_projection_factor*finala)
    loyt.append(light_projection_factor*finalb)
    lozt.append(light_projection_factor*finalret0)
    if filled_once==False:
        lotuples.append([a,b,ret0,ret1])

    a=-edge_half_length
    b=edge_half_length
    ret0=edge_half_length
    ret1=-edge_half_length
    finala=a
    finalb=b
    finalret0=(ret0*cos(angle))+(ret1*sin(angle))
    finalret1=(ret0*(-sin(angle)))+(ret1*cos(angle))  

    light_projection_factor = ((edge_length*3)/((edge_length*3)-(finalret1)))
    loxt.append(light_projection_factor*finala)
    loyt.append(light_projection_factor*finalb)
    lozt.append(light_projection_factor*finalret0)
    if filled_once==False:
        lotuples.append([a,b,ret0,ret1])

    a=-edge_half_length
    b=-edge_half_length
    ret0=edge_half_length
    ret1=-edge_half_length
    finala=a
    finalb=b
    finalret0=(ret0*cos(angle))+(ret1*sin(angle))
    finalret1=(ret0*(-sin(angle)))+(ret1*cos(angle))  

    light_projection_factor = ((edge_length*3)/((edge_length*3)-(finalret1)))
    loxt.append(light_projection_factor*finala)
    loyt.append(light_projection_factor*finalb)
    lozt.append(light_projection_factor*finalret0)
    if filled_once==False:
        lotuples.append([a,b,ret0,ret1])
        filled_once=True

    list_of_loxt_lists.append(loxt)
    list_of_loyt_lists.append(loyt)
    list_of_lozt_lists.append(lozt)

list_of_loxi_lists=[]
list_of_loyi_lists=[]
list_of_lozi_lists=[]

list_of_loxi_lists_axis=[]
list_of_loyi_lists_axis=[]
list_of_lozi_lists_axis=[]

for ratio in range(0,rotation_accuracy):
    angle= ((2*pi)*ratio)/rotation_accuracy
    loxi=[]
    loyi=[]
    lozi=[]

    finala=int((3/4)*edge_half_length)
    finalb=int((3/4)*edge_half_length)
    finalret0=(int((3/4)*edge_half_length)*cos(angle))+(int((3/4)*edge_half_length)*sin(angle))
    finalret1=(int((3/4)*edge_half_length)*(-sin(angle)))+(int((3/4)*edge_half_length)*cos(angle))

    light_projection_factor = ((edge_length*3)/((edge_length*3)-finalret1))
    loxi.append(light_projection_factor*finala)
    loyi.append(light_projection_factor*finalb)
    lozi.append(light_projection_factor*finalret0)

    list_of_loxi_lists.append(loxi)
    list_of_loyi_lists.append(loyi)
    list_of_lozi_lists.append(lozi)
  
    # Show projection of refence axes BEGIN
    loxi=[]
    loyi=[]
    lozi=[]

    finala_axis=finala
    finalb_axis=0
    finalret0_axis=0
    finalret1_axis=0
  
    light_projection_factor = ((edge_length*3)/((edge_length*3)-finalret1_axis))
    loxi.append(light_projection_factor*finala_axis)
    loyi.append(light_projection_factor*finalb_axis)
    lozi.append(light_projection_factor*finalret0_axis)
  
    finala_axis=0
    finalb_axis=finalb
    finalret0_axis=0
    finalret1_axis=0
  
    light_projection_factor = ((edge_length*3)/((edge_length*3)-finalret1_axis))
    loxi.append(light_projection_factor*finala_axis)
    loyi.append(light_projection_factor*finalb_axis)
    lozi.append(light_projection_factor*finalret0_axis)
  
    finala_axis=0
    finalb_axis=0
    finalret0_axis=finalret0
    finalret1_axis=0
  
    light_projection_factor = ((edge_length*3)/((edge_length*3)-finalret1_axis))
    loxi.append(light_projection_factor*finala_axis)
    loyi.append(light_projection_factor*finalb_axis)
    lozi.append(light_projection_factor*finalret0_axis)
  
    finala_axis=0
    finalb_axis=0
    finalret0_axis=0
    finalret1_axis=finalret1
  
    light_projection_factor = ((edge_length*3)/((edge_length*3)-finalret1_axis))
    loxi.append(light_projection_factor*finala_axis)
    loyi.append(light_projection_factor*finalb_axis)
    lozi.append(light_projection_factor*finalret0_axis)
  
    list_of_loxi_lists_axis.append(loxi)
    list_of_loyi_lists_axis.append(loyi)
    list_of_lozi_lists_axis.append(lozi)
    # Show projection of refence axes END

for ratio in range(0,rotation_accuracy):  
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    ax.view_init(elev=17., azim=-152) 

    for i in range(0,len(lotuples)):
        for j in range(i+1,len(lotuples)):
            distance = int(sqrt(((lotuples[i][0]-lotuples[j][0])**2)+((lotuples[i][1]-lotuples[j][1])**2)+((lotuples[i][2]-lotuples[j][2])**2)+((lotuples[i][3]-lotuples[j][3])**2)))
            if distance<=edge_length:
                ax.plot([list_of_loxt_lists[ratio][i],list_of_loxt_lists[ratio][j]],[list_of_loyt_lists[ratio][i],list_of_loyt_lists[ratio][j]],[list_of_lozt_lists[ratio][i],list_of_lozt_lists[ratio][j]],"r")


        ax.plot([-edge_length],[edge_length],[-edge_length],"w")
        ax.plot([-edge_length],[-edge_length],[-edge_length],"w")
        ax.plot([edge_length],[edge_length],[-edge_length],"w")
        ax.plot([edge_length],[-edge_length],[-edge_length],"w")
        ax.plot([edge_length],[edge_length],[edge_length],"w")
        ax.plot([edge_length],[-edge_length],[edge_length],"w")
        ax.plot([-edge_length],[edge_length],[edge_length],"w")
        ax.plot([-edge_length],[-edge_length],[edge_length],"w")

    ax.plot([list_of_loxi_lists[ratio][0]], [list_of_loyi_lists[ratio][0]], [list_of_lozi_lists[ratio][0]], "r*")

    #projection of refence axes around the point  
    ax.plot([0,list_of_loxi_lists_axis[ratio][0]],[0,list_of_loyi_lists_axis[ratio][0]],[0,list_of_lozi_lists[ratio][0]],"b")
    ax.plot([0,list_of_loxi_lists_axis[ratio][1]],[0,list_of_loyi_lists_axis[ratio][1]],[0,list_of_lozi_lists_axis[ratio][1]],"b")
    ax.plot([0,list_of_loxi_lists_axis[ratio][2]],[0,list_of_loyi_lists_axis[ratio][2]],[0,list_of_lozi_lists_axis[ratio][2]],"b")
    ax.plot([0,list_of_loxi_lists_axis[ratio][3]],[0,list_of_loyi_lists_axis[ratio][3]],[0,list_of_lozi_lists_axis[ratio][3]],"b")
  
    ax.dist=8
    mpl.pyplot.savefig("tesseract_movie_"+str(ratio)+".png")
    print("End ratio "+str(ratio))


The animated gif was generated joining the jpg files with VirtualDub.

Wednesday, October 26, 2016

Testing MathJax in blogger!

Added MathJax libraries to the template of my blog (thanks to the explanation available on this blog), so this is a test: $$\lim_{x \to \infty}x^2$$ I will try to reformat my previous posts to update the contents!

Tuesday, September 6, 2016

A method to visualize the underlying patterns of fractal sequences

As Wikipedia says, a fractal sequence is one that contains itself as a proper subsequence. An example is $\{1, 1, 2, 1, 2, 3, 1, 2, 3, 4, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 6, ...\}$ If the first occurrence of each n is deleted, the remaining sequence is identical to the original. The process can be repeated indefinitely, so that actually, the original sequence contains not only one copy of itself, but rather, infinitely many.

Trying to find ways of visualizing the underlying patterns, I managed to convert a sequence to an elementary cellular automaton in which the first status (or step) of the automaton corresponds with the initial sequence, and the next status is a sequence generated by a replacement rule in which the closest element to the right of an element that is marked to be removed (following the definition of the fractal sequence) will invade (or "phagocyte") that position with its value (the invader value will increase its visual size). Applying this rule repeatedly, the evolution of the status of the automaton in time is visualized, and a fractal pattern arises.

Besides, if each element of the sequence of each status represents a binary bit (marked to be removed = $0$, not marked to be removed = $1$), it is possible to obtain a representation of the original sequence in terms of an equivalent sequence obtained by a binary manipulation of the status. Below there is an example of the algorithm and the questions are at the end:

1. For instance the following fractal sequence, (OEIS A000265), when the first appearance of each odd integer is removed -$1,3,5,7,9..$.-, the resulting sequence is the same sequence than the original one:

$$\{1,1,3,1,5,3,7,1,9,5,11,3,13,7,15,1,17,9,19,5,21...\}$$

2. The following image shows the first status on the top, which is the initial sequence marking the elements that will be removed (the first $1$, first $3$, first $5$, etc. in blue color; and the second status below the first status: the elements marked in the first status were "invaded" (or "phagocyted") by the right side elements that are not removed in the former step, so those elements increased their size. And then again we mark the elements that will be removed (the first $1$, first $3$, first $5$, etc.:


3. Repeating the process we will arrive to a fractal pattern like this one (e.g. in 5 steps):


4. Now if we check the vertical values on each element of the status, we can create an unique sequence equivalent to the original fractal sequence. Starting with the most left side column and assuming that the first status (top) represents $2^0$, the second status represents $2^1$, and the color represents the value of the multiplying term (*0 if blue or *1 if not blue), we obtain the sequence:

$$\{0,1,2,3,4,5,6,7,8,9...\}$$

For instance the first column is $0 \cdot 2^0 + 0 \cdot 2^1 + 0 \cdot 2^2 ... = 0$, the second column is $1 \cdot 2^0 + 0 \cdot  2^1 + 0 \cdot 2^2 ... = 1$, etc.

So, by this method, the equivalent sequence of $\{1,1,3,1,5,3,7,1,9,5,11,3,13,7,15,1,17,9,19,5,21...\}$ is $\{1,2,3,4,5,6,7,8,9,10,11,12...\}$.

The equivalent sequence vary depending on the original fractal sequence, these are some examples (click to expand the image):


After thinking a little bit more about the options, this is a possible way of showing the underlying patterns: for the same example as above, OEIS A000265, each initial number of the sequence (or first status of the automaton) is represented by a radius $1$ circle (yellow).

In the second step, the elements marked to be removed were "invaded" by the closest elements at their right side. The invader element grew. We will show that growth by adding a new circle with a radius that covers both the invaded element (represented by its former step circle) and the invader (also represented by its former step circle).

That new circle is e.g. shown in red color. When we repeat the algorithm, or in other words, we continue evolving the automaton shown in the question some more steps, finally the pattern starts to arise:


I have asked at MSE if there are other methods to show the underlying patterns of the fractal sequences.

Tuesday, July 12, 2016

Euler's "Lucky" constant

Based on my results in the previous post regarding a prime-representing function based on Bertrand's postulate I was able to produce a new constant related with Euler's "Lucky" numbers. Initially I did not find this constant in any bibliography or reference, and it seems a nice use of a Mills-like constant. So here are the results:

Euler's "Lucky" constant $E = 2.893392257682316134127494663$

This constant is such that $\lfloor{E^{2^n}}\rfloor-(\lfloor{E^{2^{n-1}}\rfloor})^2+\frac{\lvert n-(\frac{1}{2}) \rvert}{(\frac{1}{2})-n}$ for $n=[0..5]$  provides the sequence of Euler's "Lucky" prime numbers in growing order


Euler's lucky constant E is a Mill's like constant obtained by the encapsulation of Euler's "Lucky" numbers into integers E(n) of [N,(N+1)^2] intervals. This provides a representing function of Euler's "Lucky" primes for n=[0..5].

The calculation of the sequence of elements that will be the lower bounds of the [N,(N+1)^2] intervals is as follows:

E(0)=(2)
E(1)=3+1+(E1^2) = 8
E(2)=5+1+(E2^2) = 70
E(3)=11+1+(E3^2) = 4912
E(4)=17+1+(E4^2) = 24127762
E(5)=41+1+(E5^2) = 582148899128686

And finally Euler's "Lucky" constant (E) has the following value:

E=E(5)^(1/(2^5))=2.893392257682316134127494663





The above manipulation encapsulates the "Lucky" primes into a different upper interval, so the following calculation returns the "Lucky" primes back from the representing function:

Lucky(n+1)=floor(E^(2^n))-(floor(E^(2^(n-1))))^2+(abs(n-(1/2))/((1/2)-n))
for n=0..5

Thus:

Lucky(1)=floor(E^(2^0))-(floor(E^(2^(0-1))))^2+(abs(0-(1/2))/((1/2)-0))=2

Lucky(2)=floor(E^(2^1))-(floor(E^(2^(1-1))))^2+(abs(1-(1/2))/((1/2)-1))=3

Lucky(3)=floor(E^(2^2))-(floor(E^(2^(2-1))))^2+(abs(2-(1/2))/((1/2)-2))=5

Lucky(4)=floor(E^(2^3))-(floor(E^(2^(3-1))))^2+(abs(3-(1/2))/((1/2)-3))=11

Lucky(5)=floor(E^(2^4))-(floor(E^(2^(4-1))))^2+(abs(4-(1/2))/((1/2)-4))=17

Lucky(6)=floor(E^(2^5))-(floor(E^(2^(5-1))))^2+(abs(5-(1/2))/((1/2)-5))=41

Q.E.D.


The term (abs(n-(1/2))/((1/2)-n)) is just a correction to calculate properly the initial case Lucky(1) associated with n=0 and join it to the rest of cases. The value E(0) is an special case because there is not previous term. (abs(n-(1/2))/((1/2)-n))=+1 for n=0 and =-1 for the rest of cases.

Finally, the following polynomial in two variables provides the complete set Euler primes (A196230) for n=[0..5] and x = [1..(floor(E^(2^n))-(floor(E^(2^(n-1))))^2+(abs(n-(1/2))/((1/2)-n)))-1]:

x^2-x+[floor(E^(2^n))-(floor(E^(2^(n-1))))^2+(abs(n-(1/2))/((1/2)-n))] 

I think that this is a very nice application of a Mills-like constant, we can obtain the whole set of Euler's "Lucky" numbers just using a single expression.

Thursday, July 7, 2016

A prime-representing function based on Bertrand's postulate and Mill's constant

P.S. This is a backup of my question at MSE.

Wikipedia explains in number theory, Mills' constant is defined as:

"The smallest positive real number A such that the floor function of the double exponential function floor(A^(3^(n)) is a prime number, for all natural numbers n. This constant is named after William H. Mills who proved in 1947 the existence of A based on results of Guido Hoheisel and Albert Ingham on the prime gaps. Its value is unknown, but if the Riemann hypothesis is true, it is approximately 1.3063778838630806904686144926... (sequence A051021 in OEIS)."

In other words floor(A^(3^(n))) is said to be a prime-representing function, because f(x) is a prime number for all positive integral values of x.

The demonstration made by Mills in 1947 (pdf here) is very easy to follow (one must be careful about the already known typos of the paper to follow properly the explanation).

It is also known that floor(Q^k^n) always works if k is at least 3, so there are other constants not defined yet that will work just if k is equal to or more than 3. It is also kown that k=2 might not work because it depends on the Legendre's conjecture: is there always a prime between N2 and (N+1)^2. It it's thought to be extremely difficult.

Based on the original demonstration and using the Bertrand's postulate as key for some manipulations I was able to find a constant L that for floor(L^2^n) provides a sequence of integers (the difference is that they are not primes directly) such as there is an associated sequence of primes obtained from them. Here is how it works:

According to Bertrand's postulate for a given prime p_i it holds:

p_i < p_j < 2p_i for some existing prime p_j

Now adding  +p_i^2+1:

p_i+p_i^2+1 < p_j+p_i^2+1 < 2p_i+p_i^2+1

Which is also true for:

p_i^2 < p_j+p_i^2+1 < (p_i+1)^2

Calling e_j = p_j+p_i^2+1 it is possible to build a sequence of integers E_0, E_1,...E_n such as:

E_0=P_0=2

E_0^2 < E_1 < (E_0+1)^2

E_1^2 < E_2 < (E_1+1)^2

...

E_(n-1)^2 < E_n < (E_(n-1)+1)^2
...

The same conditions for the sequence shown in Mill's demonstration would hold for this expression, and in this case the exponent is 2 instead 3. It is possible because the sequence is not directly a sequence of primes, but a sequence of integers attached to primes by Bertrand's postulate according to the manipulation:

P_n = E_n-(E(n-1)^2-1)

An in the same fashion as Mill's constant:

E_n=floor(L^(2^(n)))

Where L is obtained after some manual calculations (up to n=11) as

L_11=1.71979778440410781908547548439826963...

The way of manually calculating the constant is as follows, in the same fashion that Mill's demonstration, the relationship between L_n (the constant calculated after knowing the nth element of the sequence E_n) and E_n is:

L_n=E_n^(1/(2^n))

These are the calculations for the first 5 elements in Excel (the accuracy of the decimals of L is not so good, but it works properly).

L slightly grows on every iteration but the growth is each time smaller, so it clearly tends to a fixed value when n tends to infinity.

(graph here)

For instance for n=1..5 floor(L^(2^(n)) is:

E_1=2

E_2=8 and P_1=E_2-(E_1)^2-1=8-2^2-1=3

E_3=76 and P_2=E_3-(E_2)^2-1=76-8^2-1=76-64-1=11

E_4=5856 and P_3=E_4-(E_3)^2-1=5856-76^2-1=5856-5776-1=79

E_5=34298594 and P_4=E_5-(E_4)^2-1=34298594-5856^2-1=5857

etc.

So in essence the prime-representing function would be:

floor(L^(2^(n))-E_(n-1) for n>=2

being the starting element of the sequence E_1=2 and L the value of L_n when n tends to infinity.

This is the Python code to calculate and test L:

def mb():
    from sympy import nextprime
    from gmpy2 import is_prime
    import matplotlib.pyplot as plt
    import matplotlib as mpl
    import decimal
   
    print("L constant test. Equivalent to Mill's constant applying Bertrand's postulate")
    print("----------------------------------------------------------------------------")
    print()
   
    # Decimal precision is required, change according to test_limit
    decimal.getcontext().prec = 2000
    # List of elements E_n generated by L^(2^(n))
    E_n = []
    # List of progressive calculations of L_n, the last one of the list is the most accurate
    # capable of calculate all the associated primes
    L_n=[]
    # List of primes obtained from L^(2^(n))-E_(n-1) associated to the Bertrand's postulate
    P_n=[]
    # depth of L: L will be able to calculate the primes L^(2^(n))-E_(n-1) for n=1 to test_limit-1
    test_limit=12
    for n in range(1,test_limit):
        if n==1:
            # E_1=2
            E_n.append(2)
        else:
            # E_n=(E_(n-1)^2)+1
            # Be aware that the Python list starts in index 0 and the last current element is index n-1:
            # that is why n-2 appears in the calculation below
            E_n.append((E_n[n-2]**(decimal.Decimal(2)))+P_n[n-2]+1)
        # Next prime greater than E_n: it will be in the interval [E_(n-1),2*E_(n-1)] (Bertrand's postulate)
        P_n.append(nextprime(E_n[n-1]))
        # Calculation of L_n
        L_n.append(E_n[n-1]**(decimal.Decimal(1)/(decimal.Decimal(2)**(decimal.Decimal(n)))))
   
    print("List of Elements of L^(2^(n)) for n = 1 to " + str(test_limit-1) + ":")
    mystr = ""
    for i in E_n:
        mystr=mystr+str(i)+"\t"
    print(mystr)
   
    print()
    print("List of Primes obtained from L constant: L^(2^(n))-E_(n-1)-1 for n = 1 to :" + str(test_limit-1) + ":")
    mystr = ""
    for i in P_n:
        mystr=mystr+str(i)+"\t"
    print(mystr)
   
    print()
    print("List of calculations of L_n (the most accurate one is the last one) for n = 1 to :" + str(test_limit-1) + ":")
    mystr = ""
    ax = plt.gca()
    ax.set_axis_bgcolor((0, 0, 0))
    figure = plt.gcf()
    figure.set_size_inches(18, 16)
    n=1
    for i in L_n:
        mystr=mystr+str(i)+"\t"
        plt.plot(n,i,"w*")
        n=n+1
    print(mystr)
   
    # Print the graph of the evolution of L_n
    # Clearly it tends to one specific value when n tend to infinity
    plt.show()
   
    #Testing the constant
    print()
    print("L Accuracy Test: using the most accurate value of L will calculate the primes associated")
    print("to the constant and will compare them to the original primes used to create the constant.")
    print("If they are the same ones, the constant is able to regenerate them and the theoretical background")
    print("applied would be correct.")
    #Using the most accurate value of L available
    L = L_n[len(L_n)-1]
    print()
    print("L most accurated value is:")
    print(L)
    print()
    for i in range(1,test_limit):
        print()
        tester = int(L**(decimal.Decimal(2)**(decimal.Decimal(i))))
        if i==1:
            tester_prime = 2
        else:
            tester_prime = decimal.Decimal(tester) - (E_n[i-2]*E_n[i-2]) - 1
        print("Current calculated E:\t" + str(tester) + " for n = " + str(i))
        print("Original value of E:\t" +str(E_n[i-1]) +  " for n = " + str(i))
        if tester == E_n[i-1]:
            print("Current calculated prime:\t" + str(tester_prime) +  " for n = " + str(i))
            if i==1:
                print("Original value of prime:\t2 for n = " + str(i))
            else:
                print("Original value of prime:\t" + str(P_n[i-2]) + " for n = " + str(i))
        else:
            # If we arrive to this point, then the constant and the theory would not be correct
            print("ERROR L does not generate the correct original E")
            return
    print()
    print("TEST FINISHED CORRECTLY: L generates properly the primes associated to the Bertrand's postulate")
mb()


Initially L=1.71979778... seems not have been defined before, but who knows. At least I was not able to find such constant Mill's constant-related papers on Internet. I find it interesting. I am trying to calculate a more accurate value of the constant (will try to update this post as soon as possible). The reality is that it is able to obtain only prime numbers if the precision is correct (in the same way as the Mill's constant).

The advantage of L versus A (the original Mill's constant) is:

1. The growing rate of the double exponential is lower due to the use of powers of 2 instead of powers of 3 or more, so for the same quantity of n values calculated L provides smaller primes than A in less time.

2. The use of L^(2^(n)) does not require the validation of Legendre's conjecture because it depend on Bertrand's postulate, so the theory would hold.

I have made a question at MSE about this calculations. According to the feedback, initially it seems that the manipulations I did for the theory are correct and could open an alternative way of calculating prime-representing functions with lower rate of growth. I am wrapping-up the ideas and preparing a paper. I will try to update this post as soon as I finish it.