libzahl

big integer library
git clone git://git.suckless.org/libzahl
Log | Files | Refs | README | LICENSE

commit 16a001c8fe4e5ca99d5aafdd8ed02a35f09b6caa
parent c98c0af42f56d0859d6f3f5e3743945906c681c8
Author: Mattias Andrée <maandree@kth.se>
Date:   Fri, 13 May 2016 04:38:09 +0200

Miscellaneous stuff

Signed-off-by: Mattias Andrée <maandree@kth.se>

Diffstat:
Makefile | 6+++++-
TODO | 4++++
doc/arithmetic.tex | 5+++--
doc/bit-operations.tex | 56++++++++++++++++++++++++++++++++++++++++++++++++++++++++
doc/libzahl.tex | 4++++
doc/not-implemented.tex | 554+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
doc/number-theory.tex | 35+++++++++++++++++++++++++++++++++++
doc/random-numbers.tex | 28++++++++++++++++++++++++++++
src/zmodmul.c | 1+
src/zmodpow.c | 2++
src/zmodpowu.c | 5+++--
src/zpow.c | 2++
src/zpowu.c | 5+++--
13 files changed, 700 insertions(+), 7 deletions(-)

diff --git a/Makefile b/Makefile @@ -78,7 +78,11 @@ TEXSRC =\ doc/libzahls-design.tex\ doc/get-started.tex\ doc/miscellaneous.tex\ - doc/arithmetic.tex + doc/arithmetic.tex\ + doc/bit-operations.tex\ + doc/number-theory.tex\ + doc/random-numbers.tex\ + doc/not-implemented.tex HDR_PUBLIC = zahl.h $(HDR_SEMIPUBLIC) HDR = $(HDR_PUBLIC) $(HDR_PRIVATE) diff --git a/TODO b/TODO @@ -4,6 +4,10 @@ It uses optimised division algorithm that requires that d|n. Add zsets_radix Add zstr_radix +Can zmodpowu and zmodpow be improved using some other algorithm? +Is it worth implementing precomputed optimal + addition-chain exponentiation in zpowu? + Test big endian Test always having .used > 0 for zero Test negative/non-negative instead of sign diff --git a/doc/arithmetic.tex b/doc/arithmetic.tex @@ -186,8 +186,9 @@ can be expressed as a simple formula \vspace{-1em} \[ \hspace*{-0.4cm} - a^b = \prod_{i = 0}^{\lceil \log_2 b \rceil} - \left ( a^{2^i} \right )^{\left \lfloor {\displaystyle{b \over 2^i}} \hspace*{-1ex} \mod 2 \right \rfloor} + a^b = + \prod_{k \in \textbf{Z}_{+} ~:~ \left \lfloor {b \over 2^k} \hspace*{-1ex} \mod 2 \right \rfloor = 1} + a^{2^k} \] \noindent diff --git a/doc/bit-operations.tex b/doc/bit-operations.tex @@ -0,0 +1,56 @@ +\chapter{Bit operations} +\label{chap:Bit operations} + +TODO + +\vspace{1cm} +\minitoc + + +\newpage +\section{Boundary} +\label{sec:Boundary} + +TODO % zbits zlsb + + +\newpage +\section{Logic} +\label{sec:Logic} + +TODO % zand zor zxor znot + + +\newpage +\section{Shift} +\label{sec:Shift} + +TODO % zlsh zrsh + + +\newpage +\section{Truncation} +\label{sec:Truncation} + +TODO % ztrunc + + +\newpage +\section{Split} +\label{sec:Split} + +TODO % zsplit + + +\newpage +\section{Bit manipulation} +\label{sec:Bit manipulation} + +TODO % zbset + + +\newpage +\section{Bit test} +\label{sec:Bit test} + +TODO % zbtest diff --git a/doc/libzahl.tex b/doc/libzahl.tex @@ -82,6 +82,10 @@ all copies or substantial portions of the Document. \input doc/get-started.tex \input doc/miscellaneous.tex \input doc/arithmetic.tex +\input doc/bit-operations.tex +\input doc/number-theory.tex +\input doc/random-numbers.tex +\input doc/not-implemented.tex \appendix diff --git a/doc/not-implemented.tex b/doc/not-implemented.tex @@ -0,0 +1,554 @@ +\chapter{Not implemented} +\label{chap:Not implemented} + +In this chapter we maintain a list of +features we have choosen not to implement, +but would fit into libzahl had we not have +our priorities straight. Functions listed +herein will only be implemented if there +is shown that it would be overwhelmingly +advantageous. + +\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} + + +\newpage +\section{Least common multiple} +\label{sec:Least common multiple} + +\( \displaystyle{ + \mbox{lcm}(a, b) = {\lvert a \cdot b \rvert \over \mbox{gcd}(a, b)} +}\) + + +\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, mabs, apos); + \} + 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} + +TODO + + +\subsection{Jacobi symbol} +\label{sec:Jacobi symbol} + +TODO + + +\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 = 0}^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. + + +\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 {(-1)^i \over 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)} = {(x + n - 1)! \over (x - 1)!} +}\), undefined for $n < 0$. + + +\subsection{Failing factorial} +\label{sec:Failing factorial} + +\( \displaystyle{ + (x)_n = {x! \over (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{ + {n \choose k} = {n! \over k!(n - k)!} + = {1 \over (n - k)!} \prod_{i = k + 1}^n i + = {1 \over k!} \prod_{i = n - k + 1}^n i +}\) + + +\subsection{Catalan number} +\label{sec:Catalan number} + +\( \displaystyle{ + C_n = \left . {2n \choose n} \middle / (n + 1) \right . +}\) + + +\subsection{Fuss–Catalan number} +\label{sec:Fuss-Catalan number} % not en dash + +\( \displaystyle{ + A_m(p, r) = {r \over mp + r} {mp + r \choose 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) <= 1) \{ + 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); + \} + + 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 speed 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 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 number +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 improve by comparing character by +character manually with using {\tt zxor}. + + +\subsection{Character retrieval} +\label{sec:Character retrieval} + +\begin{alltt} + uint64_t + getu(z_t a) + \{ + return zzero(a) ? 0 : a->chars[0]; + \} +\end{alltt} diff --git a/doc/number-theory.tex b/doc/number-theory.tex @@ -0,0 +1,35 @@ +\chapter{Number theory} +\label{chap:Number theory} + +TODO + +\vspace{1cm} +\minitoc + + +\newpage +\section{Odd or even} +\label{sec:Odd or even} + +TODO % zodd zeven zodd_nonzero zeven_nonzero + + +\newpage +\section{Signum} +\label{sec:Signum} + +TODO % zsignum zzero + + +\newpage +\section{Greatest common divisor} +\label{sec:Greatest common divisor} + +TODO % zgcd + + +\newpage +\section{Primality test} +\label{sec:Primality test} + +TODO % zptest diff --git a/doc/random-numbers.tex b/doc/random-numbers.tex @@ -0,0 +1,28 @@ +\chapter{Random numbers} +\label{chap:Random numbers} + +TODO + +\vspace{1cm} +\minitoc + + +\newpage +\section{Generation} +\label{sec:Generation} + +TODO + + +\newpage +\section{Devices} +\label{sec:Devices} + +TODO + + +\newpage +\section{Distributions} +\label{sec:Distributions} + +TODO diff --git a/src/zmodmul.c b/src/zmodmul.c @@ -6,6 +6,7 @@ void zmodmul(z_t a, z_t b, z_t c, z_t d) { /* TODO Montgomery modular multiplication */ + /* TODO Kochanski multiplication */ if (unlikely(a == d)) { zset(libzahl_tmp_modmul, d); zmul(a, b, c); diff --git a/src/zmodpow.c b/src/zmodpow.c @@ -12,6 +12,8 @@ zmodpow(z_t a, z_t b, z_t c, z_t d) size_t i, j, n, bits; zahl_char_t x; + /* TODO use zmodpowu when possible */ + if (unlikely(zsignum(c) <= 0)) { if (zzero(c)) { if (check(zzero(b))) diff --git a/src/zmodpowu.c b/src/zmodpowu.c @@ -25,10 +25,11 @@ zmodpowu(z_t a, z_t b, unsigned long long int c, z_t d) zmod(tb, b, d); zset(td, d); - zsetu(a, 1); if (c & 1) - zmodmul(a, a, tb, td); + zset(a, tb); + else + zsetu(a, 1); while (c >>= 1) { zmodsqr(tb, tb, td); if (c & 1) diff --git a/src/zpow.c b/src/zpow.c @@ -14,6 +14,8 @@ zpow(z_t a, z_t b, z_t c) * 7↑19 = 7↑10011₂ = 7↑2⁰ ⋅ 7↑2¹ ⋅ 7↑2⁴ where a↑2↑(n + 1) = (a↑2↑n)². */ + /* TODO use zpowu when possible */ + size_t i, j, n, bits; zahl_char_t x; int neg; diff --git a/src/zpowu.c b/src/zpowu.c @@ -21,10 +21,11 @@ zpowu(z_t a, z_t b, unsigned long long int c) neg = znegative(b) && (c & 1); zabs(tb, b); - zsetu(a, 1); if (c & 1) - zmul_ll(a, a, tb); + zset(a, tb); + else + zsetu(a, 1); while (c >>= 1) { zsqr_ll(tb, tb); if (c & 1)