\chapter{Not implemented}
\label{chap:Not implemented}

In this chapter we maintain a list of
features we have chosen not to implement,
but would fit into libzahl had we not have
our priorities straight. Functions listed
herein will only be implemented if it is
shown that it would be overwhelmingly
advantageous. For each feature, a sample
implementation or a mathematical expression
on which you can base your implemention
is included. The sample implementations create
temporary integer references, this is to
simplify the examples. You should try to
use dedicated variables; in case of recursion,
a robust program should store temporary
variables on a stack, so they can be
clean up if something happens.

Research problems, like prime-factorisation
and discrete logarithms do not fit in the
scope of bignum libraries. % Unless they are extraordinarily bloated with vague mission-scope, like systemd.
And therefore do not fit into libzahl,
and will not be included in this chapter.
Operators and functions that grow so
ridiculously fast that a tiny lookup table
constructed to cover all practical input
will also not be included in this chapter,
nor in libzahl.

\vspace{1cm}
\minitoc


\newpage
\section{Extended greatest common divisor}
\label{sec:Extended greatest common divisor}

\begin{alltt}
void
extgcd(z_t bézout_coeff_1, z_t bézout_coeff_2, z_t gcd
       z_t quotient_1, z_t quotient_2, z_t a, z_t b)
\{
#define old_r gcd
#define old_s bézout_coeff_1
#define old_t bézout_coeff_2
#define s quotient_2
#define t quotient_1
    z_t r, q, qs, qt;
    int odd = 0;
    zinit(r), zinit(q), zinit(qs), zinit(qt);
    zset(r, b), zset(old_r, a);
    zseti(s, 0), zseti(old_s, 1);
    zseti(t, 1), zseti(old_t, 0);
    while (!zzero(r)) \{
        odd ^= 1;
        zdivmod(q, old_r, old_r, r), zswap(old_r, r);
        zmul(qs, q, s), zsub(old_s, old_s, qs);
        zmul(qt, q, t), zsub(old_t, old_t, qt);
        zswap(old_s, s), zswap(old_t, t);
    \}
    odd ? abs(s, s) : abs(t, t);
    zfree(r), zfree(q), zfree(qs), zfree(qt);
\}
\end{alltt}

Perhaps you are asking yourself ``wait a minute,
doesn't the extended Euclidean algorithm only
have three outputs if you include the greatest
common divisor, what is this shenanigans?''
No\footnote{Well, technically yes, but it calculates
two values for free in the same ways as division
calculates the remainder for free.}, it has five
outputs, most implementations just ignore two of
them. If this confuses you, or you want to know
more about this, I refer you to Wikipeida.


\newpage
\section{Least common multiple}
\label{sec:Least common multiple}

\( \displaystyle{
    \mbox{lcm}(a, b) = \frac{\lvert a \cdot b \rvert}{\mbox{gcd}(a, b)}
}\)

Be aware, {\tt zgcd} can return zero.


\newpage
\section{Modular multiplicative inverse}
\label{sec:Modular multiplicative inverse}

\begin{alltt}
int
modinv(z_t inv, z_t a, z_t m)
\{
    z_t x, _1, _2, _3, gcd, mabs, apos;
    int invertible, aneg = zsignum(a) < 0;
    zinit(x), zinit(_1), zinit(_2), zinit(_3), zinit(gcd);
    *mabs = *m;
    zabs(mabs, mabs);
    if (aneg) \{
        zinit(apos);
        zset(apos, a);
        if (zcmpmag(apos, mabs))
            zmod(apos, apos, mabs);
        zadd(apos, apos, mabs);
    \}
    extgcd(inv, _1, _2, _3, gcd, apos, mabs);
    if ((invertible = !zcmpi(gcd, 1))) \{
        if (zsignum(inv) < 0)
            (zsignum(m) < 0 ? zsub : zadd)(x, x, m);
        zswap(x, inv);
    \}
    if (aneg)
        zfree(apos);
    zfree(x), zfree(_1), zfree(_2), zfree(_3), zfree(gcd);
    return invertible;
\}
\end{alltt}


\newpage
\section{Random prime number generation}
\label{sec:Random prime number generation}

TODO


\newpage
\section{Symbols}
\label{sec:Symbols}

\subsection{Legendre symbol}
\label{sec:Legendre symbol}

\( \displaystyle{
  \left ( \frac{a}{p} \right ) \equiv a^{\frac{p - 1}{2}} ~(\text{Mod}~p),~
  \left ( \frac{a}{p} \right ) \in \{-1,~0,~1\},~
  p \in \textbf{P},~ p > 2
}\)
\vspace{1em}

