# libzahl

big integer library
git clone git://git.suckless.org/libzahl

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
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);
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
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}