First, we choose the value of $d$ and look for primes congruent to $1$ modulo $d$ in an interval $[A,B]$.

Note : you can plug in the value of $d$, and the bounds $A$ and $B$, and then use MAJ+Enter to run the cell (or simply click on "Run" or "Exécuter")

In [1]:
d = 5;
A = 400;
B = 500;
primelist = prime_range(A,B)
dadm = []
for p in primelist:
    if p%d == 1:
        dadm += [p]
dadm
Out[1]:
[401, 421, 431, 461, 491]

Once we have a suitable $p \equiv 1$ mod $d$, we choose $\alpha$ and put $q = p^\alpha$. Then we define the ring $R = \mathbf Z/ q \mathbf Z$ and compute once and for all the number of multiplicative units $\varphi(q)$.

In [2]:
p = 401;
alpha = 1;
q = p**alpha;
R = Integers(q);
phiq = euler_phi(q);

We find a generator $g$ of $(\mathbf Z/ q \mathbf Z)^{\times}$ and list all the units modulo $q$ as the powers of this generator.

In [3]:
gens = R.unit_gens();
g = gens[0]; # g is a generator of the cyclic group (Z/qZ)*
units = [g**i for i in range(phiq)];

By raising $g$ to the power $\varphi(q)/d$ we obtain $u$ : a generator of the unique subgroup of order $d$ of $(\mathbf Z/ q \mathbf Z)^{\times}$. Then we list the elements of the subgroup as the successive powers of $u$. $$ \left\{ x^d = 1 \right\} = \{ u^0, u^1, \dots, u^{d-1} \}$$

In [4]:
u = g**(phiq/d);
subgroupoforderd = [u**j for j in range(d)];

The vector $\textbf{m} = (m_1,m_2)$ below is the list of the exponents of the Laurent polynomials we put inside the exponentials. For instance, we will take $\textbf{m} =(1,-1)$ for Kloosterman sums, because we are summing $e\left( \frac{ax + b x^{-1}}{q} \right)$ (the powers of $x$ which appear are $1$ and $-1$).

In [5]:
m = [0,1,-1] #the 0 in the beginning is just so that m[i] is exactly what is denoted m_i in the article.
# we only deal with vectors of length 2 for this example.

Given a vector $\textbf{m} = (m_1, m_2)$ as defined above, the following function returns the list of (the real and imaginary parts of) the sums: $$ \sum_{\substack{x \in (\mathbf Z/ q \mathbf Z)^{\times} \\ x^d = 1}}^{} e\left( \frac{ax^{m_1} + b x^{m_2}}{q} \right)$$ for all the values of $a$ and $b$ inside $\mathbf Z/ q \mathbf Z$. Roughly speaking: $$\mathcal S(\textbf{m},q,d) = \left\{ \sum_{\substack{x \in (\mathbf Z/ q \mathbf Z)^{\times} \\ x^d = 1}}^{} e\left( \frac{ax^{m_1} + b x^{m_2}}{q} \right); \ a,b \in \mathbf Z/ q \mathbf Z \right\}$$ except that we return the list of the real parts and the list of the imaginary parts in order to use plt.scatter for the graphic.

In [6]:
def S(m,q,d):
    re = []
    im = []
    la = []
    lb = []
    for x in subgroupoforderd :
        la += [CDF(exp(2*I*pi*float(x**m[1])/q))] # CDF stands for complexdoublefield, it is used to turn the exact values
        #of our exp(2ipi...) to floats on which the computations will be faster.
        lb += [CDF(exp(2*I*pi*float(x**m[2])/q))]
    for a in range(q):
        for b in range(q):
            s = 0
            for i in range(d) :
                s += ((la[i])**a)*((lb[i])**b)
            re += [s.real()]
            im += [s.imag()]
    return re, im

We put the outcome of the function in two lists $X$ and $Y$. This next cell is where all the computations happen, so it might take some time to run.

In [7]:
[X,Y]= S(m,q,d)

We plot the set of sums $\mathcal S(\textbf{m},q,d)$ as well as the $d$-cusps hypocycloid (this part is only relevant when we are illustrating theorem A (b) in the case where $d$ is a prime number, otherwise the region of equidistribution may be different from the region delimited by a $d$-cusps hypocycloid). In other words, it is relevant to plot the hypocycloid only when $d$ is a prime number which divides neither $m_1$ nor $m_2$.

In [8]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl

mpl.rcParams['font.size'] = 8.
mpl.rcParams['font.family'] = 'serif' #these lines just set the font of the axis labels

plt.clf()
plt.gca().set_aspect('equal', adjustable='box')

# The next four lines plot the $d$-cusps hypocycloid.
T = np.linspace(0,2*np.pi,10000) 
habs = [(d-1)*np.cos(t) + np.cos((1-d)*t) for t in T]
hord = [(d-1)*np.sin(t) + np.sin((1-d)*t) for t in T]
plt.plot(habs,hord, color = 'red', linewidth=0.5)

