arithmetic.tex - libzahl - big integer library
HTML git clone git://git.suckless.org/libzahl
DIR Log
DIR Files
DIR Refs
DIR README
DIR 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}