\noindent
That is, unless $\displaystyle{a^{\frac{p - 1}{2}} \mod p \le 1}$,
$\displaystyle{a^{\frac{p - 1}{2}} \mod p = p - 1}$, so
$\displaystyle{\left ( \frac{a}{p} \right ) = -1}$.

It should be noted that
\( \displaystyle{
  \left ( \frac{a}{p} \right ) = 
  \left ( \frac{a ~\text{Mod}~ p}{p} \right ),
}\)
so a compressed lookup table can be used for small $p$.


\subsection{Jacobi symbol}
\label{sec:Jacobi symbol}

\( \displaystyle{
  \left ( \frac{a}{n} \right ) = 
  \prod_k \left ( \frac{a}{p_k} \right )^{n_k},
}\)
where $n$ = $\displaystyle{\prod_k p_k^{n_k}}$, and $p_k \in \textbf{P}$.
\vspace{1em}

Like the Legendre symbol, the Jacobi symbol is $n$-period over $a$.
If $n$, is prime, the Jacobi symbol is the Legendre symbol, but
the Jacobi symbol is defined for all odd numbers $n$, where the
Legendre symbol is defined only for odd primes $n$.

Use the following algorithm to calculate the Jacobi symbol:

\vspace{1em}
\hspace{-2.8ex}
\begin{minipage}{\linewidth}
\begin{algorithmic}
    \STATE $a \gets a \mod n$
    \STATE $r \gets \mbox{lsb}~ a$
    \STATE $a \gets a \cdot 2^{-r}$
    \STATE \(\displaystyle{
      r \gets \left \lbrace \begin{array}{rl}
        1 &
          \text{if}~ n \equiv 1, 7 ~(\text{Mod}~ 8)
          ~\text{or}~ r \equiv 0 ~(\text{Mod}~ 2) \\
        -1 & \text{otherwise} \\
      \end{array} \right .
    }\)
    \IF{$a = 1$}
        \RETURN $r$
    \ELSIF{$\gcd(a, n) \neq 1$}
        \RETURN $0$
    \ENDIF
    \STATE $(a, n) = (n, a)$
    \STATE \textbf{start over}
\end{algorithmic}
\end{minipage}
\vspace{1em}



\subsection{Kronecker symbol}
\label{sec:Kronecker symbol}

TODO


\subsection{Power residue symbol}
\label{sec:Power residue symbol}

TODO


\subsection{Pochhammer \emph{k}-symbol}
\label{sec:Pochhammer k-symbol}

\( \displaystyle{
    (x)_{n,k} = \prod_{i = 1}^n (x + (i - 1)k)
}\)


\newpage
\section{Logarithm}
\label{sec:Logarithm}

TODO


\newpage
\section{Roots}
\label{sec:Roots}

TODO


\newpage
\section{Modular roots}
\label{sec:Modular roots}

TODO % Square: Cipolla's algorithm, Pocklington's algorithm, Tonelli–Shanks algorithm


\newpage
\section{Combinatorial}
\label{sec:Combinatorial}

\subsection{Factorial}
\label{sec:Factorial}

\( \displaystyle{
    n! = \left \lbrace \begin{array}{ll}
        \displaystyle{\prod_{i = 1}^n i} & \textrm{if}~ n \ge 0 \\
        \textrm{undefined} & \textrm{otherwise}
    \end{array} \right .
}\)
\vspace{1em}

This can be implemented much more efficently
than using the naïve method, and is a very
important function for many combinatorial
applications, therefore it may be implemented
in the future if the demand is high enough.

An efficient, yet not optimal, implementation
of factorials that about halves the number of
required multiplications compared to the naïve
method can be derived from the observation

\vspace{1em}
\( \displaystyle{
    n! = n!! ~ \lfloor n / 2 \rfloor! ~ 2^{\lfloor n / 2 \rfloor}
}\), $n$ odd.
\vspace{1em}

\noindent
The resulting algorithm can be expressed

