arithmetic.tex - 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*… | |
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 = … | |
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} … | |
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} |