path: root/doc/exercises.tex

\chapter{Exercises}
\label{chap:Exercises}

% TODO
% 
% All exercises should be group with a chapter
% 
% ▶   Recommended
% M   Matematically oriented
% HM  Higher matematical education required
% 
% 00  Immediate
% 10  Simple
% 20  Medium
% 30  Moderately hard
% 40  Term project
% 50  Research project
% 
% ⁺   High risk of undershoot difficulty


\begin{enumerate}[label=\textbf{\arabic*}.]



\item {[$\RHD$\textit{02}]} \textbf{Saturated subtraction}

Implement the function

\vspace{-1em}
\begin{alltt}
   void monus(z_t r, z_t a, z_t b);
\end{alltt}
\vspace{-1em}

\noindent
which calculates $r = a \dotminus b = \max \{ 0,~ a - b \}$.



\item {[\textit{M10}]} \textbf{Convergence of the Lucas Number ratios}

Find an approximation for
$\displaystyle{ \lim_{n \to \infty} \frac{L_{n + 1}}{L_n}}$,
where $L_n$ is the $n^{\text{th}}$
Lucas Number \psecref{sec:Lucas numbers}.

\( \displaystyle{
    L_n \stackrel{\text{\tiny{def}}}{\text{=}} \left \{ \begin{array}{ll}
      2 & \text{if} ~ n = 0 \\
      1 & \text{if} ~ n = 1 \\
      L_{n - 1} + L_{n + 1} & \text{otherwise}
    \end{array} \right .
}\)



\item {[\textit{M12${}^+$}]} \textbf{Factorisation of factorials}

Implement the function

\vspace{-1em}
\begin{alltt}
   void factor_fact(z_t n);
\end{alltt}
\vspace{-1em}

\noindent
which prints the prime factorisation of $n!$
(the $n^{\text{th}}$ factorial). The function shall
be efficient for all $n$ where all primes $p \le n$
can be found efficiently. You can assume that
$n \ge 2$. You should not evaluate $n!$.



\item {[\textit{M20}]} \textbf{Reverse factorisation of factorials}

{\small\textit{You should already have solved
``Factorisation of factorials'' before you solve
this problem.}}

Implement the function

\vspace{-1em}
\begin{alltt}
   void unfactor_fact(z_t x, z_t *P,
        unsigned long long int *K, size_t n);
\end{alltt}
\vspace{-1em}

\noindent
which given the factorsation of $x!$ determines $x$.
The factorisation of $x!$ is
$\displaystyle{\prod_{i = 1}^{n} P_i^{K_i}}$, where
$P_i$ is \texttt{P[i - 1]} and $K_i$ is \texttt{K[i - 1]}.



\item {[$\RHD$\textit{M17}]} \textbf{Factorials inverted}

Implement the function

\vspace{-1em}
\begin{alltt}
   void unfact(z_t x, z_t n);
\end{alltt}
\vspace{-1em}

\noindent
which given a factorial number $n$, i.e. on the form
$x! = 1 \cdot 2 \cdot 3 \cdot \ldots \cdot x$,
calculates $x = n!^{-1}$. You can assume that
$n$ is a perfect factorial number and that $x \ge 1$.
Extra credit if you can detect when the input, $n$,
is not a factorial number. Such function would of
course return an \texttt{int} 1 if the input is a
factorial and 0 otherwise, or alternatively 0
on success and $-1$ with \texttt{errno} set to
\texttt{EDOM} if the input is not a factorial.



\item {[\textit{05}]} \textbf{Fast primality test}

$(x + y)^p \equiv x^p + y^p ~(\text{Mod}~p)$
for all primes $p$ and for a few composites $p$,
which are know as pseudoprimes. Use this to implement
a fast primality tester.



\item {[\textit{10}]} \textbf{Fermat primality test}

$a^{p - 1} \equiv 1 ~(\text{Mod}~p) ~\forall~ 1 < a < p$
for all primes $p$ and for a few composites $p$,
which are know as pseudoprimes\footnote{If $p$ is composite
but passes the test for all $a$, $p$ is a Carmichael
number.}. Use this to implement a heuristic primality
tester. Try to mimic \texttt{zptest} as much as possible.
GNU~MP uses $a = 210$, but you don't have to. ($a$ is
called a base.)