\begin{alltt}
   void
   fact(z_t r, uint64_t n)
   \{
       z_t p, f, two;
       uint64_t *ns, s = 1, i = 1;
       zinit(p), zinit(f), zinit(two);
       zseti(r, 1), zseti(p, 1), zseti(f, n), zseti(two, 2);
       ns = alloca(zbits(f) * sizeof(*ns));
       while (n > 1) \{
           if (n & 1) \{
               ns[i++] = n;
               s += n >>= 1;
           \} else \{
               zmul(r, r, (zsetu(f, n), f));
               n -= 1;
           \}
       \}
       for (zseti(f, 1); i-- > 0; zmul(r, r, p);)
           for (n = ns[i]; zcmpu(f, n); zadd(f, f, two))
               zmul(p, p, f);
       zlsh(r, r, s);
       zfree(two), zfree(f), zfree(p);
   \}
\end{alltt}


\subsection{Subfactorial}
\label{sec:Subfactorial}

\( \displaystyle{
    !n = \left \lbrace \begin{array}{ll}
      n(!(n - 1)) + (-1)^n & \textrm{if}~ n > 0 \\
      1 & \textrm{if}~ n = 0 \\
      \textrm{undefined} & \textrm{otherwise}
    \end{array} \right . =
    n! \sum_{i = 0}^n \frac{(-1)^i}{i!}
}\)


\subsection{Alternating factorial}
\label{sec:Alternating factorial}

\( \displaystyle{
    \mbox{af}(n) = \sum_{i = 1}^n {(-1)^{n - i} i!}
}\)


\subsection{Multifactorial}
\label{sec:Multifactorial}

\( \displaystyle{
    n!^{(k)} = \left \lbrace \begin{array}{ll}
      1 & \textrm{if}~ n = 0 \\
      n & \textrm{if}~ 0 < n \le k \\
      n((n - k)!^{(k)}) & \textrm{if}~ n > k \\
      \textrm{undefined} & \textrm{otherwise}
    \end{array} \right .
}\)


\subsection{Quadruple factorial}
\label{sec:Quadruple factorial}

\( \displaystyle{
    (4n - 2)!^{(4)}
}\)


\subsection{Superfactorial}
\label{sec:Superfactorial}

\( \displaystyle{
    \mbox{sf}(n) = \prod_{k = 1}^n k^{1 + n - k}
}\), undefined for $n < 0$.


\subsection{Hyperfactorial}
\label{sec:Hyperfactorial}

\( \displaystyle{
    H(n) = \prod_{k = 1}^n k^k
}\), undefined for $n < 0$.


\subsection{Raising factorial}
\label{sec:Raising factorial}

\( \displaystyle{
    x^{(n)} = \frac{(x + n - 1)!}{(x - 1)!}
}\), undefined for $n < 0$.


\subsection{Falling factorial}
\label{sec:Falling factorial}

\( \displaystyle{
    (x)_n = \frac{x!}{(x - n)!}
}\), undefined for $n < 0$.


\subsection{Primorial}
\label{sec:Primorial}

\( \displaystyle{
    n\# = \prod_{\lbrace i \in \textbf{P} ~:~ i \le n \rbrace} i
}\)
\vspace{1em}

\noindent
\( \displaystyle{
    p_n\# = \prod_{i \in \textbf{P}_{\pi(n)}} i
}\)


\subsection{Gamma function}
\label{sec:Gamma function}

$\Gamma(n) = (n - 1)!$, undefined for $n \le 0$.


\subsection{K-function}
\label{sec:K-function}

\( \displaystyle{
    K(n) = \left \lbrace \begin{array}{ll}
      \displaystyle{\prod_{i = 1}^{n - 1} i^i}  & \textrm{if}~ n \ge 0 \\
      1 & \textrm{if}~ n = -1 \\
      0 & \textrm{otherwise (result is truncated)}
    \end{array} \right .
}\)


\subsection{Binomial coefficient}
\label{sec:Binomial coefficient}

\( \displaystyle{
    \binom{n}{k} = \frac{n!}{k!(n - k)!}
    = \frac{1}{(n - k)!} \prod_{i = k + 1}^n i
    = \frac{1}{k!} \prod_{i = n - k + 1}^n i
}\)


\subsection{Catalan number}
\label{sec:Catalan number}

\( \displaystyle{
    C_n = \left . \binom{2n}{n} \middle / (n + 1) \right .
}\)


\subsection{Fuss–Catalan number}
\label{sec:Fuss-Catalan number} % not en dash

\( \displaystyle{
    A_m(p, r) = \frac{r}{mp + r} \binom{mp + r}{m}
}\)


\newpage
\section{Fibonacci numbers}
\label{sec:Fibonacci numbers}

Fibonacci numbers can be computed efficiently
using the following algorithm:

