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