\item {[\textit{11}]} \textbf{Lucas–Lehmer primality test}

The Lucas–Lehmer primality test can be used to determine
whether a Mersenne numbers $M_n = 2^n - 1$ is a prime (a
Mersenne prime). $M_n$ is a prime if and only if
$s_{n - 1} \equiv 0 ~(\text{Mod}~M_n)$, where

\( \displaystyle{
    s_i = \left \{ \begin{array}{ll}
      4 & \text{if} ~ i = 0 \\
      s_{i - 1}^2 - 2 & \text{otherwise}.
    \end{array} \right .
}\)

\noindent
The Lucas–Lehmer primality test requires that $n \ge 3$,
however $M_2 = 2^2 - 1 = 3$ is a prime. Implement a version
of the Lucas–Lehmer primality test that takes $n$ as the
input. For some more fun, when you are done, you can
implement a version that takes $M_n$ as the input.

For improved performance, instead of using \texttt{zmod},
you can use the recursive function
%
\( \displaystyle{
    k \mod (2^n - 1) =
    \left (
      (k \mod 2^n) + \lfloor k \div 2^n \rfloor
    \right ) \mod (2^n - 1),
}\)
%
where $k \mod 2^n$ is efficiently calculated
using \texttt{zand($k$, $2^n - 1$)}. (This optimisation
is not part of the difficulty rating of this problem.)



\item {[\textit{20}]} \textbf{Fast primality test with bounded perfection}

Implement a primality test that is both very fast and
never returns \texttt{PROBABLY\_PRIME} for input less
than or equal to a preselected number.



\item {[\textit{30}]} \textbf{Powers of the golden ratio}

Implement function that returns $\varphi^n$ rounded
to the nearest integer, where $n$ is the input and
$\varphi$ is the golden ratio.



\end{enumerate}



\chapter{Solutions}
\label{chap:Solutions}


\begin{enumerate}[label=\textbf{\arabic*}.]

\item \textbf{Saturated subtraction}

\vspace{-1em}
\begin{alltt}
void monus(z_t r, z_t a, z_t b)
\{
    zsub(r, a, b);
    if (zsignum(r) < 0)
        zsetu(r, 0);
\}
\end{alltt}



\item \textbf{Convergence of the Lucas Number ratios}

It would be a mistake to use bignum, and bigint in particular,
to solve this problem. Good old mathematics is a much better solution.

$$ \lim_{n \to \infty} \frac{L_{n + 1}}{L_n} = \lim_{n \to \infty} \frac{L_{n}}{L_{n - 1}} = \lim_{n \to \infty} \frac{L_{n - 1}}{L_{n - 2}} $$

$$ \frac{L_{n}}{L_{n - 1}} = \frac{L_{n - 1}}{L_{n - 2}} $$

$$ \frac{L_{n - 1} + L_{n - 2}}{L_{n - 1}} = \frac{L_{n - 1}}{L_{n - 2}} $$

$$ 1 + \frac{L_{n - 2}}{L_{n - 1}} = \frac{L_{n - 1}}{L_{n - 2}} $$

$$ 1 + \varphi = \frac{1}{\varphi} $$

So the ratio tends toward the golden ratio.



\item \textbf{Factorisation of factorials}

Base your implementation on

\( \displaystyle{
    n! = \prod_{p~\in~\textbf{P}}^{n} p^{k_p}, ~\text{where}~
    k_p = \sum_{a = 1}^{\lfloor \log_p n \rfloor} \lfloor np^{-a} \rfloor.
}\)

There is no need to calculate $\lfloor \log_p n \rfloor$,
you will see when $p^a > n$.



\item \textbf{Reverse factorisation of factorials}

$\displaystyle{x = \max_{p ~\in~ P} ~ p \cdot f(p, k_p)}$,
where $k_p$ is the power of $p$ in the factorisation
of $x!$. $f(p, k)$ is defined as:

\vspace{1em}
\hspace{-2.8ex}
\begin{minipage}{\linewidth}
\begin{algorithmic}
    \STATE $k^\prime \gets 0$
    \WHILE{$k > 0$}
      \STATE $a \gets 0$
      \WHILE{$p^a \le k$}
        \STATE $k \gets k - p^a$
        \STATE $a \gets a + 1$
      \ENDWHILE
      \STATE $k^\prime \gets k^\prime + p^{a - 1}$
    \ENDWHILE
    \RETURN $k^\prime$
\end{algorithmic}
\end{minipage}
\vspace{1em}



\item \textbf{Factorials inverted}

Use \texttt{zlsb} to get the power of 2 in the
prime factorisation of $n$, that is, the number
of times $n$ is divisible by 2. If we write $n$ on
the form $1 \cdot 2 \cdot 3 \cdot \ldots \cdot x$,
every $2^\text{nd}$ factor is divisible by 2, every
$4^\text{th}$ factor is divisible by $2^2$, and so on.
From calling \texttt{zlsb} we know how many times,
$n$ is divisible by 2, but know how many of the factors
are divisible by 2, but this can be calculated with
the following algorithm, where $k$ is the number
times $n$ is divisible by 2:

\vspace{1em}
\hspace{-2.8ex}
\begin{minipage}{\linewidth}
\begin{algorithmic}
    \STATE $k^\prime \gets 0$
    \WHILE{$k > 0$}
      \STATE $a \gets 0$
      \WHILE{$2^a \le k$}
        \STATE $k \gets k - 2^a$
        \STATE $a \gets a + 1$
      \ENDWHILE
      \STATE $k^\prime \gets k^\prime + 2^{a - 1}$
    \ENDWHILE
    \RETURN $k^\prime$
\end{algorithmic}
\end{minipage}
\vspace{1em}

\noindent
Note that $2^a$ is efficiently calculated with,
\texttt{zlsh}, but it is more efficient to use
\texttt{zbset}.

Now that we know $k^\prime$, the number of
factors ni $1 \cdot \ldots \cdot x$ that are
divisible by 2, we have two choices for $x$:
$k^\prime$ and $k^\prime + 1$. To check which, we
calculate $(k^\prime - 1)!!$ (the semifactoral, i.e.
$1 \cdot 3 \cdot 5 \cdot \ldots \cdot (k^\prime - 1)$)
naïvely and shift the result $k$ steps to the left.
This gives us $k^\prime!$. If $x! \neq k^\prime!$, then
$x = k^\prime + 1$ unless $n$ is not factorial number.
Of course, if $x! = k^\prime!$, then $x = k^\prime$.



\item \textbf{Fast primality test}

If we select $x = y = 1$ we get
$2^p \equiv 2 ~(\text{Mod}~p)$. This gives us

\vspace{-1em}
\begin{alltt}
enum zprimality ptest_fast(z_t p)
\{
    z_t a;
    int c = zcmpu(p, 2);
    if (c <= 0)
        return c ? NONPRIME : PRIME;
    zinit(a);
    zsetu(a, 1);
    zlsh(a, a, p);
    zmod(a, a, p);
    c = zcmpu(a, 2);
    zfree(a);
    return c ? NONPRIME : PROBABLY_PRIME;
\}
\end{alltt}



\item \textbf{Fermat primality test}

Below is a simple implementation. It can be improved by
just testing against a fix base, such as $a = 210$, it
$t = 0$. It could also do an exhaustive test (all $a$
such that $1 < a < p$) if $t < 0$.

\vspace{-1em}
\begin{alltt}
enum zprimality ptest_fermat(z_t witness, z_t p, int t)
\{
    enum zprimality rc = PROBABLY_PRIME;
    z_t a, p1, p3, temp;
    int c;

    if ((c = zcmpu(p, 2)) <= 0) \{
        if (!c)
            return PRIME;
        if (witness && witness != p)
            zset(witness, p);
        return NONPRIME;
    \}

    zinit(a), zinit(p1), zinit(p3), zinit(temp);
    zsetu(temp, 3), zsub(p3, p, temp);
    zsetu(temp, 1), zsub(p1, p, temp);

    zsetu(temp, 2);
    while (t--) \{
        zrand(a, DEFAULT_RANDOM, UNIFORM, p3);
        zadd(a, a, temp);
        zmodpow(a, a, p1, p);
        if (zcmpu(a, 1)) \{
            if (witness)
                zswap(witness, a);
            rc = NONPRIME;
            break;
        \}
    \}

    zfree(temp), zfree(p3), zfree(p1), zfree(a);
    return rc;
\}
\end{alltt}



\item \textbf{Lucas–Lehmer primality test}

\vspace{-1em}
\begin{alltt}
enum zprimality ptest_llt(z_t n)
\{
    z_t s, M;
    int c;
    size_t p;

    if ((c = zcmpu(n, 2)) <= 0)
        return c ? NONPRIME : PRIME;

    if (n->used > 1) \{
        \textcolor{c}{/* \textrm{An optimised implementation would not need this} */}
        errno = ENOMEM;
        return (enum zprimality)(-1);
    \}

    zinit(s), zinit(M), zinit(2);

    p = (size_t)(n->chars[0]);
    zsetu(s, 1), zsetu(M, 0);
    zbset(M, M, p, 1), zsub(M, M, s);
    zsetu(s, 4);
    zsetu(two, 2);

    p -= 2;
    while (p--) \{
        zsqr(s, s);
        zsub(s, s, two);
        zmod(s, s, M);
    \}
    c = zzero(s);

    zfree(two), zfree(M), zfree(s);
    return c ? PRIME : NONPRIME;
\}
\end{alltt}

$M_n$ is composite if $n$ is composite, therefore,
if you do not expect prime-only values on $n$, the
performance can be improve by using some other
primality test (or this same test if $n$ is a
Mersenne number) to first check that $n$ is prime.



\item \textbf{Fast primality test with bounded perfection}

First we select a fast primality test. We can use
$2^p \equiv 2 ~(\texttt{Mod}~ p) ~\forall~ p \in \textbf{P}$,
as describe in the solution for the problem
\textit{Fast primality test}.

Next, we use this to generate a large list of primes and
pseudoprimes. Use a perfect primality test, such as a
naïve test or the AKS primality test, to filter out all
primes and retain only the pseudoprimes. This is not in
runtime so it does not matter that this is slow, but to
speed it up, we can use a probabilistic test such the
Miller–Rabin primality test (\texttt{zptest}) before we
use the perfect test.

Now that we have a quite large — but not humongous — list
of pseudoprimes, we can incorporate it into our fast
primality test. For any input that passes the test, and
is less or equal to the largest pseudoprime we found,
binary search our list of pseudoprime for the input.

For input, larger than our limit, that passes the test,
we can run it through \texttt{zptest} to reduce the
number of false positives.

As an alternative solution, instead of comparing against
known pseudoprimes. Find a minimal set of primes that
includes divisors for all known pseudoprimes, and do
trail division with these primes when a number passes
the test. No pseudoprime need to have more than one divisor
included in the set, so this set cannot be larger than
the set of pseudoprimes.



\item \textbf{Powers of the golden ratio}

This was an information gathering exercise.
For $n \ge 1$, $L_n = [\varphi^n]$, where
$L_n$ is the $n^\text{th}$ Lucas number.

\( \displaystyle{
    L_n \stackrel{\text{\tiny{def}}}{\text{=}} \left \{ \begin{array}{ll}
      2 & \text{if} ~ n = 0 \\
      1 & \text{if} ~ n = 1 \\
      L_{n - 1} + L_{n + 1} & \text{otherwise}
    \end{array} \right .
}\)

\noindent
but for efficiency and briefness, we will use
\texttt{lucas} from \secref{sec:Lucas numbers}.

\vspace{-1em}
\begin{alltt}
void golden_pow(z_t r, z_t p)
\{
    if (zsignum(p) <= 0)
        zsetu(r, zzero(p));
    else
        lucas(r, p);
\}
\end{alltt}



\end{enumerate}