\begin{alltt}
   static void
   fib_ll(z_t f, z_t g, z_t n)
   \{
       z_t a, k;
       int odd;
       if (zcmpi(n, 1) <= 0) \{
           zseti(f, !zzero(n));
           zseti(f, zzero(n));
           return;
       \}
       zinit(a), zinit(k);
       zrsh(k, n, 1);
       if (zodd(n)) \{
           odd = zodd(k);
           fib_ll(a, g, k);
           zadd(f, a, a);
           zadd(k, f, g);
           zsub(f, f, g);
           zmul(f, f, k);
           zseti(k, odd ? -2 : +2);
           zadd(f, f, k);
           zadd(g, g, g);
           zadd(g, g, a);
           zmul(g, g, a);
       \} else \{
           fib_ll(g, a, k);
           zadd(f, a, a);
           zadd(f, f, g);
           zmul(f, f, g);
           zsqr(a, a);
           zsqr(g, g);
           zadd(g, a);
       \}
       zfree(k), zfree(a);
   \}
\end{alltt}

\newpage
\begin{alltt}
   void
   fib(z_t f, z_t n)
   \{
       z_t tmp, k;
       zinit(tmp), zinit(k);
       zset(k, n);
       fib_ll(f, tmp, k);
       zfree(k), zfree(tmp);
   \}
\end{alltt}

\noindent
This algorithm is based on the rules

\vspace{1em}
\( \displaystyle{
    F_{2k + 1} = 4F_k^2 - F_{k - 1}^2 + 2(-1)^k = (2F_k + F_{k-1})(2F_k - F_{k-1}) + 2(-1)^k
}\)
\vspace{1em}

\( \displaystyle{
    F_{2k} = F_k \cdot (F_k + 2F_{k - 1})
}\)
\vspace{1em}

\( \displaystyle{
    F_{2k - 1} = F_k^2 + F_{k - 1}^2
}\)
\vspace{1em}

\noindent
Each call to {\tt fib\_ll} returns $F_n$ and $F_{n - 1}$
for any input $n$. $F_{k}$ is only correctly returned
for $k \ge 0$. $F_n$ and $F_{n - 1}$ is used for
calculating $F_{2n}$ or $F_{2n + 1}$. The algorithm
can be sped up with a larger lookup table than one
covering just the base cases. Alternatively, a naïve
calculation could be used for sufficiently small input.


\newpage
\section{Lucas numbers}
\label{sec:Lucas numbers}

Lucas [lyk\textscripta] numbers can be calculated by utilising
{\tt fib\_ll} from \secref{sec:Fibonacci numbers}:

\begin{alltt}
   void
   lucas(z_t l, z_t n)
   \{
       z_t k;
       int odd;
       if (zcmp(n, 1) <= 0) \{
           zset(l, 1 + zzero(n));
           return;
       \}
       zinit(k);
       zrsh(k, n, 1);
       if (zeven(n)) \{
           lucas(l, k);
           zsqr(l, l);
           zseti(k, zodd(k) ? +2 : -2);
           zadd(l, k);
       \} else \{
           odd = zodd(k);
           fib_ll(l, k, k);
           zadd(l, l, l);
           zadd(l, l, k);
           zmul(l, l, k);
           zseti(k, 5);
           zmul(l, l, k);
           zseti(k, odd ? +4 : -4);
           zadd(l, l, k);
       \}
       zfree(k);
   \}
\end{alltt}

\noindent
This algorithm is based on the rules

\vspace{1em}
\( \displaystyle{
    L_{2k} = L_k^2 - 2(-1)^k
}\)
\vspace{1ex}

\( \displaystyle{
    L_{2k + 1} = 5F_{k - 1} \cdot (2F_k + F_{k - 1}) - 4(-1)^k
}\)
\vspace{1em}

\noindent
Alternatively, the function can be implemented
trivially using the rule

\vspace{1em}
\( \displaystyle{
    L_k = F_k + 2F_{k - 1}
}\)


\newpage
\section{Bit operation}
\label{sec:Bit operation unimplemented}

\subsection{Bit scanning}
\label{sec:Bit scanning}

Scanning for the next set or unset bit can be
trivially implemented using {\tt zbtest}. A
more efficient, although not optimally efficient,
implementation would be

\begin{alltt}
   size_t
   bscan(z_t a, size_t whence, int direction, int value)
   \{
       size_t ret;
       z_t t;
       zinit(t);
       value ? zset(t, a) : znot(t, a);
       ret = direction < 0
           ? (ztrunc(t, t, whence + 1), zbits(t) - 1)
           : (zrsh(t, t, whence), zlsb(t) + whence);
       zfree(t);
       return ret;
   \}
