\chapter{Dixon's factorization method\label{chap:dixon}}

~\cite{dixon} describes a class of ``probabilistic algorithms'' for finding a
factor of any composite number, at a sub-exponential cost. They basically
consists into taking random integers $r$ in $\{1, \ldots, N\}$ and look for those
where $r^2 \mod{N}$ is \emph{smooth}. If enough are found, then those integers
can somehow be assembled, and so a fatorization of N attemped.

%% that's not really academic to be stated officially, but I would have never
%% understood this section without Firas (thanks).
%% <http://blog.fkraiem.org/2013/12/08/factoring-integers-dixons-algorithm/>
%% I kept the voila` phrase, that was so lovely.
\section{A little bit of History \label{sec:dixon:history}}
During the latest century there has been a huge effort to approach the problem
formulated by Fermat ~\ref{eq:fermat_problem} from different perspecives. This
led to an entire family of algorithms, like \emph{Quadratic Sieve},
\emph{Dixon}, \ldots.

The core idea is still to find a pair of perfect squares whose difference can
factorize $N$, but maybe Fermat's hypotesis can be made weaker.

\paragraph{Kraitchick} was the first one popularizing the idea the instead of
looking for integers $\angular{x, y}$ such that $x^2 -y^2 = N$ it is sufficient
to look for \emph{multiples} of $N$:
\begin{align}
  x^2 - y^2 \equiv 0 \pmod{N}
\end{align}
and, once found, claim that $\gcd(N, x \pm y)$ are non-trial divisors of $N$
just as we did in \ref{sec:fermat:implementation}.
Kraitchick did not stop here: instead of trying $x^2 \equiv y^2 \pmod{N}$ he
kept the value of previous attempt, and tries to find \emph{a product} of such
values which is also a square. So we have a sequence
\begin{align}
  \label{eq:dixon:x_sequence}
  \angular{x_0, \ldots, x_k} \mid \forall i \leq k \quad x_i^2 - N
  \; \text{ is a perfect square}
\end{align}
and hence
\begin{align*}
  \prod_i (x_i^2 - N) = y^2
\end{align*}
that $\mod{N}$ is equivalent to:
\begin{align}
  \label{eq:dixon:fermat_revisited}
  y^2 \equiv \prod_i (x_i^2 - N) \equiv \big( \prod_i x_i \big) ^2 \pmod{N}
\end{align}
and voil\`a our congruence of squares (\cite{discretelogs} \S 4). For what
concerns the generation of $x_i$ with the property \ref{eq:dixon:x_sequence},
they can simply taken at random and tested using trial division.

\paragraph{Brillhart and Morrison} later proposed (\cite{morrison-brillhart}
p.187) a better approach than trial division to find such $x$. Their idea aims
to ease the enormous effort required by the trial division. In order to achieve
this. they introduce a \emph{factor base} $\factorBase$ and generate random $x$
such that $x^2 - N$ is $\factorBase$-smooth. Recalling what we anticipated in
~\ref{chap:preq}, $\factorBase$ is a precomputed set of primes
$p_i \in \naturalPrime$.
This way the complexity of generating a new $x$ is dominated by
\bigO{|\factorBase|}. Now that the right side of \ref{eq:dixon:fermat_revisited}
has been satisfied, we have to select a subset of those $x$ so that their
product can be seen as a square. Consider an \emph{exponent vector}
$v_i = (\alpha_0, \alpha_1, \ldots, \alpha_r)$ associated with each $x_i$, where
\begin{align}
  \label{eq:dixon:alphas}
  \alpha_j = \begin{cases}
    1 \quad \text{if $p_j$ divides $x_i$ to an odd power} \\
    0 \quad \text{otherwise}
  \end{cases}
\end{align}
for each $1 \leq j \leq r $. There is no need to restrict ourselves for positive
values of $x^2 -N$, so we are going to use $\alpha_0$ to indicate the sign. This
benefit has a neglegible cost: we have to add the non-prime $-1$ to our factor
base $\factorBase$.

Let now $\mathcal{M}$ be the rectangular matrix having per each $i$-th row the
$v_i$ associated to $x_i$: this way each element $m_{ij}$ will be $v_i$'s
$\alpha_j$. We are interested in finding set(s) of the subsequences of $x_i$
whose product always have even powers (\ref{eq:dixon:fermat_revisited}).
Turns out that this is equivalent to look for the set of vectors
$\{ w \mid wM = 0 \} = \ker(\mathcal{M})$ by definition of matrix multiplication
in $\mathbb{F}_2$.


\paragraph{Dixon} Morrison and Brillhart's ideas of \cite{morrison-brillhart}
were actually used for a slightly different factorization method, employing
continued fractions instead of the square difference polynomial. Dixon simply
ported these to the square problem, achieving a probabilistic factorization
method working at a computational cost asymptotically  best than all other ones
previously described: \bigO{\beta(\log N \log \log N)^{\rfrac{1}{2}}} for some
constant $\beta > 0$ \cite{dixon}.

\section{Reduction Procedure}

The following reduction procedure, extracted from ~\cite{morrison-brillhart}, is
a forward part of the Gauss-Jordan elimination algorithm (carried out from right
to left), and can be used to determine whether the set of exponent vectors is
linearly dependent.

For each $v_i$ described as above, associate a \emph{companion history vector}
$h_i = (\beta_0, \beta_1, \ldots, \beta_f)$, where for $0 \leq m \leq f$:
\begin{align*}
  \beta_m = \begin{cases}
    1 \quad \text{ if $m = i$} \\
    0 \quad \text{ otherwise}
    \end{cases}
\end{align*}
At this point, we have all data structures needed:
\\
\\
\\


\begin{center}
  \emph{Reduction Procedure}
\end{center}
\begin{enumerate}[(i)]
  \item Set $j=r$;
  \item find the ``pivot vector'', i.e. the first vector
    $e_i, \quad 0 \leq i \leq f$ such that $\alpha_j = 1$. If none is found, go
    to (iv);
  \item
    \begin{enumerate}[(a)]
      \item replace every following vector $e_m, \quad i < m \leq f$
        whose rightmost $1$ is the $j$-th component, by the sum $e_i \xor e_m$;
      \item whenever $e_m$ is replaced by $e_i \xor e_m$, replace also the
        associated history vector $h_m$ with $h_i \xor h_m$;
    \end{enumerate}
  \item Reduce $j$ by $1$. If $j \geq 0$, return to (ii); otherwise stop.
\end{enumerate}

Algorithm \ref{alg:dixon:kernel} formalizes concepts so far discussed, by
presenting a function \texttt{ker}, discovering linear dependencies in any
rectangular matrix $\mathcal{M} \in (\mathbb{F}_2)^{(f \times r)}$
and storing dependencies into a \emph{history matrix} $\mathcal{H}$.

\begin{remark}
  We are proceeding from right to left in order to conform with
  \cite{morrison-brillhart}.
  Instead, their choice lays on optimization reasons, which does
  not apply any more to a modern calculator.
\end{remark}

\begin{algorithm}
  \caption{Reduction Procedure  \label{alg:dixon:kernel}}
  \begin{algorithmic}[1]
    \Function{Ker}{$\mathcal{M}$}
    \State $\mathcal{H} \gets \texttt{Id}(f \times f)$
    \Comment the initial $\mathcal{H}$ is the identity matrix

    \For{$j = r \strong{ downto } 0$}
    \Comment reduce
      \For{$i=0 \strong{ to } f$}
        \If{$\mathcal{M}_{i, j} = 1$}
          \For{$i' = i \strong{ to } f$}
            \If{$\mathcal{M}_{i', k} = 1$}
              \State $\mathcal{M}_{i'} = \mathcal{M}_i \xor \mathcal{M}_{i'}$
              \State $\mathcal{H}_{i'} = \mathcal{H}_i \xor \mathcal{H}_{i'}$
            \EndIf
          \EndFor
          \State \strong{break}
        \EndIf
      \EndFor
    \EndFor

    \For{$i = 0 \strong{ to } f$}
    \Comment yield linear dependencies
      \If{$\mathcal{M}_i = (0, \ldots, 0)$}
        \strong{yield} $\{\mu  \mid \mathcal{H}_{i,\mu} = 1\}$
      \EndIf
    \EndFor
    \EndFunction
  \end{algorithmic}
\end{algorithm}


\section{Implementation}

Before gluing all toghether, we need one last building brick necessary for
Dixon's factorization algorithm: a \texttt{smooth}($x$) function. In our
specific case, we need a function that, given as input a number $x$, returns the
empty set $\emptyset$ if $x^2 -N$ is not $\factorBase$-smooth. Otherwise,
returns a vector $v = (\alpha_0, \ldots, \alpha_r)$ such that each $\alpha_j$ is
defined just as in \ref{eq:dixon:alphas}. Once we have established $\factorBase$, its
implementation is fairly straightforward:

\begin{algorithm}
  \caption{Discovering Smoothness}
  \begin{algorithmic}[1]
    \Require $\factorBase$, the factor base
    \Procedure{smooth}{$x$}
      \State $v \gets (\alpha_0 = 0, \ldots, \alpha_{|\factorBase|} = 0)$

      \If{$x < 0$} $\alpha_0 \gets 1$ \EndIf
      \For{$i = 1 \strong{ to } |\factorBase|$}
        \If{$\factorBase_i \nmid x$} \strong{continue} \EndIf
        \State $x \gets x// \factorBase_i$
        \State $\alpha_i \gets \alpha_i \xor 1$
      \EndFor
      \If{$x = 1$}
        \State \Return $v$
      \Else
        \State \Return \strong{nil}
      \EndIf
    \EndProcedure
  \end{algorithmic}
\end{algorithm}
\paragraph{How do we choose $\factorBase$?}
It's not easy to answer: if we choose $\factorBase$ small, we will rarely find
$x^2 -N$ \emph{smooth}. If we chose it large, attempting to factorize $x^2 -N$
with $\factorBase$ will pay the price of iterating through a large set.
\cite{Crandall} \S 6.1 finds a solution for this employng complex analytic
number theory. As a  result, the ideal value for $|\factorBase|$ is
$e^{\sqrt{\ln N \ln \ln N}}$.

\begin{algorithm}
  \caption{Dixon}
  \begin{algorithmic}[1]
    \Require $\factorBase$, the factor base
    \Function{dixon}{ }
    \State $i \gets 0$
    \State $r \gets |\factorBase| + 5$
    \Comment finding linearity requires redundance
    \While{$i < r$}
    \Comment search for suitable pairs
    \State $x_i \getsRandom \{0, \ldots N\}$
    \State $y_i \gets x_i^2 - N$
    \State $v_i \gets \texttt{smooth}(y_i)$
    \If{$v_i$} $i \gets i+1$ \EndIf
  \EndWhile
  \State $\mathcal{M} \gets \texttt{matrix}(v_0, \ldots, v_f)$
  \For{$\lambda = \{\mu_0, \ldots, \mu_k\}
    \strong{ in } \texttt{ker}(\mathcal{M})$}
  \Comment get relations
    \State $x \gets \prod\limits_{\mu \in \lambda} x_\mu \pmod{N}$
    \State $y, r \gets \dsqrt{\prod\limits_{\mu \in \lambda} y_\mu \pmod{N}}$
    \If{$\gcd(x+y, N) > 1$}
      \State $p \gets \gcd(x+y, N)$
      \State $q \gets \gcd(x-y, N)$
      \State \Return $p, q$
    \EndIf
  \EndFor
  \EndFunction
  \end{algorithmic}
\end{algorithm}

%%% Local Variables:
%%% mode: latex
%%% TeX-master: "question_authority"
%%% End:
