%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        REPORT OF THE PROJECT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentclass[a4paper,11pt]{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage[english]{babel}
\usepackage{amsmath,amssymb,amsthm,amsopn}
\usepackage{mathrsfs}
\usepackage{graphicx}
%\usepackage{tikz}
%\usepackage{array}
%\usepackage[top=1cm,bottom=1cm]{geometry}
%\usepackage{listings}
%\usepackage{xcolor}
\usepackage{hyperref}
\hypersetup{
    colorlinks=true,
    linkcolor=blue,
    citecolor=red,
}

% fancy headers and footers

%\usepackage{fancyhdr}
%\pagestyle{fancy}
%\fancyhead[L]{BCPST 2 - <NAME>}
%\fancyhead[R]{Rappels d'analyse}
%\pagenumbering{gobble} % no page numbering

% Création des labels Théorème, Lemme, etc...

\newtheoremstyle{break}%
{}{}%
{\itshape}{}%
{\bfseries}{}%  % Note that final punctuation is omitted.
{\newline}{}

\newtheoremstyle{sc}%
{}{}%
{}{}%
{\scshape}{}%  % Note that final punctuation is omitted.
{\newline}{}

\theoremstyle{break}
\newtheorem{thm}{Theorem}[section]
\newtheorem{lm}[thm]{Lemma}
\newtheorem{prop}[thm]{Proposition}
\newtheorem{cor}[thm]{Corollary}

\theoremstyle{sc}
\newtheorem{exo}{Exercice}

\theoremstyle{definition}
\newtheorem{defi}[thm]{Definition}
\newtheorem{ex}[thm]{Example}

\theoremstyle{remark}
\newtheorem{rem}[thm]{Remark}

% Raccourcis pour les opérateurs mathématiques (les espaces avant-après sont modifiés pour mieux rentrer dans les codes mathématiques usuels)
\DeclareMathOperator{\Ker}{Ker}
\DeclareMathOperator{\Id}{Id}
\DeclareMathOperator{\Img}{Im}
\DeclareMathOperator{\Card}{Card}
\DeclareMathOperator{\Vect}{Vect}
\DeclareMathOperator{\Tr}{Tr}
\DeclareMathOperator{\Mod}{mod}
\DeclareMathOperator{\Ord}{Ord}
\DeclareMathOperator{\lcm}{lcm}


% Nouvelles commandes
\newcommand{\ps}[2]{\left\langle#1,#2\right\rangle}
\newcommand{\ent}[2]{[\![#1,#2]\!]}
\newcommand{\diff}{\mathop{}\!\mathrm{d}}
\newcommand{\ie}{\emph{i.e. }}

% opening
\title{Normal bases generation in C}
\author{Édouard \textsc{Rousseau}\\Supervised by Michaël \textsc{Quisquater}}



\begin{document}

\maketitle

\begin{abstract}
  This is the report of a C project about the generation of normal elements in
  finite fields. Consider a field extension $\mathbb{F}_{p^d}/\mathbb{F}_p$, we
  say that $\alpha\in\mathbb{F}_{p^d}$ is normal if $\left\{
  \alpha,\alpha^p,\alpha^{p^2},\cdots,\alpha^{p^{d-1}} \right\}$ is a basis of
  $\mathbb{F}_{p^d}$ over $\mathbb{F}_p$. We first give some theory to
  characterize normal elements, then we describe three algorithms to compute
  normal elements : a randomized algorithm, Lüneburg's algorithm, and
  Lenstra's algorithm. Finally, we give experimental results about our
  implementation. All this work is based on Gao's PhD thesis~\cite{Ga93}.

\end{abstract}

\tableofcontents

\clearpage

\section{Introduction}
\subsection{Normal bases}
With the rise of electronics and computer science, areas like
cryptography and coding theory have been largely studied. In both these
domains, finite fields often have a fundamental role. Normal bases can be
used to implement efficiently the finite fields arithmetic in the hardware,
consumming less power than other bases. But normal bases can also be used as a
theoretic tools to understand field extensions, and it is probably why normal bases were
studied since the $19^\textrm{th}$ century. Gauss, for example, used normal
bases to study when regular polygons can be drawn with ruler and compass
alone, a \emph{very} old problem. The notion of normal bases is not linked
with finite fields, if $\mathbb{K}$ is a field and $\mathbb{L}$ is a finite
Galois extension of $\mathbb{K}$ of Galois group $G$, a normal basis of
$\mathbb{L}$ over $\mathbb{K}$ is a basis of the form $\left\{
  \sigma(\alpha)\;|\;\sigma\in G
\right\}$ where $\alpha\in\mathbb{L}$. In other words, it is a basis composed
of an element $\alpha$ and all its conjugates. In the case of finite fields,
the definition that we give is indeed the same as this one, since, in this case,
$G$ is generated by the Frobenius automorphism. Given a finite Galois extension
$\mathbb{L}/\mathbb{K}$, normal bases can also be used to realize the Galois
correspondence of $\mathbb{L}/\mathbb{K}$. Theory tells us that there exists a
correspondence between intermediate extensions of $\mathbb{L}/\mathbb{K}$ and subgroups
of $G$, but it does not give an \emph{effective} way of realising that
correspondence. Normal bases are a way to solve that problem. In this report, we
focus on the problem of finding normal elements in finite fields. That is
why, from now on, we do not look at cases other then finite extensions
of finite fields (\ie extensions of type
$\mathbb{F}_{q^d}/\mathbb{F}_q$).
\subsection{Recalls and notations}
In all the text, we denote by $\mathbb{F}_n$ the field with $n$ elements. We
recall that $n$ must be a \emph{prime power} (\ie $n=p^d$ where $p$ is a prime
number and $d$ is a positive number). We have
$\mathbb{F}_{p^d}\cong\mathbb{F}_p[X]/(P)$, where $P$ is an irreducible
polynomial of degree $d$ in $\mathbb{F}_p[X]$, and this is the representation
that will be used in the C code. From now, we note
$X$ both for the indeterminate $X$ and for its image $\bar X$ in the
quotient $\mathbb{F}_p[X]/(P)$. $\mathbb{F}_{p^d}$ is a vector space over
$\mathbb{F}_p$, of dimension $d$, the basis $\left\{ 1, X, X^2, \dots, 
X^{d-1} \right\}$ is called the \emph{polynomial basis} in the following. We also recall
that $\mathbb{F}_{p^d}$ is a \emph{field extension} of $\mathbb{F}_p$ and the
characteristic of $\mathbb{F}_{p^d}$ is $p$. Lastly, we denote by $\sigma$ the
\emph{Frobenius map}, defined by $\sigma(x)=x^p$ for all
$x\in\mathbb{F}_{p^d}$. This map is a $\mathbb{F}_p$-automorphism of
$\mathbb{F}_{p^d}$ (\ie $\sigma$ is a field morphism, a bijection, and $\forall
y\in\mathbb{F}_p,\;\sigma(y)=y$), as a consequence, $\sigma$ is also a linear
map.

Note that we look only
at \emph{prime extensions}, (\ie extensions of type
$\mathbb{F}_{p^d}/\mathbb{F}_p$). This choice has been made because the theory
of normal elements in extensions of type $\mathbb{F}_{q^n}/\mathbb{F}_q$ is not
different (replacing $p$ by $q$ in the following demontrations would be
sufficient), but the implementation of the algorithms is simpler in the case of
prime extensions. Speaking about implementation, we must introduce Flint
(Fast Library for Number Theory), because \emph{every} function wrote in this
project uses Flint.

\subsection{Flint}

Flint~\cite{Ha10} is a C library maintained by <NAME>, and developed
by him and many others since 2007. In Flint, we use the type
\texttt{fq\_t}, where the elements of $F_{p^d}\cong\mathbb{F}_p[X]/(P)$ are
represented as polynomials of degree less than $d$. The underlying data
structure is \texttt{fmpz\_poly\_t}, that is a \texttt{struct} containing
\begin{enumerate}
  \item a pointer to arbitrary large integers \texttt{fmpz}, representing the
    coefficients of the polynomial;
  \item the number of memory allocations for the pointer of type
    \texttt{slong}: Flint's own \texttt{unsigned long};
  \item the degree of the polynomial, a \texttt{slong} too.
\end{enumerate}
In Flint, it is possible to work with finite fields of arbitrary large
order and with polynomials or matrices over these fields. The basic functions
are already implemented, such as gcd or derivative for polynomials, or reduced
row echelon form for matrices. Using Flint save a lot of time, but it also has
its limits. It is possible to deal with arbitrary \emph{prime}
extensions. For extension of type $\mathbb{F}_{q^n}/\mathbb{F}_q$, it is still
possible, using polynomial arithmetic over $\mathbb{F}_q$. But using
polynomials
over $\mathbb{F}_{q^n}$ would have required polynomial in two
indeterminates, and everything would be more complicated. That is the reason we
look only at prime extensions.

All these recalls being made, and the true hero
(Flint) of our functions being introduced, we can begin our journey.

\section{Basics on normal bases}

In all this section, we set $p$ a prime number and $d$ a positive number. We
work with the extension $\mathbb{F}_{p^d}/\mathbb{F}_p$ and the Frobenius
morphism $\sigma$ defined above. We are now able to give a formal definition
of a normal element and a normal basis.

\begin{defi}[normal element, normal basis]
  Let $\alpha\in\mathbb{F}_{p^d}$, we say that $\alpha$ is a \emph{normal
  element} if $\left\{ \alpha, \sigma(\alpha),\dots,\sigma^{d-1}(\alpha)
  \right\}=\left\{ \alpha, \alpha^{p}, \dots, \alpha^{p^{d-1}} \right\}$ is a
  basis of $\mathbb{F}_{p^d}$. Such a basis is called a \emph{normal basis}. 
\end{defi}

In order to recognize a normal element $\alpha$, we can compute the dimension of the
linear span of $\left\{ \sigma^i(\alpha) \;|\; 0\leq i < k \right\}$. A way of
doing that is to construct the matrix $M$ which columns are the coordinates of the
elements $\sigma^i(\alpha)$ in the polynomial basis, and to check if $M$ is
non-singular. This can be done using Gauss algorithm. This method is not
efficient so we work on other characterizations of normal elements. We need
two more definitions to be able to state our results.

\begin{defi}[trace function]
  The \emph{trace function}
  $\Tr_{p^d|p}:\mathbb{F}_{p^d}\to\mathbb{F}_p$ of the extension
$\mathbb{F}_{p^d}/\mathbb{F}_p$ is the function defined by
\[
  \Tr_{p^d|p}(\alpha) = \sum_{i=0}^{d-1}\alpha^{p^i}.
\]
It takes its values in $\mathbb{F}_p$ since $(\Tr_{p^d|p}(\alpha))^p
= \Tr_{p^d|p}(\alpha)$. 
\end{defi}

The trace is a tool that permit to know when a family of elements $\alpha_1,
\dots, \alpha_d$ forms a basis of $\mathbb{F}_{p^d}$ over
$\mathbb{F}_p$, using the discrinant of these elements.

\begin{defi}[discriminant]
  Let $\alpha_1, \dots, \alpha_d$ be elements in $\mathbb{F}_{p^d}$, the
  \emph{discriminant} $\Delta(\alpha_1, \dots, \alpha_d)$ of these elements
  is the determinant
  \[
    \Delta(\alpha_1, \dots, \alpha_d)=\det
    \begin{pmatrix}
      \Tr(\alpha_1\alpha_1) & \Tr(\alpha_1\alpha_2) & \dots &
      \Tr(\alpha_1\alpha_d) \\
      \Tr(\alpha_2\alpha_1) & \Tr(\alpha_2\alpha_2) & \dots &
      \Tr(\alpha_2\alpha_d) \\
      \vdots & \vdots & & \vdots \\
      \Tr(\alpha_d\alpha_1) & \Tr(\alpha_d\alpha_2) & \dots &
      \Tr(\alpha_d\alpha_d)
    \end{pmatrix}
  \]
  where $\Tr$ is the same trace as before. We ommit the indices because
  the extension is always the same.
\end{defi}

With this last tool, we can state our first theorem.

\begin{thm}[Theorem 2.2.1, \cite{Ga93}]
  For any $n$ elements $\alpha_1, \dots, \alpha_d$ in
  $\mathbb{F}_{p^d}$, they form a basis of $\mathbb{F}_{p^d}$ over
  $\mathbb{F}_p$ if and only if $\Delta(\alpha_1, \dots,
  \alpha_d)\neq0$.
\end{thm}

\begin{proof}
  First assume that $\alpha_1, \dots, \alpha_d$ form a basis of
  $\mathbb{F}_{p^d}$ over $\mathbb{F}_p$. We prove that
  $\Delta(\alpha_1, \dots, \alpha_d)\neq0$ by showing that the row vectors
  $L_1, \dots, L_d$ of
  the matrix in the definition of $\Delta(\alpha_1, \dots, \alpha_d)$ are
  linearly independent over $\mathbb{F}_p$, and thus the matrix is
  nonsingular. Suppose that there exists $c_1, \dots,
  c_d\in\mathbb{F}_p$ with $\sum_i c_iL_i = 0$, it means
  \[
    c_1\Tr(\alpha_1\alpha_j) + \dots + c_n\Tr(\alpha_d\alpha_j) =
    0\textrm{ for }1\leq j\leq d.
  \]
  Then with $\beta=c_1\alpha_1 + \dots + c_d\alpha_d$, by linearity of the
  trace function, we get $\Tr(\beta\alpha_j) = 0$ for $1\leq j\leq d$. Since
  $\alpha_1, \dots, \alpha_d$ is a basis of $\mathbb{F}_{p^d}$, it follows
  that $\Tr(\beta\alpha) = 0$ for all $\alpha\in\mathbb{F}_{p^d}$. If
  $\beta\neq0$, it means that $\Tr(\gamma)=0$ for all
  $\gamma\in\mathbb{F}_{p^d}$ : we can see that this is impossible by thinking
  of $\Tr(\gamma)$ as the trace of the multiplication-by-$\gamma$ linear map. So we have $\beta=0$, and then
  $c_1\alpha_1 + \dots + c_d\alpha_d = 0$ implies $c_1=\dots=c_d=0$.

  Conversely, assume that $\Delta(\alpha_1, \dots, \alpha_d)\neq0$ and
  $c_1\alpha_1+\dots+c_d\alpha_d=0$ for some $c_1, \dots,
  c_d\in\mathbb{F}_p$. By multiplying by $\alpha_j$, we have
  \[
    c_1\alpha_1\alpha_j+\dots+c_d\alpha_d\alpha_j\textrm{ for }1\leq j \leq d,
  \]
  and by applying the trace function, we get
  \[
    c_1\Tr(\alpha_1\alpha_j)+\dots+c_n\Tr(\alpha_d\alpha_j)\textrm{ for }1\leq
    j\leq d.
  \]
  But this is the same as $\sum c_iL_i=0$ where the $L_i$ are the rows of the
  matrix in $\Delta(\alpha_1, \dots, \alpha_d)$. Since this matrix is
  nonsingular by hypothesis, its rows are linearly independent and it follow
  that $c_1=\dots=c_d$. Therefore $\alpha_1, \dots, \alpha_d$ are linearly
  independant over $\mathbb{F}_p$, and form a family of $d$ vectors in a
  $d$-dimensionnal vector space, so $\alpha_1, \dots, \alpha_d$ form  basis of
  $\mathbb{F}_{p^d}$ over $\mathbb{F}_p$.
\end{proof}

This characterization is interesting in itself, but the matrix defined in
$\Delta(\alpha_1, \dots, \alpha_d)$ is quite complicated. Fortunately,
we can use a much simpler matrix.

\begin{cor}
  \label{circulantCor}
 The elements $\alpha_1, \dots, \alpha_d$ form a basis of
 $\mathbb{F}_{p^d}$ over $\mathbb{F}_p$ if and only if the matrix $A$ is
 nonsigular, where
 \[
   A = 
   \begin{pmatrix}
     \alpha_1 & \alpha_2 & \dots & \alpha_d \\
     \alpha_1^p & \alpha_2^p & \dots & \alpha_d^p \\
     \vdots & \vdots & & \vdots \\
     \alpha_1^{p^{d-1}} & \alpha_2^{p^{d-1}} & \dots &
     \alpha_d^{p^{d-1}}.
   \end{pmatrix}
 \]
\end{cor}
\begin{proof}
  We have $\Delta(\alpha_1, \dots, \alpha_d)=\det(A^tA)=(\det A)^2$.   
\end{proof}

We now have a simpler matrix to study, but it is still a matrix, so the nature
of the problem is the same as before. Next lemma allow us to change the nature
of the problem, by working with polynomials.

\begin{lm}
  \label{henselLm}
  For any $n$ elements $a_0, a_1, \dots, a_{d-1}\in\mathbb{F}_{p^d}$, the
  $d\times d$ circulant matrix
  \[
    c[a_0, a_1, \dots, a_{d-1}] =
    \begin{pmatrix}
      a_0 & a_1 & a_2 & \dots & a_{d-1} \\
      a_{d-1} & a_0 & a_1 & \dots & a_{d-2} \\
      a_{d-2} & a_{d-1} & a_0 & \dots & a_{d-3} \\
      \vdots & \vdots & \vdots & & \vdots \\
      a_1 & a_2 & a_3 & \dots & a_0
    \end{pmatrix}
  \]
  is nonsigular if and only if the polynomial $\sum a_iX^i$ is relatively
  prime to $X^d-1$.
\end{lm}
\begin{proof}
  Let $A$ be the following $d\times d$ permutation matrix 
  \[
    \begin{pmatrix}
      0 & 1 & 0 & \dots & 0 & 0 \\
      0 & 0 & 1 & \dots & 0 & 0 \\
      \vdots & \vdots & \vdots & & \vdots & \vdots \\
      0 & 0 & 0 & \dots & 0 & 1 \\
      1 & 0 & 0 & \dots & 0 & 0
    \end{pmatrix}.
  \]
  Then we see that $c[a_0, a_1, \dots, a_{d-1}] =
  \sum_{i=0}^{d-1}a_iA^i=f(A)$, where $f(X)=\sum_{i=0}^{d-1}a_iX^i$. On the
  first line of $A^i$, there is only one nonzero value, it is a $1$ in
  position $i+1$. Hence the matrices $A^0, A, \dots, A^{d-1}$, are linearly
  independant, and since $A^d=I_d$, we know that the minimal polynomial
  of $A$ is $X^d-1$. Assume that $f(X)$ is relatively prime to $X^d-1$. Then
  there are polynomials $a(X), b(X)$ such that
  \[ 
  a(X)f(X) + b(X)(X^d-1)=1,
\]
and so
\[
  a(A)f(A) = I_d,
\]
as $A^d-I_d=0$. This implies that $f(A)$ is invertible and so
nonsigular. Now assume that $f(X)$ and $X^d-1$ are not relatively prime, and
note $d(X)\neq1$ their gcd. Let $f(X)=f_1(X)d(X)$ and $X^d-1=h(X)d(X)$.
Since $\deg h < n$, we have $h(A)\neq0$. But $h(A)d(A)=0$, so we see that
$d(A)$ is singular. Therefore $f(A)=f_1(A)d(A)$ is also singular. We have shown
that $f(A)$ is singular if and only if $f(X)$ is relatively prime to $X^d-1$.
\end{proof}

We are now able to state the last result oh this section.

\begin{thm}[Hensel, 1888]
  \label{hensel}
  Let $\alpha\in\mathbb{F}_{p^d}$, $\alpha$ is a normal element of the
  extension $\mathbb{F}_{p^d}/\mathbb{F}_p$ if and only if the polynomial
  $\alpha^{p^{d-1}}X^{d-1}+\dots+\alpha^pX+\alpha\in\mathbb{F}_{p^d}[X]$ 
  is relatively prime to $X^d-1$.
\end{thm}
\begin{proof}
  $\alpha$ is a normal element, if, by definition, the elements $\alpha,
  \alpha^p, \dots, \alpha^{p^{d-1}}$ form a basis of
  $\mathbb{F}_{p^d}$ over $\mathbb{F}_p$. By
  Corollary~\ref{circulantCor}, this is
  true if and only if the matrix
  \[
    A =
    \begin{pmatrix}
      \alpha & \alpha^p & \alpha^{p^2} & \dots & \alpha^{p^{d-1}} \\
      \alpha^p & \alpha^{p^2} & \alpha^{p^3} & \dots & \alpha \\
      \vdots & \vdots & \vdots & & \vdots \\
      \alpha^{p^{d-1}} & \alpha & \alpha^p & \dots & \alpha^{p^{d-2}}
    \end{pmatrix}
  \]
  is nonsingular. But, reversing the order of the rows in $A$, from the second
  row to the last, we get the matrix $c[\alpha, \alpha^p, \dots,
  \alpha^{p^{d-1}}]$, that is nonsingular if and only if $A$ if
  nonsingular. By Lemma~\ref{henselLm}, $c[\alpha, \alpha^p, \dots,
  \alpha^{d-1}]$ is nonsingular if and only if $X^d-1$ and
  $\alpha^{p^{d-1}}X^{d-1}+\dots+\alpha^pX+\alpha$ are relatively prime.
  This proves the theorem.
\end{proof}
We now have a new characterization of normal elements. We can decide if an
element is normal or not by computing a gcd. It is more efficient that the
naive matrix method because we have efficient algorithms to compute the gcd.
The function \texttt{is\_normal} is exactly the implementation of this
method. The code of this function, and of all the others, is available at
\url{github.com/erou/normalBases}.

\section{Computation of normal bases}
Before trying to compute normal elements, a natural question is to wander if
normal elements always exist in extensions of type
$\mathbb{F}_{p^{d}}/\mathbb{F}_p$. We do not prove this result here (it is done
in~\cite{Ga93}), but the answer is yes. Now that this natural fear has been
eliminated, we give some ways to compute normal elements. We use the same
notations as in the previous section.

\subsection{Randomized algorithm}
First, we give a randomized algorithm. Such algorithms are often based on the
same strategy.
\begin{enumerate}
  \item Take a random element, following a certain protocol.
  \item Check if this element verify the wanted property.
\end{enumerate}
Our algorithm also uses this strategy, that is why our first function
\texttt{is\_normal} is very important. But the protocol is also essential.
For example, if we just take an element completly at random in
$\mathbb{F}_{p^d}$, our algorithm will not be very efficient. The following
theorem helps us defining a better protocol.
\begin{thm}
  Let $f(X)$ be an irreducible polynomial of degree $d$ over
  $\mathbb{F}_p$ and $\alpha$ a root of $f(X)$. Let
  \[
    g(X)=\cfrac{f(X)}{(X-\alpha)f'(\alpha)}.
  \]
  Then there are at least $p -d(d-1)$ elements $u$ in $\mathbb{F}_p$ such that
  $g(u)$ is a normal element of $\mathbb{F}_{p^d}$ over $\mathbb{F}_p$.
\end{thm}
\begin{proof}
  Let $\sigma_i=\sigma^i$ be the automorphism
  $\theta\to\theta^{p^i},\theta\in\mathbb{F}_{p^d}$, for $0\leq i <d$. Then
  $\alpha_i=\sigma_i(\alpha)$ is also a root of $f(X)$, for $0\leq i <d$. The
  automorphism $\sigma_i:\mathbb{F}_{p^d}\to\mathbb{F}_{p^d}$ can be extended
  into a ring morphism
  $\sigma_i:\mathbb{F}_{p^d}[X]\to\mathbb{F}_{p^d}[X]$ by setting
  $\sigma_i(X)=X$. We get
  \[
    g_i(X):=\sigma_i(g(X)) = \cfrac{f(X)}{(X-\alpha_i)f'(\alpha_i)},
  \]
  and we note that $\sigma_i\sigma_i(g(X)) = \sigma_{i+j}(g(X))$. Then
  $g_i(X)$ is a polynomial in $\mathbb{F}_{p^d}[X]$ having $\alpha_k$ as a root
  for $k\neq i$ and $g_i(\alpha_i)=1$. Hence, for $i\neq k$, we have that every
  root of $f(X)$ is also a root of $g_i(X)g_k(X)$, so
  \begin{equation}
    g_i(X)g_k(X) \equiv 0 \;(\Mod f(X)), \textrm{ for }i\neq k.
    \label{3.1}
  \end{equation}
  Note that 
  \begin{equation}
    g_1(X)+g_2(X)+\dots+g_d(X)-1=0,
    \label{3.2}
  \end{equation}
  since the left side is a polynomial of degree at most $d-1$ (all the
  $g_i$
  are of degree $d-1$) and has $\alpha_0,
  \alpha_1, \dots, \alpha_{d-1}$ as roots. Multiplying (\ref{3.2}) by $g_i(X)$
  and using (\ref{3.1}) to simplify, we have
  \begin{equation}
    g_i(X)^2 \equiv g_i(X)\;(\Mod f(x)).
    \label{3.3}
  \end{equation}
  We next set the matrix
  \[
    D=(\sigma_{i+j}(g(X)))_{0\leq i,j <0}
  \]
  and we study its determinant, $\Delta(X)$. From the equations
  (\ref{3.1}), (\ref{3.2}) and (\ref{3.3}), we wee that the entries of $D^tD$
  molulo $f(X)$ are all $0$, except on the main diagonal where they are all
  $1$. It follows that
  \[
    \Delta(X)^2 = \det(D^tD) \equiv 1\;(\Mod f(x)).
  \]
  This proves that $\Delta(X)$ is a nonzero polynomial of degree at most
  $d(d-1)$, since it is a sum of $d$ products of polynomials of degree
  $d-1$. Therefore $\Delta(X)$ has at most $d(d-1)$ roots in
  $\mathbb{F}_p$. The result follows from Corollary~\ref{circulantCor}, since the matrix
  $D(u)$ is exactly the one defined in the corollary applied to the elements
  $g(u), \sigma(g(u)),\dots,\sigma^{d-1}(g(u))$.
\end{proof}
Now we have our protocol. We take $u\in\mathbb{F}_p$ at random, then we check
if $g(u)$ is normal. If $p$ is large enough, for example $p>2d(d-1)$, then
$g(u)$ is normal with probability at least $1/2$. We will not discuss it, but
the entire computation takes $O((n+\log q)(n\log q)^2)$ bit operations.
A reference for this result is given in~\cite{Ga93}. As always, we are not
\emph{completely} satisfied with randomized algorithm, so we give two
deterministic algorithms.

\subsection{Deterministic algorithms}

First, we define the $\sigma$-Order of an element
$\theta\in\mathbb{F}_{p^d}$, both Lenstra's and Lüneburg's algorithms use it.

\subsubsection{The $\sigma$-Order polynomials}

\begin{defi}[$\sigma$-Order]
  Let $\theta\in\mathbb{F}_{p^d}$ be an arbitrary element. Let $k$ be the least
  positive integer such that $\sigma^k(\theta)=\theta^{p^k}$ belongs to
  $\mathbb{F}_p$-linear span of $\left\{ \sigma^i\theta\;|\;0\leq i < k
  \right\}$. If $\sigma^k\theta=\sum_{i=0}^{k-1}c_i\sigma^i\theta$ for that
  $k$, then the \emph{$\sigma$-Order} of $\theta$ is the polynomial
  \[
    \Ord_\theta(X)=X^k-\sum_{i=0}^{k-1}c_iX^i.
  \]
\end{defi}
The $\sigma$-Order polynomials are widely used in our functions. We compute
them using Gauss algorithm and the row reduced echelon form. To know if
$\sigma^j(\theta)$ is in the linear span of $\left\{
\sigma^i(\theta)\;|\; 0\leq i <j \right\}$, we construct the matrix $M$ whose columns are the
vectors $\theta, \sigma(\theta), \dots, \sigma^{j}(\theta)$, we compute
the rank using Gauss algorithm. If the rank is $j+1$ then $\theta,
\sigma(\theta), \dots, \sigma^j(\theta)$ are linearly independant. In order to
compute $k$, we do this operation for $j=1$ and while the vectors are linearly
independant, we increase the value of $j$ by $1$. When we get that the
vectors are dependant, we use the reduced row echelon form to get a
$d\times(k+1)$ matrix of type 
\[
\begin{pmatrix}
  1 & 0 & 0 & \dots & 0 & c_0 \\
  0 & 1 & 0 & \dots & 0 & c_1 \\
  0 & 0 & 1 & \dots & 0 & c_2 \\
  \vdots & \vdots & & \ddots & & \vdots \\
  0 & 0 & 0 & \dots & 1 & c_{k-1} \\
  0 & 0 & 0 & \dots & 0 & 0 \\
  \vdots & \vdots & \vdots & & \vdots & \vdots \\
  0 & 0 & 0 & \dots & 0 & 0 \\
\end{pmatrix}
\]
where the coefficients $c_i$ are in the last column. This strategy is
implemented in the function \texttt{sigma\_order}.
  
  
We also see that an element
$\alpha\in\mathbb{F}_{p^d}$ is normal is and only if
$\Ord_\alpha(X)=X^d-1$. This characterization is the same as the
definition, thus it is not very interesting. The following property is used in
further demonstrations.

\begin{prop}
  \label{orderDiv}
  Let $P\in\mathbb{F}_p[X]$, if $P(\sigma)\zeta=0$, then $P$ is divisible
  by $\Ord_{\zeta}$.
\end{prop}
\begin{proof}
  Let $I_\zeta$ be the set defined by $I_\zeta=\left\{P\in\mathbb{F}_p[X]\;|\; P(\sigma)\zeta
  = 0
\right\}$. This is an ideal of $\mathbb{F}_p[X]$, so $I_\zeta$ is principal. The
polynomial 
$\Ord_\zeta$ belongs to $I_\zeta$ and has minimum degree in $I_\zeta$ by definition.
Hence $I_\zeta=(\Ord_\zeta)$.
\end{proof}
We always have $(X^d-1)(\sigma)\theta=0$, so for any
$\theta\in\mathbb{F}_{p^d}$, we have that $X^d-1$ is divisible by
$\Ord_\theta$. The result below is also used both in Lenstra's and Lüneburg's
algorithm.
\begin{prop}
  \label{orderLcm}
  Let $\alpha$ and $\beta$ be two elements in $\mathbb{F}_{p^d}$. Then we
  have \[
    \Ord_{\alpha+\beta}\;|\;\lcm(\Ord_\alpha,\Ord_\beta).
  \]
\end{prop}
\begin{proof}
 Let $P=\lcm(\Ord_\alpha,\Ord_\beta)$. By linearity, we have
  \[
  P(\sigma)(\alpha+\beta)=P(\sigma)\alpha+P(\sigma)\beta=0+0=0.
\]
It follows from Proposition~\ref{orderDiv} that $\Ord_{\alpha+\beta}$ divides
$P$.
\end{proof}
\begin{cor}
  \label{orderMul}
  Let $\alpha$ and $\beta$ be two elements in $\mathbb{F}_{p^d}$ such that
  $\Ord_\alpha$ and $\Ord_\beta$ are relatively prime. Then we have
  $\Ord_{\alpha+\beta}=\Ord_\alpha\Ord_\beta$.
\end{cor}
\begin{proof}
  Since $\Ord_\alpha$ and $\Ord_\beta$ are relatively prime, we have
  $\lcm(\Ord_\alpha, \Ord_\beta)=\Ord_\alpha\Ord_\beta$. It follows from
  Proposition~\ref{orderLcm} that $\Ord_{\alpha+\beta}$ divides
  $\Ord_\alpha\Ord_\beta$. Now, note that for any
  $\theta\in\mathbb{F}_{p^d}$, $\Ord_{-\theta}=\Ord_\theta$, then we also have 
  \[
    \Ord_{\alpha}\;|\;\lcm(\Ord_{\alpha+\beta}, \Ord_\beta)
  \]
  \[
    \Ord_{\beta}\;|\;\lcm(\Ord_{\alpha+\beta}, \Ord_\alpha).
  \]
  Since $\Ord_\alpha$ and $\Ord_\beta$ are relatively prime, it follows
  that
  \[
    \Ord_{\alpha}\;|\;\Ord_{\alpha+\beta}
  \]
  \[
    \Ord_{\beta}\;|\;\Ord_{\alpha+\beta}.
  \]
  Therefore, $\Ord_\alpha\Ord_\beta$ divides $\Ord_{\alpha+\beta}$. Both
  polynomial are monic so we have the wanted equality.
\end{proof}

\subsubsection{Lenstra's algorithm}

Lenstra's algorithm is based on linear algebra, so our implementation is based
on the matrices module \texttt{fq\_mat} of Flint. Before describing the
algorithm, we need two lemmas. 
\begin{lm}
  \label{lmExist}
  Let $\theta\in\mathbb{F}_{p^d}$ with $\Ord_\theta(X)\neq X^d-1$. Let
  $g(X)=(X^d-1)/\Ord_\theta(X)$. Then there exists
  $\beta\in\mathbb{F}_{p^d}$ such that
  \begin{equation}
    g(\sigma)\beta=\theta.
    \label{3.4}
  \end{equation}
\end{lm}
\begin{proof}
  Let $\gamma$ be a normal element of $\mathbb{F}_{p^d}$ over
  $\mathbb{F}_p$. Then, by definition of being a normal element, there
  exists a polynomial $f(X)\in\mathbb{F}_{p}[X]$ such that
  $f(\sigma)\gamma=\theta$. Since $\Ord_\theta(\sigma)\theta=0$, we have
  $(\Ord_\theta(\sigma)f(\sigma))\gamma=0$. So, by
  Proposition~\ref{orderDiv}, $\Ord_\theta(X)f(X)$ is
  divisible by $X^d-1$. We have $g(X)=(X^d-1)/\Ord_\theta(X)$ and $X^d-1|\Ord_\theta(X)f(X)$, so
$g(X)|f(X)$. Let $f(X)=g(X)h(X)$. Then
\[
  g(\sigma)(h(\sigma)\gamma)=\theta.
\]
This proves that $\beta=h(\sigma)\gamma$ is a solution of (\ref{3.4}).
\end{proof}
\begin{lm}
  \label{lmDegree}
  Let $\theta\in\mathbb{F}_{d}$ with $\Ord_\theta(X)\neq X^d-1$. Assume
  that there exists a solution $\beta$ of (\ref{3.4}) such that
  $\deg\Ord_{\beta}\leq\deg\Ord_{\theta}$. Then there exists a nonzero
  element $\zeta\in\mathbb{F}_{p^d}$ such that 
  \begin{equation}
    g(\sigma)\zeta=0,
    \label{3.5}
  \end{equation}
  where $g(X)=(X^d-1)/\Ord_\theta(X)$. Moreover any such $\zeta$ has the
  property that
  \begin{equation}
    \deg\Ord_{\theta+\zeta}>\deg\Ord_\theta
    \label{3.6}
  \end{equation}
\end{lm}
\begin{proof}
  Let $\gamma$ be a normal element in $\mathbb{F}_{p^d}$ over
  $\mathbb{F}_p$. We can see that
  $\zeta=\Ord_\theta(\sigma)\gamma\neq0$ is a solution of~(\ref{3.5}). In
  fact, we prove that~(\ref{3.6}) is true for any solution $\zeta$
  of~(\ref{3.5}). From~(\ref{3.4}) and Proposition~\ref{orderDiv}, it follows
  that $\Ord_\theta$ divides $\Ord_\beta$, so the hypothesis that
  $\deg\Ord_\beta\leq\deg\Ord_\theta$ implies that
  $\Ord_\beta=\Ord_\theta$. It follows that $g$ and $\Ord_\theta$ are relatively
  prime: suppose $\gcd(g, \Ord_\theta)=d$ and $\deg d>0$, let
  $\Ord_\theta=\Ord_\theta'd=\Ord_\beta$ and $g=g'd$, then
  $\Ord_\theta'(\sigma)\theta=(\Ord_\theta'
  g)(\sigma)\beta=(\Ord_\beta g')(\sigma)\beta=0$ and
  $\deg\Ord_\theta'<\deg\Ord_\beta$,
  that is a contradiction. Since $\Ord_\zeta$ is a divisor of $g$,
  $\Ord_\zeta$ and $\Ord_\theta$ are also relatively prime. Therefore,
  by Corollary~\ref{orderMul}
  \[
    \Ord_{\theta+\zeta}=\Ord_\theta\Ord_\zeta,
  \]
  and then (\ref{3.6}) follows from the fact that $\zeta\neq0$ and then
  $\deg\Ord_\zeta\geq1$.
\end{proof}
We now have a deterministic algorithm to compute a normal element of
$\mathbb{F}_{p^d}$ over $\mathbb{F}_p$.
\begin{enumerate}
  \item Take an element $\theta\in\mathbb{F}_{p^d}$ and compute
    $\Ord_\theta$.
  \item If $\Ord_\theta(X)=X^d-1$ then the algorithm stops ($\theta$ is
    normal).
  \item Calculate $g(X)=(X^d-1)/\Ord_\theta$ and then find
    $\beta\in\mathbb{F}_{p^d}$ such that $g(\sigma)\beta=\theta$.
  \item Determine $\Ord_\beta$. If $\deg\Ord_\beta>\deg\Ord_\theta$ then
    replace $\theta$ by $\beta$ and go to $2$; otherwise if
    $\deg\Ord_\beta\leq\deg\Ord_\theta$ then find a nonzero element $\zeta$
    such that $g(\sigma)\zeta=0$, replace $\theta$ by $\theta+\zeta$ and
    determine the $\sigma$-Order of the new $\theta$, then go to $2$.
\end{enumerate}
This algorithm ends because with each replacement of $\theta$, the degree of
$\Ord_\theta$ increases by at least one, by Lemma~\ref{lmDegree}, and because
this degree is bounded by $d$ since the only possible $\sigma$-Order of degree
$d$ is $X^d-1$. The existence of the elements $\beta$ of step $3$ is guaranteed
by Lemma~\ref{lmExist} and the existence of $\zeta$ of step $4$ by
Lemma~\ref{lmDegree}. Since $g(\sigma)$ is a linear map, we compute $\beta$ and
$\zeta$, using Gauss algorithm and reduced row echelon form again. Computing
elements of type $g(\sigma)\theta$ is not just polynomial evaluation,
and is widely used in our algorithms, so we have a special function
\texttt{sigma\_composition} to do it. The rest of the algorithm is
essentially linear algebra, it uses Flint matrices and is in the function
\texttt{lenstra}.

\subsubsection{Lüneburg's algorithm}
Lüneburg algorithm is more elementary, it is not based on a powerfull theory
like linear algebra. Let $f$ be an irreducible polynomial of degree $d$ over
$\mathbb{F}_p$ and $\alpha$ a root of $f$. Then $\left\{ 1, \alpha, \dots,
  \alpha^{d-1}
\right\}$ is a basis of $\mathbb{F}_{p^d}$ over $\mathbb{F}_p$. Let
$f_i=\Ord_{\alpha^i}$, for $0\leq i < d$. Let
$\gamma=\sum_{i=0}^{d-1}a_i\alpha^i$ be a normal element of
$\mathbb{F}_{p^d}$ over $\mathbb{F}_p$. Note that for any
$a_i\in\mathbb{F}_p$, $\Ord_{a_i\alpha^i}=\Ord_{\alpha^i}=f_i$, and as seen in the
demonstration of Corollary~\ref{orderMul}, $\Ord_\gamma=X^d-1$ divides
$\lcm(f_0, \dots, f_{d-1})$. Since for all $i$, $f_i$ divides $X^d-1$, we have that
$\lcm(f_0, \dots, f_{d-1})$ divides $X^d-1$. Therefore $\lcm(f_0, \dots,
f_{d-1})=X^d-1$. Then we apply the \emph{factor
refinement}~\cite{BaDrSh93}. It is an algorithm that, given a list of
polynomials $f_0, \dots, f_{d-1}$, compute a new list $g_1, \dots, g_r$ of
polynomials, pairwise relatively prime, such that
\[ \prod_{i=0}^{d-1} f_i=\prod_{j=1}^{r} g_j^{e_j}. \]
The factor refinement works that way: we first set $h_i=f_i$ and $e_i=1$ for $0\leq i
<d$ and we consider the list $(h_i, e_i)_{i}$. While there exists $i\neq j$ with $p:=\gcd(h_i, h_j)\neq1$, we delete
$(h_i,e_i)$ and $(h_j,e_j)$ from the list and we add $(p, e_i+e_j),( h_i/p,
e_i), (h_j/p, e_j)$, except for the cases where the first entry is $1$. This
algorithm stops and compute what we expect. In our case, we use it to compute
the pairwise relatively prime polynomials $g_j$ $(1\leq j \leq r)$, and then
we compute the integers $e_{ij}$ such that 
\[
  f_i=\prod_{1\leq j \leq r}g_j^{e_{ij}}\textrm{, for all }0\leq i<d.
\]
This algorithm is implemented in the function
\texttt{factor\_refinement}, it uses Flint's type
\texttt{fq\_poly\_factor} to represent the list $(h_i, e_i)_i$ and to let Flint
deal with the memory allocation.

Next, for each $j$, $1\leq j \leq r$, we find an index $i(j)$ such that
$e_{ij}$ is maximized. Let 
\[
  h_j=f_{i(j)}/g_j^{e_{i(j)j}}
\] 
and take $\beta_j=h_j(\sigma)\alpha^{i(j)}$. Then
\[
  \beta=\sum_{j=1}^r\beta_j
\] 
is a normal element of $\mathbb{F}_{p^d}$ over $\mathbb{F}_p$. In fact, the
$\sigma$-Order of $\beta_j$ is $g_j^{e_{i(j)j}}$ ($\beta_j$ has been
constructed in that purpose) for $1\leq j \leq r$. As $g_1, \dots, g_r$ are
pairwise relatively prime, Corollary~\ref{orderMul} states that the
$\sigma$-Order of $\beta$ must 
\[
  \prod_{j=1}^rg_j^{e_{i(j)j}}=\lcm(f_0, \dots, f_{d-1})=X^d-1,
\]
meaning that $\beta$ is a normal element. 

To implement this algorithm, we use our functions
\texttt{sigma\_composition},
\texttt{sigma\_order}, \texttt{factor\_refinement}, and Flint's basic functions to deal with
polynomials. The implementation is a translation of the steps we followed
to obtain $\beta$.

\section{Experimental results}
In the source code, one can find two directories named \texttt{fq} and
\texttt{fq\_nmod} with very little differences between the files inside these
directories. The name of the file determines the type used in the functions. The
type \texttt{fq} supports arbitrary large integers $p$ for the characteristic of
the fields $\mathbb{F}_{p^d}$, whereas \texttt{fq\_nmod} supports only small
integers for $p$ but is faster.

\subsection{Random algorithms}
We implemented the naive algorithm for computing normal elements, in order to
compare this algorithm to the more elaborate one. Surprinsigly, as we see in
Figure~\ref{NaiveElab}, the speed is
often the same, because when we chose an element randomly in
$\mathbb{F}_{p^d}$, it is often a normal element. If the naive algorithm finds a
normal element with its first pick, there is no way that the the elaborate one is
faster, since it does more calculus. There are still a few cases where the naive
algorithm does not find a normal element easily, and where the elaborate one is
faster.

The time spent in those algorithms is essentially to time spent to perform the
tests. 

\begin{figure}
  \begin{center}
\includegraphics[scale=0.6]{../benchmarks/benchNaiveElab.pdf}
\end{center}
\caption{Benchmark of naive and elaborate algorithms with $p=90001$.}
\label{NaiveElab}
\end{figure}

\subsection{Lüneburg's and Lenstra's algorithms}

Lenstra's algorithm is faster than Lüneburg algorithm, though the complexity
obtained seems to match with the theory in both algorithms. Also, Lenstra's
algotihm is mush more stable than Lüneburg's, the latter can have really
different behaviours for entries of the same sizes. The core of
both Lüneburg's and Lenstra's algorithms is the computation of the
$\sigma$-order polynomials.

  \begin{figure}
\begin{center}
    \begin{tabular}{ccc}
      \textbf{d} & \textbf{Lüneburg} & \textbf{Lenstra} \\ 
     $50$ & $1,87$ s & $0,65$ s \\
     $114$ & $93,0$ s  & $11,7$ s \\
     $205$ & $290$ s & $187$ s \\
     $249$ & $1960$ s & $253$ s 
    \end{tabular}
    \caption{Benchmark of Lüneburg's and Lenstra's algorithms with $p=90001$.}
    \label{LuneburgLenstra}
\end{center}
  \end{figure}


\clearpage
\bibliographystyle{unsrt}
\bibliography{biblio}
\end{document}