\end{alltt}


\subsection{Population count}
\label{sec:Population count}

The following function can be used to compute
the population count, the number of set bits,
in an integer, counting the sign bit:

\begin{alltt}
   size_t
   popcount(z_t a)
   \{
       size_t i, ret = zsignum(a) < 0;
       for (i = 0; i < a->used; i++) \{
           ret += __builtin_popcountll(a->chars[i]);
       \}
       return ret;
   \}
\end{alltt}

\noindent
It requires a compiler extension, if missing,
there are other ways to computer the population
count for a word: manually bit-by-bit, or with
a fully unrolled

\begin{alltt}
   int s;
   for (s = 1; s < 64; s <<= 1)
       w = (w >> s) + w;
\end{alltt}


\subsection{Hamming distance}
\label{sec:Hamming distance}

A simple way to compute the Hamming distance,
the number of differing bits, between two
numbers is with the function

\begin{alltt}
   size_t
   hammdist(z_t a, z_t b)
   \{
       size_t ret;
       z_t t;
       zinit(t);
       zxor(t, a, b);
       ret = popcount(t);
       zfree(t);
       return ret;
   \}
\end{alltt}

\noindent
The performance of this function could
be improved by comparing character by
character manually using {\tt zxor}.


\newpage
\section{Miscellaneous}
\label{sec:Miscellaneous}


\subsection{Character retrieval}
\label{sec:Character retrieval}

\begin{alltt}
uint64_t
getu(z_t a)
\{
    return zzero(a) ? 0 : a->chars[0];
\}
\end{alltt}

\subsection{Fit test}
\label{sec:Fit test}

Some libraries have functions for testing
whether a big integer is small enough to
fit into an intrinsic type. Since libzahl
does not provide conversion to intrinsic
types this is irrelevant. But additionally,
it can be implemented with a single
one-line macro that does not have any
side-effects.

\begin{alltt}
   #define fits_in(a, type)  (zbits(a) <= 8 * sizeof(type))
   \textcolor{c}{/* \textrm{Just be sure the type is integral.} */}
\end{alltt}


\subsection{Reference duplication}
\label{sec:Reference duplication}

This could be useful for creating duplicates
with modified sign. But only if neither
{\tt r} nor {\tt a} will be modified whilst
both are in use. Because it is unsafe,
fairly simple to create an implementation
with acceptable performance — {\tt *r = *a},
— and probably seldom useful, this has not
be implemented.

\begin{alltt}
   int
   refdup(z_t r, z_t a)
   \{
       \textcolor{c}{/* \textrm{Almost fully optimised, but perfectly portable} *r = *a; */}
       r->sign    = a->sign;
       r->used    = a->used;
       r->alloced = a->alloced;
       r->chars   = a->chars;
   \}
\end{alltt}


\subsection{Variadic initialisation}
\label{sec:Variadic initialisation}

Most bignum libraries have variadic functions
for initialisation and uninitialisation. This
is not available in libzahl, because it is
not useful enough and has performance overhead.
And what's next, support {\tt va\_list},
variadic addition, variadic multiplication,
power towers, set manipulation? Anyone can
implement variadic wrapper for {\tt zinit} and
{\tt zfree} if they really need it. But if
you want to avoid the overhead, you can use
something like this:

\begin{alltt}
   /* \textrm{Call like this:} MANY(zinit, (a), (b), (c)) */
   #define MANY(f, ...)  (_MANY1(f, __VA_ARGS__,,,,,,,,,))
   
   #define _MANY1(f, a, ...)  (void)f a, _MANY2(f, __VA_ARGS__)
   #define _MANY2(f, a, ...)  (void)f a, _MANY3(f, __VA_ARGS__)
   #define _MANY3(f, a, ...)  (void)f a, _MANY4(f, __VA_ARGS__)
   #define _MANY4(f, a, ...)  (void)f a, _MANY5(f, __VA_ARGS__)
   #define _MANY5(f, a, ...)  (void)f a, _MANY6(f, __VA_ARGS__)
   #define _MANY6(f, a, ...)  (void)f a, _MANY7(f, __VA_ARGS__)
   #define _MANY7(f, a, ...)  (void)f a, _MANY8(f, __VA_ARGS__)
   #define _MANY8(f, a, ...)  (void)f a, _MANY9(f, __VA_ARGS__)
   #define _MANY9(f, a, ...)  (void)f a
\end{alltt}