 # libzahl

big integer library
git clone git://git.suckless.org/libzahl

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
17
18 To calculate the sum of two terms, we perform
20
21 \vspace{1em}
22 $r \gets a + b$
23 \vspace{1em}
24
25 \noindent
26 is written as
27
28 \begin{alltt}
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}
45 \end{alltt}
46
47 \noindent
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
91 in-place operation:
92
93 \begin{alltt}
95    zadd(b, a, b);           \textcolor{c}{/* \textrm{should be avoided} */}
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
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}
`