libzahl

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

arithmetic.tex (15166B)


      1 \chapter{Arithmetic}
      2 \label{chap:Arithmetic}
      3 
      4 In this chapter, we will learn how to perform basic
      5 arithmetic with libzahl: addition, subtraction,
      6 multiplication, division, modulus, exponentiation,
      7 and sign manipulation. \secref{sec:Division} is of
      8 special importance.
      9 
     10 \vspace{1cm}
     11 \minitoc
     12 
     13 
     14 \newpage
     15 \section{Addition}
     16 \label{sec:Addition}
     17 
     18 To calculate the sum of two terms, we perform
     19 addition using {\tt zadd}.
     20 
     21 \vspace{1em}
     22 $r \gets a + b$
     23 \vspace{1em}
     24 
     25 \noindent
     26 is written as
     27 
     28 \begin{alltt}
     29    zadd(r, a, b);
     30 \end{alltt}
     31 
     32 libzahl also provides {\tt zadd\_unsigned} which
     33 has slightly lower overhead. The calculates the
     34 sum of the absolute values of two integers.
     35 
     36 \vspace{1em}
     37 $r \gets \lvert a \rvert + \lvert b \rvert$
     38 \vspace{1em}
     39 
     40 \noindent
     41 is written as
     42 
     43 \begin{alltt}
     44    zadd_unsigned(r, a, b);
     45 \end{alltt}
     46 
     47 \noindent
     48 {\tt zadd\_unsigned} has lower overhead than
     49 {\tt zadd} because it does not need to inspect
     50 or change the sign of the input, the low-level
     51 function that performs the addition inherently
     52 calculates the sum of the absolute values of
     53 the input.
     54 
     55 In libzahl, addition is implemented using a
     56 technique called ripple-carry. It is derived
     57 from that observation that
     58 
     59 \vspace{1em}
     60 $f : \textbf{Z}_n, \textbf{Z}_n \rightarrow \textbf{Z}_n$
     61 \\ \indent
     62 $f : a, b \mapsto a + b + 1$
     63 \vspace{1em}
     64 
     65 \noindent
     66 only wraps at most once, that is, the
     67 carry cannot exceed 1. CPU:s provide an
     68 instruction specifically for performing
     69 addition with ripple-carry over multiple words,
     70 adds twos numbers plus the carry from the
     71 last addition. libzahl uses assembly to
     72 implement this efficiently. If, however, an
     73 assembly implementation is not available for
     74 the on which machine it is running, libzahl
     75 implements ripple-carry less efficiently
     76 using compiler extensions that check for
     77 overflow. In the event that neither an
     78 assembly implementation is available nor
     79 the compiler is known to support this
     80 extension, it is implemented using inefficient
     81 pure C code. This last resort manually
     82 predicts whether an addition will overflow;
     83 this could be made more efficient, by never
     84 using the highest bit in each character,
     85 except to detect overflow. This optimisation
     86 is however not implemented because it is
     87 not deemed important enough and would
     88 be detrimental to libzahl's simplicity.
     89 
     90 {\tt zadd} and {\tt zadd\_unsigned} support
     91 in-place operation:
     92 
     93 \begin{alltt}
     94    zadd(a, a, b);
     95    zadd(b, a, b);           \textcolor{c}{/* \textrm{should be avoided} */}
     96    zadd_unsigned(a, a, b);
     97    zadd_unsigned(b, a, b);  \textcolor{c}{/* \textrm{should be avoided} */}
     98 \end{alltt}
     99 
    100 \noindent
    101 Use this whenever possible, it will improve
    102 your performance, as it will involve less
    103 CPU instructions for each character addition
    104 and it may be possible to eliminate some
    105 character additions.
    106 
    107 
    108 \newpage
    109 \section{Subtraction}
    110 \label{sec:Subtraction}
    111 
    112 TODO % zsub zsub_unsigned
    113 
    114 
    115 \newpage
    116 \section{Multiplication}
    117 \label{sec:Multiplication}
    118 
    119 TODO % zmul zmodmul
    120 
    121 
    122 \newpage
    123 \section{Division}
    124 \label{sec:Division}
    125 
    126 To calculate the quotient or modulus of two integers,
    127 use either of
    128 
    129 \begin{alltt}
    130    void zdiv(z_t quotient, z_t dividend, z_t divisor);
    131    void zmod(z_t remainder, z_t dividend, z_t divisor);
    132    void zdivmod(z_t quotient, z_t remainder,
    133                 z_t dividend, z_t divisor);
    134 \end{alltt}
    135 
    136 \noindent
    137 These function \emph{do not} allow {\tt NULL}
    138 for the output parameters: {\tt quotient} and
    139 {\tt remainder}. The quotient and remainder are
    140 calculated simultaneously and indivisibly, hence
    141 {\tt zdivmod} is provided to calculated both; if
    142 you are only interested in the quotient or only
    143 interested in the remainder, use {\tt zdiv} or
    144 {\tt zmod}, respectively.
    145 
    146 These functions calculate a truncated quotient.
    147 That is, the result is rounded towards zero. This
    148 means for example that if the quotient is in
    149 $(-1,~1)$, {\tt quotient} gets 0. That is, this % TODO try to clarify
    150 would not be the case for one of the sides of zero.
    151 For example, if the quotient would have been
    152 floored, negative quotients would have been rounded
    153 away from zero. libzahl only provides truncated
    154 division.
    155 
    156 The remainder is defined such that $n = qd + r$ after
    157 calling {\tt zdivmod(q, r, n, d)}. There is no
    158 difference in the remainer between {\tt zdivmod}
    159 and {\tt zmod}. The sign of {\tt d} has no affect
    160 on {\tt r}, {\tt r} will always, unless it is zero,
    161 have the same sign as {\tt n}.
    162 
    163 There are of course other ways to define integer
    164 division (that is, \textbf{Z} being the codomain)
    165 than as truncated division. For example integer
    166 divison in Python is floored — yes, you did just
    167 read `integer divison in Python is floored,' and
    168 you are correct, that is not the case in for
    169 example C. Users that want another definition
    170 for division than truncated division are required
    171 to implement that themselves. We will however
    172 lend you a hand.
    173 
    174 \begin{alltt}
    175    #define isneg(x) (zsignum(x) < 0)
    176    static z_t one;
    177    \textcolor{c}{__attribute__((constructor)) static}
    178    void init(void) \{ zinit(one), zseti(one, 1); \}
    179 
    180    static int
    181    cmpmag_2a_b(z_t a, z_b b)
    182    \{
    183        int r;
    184        zadd(a, a, a), r = zcmpmag(a, b), zrsh(a, a, 1);
    185        return r;
    186    \}
    187 \end{alltt}
    188 
    189 % Floored division
    190 \begin{alltt}
    191    void \textcolor{c}{/* \textrm{All arguments must be unique.} */}
    192    divmod_floor(z_t q, z_t r, z_t n, z_t d)
    193    \{
    194        zdivmod(q, r, n, d);
    195        if (!zzero(r) && isneg(n) != isneg(d))
    196            zsub(q, q, one), zadd(r, r, d);
    197    \}
    198 \end{alltt}
    199 
    200 % Ceiled division
    201 \begin{alltt}
    202    void \textcolor{c}{/* \textrm{All arguments must be unique.} */}
    203    divmod_ceiling(z_t q, z_t r, z_t n, z_t d)
    204    \{
    205        zdivmod(q, r, n, d);
    206        if (!zzero(r) && isneg(n) == isneg(d))
    207            zadd(q, q, one), zsub(r, r, d);
    208    \}
    209 \end{alltt}
    210 
    211 % Division with round half aways from zero
    212 % This rounding method is also called:
    213 %    round half toward infinity
    214 %    commercial rounding
    215 \begin{alltt}
    216    /* \textrm{This is how we normally round numbers.} */
    217    void \textcolor{c}{/* \textrm{All arguments must be unique.} */}
    218    divmod_half_from_zero(z_t q, z_t r, z_t n, z_t d)
    219    \{
    220        zdivmod(q, r, n, d);
    221        if (!zzero(r) && cmpmag_2a_b(r, d) >= 0) \{
    222            if (isneg(n) == isneg(d))
    223                zadd(q, q, one), zsub(r, r, d);
    224            else
    225                zsub(q, q, one), zadd(r, r, d);
    226        \}
    227    \}
    228 \end{alltt}
    229 
    230 \noindent
    231 Now to the weird ones that will more often than
    232 not award you a face-slap. % Had positive punishment
    233 % been legal or even mildly pedagogical. But I would
    234 % not put it past Coca-Cola.
    235 
    236 % Division with round half toward zero
    237 % This rounding method is also called:
    238 %     round half away from infinity
    239 \begin{alltt}
    240    void \textcolor{c}{/* \textrm{All arguments must be unique.} */}
    241    divmod_half_to_zero(z_t q, z_t r, z_t n, z_t d)
    242    \{
    243        zdivmod(q, r, n, d);
    244        if (!zzero(r) && cmpmag_2a_b(r, d) > 0) \{
    245            if (isneg(n) == isneg(d))
    246                zadd(q, q, one), zsub(r, r, d);
    247            else
    248                zsub(q, q, one), zadd(r, r, d);
    249        \}
    250    \}
    251 \end{alltt}
    252 
    253 % Division with round half up
    254 % This rounding method is also called:
    255 %     round half towards positive infinity
    256 \begin{alltt}
    257    void \textcolor{c}{/* \textrm{All arguments must be unique.} */}
    258    divmod_half_up(z_t q, z_t r, z_t n, z_t d)
    259    \{
    260        int cmp;
    261        zdivmod(q, r, n, d);
    262        if (!zzero(r) && (cmp = cmpmag_2a_b(r, d)) >= 0) \{
    263            if (isneg(n) == isneg(d))
    264                zadd(q, q, one), zsub(r, r, d);
    265            else if (cmp)
    266                zsub(q, q, one), zadd(r, r, d);
    267        \}
    268    \}
    269 \end{alltt}
    270 
    271 % Division with round half down
    272 % This rounding method is also called:
    273 %     round half towards negative infinity
    274 \begin{alltt}
    275    void \textcolor{c}{/* \textrm{All arguments must be unique.} */}
    276    divmod_half_down(z_t q, z_t r, z_t n, z_t d)
    277    \{
    278        int cmp;
    279        zdivmod(q, r, n, d);
    280        if (!zzero(r) && (cmp = cmpmag_2a_b(r, d)) >= 0) \{
    281            if (isneg(n) != isneg(d))
    282                zsub(q, q, one), zadd(r, r, d);
    283            else if (cmp)
    284                zadd(q, q, one), zsub(r, r, d);
    285        \}
    286    \}
    287 \end{alltt}
    288 
    289 % Division with round half to even
    290 % This rounding method is also called:
    291 %     unbiased rounding        (really stupid name)
    292 %     convergent rounding      (also quite stupid name)
    293 %     statistician's rounding
    294 %     Dutch rounding
    295 %     Gaussian rounding
    296 %     odd–even rounding
    297 %     bankers' rounding
    298 % It is the default rounding method used in IEEE 754.
    299 \begin{alltt}
    300    void \textcolor{c}{/* \textrm{All arguments must be unique.} */}
    301    divmod_half_to_even(z_t q, z_t r, z_t n, z_t d)
    302    \{
    303        int cmp;
    304        zdivmod(q, r, n, d);
    305        if (!zzero(r) && (cmp = cmpmag_2a_b(r, d)) >= 0) \{
    306            if (cmp || zodd(q)) \{
    307                if (isneg(n) != isneg(d))
    308                    zsub(q, q, one), zadd(r, r, d);
    309                else
    310                    zadd(q, q, one), zsub(r, r, d);
    311            \}
    312        \}
    313    \}
    314 \end{alltt}
    315 
    316 % Division with round half to odd
    317 \newpage
    318 \begin{alltt}
    319    void \textcolor{c}{/* \textrm{All arguments must be unique.} */}
    320    divmod_half_to_odd(z_t q, z_t r, z_t n, z_t d)
    321    \{
    322        int cmp;
    323        zdivmod(q, r, n, d);
    324        if (!zzero(r) && (cmp = cmpmag_2a_b(r, d)) >= 0) \{
    325            if (cmp || zeven(q)) \{
    326                if (isneg(n) != isneg(d))
    327                    zsub(q, q, one), zadd(r, r, d);
    328                else
    329                    zadd(q, q, one), zsub(r, r, d);
    330            \}
    331        \}
    332    \}
    333 \end{alltt}
    334 
    335 % Other standard methods include stochastic rounding
    336 % and round half alternatingly, and what is is
    337 % New Zealand called “Swedish rounding”, which is
    338 % no longer used in Sweden, and is just normal round
    339 % half aways from zero but with 0.5 rather than
    340 % 1 as the integral unit, and is just a special case
    341 % of a more general rounding method.
    342 
    343 Currently, libzahl uses an almost trivial division
    344 algorithm. It operates on positive numbers. It begins
    345 by left-shifting the divisor as much as possible with
    346 letting it exceed the dividend. Then, it subtracts
    347 the shifted divisor from the dividend and add 1,
    348 left-shifted as much as the divisor, to the quotient.
    349 The quotient begins at 0. It then right-shifts
    350 the shifted divisor as little as possible until
    351 it no longer exceeds the diminished dividend and
    352 marks the shift in the quotient. This process is
    353 repeated until the unshifted divisor is greater
    354 than the diminished dividend. The final diminished
    355 dividend is the remainder.
    356 
    357 
    358 
    359 \newpage
    360 \section{Exponentiation}
    361 \label{sec:Exponentiation}
    362 
    363 Exponentiation refers to raising a number to
    364 a power. libzahl provides two functions for
    365 regular exponentiation, and two functions for
    366 modular exponentiation. libzahl also provides
    367 a function for raising a number to the second
    368 power, see \secref{sec:Multiplication} for
    369 more details on this. The functions for regular
    370 exponentiation are
    371 
    372 \begin{alltt}
    373    void zpow(z_t power, z_t base, z_t exponent);
    374    void zpowu(z_t, z_t, unsigned long long int);
    375 \end{alltt}
    376 
    377 \noindent
    378 They are identical, except {\tt zpowu} expects
    379 an intrinsic type as the exponent. Both functions
    380 calculate
    381 
    382 \vspace{1em}
    383 $power \gets base^{exponent}$
    384 \vspace{1em}
    385 
    386 \noindent
    387 The functions for modular exponentiation are
    388 \begin{alltt}
    389    void zmodpow(z_t, z_t, z_t, z_t modulator);
    390    void zmodpowu(z_t, z_t, unsigned long long int, z_t);
    391 \end{alltt}
    392 
    393 \noindent
    394 They are identical, except {\tt zmodpowu} expects
    395 and intrinsic type as the exponent. Both functions
    396 calculate
    397 
    398 \vspace{1em}
    399 $power \gets base^{exponent} \mod modulator$
    400 \vspace{1em}
    401 
    402 The sign of {\tt modulator} does not affect the
    403 result, {\tt power} will be negative if and only
    404 if {\tt base} is negative and {\tt exponent} is
    405 odd, that is, under the same circumstances as for
    406 {\tt zpow} and {\tt zpowu}.
    407 
    408 These four functions are implemented using
    409 exponentiation by squaring. {\tt zmodpow} and
    410 {\tt zmodpowu} are optimised, they modulate
    411 results for multiplication and squaring at
    412 every multiplication and squaring, rather than
    413 modulating the result at the end. Exponentiation
    414 by modulation is a very simple algorithm which
    415 can be expressed as a simple formula
    416 
    417 \vspace{-1em}
    418 \[ \hspace*{-0.4cm}
    419     a^b =
    420     \prod_{k \in \textbf{Z}_{+} ~:~ \left \lfloor \frac{b}{2^k} \hspace*{-1ex} \mod 2 \right \rfloor = 1}
    421     a^{2^k}
    422 \]
    423 
    424 \noindent
    425 This is a natural extension to the
    426 observations\footnote{The first of course being
    427 that any non-negative number can be expressed
    428 with the binary positional system. The latter
    429 should be fairly self-explanatory.}
    430 
    431 \vspace{-1em}
    432 \[ \hspace*{-0.4cm}
    433     \forall b \in \textbf{Z}_{+} \exists B \subset \textbf{Z}_{+} : b = \sum_{i \in B} 2^i
    434     ~~~~ \textrm{and} ~~~~
    435     a^{\sum x} = \prod a^x.
    436 \]
    437 
    438 \noindent
    439 The algorithm can be expressed in psuedocode as
    440 
    441 \vspace{1em}
    442 \hspace{-2.8ex}
    443 \begin{minipage}{\linewidth}
    444 \begin{algorithmic}
    445     \STATE $r, f \gets 1, a$
    446     \WHILE{$b \neq 0$}
    447       \STATE $r \gets r \cdot f$ {\bf unless} $2 \vert b$
    448       \STATE $f \gets f^2$ \textcolor{c}{\{$f \gets f \cdot f$\}}
    449       \STATE $b \gets \lfloor b / 2 \rfloor$
    450     \ENDWHILE
    451     \RETURN $r$ 
    452 \end{algorithmic}
    453 \end{minipage}
    454 \vspace{1em}
    455 
    456 \noindent
    457 Modular exponentiation ($a^b \mod m$) by squaring can be
    458 expressed as
    459 
    460 \vspace{1em}
    461 \hspace{-2.8ex}
    462 \begin{minipage}{\linewidth}
    463 \begin{algorithmic}
    464     \STATE $r, f \gets 1, a$
    465     \WHILE{$b \neq 0$}
    466       \STATE $r \gets r \cdot f \hspace*{-1ex}~ \mod m$ \textbf{unless} $2 \vert b$
    467       \STATE $f \gets f^2 \hspace*{-1ex}~ \mod m$
    468       \STATE $b \gets \lfloor b / 2 \rfloor$
    469     \ENDWHILE
    470     \RETURN $r$ 
    471 \end{algorithmic}
    472 \end{minipage}
    473 \vspace{1em}
    474 
    475 {\tt zmodpow} does \emph{not} calculate the
    476 modular inverse if the exponent is negative,
    477 rather, you should expect the result to be
    478 1 and 0 depending of whether the base is 1
    479 or not 1.
    480 
    481 
    482 \newpage
    483 \section{Sign manipulation}
    484 \label{sec:Sign manipulation}
    485 
    486 libzahl provides two functions for manipulating
    487 the sign of integers:
    488 
    489 \begin{alltt}
    490    void zabs(z_t r, z_t a);
    491    void zneg(z_t r, z_t a);
    492 \end{alltt}
    493 
    494 {\tt zabs} stores the absolute value of {\tt a}
    495 in {\tt r}, that is, it creates a copy of
    496 {\tt a} to {\tt r}, unless {\tt a} and {\tt r}
    497 are the same reference, and then removes its sign;
    498 if the value is negative, it becomes positive.
    499 
    500 \vspace{1em}
    501 \(
    502     r \gets \lvert a \rvert =
    503     \left \lbrace \begin{array}{rl}
    504         -a & \quad \textrm{if}~a \le 0 \\
    505         +a & \quad \textrm{if}~a \ge 0 \\
    506     \end{array} \right .
    507 \)
    508 \vspace{1em}
    509 
    510 {\tt zneg} stores the negated of {\tt a}
    511 in {\tt r}, that is, it creates a copy of
    512 {\tt a} to {\tt r}, unless {\tt a} and {\tt r}
    513 are the same reference, and then flips sign;
    514 if the value is negative, it becomes positive,
    515 if the value is positive, it becomes negative.
    516 
    517 \vspace{1em}
    518 \(
    519     r \gets -a
    520 \)
    521 \vspace{1em}
    522 
    523 Note that there is no function for
    524 
    525 \vspace{1em}
    526 \(
    527     r \gets -\lvert a \rvert =
    528     \left \lbrace \begin{array}{rl}
    529          a & \quad \textrm{if}~a \le 0 \\
    530         -a & \quad \textrm{if}~a \ge 0 \\
    531     \end{array} \right .
    532 \)
    533 \vspace{1em}
    534 
    535 \noindent
    536 calling {\tt zabs} followed by {\tt zneg}
    537 should be sufficient for most users:
    538 
    539 \begin{alltt}
    540    #define my_negabs(r, a)  (zabs(r, a), zneg(r, r))
    541 \end{alltt}