plt.scatter(X,Y,s=1, color='blue', marker = '+', linewidths=1.0)
plt.show()
#plt.savefig('Kd'+str(d)+'p'+str(p)+'alpha'+str(alpha)+'.png',bbox_inches ='tight', dpi=1000) #dpi determines the output
# resolution.

In order to force the parameters $a$ and $b$ to live only in sufficiently large subgroups, we list the divisors of $\varphi(q)$ which are $> \sqrt q$. The following cell returns $\texttt{SGadm}$: the list of the "admissible sizes of subgroups".

In [9]:
div = divisors(phiq)
SGadm = []
for k in div:
    if k > sqrt(q):
        SGadm += [k]
SGadm
Out[9]:
[25, 40, 50, 80, 100, 200, 400]

Then, we choose two admissible sizes of subgroups $\texttt{sizehq1}$ and $\texttt{sizehq2}$ from the output of the previous cell, and build the corresponding subgroups $H_q^{(1)}$ and $H_q^{(2)}$ in which we will let $a$ and $b$ vary. For $H_q^{(1)}$ for instance, we raise the generator $g$ of $(\mathbf Z/ q \mathbf Z)^{\times}$ to the power $\varphi(q)/ \texttt{sizehq1}$: this gives a generator $\texttt{genhq1}$ of the unique subgroup of size $\texttt{sizehq1}$ in $(\mathbf Z/ q \mathbf Z)^{\times}$, then we list the elements of this subgroup as the successive powers of $\texttt{genhq1}$.

In [10]:
sizehq1 = 200
sizehq2 = 80
genhq1 = g**(phiq/sizehq1)
genhq2 = g**(phiq/sizehq2)
hq1 = [(genhq1)**i for i in range(sizehq1)]
hq2 = [(genhq2)**i for i in range(sizehq2)]

Finally, we almost copy-paste the function for ploting the sums of the set $\mathcal S(\textbf{m},q,d)$, except that we only let $a$ and $b$ vary in $H_q^{(1)}$ and $H_q^{(2)}$ instead of all $\mathbf Z/ q \mathbf Z$. The output of the function $\texttt{SSG}(\mathbf{m},q,d,\texttt{hq1},\texttt{hq1})$ is roughly the list of points: $$ \left\{ \sum_{\substack{x \in (\mathbf Z/ q \mathbf Z)^{\times} \\ x^d = 1}}^{} e\left( \frac{ax^{m_1} + b x^{m_2}}{q} \right); \ a \in H_q^{(1)}, \ b \in H_q^{(2)} \right\}$$ except that we return the list of the real parts and the list of the imaginary parts in order to use plt.scatter for the graphic.

In [11]:
def SSG(m,q,d,hq1,hq2):
    re = []
    im = []
    la = []
    lb = []
    for x in subgroupoforderd :
        la += [CDF(exp(2*I*pi*float(x**m[1])/q))] 
        lb += [CDF(exp(2*I*pi*float(x**m[2])/q))]
    for a in hq1:
        fa = float(a)
        for b in hq2:
            fb = float(b)
            s = 0
            for i in range(d) :
                s += ((la[i])**fa)*((lb[i])**fb)
            re += [s.real()]
            im += [s.imag()]
    return re, im
In [12]:
[V,W]= SSG(m,q,d,hq1,hq2)

Finally, we plot the sets $$ \left\{ \sum_{\substack{x \in (\mathbf Z/ q \mathbf Z)^{\times} \\ x^d = 1}}^{} e\left( \frac{ax^{m_1} + b x^{m_2}}{q} \right); \ a \in H_q^{(1)}, \ b \in H_q^{(2)} \right\}$$ and (optional) the $d$-cusps hypocycloid.

In [13]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl

mpl.rcParams['font.size'] = 8.
mpl.rcParams['font.family'] = 'serif' 

plt.clf()
plt.gca().set_aspect('equal', adjustable='box')

#The next four lines plot the $d$-cusps hypocycloid.
T = np.linspace(0,2*np.pi,10000) 
habs = [(d-1)*np.cos(t) + np.cos((1-d)*t) for t in T]
hord = [(d-1)*np.sin(t) + np.sin((1-d)*t) for t in T]
plt.plot(habs,hord, color = 'red', linewidth=0.5)

plt.scatter(V,W,s=1, color='blue', marker = '+', linewidths=1.0)
plt.show()

#plt.savefig('KSGd'+str(d)+'p'+str(p)+'alpha'+str(alpha)+'hq1_'+str(sizehq1)+'hq2_'+str(sizehq2)+'.png',bbox_inches ='tight', dpi=1000)
In [ ]: