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}