exercises.tex (21757B)
1 \chapter{Exercises} 2 \label{chap:Exercises} 3 4 % TODO 5 % 6 % All exercises should be group with a chapter 7 % 8 % ▶ Recommended 9 % M Matematically oriented 10 % HM Higher matematical education required 11 % MP Matematical problem solving skills required, 12 % double the rating otherwise; difficult is 13 % very personal 14 % 15 % 00 Immediate 16 % 10 Simple 17 % 20 Medium 18 % 30 Moderately hard 19 % 40 Term project 20 % 50 Research project 21 % 22 % ⁺ High risk of undershoot difficulty 23 24 25 \begin{enumerate}[label=\textbf{\arabic*}.] 26 27 28 29 \item {[$\RHD$\textit{02}]} \textbf{Saturated subtraction} 30 31 Implement the function 32 33 \vspace{-1em} 34 \begin{alltt} 35 void monus(z_t r, z_t a, z_t b); 36 \end{alltt} 37 \vspace{-1em} 38 39 \noindent 40 which calculates $r = a \dotminus b = \max \{ 0,~ a - b \}$. 41 42 43 44 \item {[$\RHD$\textit{10}]} \textbf{Modular powers of 2} 45 46 What is the advantage of using \texttt{zmodpow} 47 over \texttt{zbset} or \texttt{zlsh} in combination 48 with \texttt{zmod}? 49 50 51 52 \item {[\textit{M15}]} \textbf{Convergence of the Lucas Number ratios} 53 54 Find an approximation for 55 $\displaystyle{ \lim_{n \to \infty} \frac{L_{n + 1}}{L_n}}$, 56 where $L_n$ is the $n^{\text{th}}$ 57 Lucas Number \psecref{sec:Lucas numbers}. 58 59 \( \displaystyle{ 60 L_n \stackrel{\text{\tiny{def}}}{\text{=}} \left \{ \begin{array}{ll} 61 2 & \text{if} ~ n = 0 \\ 62 1 & \text{if} ~ n = 1 \\ 63 L_{n - 1} + L_{n + 1} & \text{otherwise} 64 \end{array} \right . 65 }\) 66 67 68 69 \item {[\textit{MP12}]} \textbf{Factorisation of factorials} 70 71 Implement the function 72 73 \vspace{-1em} 74 \begin{alltt} 75 void factor_fact(z_t n); 76 \end{alltt} 77 \vspace{-1em} 78 79 \noindent 80 which prints the prime factorisation of $n!$ 81 (the $n^{\text{th}}$ factorial). The function shall 82 be efficient for all $n$ where all primes $p \le n$ 83 can be found efficiently. You can assume that 84 $n \ge 2$. You should not evaluate $n!$. 85 86 87 88 \item {[\textit{M20}]} \textbf{Reverse factorisation of factorials} 89 90 {\small\textit{You should already have solved 91 ``Factorisation of factorials'' before you solve 92 this problem.}} 93 94 Implement the function 95 96 \vspace{-1em} 97 \begin{alltt} 98 void unfactor_fact(z_t x, z_t *P, 99 unsigned long long int *K, size_t n); 100 \end{alltt} 101 \vspace{-1em} 102 103 \noindent 104 which given the factorsation of $x!$ determines $x$. 105 The factorisation of $x!$ is 106 $\displaystyle{\prod_{i = 1}^{n} P_i^{K_i}}$, where 107 $P_i$ is \texttt{P[i - 1]} and $K_i$ is \texttt{K[i - 1]}. 108 109 110 111 \item {[$\RHD$\textit{MP17}]} \textbf{Factorials inverted} 112 113 Implement the function 114 115 \vspace{-1em} 116 \begin{alltt} 117 void unfact(z_t x, z_t n); 118 \end{alltt} 119 \vspace{-1em} 120 121 \noindent 122 which given a factorial number $n$, i.e. on the form 123 $x! = 1 \cdot 2 \cdot 3 \cdot \ldots \cdot x$, 124 calculates $x = n!^{-1}$. You can assume that 125 $n$ is a perfect factorial number and that $x \ge 1$. 126 Extra credit if you can detect when the input, $n$, 127 is not a factorial number. Such function would of 128 course return an \texttt{int} 1 if the input is a 129 factorial and 0 otherwise, or alternatively 0 130 on success and $-1$ with \texttt{errno} set to 131 \texttt{EDOM} if the input is not a factorial. 132 133 134 135 \item {[\textit{05}]} \textbf{Fast primality test} 136 137 $(x + y)^p \equiv x^p + y^p ~(\text{Mod}~p)$ 138 for all primes $p$ and for a few composites $p$, 139 which are know as pseudoprimes. Use this to implement 140 a fast primality tester. 141 142 143 144 \item {[\textit{10}]} \textbf{Fermat primality test} 145 146 $a^{p - 1} \equiv 1 ~(\text{Mod}~p) ~\forall~ 1 < a < p$ 147 for all primes $p$ and for a few composites $p$, 148 which are know as pseudoprimes\footnote{If $p$ is composite 149 but passes the test for all $a$, $p$ is a Carmichael 150 number.}. Use this to implement a heuristic primality 151 tester. Try to mimic \texttt{zptest} as much as possible. 152 GNU~MP uses $a = 210$, but you don't have to. ($a$ is 153 called a base.) 154 155 156 157 \item {[\textit{11}]} \textbf{Lucas–Lehmer primality test} 158 159 The Lucas–Lehmer primality test can be used to determine 160 whether a Mersenne numbers $M_n = 2^n - 1$ is a prime (a 161 Mersenne prime). $M_n$ is a prime if and only if 162 $s_{n - 1} \equiv 0 ~(\text{Mod}~M_n)$, where 163 164 \( \displaystyle{ 165 s_i = \left \{ \begin{array}{ll} 166 4 & \text{if} ~ i = 0 \\ 167 s_{i - 1}^2 - 2 & \text{otherwise}. 168 \end{array} \right . 169 }\) 170 171 \noindent 172 The Lucas–Lehmer primality test requires that $n \ge 3$, 173 however $M_2 = 2^2 - 1 = 3$ is a prime. Implement a version 174 of the Lucas–Lehmer primality test that takes $n$ as the 175 input. For some more fun, when you are done, you can 176 implement a version that takes $M_n$ as the input. 177 178 For improved performance, instead of using \texttt{zmod}, 179 you can use the recursive function 180 % 181 \( \displaystyle{ 182 k \text{ mod } (2^n - 1) = 183 \left ( 184 (k \text{ mod } 2^n) + \lfloor k \div 2^n \rfloor 185 \right ) \text{ mod } (2^n - 1), 186 }\) 187 % 188 where $k \mod 2^n$ is efficiently calculated 189 using \texttt{zand($k$, $2^n - 1$)}. (This optimisation 190 is not part of the difficulty rating of this problem.) 191 192 193 194 \item {[\textit{20}]} \textbf{Fast primality test with bounded perfection} 195 196 Implement a primality test that is both very fast and 197 never returns \texttt{PROBABLY\_PRIME} for input less 198 than or equal to a preselected number. 199 200 201 202 \item {[\textit{30}]} \textbf{Powers of the golden ratio} 203 204 Implement function that returns $\varphi^n$ rounded 205 to the nearest integer, where $n$ is the input and 206 $\varphi$ is the golden ratio. 207 208 209 210 \item {[\textit{$\RHD$05}]} \textbf{zlshu and zrshu} 211 212 Why does libzahl have 213 214 \vspace{-1em} 215 \begin{alltt} 216 void zlsh(z_t, z_t, size_t); 217 void zrsh(z_t, z_t, size_t); 218 \end{alltt} 219 \vspace{-1em} 220 221 \noindent 222 rather than 223 224 \vspace{-1em} 225 \begin{alltt} 226 void zlsh(z_t, z_t, z_t); 227 void zrsh(z_t, z_t, z_t); 228 void zlshu(z_t, z_t, size_t); 229 void zrshu(z_t, z_t, size_t); 230 \end{alltt} 231 \vspace{-1em} 232 233 234 235 \item {[\textit{$\RHD$M15}]} \textbf{Modular left-shift} 236 237 Implement a function that calculates 238 $2^a \text{ mod } b$, using \texttt{zmod} and 239 only cheap functions. You can assume $a \ge 0$, 240 $b \ge 1$. You can also assume that all 241 parameters are unique pointers. 242 243 244 245 \item {[\textit{$\RHD$08}]} \textbf{Modular left-shift, extended} 246 247 {\small\textit{You should already have solved 248 ``Modular left-shift'' before you solve this 249 problem.}} 250 251 Extend the function you wrote in ``Modular left-shift'' 252 to accept negative $b$ and non-unique pointers. 253 254 255 256 \item {[\textit{13}]} \textbf{The totient} 257 258 The totient of $n$ is the number of integer $a$, 259 $0 < a < n$ that are relatively prime to $n$. 260 Implement Euler's totient function $\varphi(n)$ 261 which calculates the totient of $n$. Its 262 formula is 263 264 \( \displaystyle{ 265 \varphi(n) = |n| \prod_{p \in \textbf{P} : p | n} 266 \left ( 1 - \frac{1}{p} \right ). 267 }\) 268 269 Note that $\varphi(-n) = \varphi(n)$, $\varphi(0) = 0$, 270 and $\varphi(1) = 1$. 271 272 273 274 \item {[\textit{M13}]} \textbf{The totient from factorisation} 275 276 Implement the function 277 278 \vspace{-1em} 279 \begin{alltt} 280 void totient_fact(z_t t, z_t *P, 281 unsigned long long int *K, size_t n); 282 \end{alltt} 283 \vspace{-1em} 284 285 \noindent 286 which calculates the totient $t = \varphi(n)$, where 287 $n = \displaystyle{\prod_{i = 1}^n P_i^{K_i}} > 0$, 288 and $P_i = \texttt{P[i - 1]} \in \textbf{P}$, 289 $K_i = \texttt{K[i - 1]} \ge 1$. All values \texttt{P} 290 are mutually unique. \texttt{P} and \texttt{K} make up 291 the prime factorisation of $n$. 292 293 You can use the following rules: 294 295 \( \displaystyle{ 296 \begin{array}{ll} 297 \varphi(1) = 1 & \\ 298 \varphi(p) = p - 1 & \text{if } p \in \textbf{P} \\ 299 \varphi(nm) = \varphi(n)\varphi(m) & \text{if } \gcd(n, m) = 1 \\ 300 n^a\varphi(n) = \varphi(n^{a + 1}) & 301 \end{array} 302 }\) 303 304 305 306 \item {[\textit{HMP32}]} \textbf{Modular tetration} 307 308 Implement the function 309 310 \vspace{-1em} 311 \begin{alltt} 312 void modtetra(z_t r, z_t b, unsigned long n, z_t m); 313 \end{alltt} 314 \vspace{-1em} 315 316 \noindent 317 which calculates $r = {}^n{}b \text{ mod } m$, where 318 ${}^0{}b = 1$, ${}^1{}b = b$, ${}^2{}b = b^b$, 319 ${}^3{}b = b^{b^b}$, ${}^4{}b = b^{b^{b^b}}$, and so on. 320 You can assume $b > 0$ and $m > 0$. You can also assume 321 \texttt{r}, \texttt{b}, and \texttt{m} are mutually 322 unique pointers. 323 324 325 326 \item {[\textit{13}]} \textbf{Modular generalised power towers} 327 328 {\small\textit{This problem requires a working 329 solution for ``Modular tetration''.}} 330 331 Modify your solution for ``Modular tetration'' to 332 evaluate any expression on the forms 333 $a^b,~a^{b^c},~a^{b^{c^d}},~\ldots \text{ mod } m$. 334 335 336 337 \end{enumerate} 338 339 340 341 \chapter{Solutions} 342 \label{chap:Solutions} 343 344 345 \begin{enumerate}[label=\textbf{\arabic*}.] 346 347 \item \textbf{Saturated subtraction} 348 349 \vspace{-1em} 350 \begin{alltt} 351 void monus(z_t r, z_t a, z_t b) 352 \{ 353 zsub(r, a, b); 354 if (zsignum(r) < 0) 355 zsetu(r, 0); 356 \} 357 \end{alltt} 358 \vspace{-1em} 359 360 361 \item \textbf{Modular powers of 2} 362 363 \texttt{zbset} and \texttt{zbit} requires $\Theta(n)$ 364 memory to calculate $2^n$. \texttt{zmodpow} only 365 requires $\mathcal{O}(\min \{n, \log m\})$ memory 366 to calculate $2^n \text{ mod } m$. $\Theta(n)$ 367 memory complexity becomes problematic for very 368 large $n$. 369 370 371 \item \textbf{Convergence of the Lucas Number ratios} 372 373 It would be a mistake to use bignum, and bigint in particular, 374 to solve this problem. Good old mathematics is a much better solution. 375 376 $$ \lim_{n \to \infty} \frac{L_{n + 1}}{L_n} = \lim_{n \to \infty} \frac{L_{n}}{L_{n - 1}} = \lim_{n \to \infty} \frac{L_{n - 1}}{L_{n - 2}} $$ 377 378 $$ \frac{L_{n}}{L_{n - 1}} = \frac{L_{n - 1}}{L_{n - 2}} $$ 379 380 $$ \frac{L_{n - 1} + L_{n - 2}}{L_{n - 1}} = \frac{L_{n - 1}}{L_{n - 2}} $$ 381 382 $$ 1 + \frac{L_{n - 2}}{L_{n - 1}} = \frac{L_{n - 1}}{L_{n - 2}} $$ 383 384 $$ 1 + \varphi = \frac{1}{\varphi} $$ 385 386 So the ratio tends toward the golden ratio. 387 388 389 390 \item \textbf{Factorisation of factorials} 391 392 Base your implementation on 393 394 \( \displaystyle{ 395 n! = \prod_{p~\in~\textbf{P}}^{n} p^{k_p}, ~\text{where}~ 396 k_p = \sum_{a = 1}^{\lfloor \log_p n \rfloor} \lfloor np^{-a} \rfloor. 397 }\) 398 399 There is no need to calculate $\lfloor \log_p n \rfloor$, 400 you will see when $p^a > n$. 401 402 403 404 \item \textbf{Reverse factorisation of factorials} 405 406 $\displaystyle{x = \max_{p ~\in~ P} ~ p \cdot f(p, k_p)}$, 407 where $k_p$ is the power of $p$ in the factorisation 408 of $x!$. $f(p, k)$ is defined as: 409 410 \vspace{1em} 411 \hspace{-2.8ex} 412 \begin{minipage}{\linewidth} 413 \begin{algorithmic} 414 \STATE $k^\prime \gets 0$ 415 \WHILE{$k > 0$} 416 \STATE $a \gets 0$ 417 \WHILE{$p^a \le k$} 418 \STATE $k \gets k - p^a$ 419 \STATE $a \gets a + 1$ 420 \ENDWHILE 421 \STATE $k^\prime \gets k^\prime + p^{a - 1}$ 422 \ENDWHILE 423 \RETURN $k^\prime$ 424 \end{algorithmic} 425 \end{minipage} 426 \vspace{1em} 427 428 429 430 \item \textbf{Factorials inverted} 431 432 Use \texttt{zlsb} to get the power of 2 in the 433 prime factorisation of $n$, that is, the number 434 of times $n$ is divisible by 2. If we write $n$ on 435 the form $1 \cdot 2 \cdot 3 \cdot \ldots \cdot x$, 436 every $2^\text{nd}$ factor is divisible by 2, every 437 $4^\text{th}$ factor is divisible by $2^2$, and so on. 438 From calling \texttt{zlsb} we know how many times, 439 $n$ is divisible by 2, but know how many of the factors 440 are divisible by 2, but this can be calculated with 441 the following algorithm, where $k$ is the number 442 times $n$ is divisible by 2: 443 444 \vspace{1em} 445 \hspace{-2.8ex} 446 \begin{minipage}{\linewidth} 447 \begin{algorithmic} 448 \STATE $k^\prime \gets 0$ 449 \WHILE{$k > 0$} 450 \STATE $a \gets 0$ 451 \WHILE{$2^a \le k$} 452 \STATE $k \gets k - 2^a$ 453 \STATE $a \gets a + 1$ 454 \ENDWHILE 455 \STATE $k^\prime \gets k^\prime + 2^{a - 1}$ 456 \ENDWHILE 457 \RETURN $k^\prime$ 458 \end{algorithmic} 459 \end{minipage} 460 \vspace{1em} 461 462 \noindent 463 Note that $2^a$ is efficiently calculated with, 464 \texttt{zlsh}, but it is more efficient to use 465 \texttt{zbset}. 466 467 Now that we know $k^\prime$, the number of 468 factors ni $1 \cdot \ldots \cdot x$ that are 469 divisible by 2, we have two choices for $x$: 470 $k^\prime$ and $k^\prime + 1$. To check which, we 471 calculate $(k^\prime - 1)!!$ (the semifactoral, i.e. 472 $1 \cdot 3 \cdot 5 \cdot \ldots \cdot (k^\prime - 1)$) 473 naïvely and shift the result $k$ steps to the left. 474 This gives us $k^\prime!$. If $x! \neq k^\prime!$, then 475 $x = k^\prime + 1$ unless $n$ is not factorial number. 476 Of course, if $x! = k^\prime!$, then $x = k^\prime$. 477 478 479 480 \item \textbf{Fast primality test} 481 482 If we select $x = y = 1$ we get 483 $2^p \equiv 2 ~(\text{Mod}~p)$. This gives us 484 485 \vspace{-1em} 486 \begin{alltt} 487 enum zprimality 488 ptest_fast(z_t p) 489 \{ 490 z_t a; 491 int c = zcmpu(p, 2); 492 if (c <= 0) 493 return c ? NONPRIME : PRIME; 494 zinit(a); 495 zsetu(a, 1); 496 zlsh(a, a, p); 497 zmod(a, a, p); 498 c = zcmpu(a, 2); 499 zfree(a); 500 return c ? NONPRIME : PROBABLY_PRIME; 501 \} 502 \end{alltt} 503 \vspace{-1em} 504 505 506 507 \item \textbf{Fermat primality test} 508 509 Below is a simple implementation. It can be improved by 510 just testing against a fix base, such as $a = 210$, it 511 $t = 0$. It could also do an exhaustive test (all $a$ 512 such that $1 < a < p$) if $t < 0$. 513 514 \vspace{-1em} 515 \begin{alltt} 516 enum zprimality 517 ptest_fermat(z_t witness, z_t p, int t) 518 \{ 519 enum zprimality rc = PROBABLY_PRIME; 520 z_t a, p1, p3, temp; 521 int c; 522 523 if ((c = zcmpu(p, 2)) <= 0) \{ 524 if (!c) 525 return PRIME; 526 if (witness && witness != p) 527 zset(witness, p); 528 return NONPRIME; 529 \} 530 531 zinit(a), zinit(p1), zinit(p3), zinit(temp); 532 zsetu(temp, 3), zsub(p3, p, temp); 533 zsetu(temp, 1), zsub(p1, p, temp); 534 535 zsetu(temp, 2); 536 while (t--) \{ 537 zrand(a, DEFAULT_RANDOM, UNIFORM, p3); 538 zadd(a, a, temp); 539 zmodpow(a, a, p1, p); 540 if (zcmpu(a, 1)) \{ 541 if (witness) 542 zswap(witness, a); 543 rc = NONPRIME; 544 break; 545 \} 546 \} 547 548 zfree(temp), zfree(p3), zfree(p1), zfree(a); 549 return rc; 550 \} 551 \end{alltt} 552 \vspace{-1em} 553 554 555 556 \item \textbf{Lucas–Lehmer primality test} 557 558 \vspace{-1em} 559 \begin{alltt} 560 enum zprimality 561 ptest_llt(z_t n) 562 \{ 563 z_t s, M; 564 int c; 565 size_t p; 566 567 if ((c = zcmpu(n, 2)) <= 0) 568 return c ? NONPRIME : PRIME; 569 570 if (n->used > 1) \{ 571 \textcolor{c}{/* \textrm{An optimised implementation would not need this} */} 572 errno = ENOMEM; 573 return (enum zprimality)(-1); 574 \} 575 576 zinit(s), zinit(M), zinit(2); 577 578 p = (size_t)(n->chars[0]); 579 zsetu(s, 1), zsetu(M, 0); 580 zbset(M, M, p, 1), zsub(M, M, s); 581 zsetu(s, 4); 582 zsetu(two, 2); 583 584 p -= 2; 585 while (p--) \{ 586 zsqr(s, s); 587 zsub(s, s, two); 588 zmod(s, s, M); 589 \} 590 c = zzero(s); 591 592 zfree(two), zfree(M), zfree(s); 593 return c ? PRIME : NONPRIME; 594 \} 595 \end{alltt} 596 597 $M_n$ is composite if $n$ is composite, therefore, 598 if you do not expect prime-only values on $n$, the 599 performance can be improved by using some other 600 primality test (or this same test if $n$ is a 601 Mersenne number) to first check that $n$ is prime. 602 603 604 605 \item \textbf{Fast primality test with bounded perfection} 606 607 First we select a fast primality test. We can use 608 $2^p \equiv 2 ~(\texttt{Mod}~ p) ~\forall~ p \in \textbf{P}$, 609 as describe in the solution for the problem 610 \textit{Fast primality test}. 611 612 Next, we use this to generate a large list of primes and 613 pseudoprimes. Use a perfect primality test, such as a 614 naïve test or the AKS primality test, to filter out all 615 primes and retain only the pseudoprimes. This is not in 616 runtime so it does not matter that this is slow, but to 617 speed it up, we can use a probabilistic test such the 618 Miller–Rabin primality test (\texttt{zptest}) before we 619 use the perfect test. 620 621 Now that we have a quite large — but not humongous — list 622 of pseudoprimes, we can incorporate it into our fast 623 primality test. For any input that passes the test, and 624 is less or equal to the largest pseudoprime we found, 625 binary search our list of pseudoprime for the input. 626 627 For input, larger than our limit, that passes the test, 628 we can run it through \texttt{zptest} to reduce the 629 number of false positives. 630 631 As an alternative solution, instead of comparing against 632 known pseudoprimes. Find a minimal set of primes that 633 includes divisors for all known pseudoprimes, and do 634 trail division with these primes when a number passes 635 the test. No pseudoprime need to have more than one divisor 636 included in the set, so this set cannot be larger than 637 the set of pseudoprimes. 638 639 640 641 \item \textbf{Powers of the golden ratio} 642 643 This was an information gathering exercise. 644 For $n \ge 2$, $L_n = [\varphi^n]$, where 645 $L_n$ is the $n^\text{th}$ Lucas number. 646 647 \( \displaystyle{ 648 L_n \stackrel{\text{\tiny{def}}}{\text{=}} \left \{ \begin{array}{ll} 649 2 & \text{if} ~ n = 0 \\ 650 1 & \text{if} ~ n = 1 \\ 651 L_{n - 1} + L_{n + 1} & \text{otherwise} 652 \end{array} \right . 653 }\) 654 655 \noindent 656 but for efficiency and briefness, we will use 657 \texttt{lucas} from \secref{sec:Lucas numbers}. 658 659 \vspace{-1em} 660 \begin{alltt} 661 void 662 golden_pow(z_t r, z_t n) 663 \{ 664 if (zsignum(n) <= 0) 665 zsetu(r, zcmpi(n, -1) >= 0); 666 else if (!zcmpu(n, 1)) 667 zsetu(r, 2); 668 else 669 lucas(r, n); 670 \} 671 \end{alltt} 672 \vspace{-1em} 673 674 675 676 \item \textbf{zlshu and zrshu} 677 678 You are in big trouble, memorywise, of you 679 need to evaluate $2^{2^{64}}$. 680 681 682 683 \item \textbf{Modular left-shift} 684 685 \vspace{-1em} 686 \begin{alltt} 687 void 688 modlsh(z_t r, z_t a, z_t b) 689 \{ 690 z_t t, at; 691 size_t s = zbits(b); 692 693 zinit(t), zinit(at); 694 zset(at, a); 695 zsetu(r, 1); 696 zsetu(t, s); 697 698 while (zcmp(at, t) > 0) \{ 699 zsub(at, at, t); 700 zlsh(r, r, t); 701 zmod(r, r, b); 702 if (zzero(r)) 703 break; 704 \} 705 if (!zzero(a) && !zzero(b)) \{ 706 zlsh(r, r, a); 707 zmod(r, r, b); 708 \} 709 710 zfree(at), zfree(t); 711 \} 712 \end{alltt} 713 \vspace{-1em} 714 715 It is worth noting that this function is 716 not necessarily faster than \texttt{zmodpow}. 717 718 719 720 \item \textbf{Modular left-shift, extended} 721 722 The sign of \texttt{b} shall not effect the 723 result. Use \texttt{zabs} to create a copy 724 of \texttt{b} with its absolute value. 725 726 727 728 \item \textbf{The totient} 729 730 \( \displaystyle{ 731 \varphi(n) 732 = n \prod_{p \in \textbf{P} : p | n} \left ( 1 - \frac{1}{p} \right ) 733 = n \prod_{p \in \textbf{P} : p | n} \left ( \frac{p - 1}{p} \right ) 734 }\) 735 736 \noindent 737 So, if we set $a = n$ and $b = 1$, then we iterate 738 of all integers $p$, $2 \le p \le n$. For which $p$ 739 that is prime, we set $a \gets a \cdot (p - 1)$ and 740 $b \gets b \cdot p$. After the iteration, $b | a$, 741 and $\varphi(n) = \frac{a}{b}$. However, if $n < 0$, 742 then, $\varphi(n) = \varphi|n|$. 743 744 745 746 \item \textbf{The totient from factorisation} 747 748 \vspace{-1em} 749 \begin{alltt} 750 void 751 totient_fact(z_t t, z_t *P, 752 unsigned long long *K, size_t n) 753 \{ 754 z_t a, one; 755 zinit(a), zinit(one); 756 zseti(t, 1); 757 zseti(one, 1); 758 while (n--) \{ 759 zpowu(a, P[n], K[n] - 1); 760 zmul(t, t, a); 761 zsub(a, P[n], one); 762 zmul(t, t, a); 763 \} 764 zfree(a), zfree(one); 765 \} 766 \end{alltt} 767 \vspace{-1em} 768 769 770 771 \item \textbf{Modular tetration} 772 773 Let \texttt{totient} be Euler's totient function. 774 It is described in the problem ``The totient''. 775 776 We need two help function: \texttt{tetra(r, b, n)} 777 which calculated $r = {}^n{}b$, and \texttt{cmp\_tetra(a, b, n)} 778 which is compares $a$ to ${}^n{}b$. 779 780 \vspace{-1em} 781 \begin{alltt} 782 void 783 tetra(z_t r, z_t b, unsigned long n) 784 \{ 785 zsetu(r, 1); 786 while (n--) 787 zpow(r, b, r); 788 \} 789 \end{alltt} 790 \vspace{-1em} 791 792 \vspace{-1em} 793 \begin{alltt} 794 int 795 cmp_tetra(z_t a, z_t b, unsigned long n) 796 \{ 797 z_t B; 798 int cmp; 799 800 if (!n || !zcmpu(b, 1)) 801 return zcmpu(a, 1); 802 if (n == 1) 803 return zcmp(a, b); 804 if (zcmp(a, b) >= 0) 805 return +1; 806 807 zinit(B); 808 zsetu(B, 1); 809 while (n) \{ 810 zpow(B, b, B); 811 if (zcmp(a, B) < 0) \{ 812 zfree(B); 813 return -1; 814 \} 815 \} 816 cmp = zcmp(a, B); 817 zfree(B); 818 return cmp; 819 820 \} 821 \end{alltt} 822 \vspace{-1em} 823 824 \texttt{tetra} can generate unmaintainably huge 825 numbers. Will however only call \texttt{tetra} 826 when this is not the case. 827 828 \vspace{-1em} 829 \begin{alltt} 830 void 831 modtetra(z_t r, z_t b, unsigned long n, z_t m) 832 \{ 833 z_t t, temp; 834 835 if (n <= 1) \{ 836 if (!n) 837 zsetu(r, zcmpu(m, 1)); 838 else 839 zmod(r, b, m); 840 return; 841 \} 842 843 zmod(r, b, m); 844 if (zcmpu(r, 1) <= 0) 845 return; 846 847 zinit(t); 848 zinit(temp); 849 850 t = totient(m); 851 zgcd(temp, b, m); 852 853 if (!zcmpu(temp, 1)) \{ 854 modtetra(temp, b, n - 1, t); 855 zmodpow(r, r, temp, m); 856 \} else if (cmp_tetra(t, b, n - 1) > 0) \{ 857 temp = tetra(b, n - 1); 858 zpowmod(r, r, temp, m); 859 \} else \{ 860 modtetra(temp, b, n - 1, t); 861 zmodpow(temp, r, temp, m); 862 zmodpow(r, r, t, m); 863 zmodmul(r, r, temp, m); 864 \} 865 866 zfree(temp); 867 zfree(t); 868 \} 869 \end{alltt} 870 \vspace{-1em} 871 872 873 874 \item \textbf{Modular generalised power towers} 875 876 Instead of the signature 877 878 \vspace{-1em} 879 \begin{alltt} 880 void modtetra(z_t r, z_t b, unsigned long n, z_t m); 881 \end{alltt} 882 \vspace{-1em} 883 884 \noindent 885 you want to use the signature 886 887 \vspace{-1em} 888 \begin{alltt} 889 void modtower_(z_t r, z_t *a, size_t n, z_t m); 890 \end{alltt} 891 \vspace{-1em} 892 893 Instead of using, \texttt{b} (in \texttt{modtetra}), 894 use \texttt{*a}. At every recursion, instead of 895 passing on \texttt{a}, pass on \texttt{a + 1}. 896 897 The function \texttt{tetra} is modified into 898 899 \vspace{-1em} 900 \begin{alltt} 901 void tower(z_t r, z_t *a, size_t n) 902 \{ 903 zsetu(r, 1); 904 while (n--); 905 zpow(r, a[n], r); 906 \} 907 \end{alltt} 908 \vspace{-1em} 909 910 \noindent 911 \texttt{cmp\_tetra} is changed analogously. 912 913 To avoid problems in the evaluation, define 914 915 \vspace{-1em} 916 \begin{alltt} 917 void modtower(z_t r, z_t *a, size_t n, z_t m); 918 \end{alltt} 919 \vspace{-1em} 920 921 \noindent 922 which cuts the power short at the first element 923 of of \texttt{a} that equals 1. For example, if 924 $a = \{2, 3, 4, 5, 1, 2, 3\}$, and $n = 7$ 925 ($n = |a|$), then \texttt{modtower}, sets $n = 4$, 926 and effectively $a = \{2, 3, 4, 5\}$. After this 927 \texttt{modtower} calls \texttt{modtower\_}. 928 929 930 931 \end{enumerate}