Remove some debugging prints and add comments
[maxima.git] / doc / info / Number.texi
blob2628dac3a52c56b9e0270828db30b2729808ce28
1 @menu
2 * Functions and Variables for Number Theory::  
3 @end menu
5 @c -----------------------------------------------------------------------------
6 @node Functions and Variables for Number Theory,  , Number Theory, Number Theory
7 @section Functions and Variables for Number Theory
8 @c -----------------------------------------------------------------------------
10 @c -----------------------------------------------------------------------------
11 @anchor{bern}
12 @deffn {Function} bern (@var{n})
14 Returns the @var{n}'th Bernoulli number for integer @var{n}.
15 @c WELL, ACTUALLY bern SIMPLIFIES, LIKE FACTORIAL -- DO WE WANT TO GET INTO THAT ???
16 @c OR JUST PRETEND IT'S "RETURNED" ???
17 Bernoulli numbers equal to zero are suppressed if @code{zerobern} is
18 @code{false}.
20 See also @mrefdot{burn}
22 @example
23 (%i1) zerobern: true$
24 (%i2) map (bern, [0, 1, 2, 3, 4, 5, 6, 7, 8]);
25                       1  1       1      1        1
26 (%o2)           [1, - -, -, 0, - --, 0, --, 0, - --]
27                       2  6       30     42       30
28 (%i3) zerobern: false$
29 (%i4) map (bern, [0, 1, 2, 3, 4, 5, 6, 7, 8]);
30                       1  1    1   1     1   5     691   7
31 (%o4)           [1, - -, -, - --, --, - --, --, - ----, -]
32                       2  6    30  42    30  66    2730  6
33 @end example
35 @opencatbox{Categories:}
36 @category{Number theory}
37 @closecatbox
38 @end deffn
40 @c -----------------------------------------------------------------------------
41 @anchor{bernpoly}
42 @deffn {Function} bernpoly (@var{x}, @var{n})
44 Returns the @var{n}'th Bernoulli polynomial in the
45 variable @var{x}.
47 @opencatbox{Categories:}
48 @category{Number theory}
49 @closecatbox
50 @end deffn
52 @c -----------------------------------------------------------------------------
53 @anchor{bfzeta}
54 @deffn {Function} bfzeta (@var{s}, @var{n})
56 Returns the Riemann zeta function for the argument @var{s}.
57 The return value is a big float (bfloat);
58 @var{n} is the number of digits in the return value.
60 @opencatbox{Categories:}
61 @category{Number theory}
62 @category{Numerical evaluation}
63 @closecatbox
64 @end deffn
66 @c -----------------------------------------------------------------------------
67 @anchor{bfhzeta}
68 @deffn {Function} bfhzeta (@var{s}, @var{h}, @var{n})
70 Returns the Hurwitz zeta function for the arguments @var{s} and @var{h}.
71 The return value is a big float (bfloat);
72 @var{n} is the number of digits in the return value.
74 The Hurwitz zeta function is defined as
76 @tex
77 $$\zeta \left(s,h\right) = \sum_{k=0}^\infty {1 \over \left(k+h\right)^{s}}$$
78 @end tex
79 @ifnottex
80 @example
81                         inf
82                         ====
83                         \        1
84          zeta (s,h)  =   >    --------
85                         /            s
86                         ====  (k + h)
87                         k = 0
88 @end example
89 @end ifnottex
91 @code{load ("bffac")} loads this function.
93 @opencatbox{Categories:}
94 @category{Number theory}
95 @category{Numerical evaluation}
96 @closecatbox
97 @end deffn
99 @c -----------------------------------------------------------------------------
100 @anchor{burn}
101 @deffn {Function} burn (@var{n})
103 Returns a rational number, which is an approximation of the @var{n}'th Bernoulli
104 number for integer @var{n}.  @code{burn} exploits the observation that
105 (rational) Bernoulli numbers can be approximated by (transcendental) zetas with
106 tolerable efficiency:
108 @example
109                    n - 1  1 - 2 n
110               (- 1)      2        zeta(2 n) (2 n)!
111      B(2 n) = ------------------------------------
112                                 2 n
113                              %pi
114 @end example
116 @code{burn} may be more efficient than @mref{bern} for large, isolated @var{n}
117 as @mref{bern} computes all the Bernoulli numbers up to index @var{n} before 
118 returning.  @code{burn} invokes the approximation for even integers @code{n >
119 255}.  For odd integers and @code{n <= 255} the function @mref{bern} is called.
121 @code{load ("bffac")} loads this function.  See also @mrefdot{bern}
123 @opencatbox{Categories:}
124 @category{Number theory}
125 @closecatbox
126 @end deffn
128 @c -----------------------------------------------------------------------------
129 @anchor{solve_congruences}
130 @deffn {Function} solve_congruences ([@var{r_1}, @dots{}, @var{r_n}], [@var{m_1}, @dots{}, @var{m_n}])
132 Solves the system of congruences @code{x = r_1 mod m_1}, @dots{}, @code{x = r_n mod m_n}.
133 The remainders @var{r_n} may be arbitrary integers while the moduli @var{m_n} have to be 
134 positive and pairwise coprime integers.
136 @example
137 (%i1) mods : [1000, 1001, 1003, 1007];
138 (%o1)                   [1000, 1001, 1003, 1007]
139 (%i2) lreduce('gcd, mods);
140 (%o2)                               1
141 (%i3) x : random(apply("*", mods));
142 (%o3)                         685124877004
143 (%i4) rems : map(lambda([z], mod(x, z)), mods);
144 (%o4)                       [4, 568, 54, 624]
145 (%i5) solve_congruences(rems, mods);
146 (%o5)                         685124877004
147 (%i6) solve_congruences([1, 2], [3, n]);
148 (%o6)               solve_congruences([1, 2], [3, n])
149 (%i7) %, n = 4;
150 (%o7)                              10
151 @end example
153 @opencatbox{Categories:}
154 @category{Number theory}
155 @closecatbox
156 @end deffn
158 @c -----------------------------------------------------------------------------
159 @anchor{cf}
160 @deffn {Function} cf (@var{expr})
162 Computes a continued fraction approximation.
163 @var{expr} is an expression comprising continued fractions,
164 square roots of integers, and literal real numbers
165 (integers, rational numbers, ordinary floats, and bigfloats).
166 @code{cf} computes exact expansions for rational numbers,
167 but expansions are truncated at @code{ratepsilon} for ordinary floats
168 and @code{10^(-fpprec)} for bigfloats.
170 Operands in the expression may be combined with arithmetic operators.
171 Maxima does not know about operations on continued fractions
172 outside of @code{cf}.
174 @code{cf} evaluates its arguments after binding @code{listarith} to
175 @code{false}.  @code{cf} returns a continued fraction, represented as a list.
177 A continued fraction @code{a + 1/(b + 1/(c + ...))} is represented by the list
178 @code{[a, b, c, ...]}.  The list elements @code{a}, @code{b}, @code{c}, @dots{}
179 must evaluate to integers.  @var{expr} may also contain @code{sqrt (n)} where
180 @code{n} is an integer.  In this case @code{cf} will give as many terms of the
181 continued fraction as the value of the variable @mref{cflength} times the
182 period.
184 A continued fraction can be evaluated to a number by evaluating the arithmetic
185 representation returned by @mrefdot{cfdisrep}  See also @mref{cfexpand} for
186 another way to evaluate a continued fraction.
188 See also @mrefcomma{cfdisrep} @mrefcomma{cfexpand} and @mrefdot{cflength}
190 Examples:
192 @itemize @bullet
193 @item
194 @var{expr} is an expression comprising continued fractions and square roots of
195 integers.
197 @example
198 (%i1) cf ([5, 3, 1]*[11, 9, 7] + [3, 7]/[4, 3, 2]);
199 (%o1)               [59, 17, 2, 1, 1, 1, 27]
200 (%i2) cf ((3/17)*[1, -2, 5]/sqrt(11) + (8/13));
201 (%o2)        [0, 1, 1, 1, 3, 2, 1, 4, 1, 9, 1, 9, 2]
202 @end example
204 @item
205 @code{cflength} controls how many periods of the continued fraction
206 are computed for algebraic, irrational numbers.
208 @example
209 (%i1) cflength: 1$
210 (%i2) cf ((1 + sqrt(5))/2);
211 (%o2)                    [1, 1, 1, 1, 2]
212 (%i3) cflength: 2$
213 (%i4) cf ((1 + sqrt(5))/2);
214 (%o4)               [1, 1, 1, 1, 1, 1, 1, 2]
215 (%i5) cflength: 3$
216 (%i6) cf ((1 + sqrt(5))/2);
217 (%o6)           [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2]
218 @end example
220 @item
221 A continued fraction can be evaluated by evaluating the arithmetic
222 representation returned by @mrefdot{cfdisrep}
224 @example
225 (%i1) cflength: 3$
226 (%i2) cfdisrep (cf (sqrt (3)))$
227 (%i3) ev (%, numer);
228 (%o3)                   1.731707317073171
229 @end example
231 @item
232 Maxima does not know about operations on continued fractions outside of
233 @code{cf}.
235 @example
236 (%i1) cf ([1,1,1,1,1,2] * 3);
237 (%o1)                     [4, 1, 5, 2]
238 (%i2) cf ([1,1,1,1,1,2]) * 3;
239 (%o2)                  [3, 3, 3, 3, 3, 6]
240 @end example
242 @end itemize
244 @opencatbox{Categories:}
245 @category{Continued fractions}
246 @closecatbox
247 @end deffn
249 @c NEEDS CLARIFICATION -- MAKE EXPLICIT HOW list IS RELATED TO a, b, c, ...
250 @c ALSO, CAN list CONTAIN ANYTHING OTHER THAN LITERAL INTEGERS ??
252 @c -----------------------------------------------------------------------------
253 @anchor{cfdisrep}
254 @deffn {Function} cfdisrep (@var{list})
256 Constructs and returns an ordinary arithmetic expression
257 of the form @code{a + 1/(b + 1/(c + ...))}
258 from the list representation of a continued fraction @code{[a, b, c, ...]}.
260 @example
261 (%i1) cf ([1, 2, -3] + [1, -2, 1]);
262 (%o1)                     [1, 1, 1, 2]
263 (%i2) cfdisrep (%);
264                                   1
265 (%o2)                     1 + ---------
266                                     1
267                               1 + -----
268                                       1
269                                   1 + -
270                                       2
271 @end example
273 @opencatbox{Categories:}
274 @category{Continued fractions}
275 @closecatbox
276 @end deffn
278 @c -----------------------------------------------------------------------------
279 @anchor{cfexpand}
280 @deffn {Function} cfexpand (@var{x})
282 Returns a matrix of the numerators and denominators of the last (column 1) and
283 next-to-last (column 2) convergents of the continued fraction @var{x}.
285 @example
286 (%i1) cf (rat (ev (%pi, numer)));
288 `rat' replaced 3.141592653589793 by 103993/33102 =3.141592653011902
289 (%o1)                  [3, 7, 15, 1, 292]
290 (%i2) cfexpand (%); 
291                          [ 103993  355 ]
292 (%o2)                    [             ]
293                          [ 33102   113 ]
294 (%i3) %[1,1]/%[2,1], numer;
295 (%o3)                   3.141592653011902
296 @end example
298 @opencatbox{Categories:}
299 @category{Continued fractions}
300 @closecatbox
301 @end deffn
303 @c -----------------------------------------------------------------------------
304 @anchor{cflength}
305 @defvr {Option variable} cflength
306 Default value: 1
308 @code{cflength} controls the number of terms of the continued fraction the
309 function @code{cf} will give, as the value @code{cflength} times the period.
310 Thus the default is to give one period.
312 @example
313 (%i1) cflength: 1$
314 (%i2) cf ((1 + sqrt(5))/2);
315 (%o2)                    [1, 1, 1, 1, 2]
316 (%i3) cflength: 2$
317 (%i4) cf ((1 + sqrt(5))/2);
318 (%o4)               [1, 1, 1, 1, 1, 1, 1, 2]
319 (%i5) cflength: 3$
320 (%i6) cf ((1 + sqrt(5))/2);
321 (%o6)           [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2]
322 @end example
324 @opencatbox{Categories:}
325 @category{Continued fractions}
326 @closecatbox
327 @end defvr
329 @c -----------------------------------------------------------------------------
330 @anchor{divsum}
331 @deffn  {Function} divsum @
332 @fname{divsum} (@var{n}, @var{k}) @
333 @fname{divsum} (@var{n})
335 @code{divsum (@var{n}, @var{k})} returns the sum of the divisors of @var{n}
336 raised to the @var{k}'th power.
338 @code{divsum (@var{n})} returns the sum of the divisors of @var{n}.
340 @example
341 (%i1) divsum (12);
342 (%o1)                          28
343 (%i2) 1 + 2 + 3 + 4 + 6 + 12;
344 (%o2)                          28
345 (%i3) divsum (12, 2);
346 (%o3)                          210
347 (%i4) 1^2 + 2^2 + 3^2 + 4^2 + 6^2 + 12^2;
348 (%o4)                          210
349 @end example
351 @opencatbox{Categories:}
352 @category{Number theory}
353 @closecatbox
354 @end deffn
356 @c -----------------------------------------------------------------------------
357 @anchor{euler}
358 @deffn {Function} euler (@var{n})
360 Returns the @var{n}'th Euler number for nonnegative integer @var{n}.
361 Euler numbers equal to zero are suppressed if @code{zerobern} is
362 @code{false}.
364 For the Euler-Mascheroni constant, see @code{%gamma}.
366 @example
367 (%i1) zerobern: true$
368 (%i2) map (euler, [0, 1, 2, 3, 4, 5, 6]);
369 (%o2)               [1, 0, - 1, 0, 5, 0, - 61]
370 (%i3) zerobern: false$
371 (%i4) map (euler, [0, 1, 2, 3, 4, 5, 6]);
372 (%o4)               [1, - 1, 5, - 61, 1385, - 50521, 2702765]
373 @end example
375 @opencatbox{Categories:}
376 @category{Number theory}
377 @closecatbox
378 @end deffn
380 @c -----------------------------------------------------------------------------
381 @anchor{factors_only}
382 @defvr {Option variable} factors_only
383 Default value: @code{false}
385 Controls the value returned by @mrefdot{ifactors} The default @code{false} 
386 causes @code{ifactors} to provide information about multiplicities of the 
387 computed prime factors. If @code{factors_only} is set to @code{true}, 
388 @code{ifactors} returns nothing more than a list of prime factors.
390 Example: See @mrefdot{ifactors}
392 @opencatbox{Categories:}
393 @category{Number theory}
394 @closecatbox
395 @end defvr
397 @c -----------------------------------------------------------------------------
398 @anchor{fib}
399 @deffn {Function} fib (@var{n})
401 Returns the @var{n}'th Fibonacci number.
402 @code{fib(0)} is equal to 0 and @code{fib(1)} equal to 1, and 
403 @code{fib (-@var{n})} equal to @code{(-1)^(@var{n} + 1) * fib(@var{n})}.
405 @example
406 (%i1) map (fib, [-4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8]);
407 (%o1)           [- 3, 2, - 1, 1, 0, 1, 1, 2, 3, 5, 8, 13, 21]
408 @end example
410 @opencatbox{Categories:}
411 @category{Number theory}
412 @closecatbox
413 @end deffn
415 @c -----------------------------------------------------------------------------
416 @anchor{fibtophi}
417 @deffn {Function} fibtophi (@var{expr})
419 Expresses Fibonacci numbers in @var{expr} in terms of the constant @code{%phi},
420 which is @code{(1 + sqrt(5))/2}, approximately 1.61803399.
422 Examples:
424 @c ===beg===
425 @c fibtophi (fib (n));
426 @c fib (n-1) + fib (n) - fib (n+1);
427 @c fibtophi (%);
428 @c ratsimp (%);
429 @c ===end===
430 @example
431 (%i1) fibtophi (fib (n));
432                            n             n
433                        %phi  - (1 - %phi)
434 (%o1)                  -------------------
435                            2 %phi - 1
436 (%i2) fib (n-1) + fib (n) - fib (n+1);
437 (%o2)          - fib(n + 1) + fib(n) + fib(n - 1)
438 (%i3) fibtophi (%);
439             n + 1             n + 1       n             n
440         %phi      - (1 - %phi)        %phi  - (1 - %phi)
441 (%o3) - --------------------------- + -------------------
442                 2 %phi - 1                2 %phi - 1
443                                           n - 1             n - 1
444                                       %phi      - (1 - %phi)
445                                     + ---------------------------
446                                               2 %phi - 1
447 (%i4) ratsimp (%);
448 (%o4)                           0
449 @end example
451 @opencatbox{Categories:}
452 @category{Number theory}
453 @closecatbox
454 @end deffn
456 @c -----------------------------------------------------------------------------
457 @anchor{ifactors}
458 @deffn {Function} ifactors (@var{n})
460 For a positive integer @var{n} returns the factorization of @var{n}.  If
461 @code{n=p1^e1..pk^nk} is the decomposition of @var{n} into prime
462 factors, ifactors returns @code{[[p1, e1], ... , [pk, ek]]}.
464 Factorization methods used are trial divisions by primes up to 9973,
465 Pollard's rho and p-1 method and elliptic curves.
467 If the variable @code{ifactor_verbose} is set to @code{true}
468 ifactor produces detailed output about what it is doing including
469 immediate feedback as soon as a factor has been found.
471 The value returned by @code{ifactors} is controlled by the option variable @mrefdot{factors_only}
472 The default @code{false} causes @code{ifactors} to provide information about 
473 the multiplicities of the computed prime factors. If @code{factors_only} 
474 is set to @code{true}, @code{ifactors} simply returns the list of 
475 prime factors.
477 @example
478 (%i1) ifactors(51575319651600);
479 (%o1)     [[2, 4], [3, 2], [5, 2], [1583, 1], [9050207, 1]]
480 (%i2) apply("*", map(lambda([u], u[1]^u[2]), %));
481 (%o2)                        51575319651600
482 (%i3) ifactors(51575319651600), factors_only : true;
483 (%o3)                   [2, 3, 5, 1583, 9050207]
484 @end example
486 @opencatbox{Categories:}
487 @category{Number theory}
488 @closecatbox
489 @end deffn
491 @c -----------------------------------------------------------------------------
492 @anchor{igcdex}
493 @deffn {Function} igcdex (@var{n}, @var{k})
495 Returns a list @code{[@var{a}, @var{b}, @var{u}]} where @var{u} is the greatest
496 common divisor of @var{n} and @var{k}, and @var{u} is equal to
497 @code{@var{a} @var{n} + @var{b} @var{k}}.  The arguments @var{n} and @var{k}
498 must be integers.
500 @code{igcdex} implements the Euclidean algorithm.  See also @mrefdot{gcdex}
502 The command @code{load("gcdex")} loads the function.
504 Examples:
506 @example
507 (%i1) load("gcdex")$
509 (%i2) igcdex(30,18);
510 (%o2)                      [- 1, 2, 6]
511 (%i3) igcdex(1526757668, 7835626735736);
512 (%o3)            [845922341123, - 164826435, 4]
513 (%i4) igcdex(fib(20), fib(21));
514 (%o4)                   [4181, - 2584, 1]
515 @end example
517 @opencatbox{Categories:}
518 @category{Number theory}
519 @closecatbox
520 @end deffn
522 @c -----------------------------------------------------------------------------
523 @anchor{inrt}
524 @deffn {Function} inrt (@var{x}, @var{n})
526 Returns the integer @var{n}'th root of the absolute value of @var{x}.
528 @example
529 (%i1) l: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]$
530 (%i2) map (lambda ([a], inrt (10^a, 3)), l);
531 (%o2) [2, 4, 10, 21, 46, 100, 215, 464, 1000, 2154, 4641, 10000]
532 @end example
534 @opencatbox{Categories:}
535 @category{Number theory}
536 @closecatbox
537 @end deffn
539 @c -----------------------------------------------------------------------------
540 @anchor{inv_mod}
541 @deffn {Function} inv_mod (@var{n}, @var{m})
543 Computes the inverse of @var{n} modulo @var{m}.
544 @code{inv_mod (n,m)} returns @code{false}, 
545 if @var{n} is a zero divisor modulo @var{m}.
547 @example
548 (%i1) inv_mod(3, 41);
549 (%o1)                           14
550 (%i2) ratsimp(3^-1), modulus = 41;
551 (%o2)                           14
552 (%i3) inv_mod(3, 42);
553 (%o3)                          false
554 @end example
556 @opencatbox{Categories:}
557 @category{Number theory}
558 @closecatbox
559 @end deffn
561 @c -----------------------------------------------------------------------------
562 @anchor{isqrt}
563 @deffn {Function} isqrt (@var{x})
565 Returns the "integer square root" of the absolute value of @var{x}, which is an
566 integer.
568 @opencatbox{Categories:}
569 @category{Mathematical functions}
570 @closecatbox
571 @end deffn
573 @c -----------------------------------------------------------------------------
574 @anchor{jacobi}
575 @deffn {Function} jacobi (@var{p}, @var{q})
577 Returns the Jacobi symbol of @var{p} and @var{q}.
579 @example
580 (%i1) l: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]$
581 (%i2) map (lambda ([a], jacobi (a, 9)), l);
582 (%o2)         [1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0]
583 @end example
585 @opencatbox{Categories:}
586 @category{Number theory}
587 @closecatbox
588 @end deffn
590 @c -----------------------------------------------------------------------------
591 @anchor{lcm}
592 @deffn {Function} lcm (@var{expr_1}, @dots{}, @var{expr_n})
594 Returns the least common multiple of its arguments.
595 The arguments may be general expressions as well as integers.
597 @code{load ("functs")} loads this function.
599 @opencatbox{Categories:}
600 @category{Number theory}
601 @closecatbox
602 @end deffn
604 @c -----------------------------------------------------------------------------
605 @anchor{lucas}
606 @deffn {Function} lucas (@var{n})
608 Returns the @var{n}'th Lucas number.
609 @code{lucas(0)} is equal to 2 and @code{lucas(1)} equal to 1, and
610 in general, @code{lucas(@var{n}) = lucas(@var{n}-1) + lucas(@var{n}-2)}.  Also 
611 @code{lucas(-@var{n})} is equal to @code{(-1)^(-@var{n}) * lucas(@var{n})}.
613 @example
614 (%i1) map (lucas, [-4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8]);
615 (%o1)    [7, - 4, 3, - 1, 2, 1, 3, 4, 7, 11, 18, 29, 47]
616 @end example
619 @opencatbox{Categories:}
620 @category{Number theory}
621 @closecatbox
622 @end deffn
624 @c -----------------------------------------------------------------------------
625 @anchor{mod}
626 @deffn {Function} mod (@var{x}, @var{y})
628 If @var{x} and @var{y} are real numbers and @var{y} is nonzero, return
629 @code{@var{x} - @var{y} * floor(@var{x} / @var{y})}.  Further for all real
630 @var{x}, we have @code{mod (@var{x}, 0) = @var{x}}.  For a discussion of the
631 definition @code{mod (@var{x}, 0) = @var{x}}, see Section 3.4, of
632 "Concrete Mathematics," by Graham, Knuth, and Patashnik.  The function
633 @code{mod (@var{x}, 1)} is a sawtooth function with period 1 with
634 @code{mod (1, 1) = 0} and @code{mod (0, 1) = 0}.
636 To find the principal argument (a number in the interval @code{(-%pi, %pi]}) of
637 a complex number, use the function
638 @code{@var{x} |-> %pi - mod (%pi - @var{x}, 2*%pi)}, where @var{x} is an
639 argument.
641 When @var{x} and @var{y} are constant expressions (@code{10 * %pi}, for 
642 example), @code{mod} uses the same big float evaluation scheme that @code{floor}
643 and @code{ceiling} uses.  Again, it's possible, although unlikely, that
644 @code{mod} could return an erroneous value in such cases.
646 For nonnumerical arguments @var{x} or @var{y}, @code{mod} knows several
647 simplification rules:
649 @c ===beg===
650 @c mod (x, 0);
651 @c mod (a*x, a*y);
652 @c mod (0, x);
653 @c ===end===
654 @example
655 (%i1) mod (x, 0);
656 (%o1)                           x
657 (%i2) mod (a*x, a*y);
658 (%o2)                      a mod(x, y)
659 (%i3) mod (0, x);
660 (%o3)                           0
661 @end example
663 @opencatbox{Categories:}
664 @category{Mathematical functions}
665 @closecatbox
666 @end deffn
668 @c -----------------------------------------------------------------------------
669 @anchor{next_prime}
670 @deffn {Function} next_prime (@var{n})
672 Returns the smallest prime bigger than @var{n}.
674 @example
675 (%i1) next_prime(27);
676 (%o1)                       29
677 @end example
679 @opencatbox{Categories:}
680 @category{Number theory}
681 @closecatbox
682 @end deffn
684 @c -----------------------------------------------------------------------------
685 @anchor{partfrac}
686 @deffn {Function} partfrac (@var{expr}, @var{var})
688 Expands the expression @var{expr} in partial fractions
689 with respect to the main variable @var{var}.  @code{partfrac} does a complete
690 partial fraction decomposition.  The algorithm employed is based on
691 the fact that the denominators of the partial fraction expansion (the
692 factors of the original denominator) are relatively prime.  The
693 numerators can be written as linear combinations of denominators, and
694 the expansion falls out.
696 @code{partfrac} ignores the value @code{true} of the option variable
697 @code{keepfloat}.
699 @example
700 (%i1) 1/(1+x)^2 - 2/(1+x) + 2/(2+x);
701                       2       2        1
702 (%o1)               ----- - ----- + --------
703                     x + 2   x + 1          2
704                                     (x + 1)
705 (%i2) ratsimp (%);
706                                  x
707 (%o2)                 - -------------------
708                          3      2
709                         x  + 4 x  + 5 x + 2
710 (%i3) partfrac (%, x);
711                       2       2        1
712 (%o3)               ----- - ----- + --------
713                     x + 2   x + 1          2
714                                     (x + 1)
715 @end example
716 @end deffn
718 @c -----------------------------------------------------------------------------
719 @anchor{power_mod}
720 @deffn {Function} power_mod (@var{a}, @var{n}, @var{m})
722 Uses a modular algorithm to compute @code{a^n mod m} 
723 where @var{a} and @var{n} are integers and @var{m} is a positive integer.
724 If @var{n} is negative, @code{inv_mod} is used to find the modular inverse.
726 @example
727 (%i1) power_mod(3, 15, 5);
728 (%o1)                          2
729 (%i2) mod(3^15,5);
730 (%o2)                          2
731 (%i3) power_mod(2, -1, 5);
732 (%o3)                          3
733 (%i4) inv_mod(2,5);
734 (%o4)                          3
735 @end example
737 @opencatbox{Categories:}
738 @category{Number theory}
739 @closecatbox
740 @end deffn
742 @c -----------------------------------------------------------------------------
743 @anchor{primep}
744 @deffn {Function} primep (@var{n})
746 Primality test.  If @code{primep (@var{n})} returns @code{false}, @var{n} is a
747 composite number and if it returns @code{true}, @var{n} is a prime number
748 with very high probability.
750 For @var{n} less than 3317044064679887385961981 a deterministic version of
751 Miller-Rabin's test is used.  If @code{primep (@var{n})} returns
752 @code{true}, then @var{n} is a prime number.
754 For @var{n} bigger than 3317044064679887385961981 @code{primep} uses
755 @code{primep_number_of_tests} Miller-Rabin's pseudo-primality tests and one 
756 Lucas pseudo-primality test.  The probability that a non-prime @var{n} will 
757 pass one Miller-Rabin test is less than 1/4.  Using the default value 25 for
758 @code{primep_number_of_tests}, the probability of @var{n} being
759 composite is much smaller that 10^-15.
761 @opencatbox{Categories:}
762 @category{Predicate functions}
763 @category{Number theory}
764 @closecatbox
765 @end deffn
767 @c -----------------------------------------------------------------------------
768 @anchor{primep_number_of_tests}
769 @defvr {Option variable} primep_number_of_tests
770 Default value: 25
772 Number of Miller-Rabin's tests used in @code{primep}.
774 @opencatbox{Categories:}
775 @category{Number theory}
776 @closecatbox
777 @end defvr
779 @c -----------------------------------------------------------------------------
780 @anchor{primes}
781 @deffn {Function} primes (@var{start}, @var{end})
783 Returns the list of all primes from @var{start} to @var{end}.
785 @example
786 (%i1) primes(3, 7);
787 (%o1)                     [3, 5, 7]
788 @end example
790 @opencatbox{Categories:}
791 @category{Number theory}
792 @closecatbox
793 @end deffn
795 @c -----------------------------------------------------------------------------
796 @anchor{prev_time}
797 @deffn {Function} prev_prime (@var{n})
799 Returns the greatest prime smaller than @var{n}.
801 @example
802 (%i1) prev_prime(27);
803 (%o1)                       23
804 @end example
806 @opencatbox{Categories:}
807 @category{Number theory}
808 @closecatbox
809 @end deffn
811 @c -----------------------------------------------------------------------------
812 @anchor{qunit}
813 @deffn {Function} qunit (@var{n})
815 Returns the principal unit of the real quadratic number field
816 @code{sqrt (@var{n})} where @var{n} is an integer,
817 i.e., the element whose norm is unity.
818 This amounts to solving Pell's equation @code{a^2 - @var{n} b^2 = 1}.
820 @example
821 (%i1) qunit (17);
822 (%o1)                     sqrt(17) + 4
823 (%i2) expand (% * (sqrt(17) - 4));
824 (%o2)                           1
825 @end example
827 @opencatbox{Categories:}
828 @category{Number theory}
829 @closecatbox
830 @end deffn
832 @c -----------------------------------------------------------------------------
833 @anchor{totient}
834 @deffn {Function} totient (@var{n})
836 Returns the number of integers less than or equal to @var{n} which
837 are relatively prime to @var{n}.
839 @opencatbox{Categories:}
840 @category{Number theory}
841 @closecatbox
842 @end deffn
844 @c -----------------------------------------------------------------------------
845 @defvr {Option variable} zerobern
846 Default value: @code{true}
848 When @code{zerobern} is @code{false}, @code{bern} excludes the Bernoulli numbers
849 and @code{euler} excludes the Euler numbers which are equal to zero.
850 See @code{bern} and @code{euler}.
852 @opencatbox{Categories:}
853 @category{Number theory}
854 @closecatbox
855 @end defvr
857 @c -----------------------------------------------------------------------------
858 @anchor{zeta}
859 @deffn {Function} zeta (@var{n})
861 Returns the Riemann zeta function.  If @var{n} is a negative integer, 0, or a
862 positive even integer, the Riemann zeta function simplifies to an exact value.
863 For a positive even integer the option variable @code{zeta%pi} has to be
864 @code{true} in addition (See @code{zeta%pi}).  For a floating point or bigfloat
865 number the Riemann zeta function is evaluated numerically.  Maxima returns a
866 noun form @code{zeta (@var{n})} for all other arguments, including rational
867 noninteger, and complex arguments, or for even integers, if @code{zeta%pi} has
868 the value @code{false}.
870 @code{zeta(1)} is undefined, but Maxima knows the limit 
871 @code{limit(zeta(x), x, 1)} from above and below.
873 The Riemann zeta function distributes over lists, matrices, and equations.
875 See also @mref{bfzeta} and @mrefdot{zeta%pi}
877 Examples:
879 @c ===beg===
880 @c zeta([-2, -1, 0, 0.5, 2, 3,1+%i]);
881 @c limit(zeta(x),x,1,plus);
882 @c limit(zeta(x),x,1,minus);
883 @c ===end===
884 @example
885 (%i1) zeta([-2, -1, 0, 0.5, 2, 3, 1+%i]);
886                                              2
887             1     1                       %pi
888 (%o1) [0, - --, - -, - 1.460354508809586, ----, zeta(3), 
889             12    2                        6
890                                                     zeta(%i + 1)]
891 (%i2) limit(zeta(x),x,1,plus);
892 (%o2)                          inf
893 (%i3) limit(zeta(x),x,1,minus);
894 (%o3)                         minf
895 @end example
897 @opencatbox{Categories:}
898 @category{Number theory}
899 @closecatbox
900 @end deffn
902 @c -----------------------------------------------------------------------------
903 @anchor{zeta%pi}
904 @defvr {Option variable} zeta%pi
905 Default value: @mref{true}
907 When @code{zeta%pi} is @code{true}, @code{zeta} returns an expression 
908 proportional to @code{%pi^n} for even integer @code{n}.  Otherwise, @code{zeta} 
909 returns a noun form @code{zeta (n)} for even integer @code{n}.
911 Examples:
913 @c ===beg===
914 @c zeta%pi: true$
915 @c zeta (4);
916 @c zeta%pi: false$
917 @c zeta (4);
918 @c ===end===
919 @example
920 (%i1) zeta%pi: true$
921 (%i2) zeta (4);
922                                  4
923                               %pi
924 (%o2)                         ----
925                                90
926 (%i3) zeta%pi: false$
927 (%i4) zeta (4);
928 (%o4)                        zeta(4)
929 @end example
931 @opencatbox{Categories:}
932 @category{Number theory}
933 @closecatbox
934 @end defvr
936 @c -----------------------------------------------------------------------------
937 @anchor{zn_add_table}
938 @deffn {Function} zn_add_table (@var{n}) 
940 Shows an addition table of all elements in (Z/@var{n}Z).
942 See also @mrefcomma{zn_mult_table}  @mrefdot{zn_power_table}
944 @opencatbox{Categories:}
945 @category{Number theory}
946 @closecatbox
947 @end deffn
949 @c -----------------------------------------------------------------------------
950 @anchor{zn_characteristic_factors}
951 @deffn {Function} zn_characteristic_factors (@var{n}) 
953 Returns a list containing the characteristic factors of the totient of @var{n}.
955 Using the characteristic factors a multiplication group modulo @var{n} 
956 can be expressed as a group direct product of cyclic subgroups.
958 In case the group itself is cyclic the list only contains the totient 
959 and using @code{zn_primroot} a generator can be computed. 
960 If the totient splits into more than one characteristic factors 
961 @code{zn_factor_generators} finds generators of the corresponding subgroups.
963 Each of the @code{r} factors in the list divides the right following factors. 
964 For the last factor @code{f_r} therefore holds @code{a^f_r = 1 (mod n)} 
965 for all @code{a} coprime to @var{n}.  
966 This factor is also known as Carmichael function or Carmichael lambda.
968 If @code{n > 2}, then @code{totient(n)/2^r} is the number of quadratic residues, 
969 and each of these has @code{2^r} square roots.
971 See also @mrefcomma{totient}  @mrefcomma{zn_primroot}  @mrefdot{zn_factor_generators}
973 Examples:
975 The multiplication group modulo @code{14} is cyclic and its @code{6} elements 
976 can be generated by a primitive root.
978 @example
979 (%i1) [zn_characteristic_factors(14), phi: totient(14)];
980 (%o1)                              [[6], 6]
981 (%i2) [zn_factor_generators(14), g: zn_primroot(14)];
982 (%o2)                              [[3], 3]
983 (%i3) M14: makelist(power_mod(g,i,14), i,0,phi-1);
984 (%o3)                         [1, 3, 9, 13, 11, 5]
985 @end example
987 The multiplication group modulo @code{15} is not cyclic and its @code{8} elements 
988 can be generated by two factor generators.
990 @example
991 (%i1) [[f1,f2]: zn_characteristic_factors(15), totient(15)];
992 (%o1)                             [[2, 4], 8]
993 (%i2) [[g1,g2]: zn_factor_generators(15), zn_primroot(15)];
994 (%o2)                           [[11, 7], false]
995 (%i3) UG1: makelist(power_mod(g1,i,15), i,0,f1-1);
996 (%o3)                               [1, 11]
997 (%i4) UG2: makelist(power_mod(g2,i,15), i,0,f2-1);
998 (%o4)                            [1, 7, 4, 13]
999 (%i5) M15: create_list(mod(i*j,15), i,UG1, j,UG2);
1000 (%o5)                      [1, 7, 4, 13, 11, 2, 14, 8]
1001 @end example
1003 For the last characteristic factor @code{4} it holds that @code{a^4 = 1 (mod 15)} 
1004 for all @code{a} in @code{M15}. 
1006 @code{M15} has two characteristic factors and therefore @code{8/2^2} quadratic residues, 
1007 and each of these has @code{2^2} square roots.
1009 @example
1010 (%i6) zn_power_table(15);
1011                                [ 1   1  1   1 ]
1012                                [              ]
1013                                [ 2   4  8   1 ]
1014                                [              ]
1015                                [ 4   1  4   1 ]
1016                                [              ]
1017                                [ 7   4  13  1 ]
1018 (%o6)                          [              ]
1019                                [ 8   4  2   1 ]
1020                                [              ]
1021                                [ 11  1  11  1 ]
1022                                [              ]
1023                                [ 13  4  7   1 ]
1024                                [              ]
1025                                [ 14  1  14  1 ]
1026 (%i7) map(lambda([i], zn_nth_root(i,2,15)), [1,4]);
1027 (%o7)                   [[1, 4, 11, 14], [2, 7, 8, 13]]
1028 @end example
1030 @opencatbox{Categories:}
1031 @category{Number theory}
1032 @closecatbox
1033 @end deffn
1035 @c -----------------------------------------------------------------------------
1036 @anchor{zn_carmichael_lambda}
1037 @deffn {Function} zn_carmichael_lambda (@var{n}) 
1039 Returns @code{1} if @var{n} is @code{1} and otherwise 
1040 the greatest characteristic factor of the totient of @var{n}.
1042 For remarks and examples see @mrefdot{zn_characteristic_factors}
1044 @opencatbox{Categories:}
1045 @category{Number theory}
1046 @closecatbox
1047 @end deffn
1049 @c -----------------------------------------------------------------------------
1050 @anchor{zn_determinant}
1051 @deffn {Function} zn_determinant (@var{matrix}, @var{p}) 
1053 Uses the technique of LU-decomposition to compute the determinant of @var{matrix} 
1054 over (Z/@var{p}Z). @var{p} must be a prime. 
1056 However if the determinant is equal to zero the LU-decomposition might fail. 
1057 In that case @code{zn_determinant} computes the determinant non-modular 
1058 and reduces thereafter.
1060 See also @mrefdot{zn_invert_by_lu}
1062 Examples:
1064 @example
1065 (%i1) m : matrix([1,3],[2,4]);
1066                                 [ 1  3 ]
1067 (%o1)                           [      ]
1068                                 [ 2  4 ]
1069 (%i2) zn_determinant(m, 5);
1070 (%o2)                               3
1071 (%i3) m : matrix([2,4,1],[3,1,4],[4,3,2]);
1072                                [ 2  4  1 ]
1073                                [         ]
1074 (%o3)                          [ 3  1  4 ]
1075                                [         ]
1076                                [ 4  3  2 ]
1077 (%i4) zn_determinant(m, 5);
1078 (%o4)                               0
1079 @end example
1081 @opencatbox{Categories:}
1082 @category{Number theory}
1083 @closecatbox
1084 @end deffn
1086 @c -----------------------------------------------------------------------------
1087 @anchor{zn_factor_generators}
1088 @deffn {Function} zn_factor_generators (@var{n}) 
1090 Returns a list containing factor generators corresponding to the 
1091 characteristic factors of the totient of @var{n}.
1093 For remarks and examples see @mrefdot{zn_characteristic_factors}
1095 @opencatbox{Categories:}
1096 @category{Number theory}
1097 @closecatbox
1098 @end deffn
1100 @c -----------------------------------------------------------------------------
1101 @anchor{zn_invert_by_lu}
1102 @deffn {Function} zn_invert_by_lu (@var{matrix}, @var{p}) 
1104 Uses the technique of LU-decomposition to compute the modular inverse of 
1105 @var{matrix} over (Z/@var{p}Z). @var{p} must be a prime and @var{matrix} 
1106 invertible. @code{zn_invert_by_lu} returns @code{false} if @var{matrix} 
1107 is not invertible.
1109 See also @mrefdot{zn_determinant}
1111 Example:
1113 @example
1114 (%i1) m : matrix([1,3],[2,4]);
1115                                 [ 1  3 ]
1116 (%o1)                           [      ]
1117                                 [ 2  4 ]
1118 (%i2) zn_determinant(m, 5);
1119 (%o2)                               3
1120 (%i3) mi : zn_invert_by_lu(m, 5);
1121                                 [ 3  4 ]
1122 (%o3)                           [      ]
1123                                 [ 1  2 ]
1124 (%i4) matrixmap(lambda([a], mod(a, 5)), m . mi);
1125                                 [ 1  0 ]
1126 (%o4)                           [      ]
1127                                 [ 0  1 ]
1128 @end example
1130 @opencatbox{Categories:}
1131 @category{Number theory}
1132 @closecatbox
1133 @end deffn
1135 @c -----------------------------------------------------------------------------
1136 @anchor{zn_log}
1137 @deffn {Function} zn_log @
1138 @fname{zn_log} (@var{a}, @var{g}, @var{n})  @
1139 @fname{zn_log} (@var{a}, @var{g}, @var{n}, [[@var{p1}, @var{e1}], @dots{}, [@var{pk}, @var{ek}]])
1141 Computes the discrete logarithm. Let (Z/@var{n}Z)* be a cyclic group, @var{g} a 
1142 primitive root modulo @var{n} or a generator of a subgroup of (Z/@var{n}Z)* 
1143 and let @var{a} be a member of this group.  
1144 @code{zn_log (a, g, n)} then solves the congruence @code{g^x = a mod n}.
1145 Please note that if @var{a} is not a power of @var{g} modulo @var{n}, 
1146 @code{zn_log} will not terminate.
1148 The applied algorithm needs a prime factorization of @code{zn_order(g)} resp. @code{totient(n)} 
1149 in case @var{g} is a primitive root modulo @var{n}. 
1150 A precomputed list of factors of @code{zn_order(g)} might be used as the optional fourth argument.
1151 This list must be of the same form as the list returned by @code{ifactors(zn_order(g))} 
1152 using the default option @code{factors_only : false}.
1153 However, compared to the running time of the logarithm algorithm 
1154 providing the list of factors has only a quite small effect.
1156 The algorithm uses a Pohlig-Hellman-reduction and Pollard's Rho-method for 
1157 discrete logarithms. The running time of @code{zn_log} primarily depends on the 
1158 bitlength of the greatest prime factor of @code{zn_order(g)}.
1160 See also @mrefcomma{zn_primroot}  @mrefcomma{zn_order}  @mrefcomma{ifactors}  @mrefdot{totient}
1162 Examples:
1164 @code{zn_log (a, g, n)} solves the congruence @code{g^x = a mod n}.
1166 @example
1167 (%i1) n : 22$
1168 (%i2) g : zn_primroot(n);
1169 (%o2)                               7
1170 (%i3) ord_7 : zn_order(7, n);
1171 (%o3)                              10
1172 (%i4) powers_7 : makelist(power_mod(g, x, n), x, 0, ord_7 - 1);
1173 (%o4)              [1, 7, 5, 13, 3, 21, 15, 17, 9, 19]
1174 (%i5) zn_log(9, g, n);
1175 (%o5)                               8
1176 (%i6) map(lambda([x], zn_log(x, g, n)), powers_7);
1177 (%o6)                [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
1178 (%i7) ord_5 : zn_order(5, n);
1179 (%o7)                               5
1180 (%i8) powers_5 : makelist(power_mod(5,x,n), x, 0, ord_5 - 1);
1181 (%o8)                       [1, 5, 3, 15, 9]
1182 (%i9) zn_log(9, 5, n);
1183 (%o9)                               4
1184 @end example
1186 The optional fourth argument must be of the same form as the list returned by 
1187 @code{ifactors(zn_order(g))}.
1188 The running time primarily depends on the bitlength of the totient's greatest prime factor.
1190 @example
1191 (%i1) (p : 2^127-1, primep(p));
1192 (%o1)                             true
1193 (%i2) ifs : ifactors(p - 1)$
1194 (%i3) g : zn_primroot(p, ifs);
1195 (%o3)                              43
1196 (%i4) a : power_mod(g, 4711, p)$
1197 (%i5) zn_log(a, g, p, ifs);
1198 (%o5)                             4711
1199 (%i6) f_max : last(ifs);  
1200 (%o6)                       [77158673929, 1]
1201 (%i7) ord_5 : zn_order(5,p,ifs)$
1202 (%i8) (p - 1)/ord_5;
1203 (%o8)                              73
1204 (%i9) ifs_5 : ifactors(ord_5)$
1205 (%i10) a : power_mod(5, 4711, p)$
1206 (%i11) zn_log(a, 5, p, ifs_5);
1207 (%o11)                            4711
1208 @end example
1210 @opencatbox{Categories:}
1211 @category{Number theory}
1212 @closecatbox
1213 @end deffn
1215 @c -----------------------------------------------------------------------------
1216 @anchor{zn_mult_table}
1217 @deffn {Function} zn_mult_table @
1218 @fname{zn_mult_table} (@var{n})  @
1219 @fname{zn_mult_table} (@var{n}, @var{gcd})
1221 Without the optional argument @var{gcd} @code{zn_mult_table(n)} shows a 
1222 multiplication table of all elements in (Z/@var{n}Z)* which are all elements 
1223 coprime to @var{n}.
1225 The optional second argument @var{gcd} allows to select a specific 
1226 subset of (Z/@var{n}Z). If @var{gcd} is an integer, a multiplication table of 
1227 all residues @code{x} with @code{gcd(x,n) = }@var{gcd} are returned.
1228 Additionally row and column headings are added for better readability. 
1229 If necessary, these can be easily removed by @code{submatrix(1, table, 1)}. 
1231 If @var{gcd} is set to @code{all}, the table is printed for all non-zero 
1232 elements in (Z/@var{n}Z).
1234 The second example shows an alternative way to create a multiplication table 
1235 for subgroups.
1237 See also @mrefcomma{zn_add_table}  @mrefdot{zn_power_table}
1239 Examples:
1241 The default table shows all elements in (Z/@var{n}Z)* and allows to demonstrate 
1242 and study basic properties of modular multiplication groups.
1243 E.g. the principal diagonal contains all quadratic residues, 
1244 each row and column contains every element, the tables are symmetric, etc..
1246 If @var{gcd} is set to @code{all}, the table is printed for all non-zero 
1247 elements in (Z/@var{n}Z).
1249 @example
1250 (%i1) zn_mult_table(8);
1251                                 [ 1  3  5  7 ]
1252                                 [            ]
1253                                 [ 3  1  7  5 ]
1254 (%o1)                           [            ]
1255                                 [ 5  7  1  3 ]
1256                                 [            ]
1257                                 [ 7  5  3  1 ]
1258 (%i2) zn_mult_table(8, all);
1259                             [ 1  2  3  4  5  6  7 ]
1260                             [                     ]
1261                             [ 2  4  6  0  2  4  6 ]
1262                             [                     ]
1263                             [ 3  6  1  4  7  2  5 ]
1264                             [                     ]
1265 (%o2)                       [ 4  0  4  0  4  0  4 ]
1266                             [                     ]
1267                             [ 5  2  7  4  1  6  3 ]
1268                             [                     ]
1269                             [ 6  4  2  0  6  4  2 ]
1270                             [                     ]
1271                             [ 7  6  5  4  3  2  1 ]
1272 @end example
1274 If @var{gcd} is an integer, row and column headings are added for better readability. 
1276 If the subset chosen by @var{gcd} is a group there is another way to create 
1277 a multiplication table. An isomorphic mapping from a group with @code{1} as 
1278 identity builds a table which is easy to read. The mapping is accomplished via CRT.
1280 In the second version of @code{T36_4} the identity, here @code{28}, is placed in 
1281 the top left corner, just like in table @code{T9}. 
1283 @example
1284 (%i1) T36_4: zn_mult_table(36,4);
1285                         [ *   4   8   16  20  28  32 ]
1286                         [                            ]
1287                         [ 4   16  32  28  8   4   20 ]
1288                         [                            ]
1289                         [ 8   32  28  20  16  8   4  ]
1290                         [                            ]
1291 (%o1)                   [ 16  28  20  4   32  16  8  ]
1292                         [                            ]
1293                         [ 20  8   16  32  4   20  28 ]
1294                         [                            ]
1295                         [ 28  4   8   16  20  28  32 ]
1296                         [                            ]
1297                         [ 32  20  4   8   28  32  16 ]
1298 (%i2) T9: zn_mult_table(36/4);
1299                              [ 1  2  4  5  7  8 ]
1300                              [                  ]
1301                              [ 2  4  8  1  5  7 ]
1302                              [                  ]
1303                              [ 4  8  7  2  1  5 ]
1304 (%o2)                        [                  ]
1305                              [ 5  1  2  7  8  4 ]
1306                              [                  ]
1307                              [ 7  5  1  8  4  2 ]
1308                              [                  ]
1309                              [ 8  7  5  4  2  1 ]
1310 (%i3) T36_4: matrixmap(lambda([x], solve_congruences([0,x],[4,9])), T9);
1311                           [ 28  20  4   32  16  8  ]
1312                           [                        ]
1313                           [ 20  4   8   28  32  16 ]
1314                           [                        ]
1315                           [ 4   8   16  20  28  32 ]
1316 (%o3)                     [                        ]
1317                           [ 32  28  20  16  8   4  ]
1318                           [                        ]
1319                           [ 16  32  28  8   4   20 ]
1320                           [                        ]
1321                           [ 8   16  32  4   20  28 ]
1322 @end example
1324 @opencatbox{Categories:}
1325 @category{Number theory}
1326 @closecatbox
1327 @end deffn
1329 @c -----------------------------------------------------------------------------
1330 @anchor{zn_nth_root}
1331 @deffn {Function} zn_nth_root @
1332 @fname{zn_nth_root} (@var{x}, @var{n}, @var{m})  @
1333 @fname{zn_nth_root} (@var{x}, @var{n}, @var{m}, [[@var{p1}, @var{e1}], @dots{}, [@var{pk}, @var{ek}]])
1335 Returns a list with all @var{n}-th roots of @var{x} from the multiplication 
1336 subgroup of (Z/@var{m}Z) which contains @var{x}, or @code{false}, if @var{x} 
1337 is no @var{n}-th power modulo @var{m} or not contained in any multiplication 
1338 subgroup of (Z/@var{m}Z).
1340 @var{x} is an element of a multiplication subgroup modulo @var{m}, if the 
1341 greatest common divisor @code{g = gcd(x,m)} is coprime to @code{m/g}.
1343 @code{zn_nth_root} is based on an algorithm by Adleman, Manders and Miller 
1344 and on theorems about modulo multiplication groups by Daniel Shanks.
1346 The algorithm needs a prime factorization of the modulus @var{m}. 
1347 So in case the factorization of @var{m} is known, the list of factors 
1348 can be passed as the fourth argument. This optional argument
1349 must be of the same form as the list returned by @code{ifactors(m)} 
1350 using the default option @code{factors_only: false}.
1353 Examples:
1355 A power table of the multiplication group modulo @code{14} 
1356 followed by a list of lists containing all @var{n}-th roots of @code{1} 
1357 with @var{n} from @code{1} to @code{6}.
1359 @example
1360 (%i1) zn_power_table(14);
1361                          [ 1   1   1   1   1   1 ]
1362                          [                       ]
1363                          [ 3   9   13  11  5   1 ]
1364                          [                       ]
1365                          [ 5   11  13  9   3   1 ]
1366 (%o1)                    [                       ]
1367                          [ 9   11  1   9   11  1 ]
1368                          [                       ]
1369                          [ 11  9   1   11  9   1 ]
1370                          [                       ]
1371                          [ 13  1   13  1   13  1 ]
1372 (%i2) makelist(zn_nth_root(1,n,14), n,1,6);
1373 (%o2)  [[1], [1, 13], [1, 9, 11], [1, 13], [1], [1, 3, 5, 9, 11, 13]]
1374 @end example
1376 In the following example @var{x} is not coprime to @var{m}, 
1377 but is a member of a multiplication subgroup of (Z/@var{m}Z) 
1378 and any @var{n}-th root is a member of the same subgroup.
1380 The residue class @code{3} is no member of any multiplication subgroup of (Z/63Z) 
1381 and is therefore not returned as a third root of @code{27}.
1383 Here @code{zn_power_table} shows all residues @code{x} in (Z/63Z) 
1384 with @code{gcd(x,63) = 9}. This subgroup is isomorphic to (Z/7Z)*  
1385 and its identity @code{36} is computed via CRT.
1387 @example
1388 (%i1) m: 7*9$
1390 (%i2) zn_power_table(m,9);
1391                          [ 9   18  36  9   18  36 ]
1392                          [                        ]
1393                          [ 18  9   36  18  9   36 ]
1394                          [                        ]
1395                          [ 27  36  27  36  27  36 ]
1396 (%o2)                    [                        ]
1397                          [ 36  36  36  36  36  36 ]
1398                          [                        ]
1399                          [ 45  9   27  18  54  36 ]
1400                          [                        ]
1401                          [ 54  18  27  9   45  36 ]
1402 (%i3) zn_nth_root(27,3,m);
1403 (%o3)                           [27, 45, 54]
1404 (%i4) id7:1$  id63_9: solve_congruences([id7,0],[7,9]);
1405 (%o5)                                36
1406 @end example
1408 In the following RSA-like example, where the modulus @code{N} is squarefree, 
1409 i.e. it splits into 
1410 exclusively first power factors, every @code{x} from @code{0} to @code{N-1} 
1411 is contained in a multiplication subgroup.
1413 The process of decryption needs the @code{e}-th root. 
1414 @code{e} is coprime to @code{totient(N)} and therefore the @code{e}-th root is unique. 
1415 In this case @code{zn_nth_root} effectively performs CRT-RSA. 
1416 (Please note that @code{flatten} removes braces but no solutions.)
1418 @example
1419 (%i1) [p,q,e]: [5,7,17]$  N: p*q$
1421 (%i3) xs: makelist(x,x,0,N-1)$
1423 (%i4) ys: map(lambda([x],power_mod(x,e,N)),xs)$
1425 (%i5) zs: flatten(map(lambda([y], zn_nth_root(y,e,N)), ys))$
1427 (%i6) is(zs = xs);
1428 (%o6)                             true
1429 @end example
1431 In the following example the factorization of the modulus is known 
1432 and passed as the fourth argument.
1434 @example
1435 (%i1) p: 2^107-1$  q: 2^127-1$  N: p*q$
1437 (%i4) ibase: obase: 16$
1439 (%i5) msg: 11223344556677889900aabbccddeeff$
1441 (%i6) enc: power_mod(msg, 10001, N);
1442 (%o6)    1a8db7892ae588bdc2be25dd5107a425001fe9c82161abc673241c8b383
1443 (%i7) zn_nth_root(enc, 10001, N, [[p,1],[q,1]]);
1444 (%o7)               [11223344556677889900aabbccddeeff]
1445 @end example
1447 @opencatbox{Categories:}
1448 @category{Number theory}
1449 @closecatbox
1450 @end deffn
1452 @c -----------------------------------------------------------------------------
1453 @anchor{zn_order}
1454 @deffn {Function} zn_order @
1455 @fname{zn_order} (@var{x}, @var{n})  @
1456 @fname{zn_order} (@var{x}, @var{n}, [[@var{p1}, @var{e1}], @dots{}, [@var{pk}, @var{ek}]])
1458 Returns the order of @var{x} if it is a unit of the finite group (Z/@var{n}Z)*
1459 or returns @code{false}.  @var{x} is a unit modulo @var{n} if it is coprime to @var{n}.
1461 The applied algorithm needs a prime factorization of @code{totient(n)}. This factorization 
1462 might be time consuming in some cases and it can be useful to factor first 
1463 and then to pass the list of factors to @code{zn_log} as the third argument. 
1464 The list must be of the same form as the list returned by @code{ifactors(totient(n))} 
1465 using the default option @code{factors_only : false}.
1467 See also @mrefcomma{zn_primroot}  @mrefcomma{ifactors}  @mrefdot{totient}
1469 Examples:
1471 @code{zn_order} computes the order of the unit @var{x} in (Z/@var{n}Z)*.
1473 @example
1474 (%i1) n: 22$
1475 (%i2) g: zn_primroot(n);
1476 (%o2)                            7
1477 (%i3) units_22: sublist(makelist(i,i,1,21), lambda([x], gcd(x,n)=1));
1478 (%o3)          [1, 3, 5, 7, 9, 13, 15, 17, 19, 21]
1479 (%i4) (ord_7: zn_order(7, n)) = totient(n);
1480 (%o4)                        10 = 10
1481 (%i5) powers_7: makelist(power_mod(g,i,n), i,0,ord_7 - 1);
1482 (%o5)          [1, 7, 5, 13, 3, 21, 15, 17, 9, 19]
1483 (%i6) map(lambda([x], zn_order(x, n)), powers_7);
1484 (%o6)          [1, 10, 5, 10, 5, 2, 5, 10, 5, 10]
1485 (%i7) map(lambda([x], ord_7/gcd(x,ord_7)), makelist(i,i,0,ord_7-1));
1486 (%o7)         [1, 10, 5, 10, 5, 2, 5, 10, 5, 10]
1487 (%i8) totient(totient(n));
1488 (%o8)                          4
1489 @end example
1491 The optional third argument must be of the same form as the list returned by 
1492 @code{ifactors(totient(n))}.
1494 @example
1495 (%i1) (p : 2^142 + 217, primep(p));
1496 (%o1)                         true
1497 (%i2) ifs: ifactors( totient(p) )$
1498 (%i3) g: zn_primroot(p, ifs);
1499 (%o3)                           3
1500 (%i4) is( (ord_3 : zn_order(g, p, ifs)) = totient(p) );
1501 (%o4)                         true
1502 (%i5) map(lambda([x], ord_3/zn_order(x,p,ifs)), makelist(i,i,2,15));
1503 (%o5)    [22, 1, 44, 10, 5, 2, 22, 2, 8, 2, 1, 1, 20, 1]
1504 @end example
1506 @opencatbox{Categories:}
1507 @category{Number theory}
1508 @closecatbox
1509 @end deffn
1511 @c -----------------------------------------------------------------------------
1512 @anchor{zn_power_table}
1513 @deffn {Function} zn_power_table @
1514 @fname{zn_power_table} (@var{n})  @
1515 @fname{zn_power_table} (@var{n}, @var{gcd})  @
1516 @fname{zn_power_table} (@var{n}, @var{gcd}, @var{max_exp})
1518 Without any optional argument @code{zn_power_table(n)} 
1519 shows a power table of all elements in (Z/@var{n}Z)* 
1520 which are all residue classes coprime to @var{n}. 
1521 The exponent loops from @code{1} to the greatest characteristic factor of 
1522 @code{totient(n)} (also known as Carmichael function or Carmichael lambda)
1523 and the table ends with a column of ones on the right side. 
1525 The optional second argument @var{gcd} allows to select powers of a specific 
1526 subset of (Z/@var{n}Z). If @var{gcd} is an integer, powers of all residue 
1527 classes @code{x} with @code{gcd(x,n) = }@var{gcd} are returned,
1528 i.e. the default value for @var{gcd} is @code{1}.   
1529 If @var{gcd} is set to @code{all}, the table contains powers of all elements 
1530 in (Z/@var{n}Z).
1532 If the optional third argument @var{max_exp} is given, the exponent loops from 
1533 @code{1} to @var{max_exp}. 
1535 See also @mrefcomma{zn_add_table}  @mrefdot{zn_mult_table}
1537 Examples:
1539 The default which is @var{gcd}@code{ = 1} allows to demonstrate and study basic 
1540 theorems of e.g. Fermat and Euler.
1542 The argument @var{gcd} allows to select subsets of (Z/@var{n}Z) and to study 
1543 multiplication subgroups and isomorphisms. 
1544 E.g. the groups @code{G10} and @code{G10_2} are under multiplication both 
1545 isomorphic to @code{G5}. @code{1} is the identity in @code{G5}. 
1546 So are @code{1} resp. @code{6} the identities in @code{G10} resp. @code{G10_2}. 
1547 There are corresponding mappings for primitive roots, n-th roots, etc..
1549 @example
1550 (%i1) zn_power_table(10);
1551                               [ 1  1  1  1 ]
1552                               [            ]
1553                               [ 3  9  7  1 ]
1554 (%o1)                         [            ]
1555                               [ 7  9  3  1 ]
1556                               [            ]
1557                               [ 9  1  9  1 ]
1558 (%i2) zn_power_table(10,2);
1559                               [ 2  4  8  6 ]
1560                               [            ]
1561                               [ 4  6  4  6 ]
1562 (%o2)                         [            ]
1563                               [ 6  6  6  6 ]
1564                               [            ]
1565                               [ 8  4  2  6 ]
1566 (%i3) zn_power_table(10,5);
1567 (%o3)                         [ 5  5  5  5 ]
1568 (%i4) zn_power_table(10,10);
1569 (%o4)                         [ 0  0  0  0 ]
1570 (%i5) G5: [1,2,3,4];
1571 (%o6)                          [1, 2, 3, 4]
1572 (%i6) G10_2: map(lambda([x], solve_congruences([0,x],[2,5])), G5);
1573 (%o6)                          [6, 2, 8, 4]
1574 (%i7) G10: map(lambda([x], power_mod(3, zn_log(x,2,5), 10)), G5);
1575 (%o7)                          [1, 3, 7, 9]
1576 @end example
1578 If @var{gcd} is set to @code{all}, the table contains powers of all elements 
1579 in (Z/@var{n}Z).
1581 The third argument @var{max_exp} allows to set the highest exponent. 
1582 The following table shows a very small example of RSA.
1584 @example
1585 (%i1) N:2*5$ phi:totient(N)$ e:7$ d:inv_mod(e,phi)$
1587 (%i5) zn_power_table(N, all, e*d);
1588       [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0 ]
1589       [                                                               ]
1590       [ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1 ]
1591       [                                                               ]
1592       [ 2  4  8  6  2  4  8  6  2  4  8  6  2  4  8  6  2  4  8  6  2 ]
1593       [                                                               ]
1594       [ 3  9  7  1  3  9  7  1  3  9  7  1  3  9  7  1  3  9  7  1  3 ]
1595       [                                                               ]
1596       [ 4  6  4  6  4  6  4  6  4  6  4  6  4  6  4  6  4  6  4  6  4 ]
1597 (%o5) [                                                               ]
1598       [ 5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5 ]
1599       [                                                               ]
1600       [ 6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6 ]
1601       [                                                               ]
1602       [ 7  9  3  1  7  9  3  1  7  9  3  1  7  9  3  1  7  9  3  1  7 ]
1603       [                                                               ]
1604       [ 8  4  2  6  8  4  2  6  8  4  2  6  8  4  2  6  8  4  2  6  8 ]
1605       [                                                               ]
1606       [ 9  1  9  1  9  1  9  1  9  1  9  1  9  1  9  1  9  1  9  1  9 ]
1607 @end example
1609 @opencatbox{Categories:}
1610 @category{Number theory}
1611 @closecatbox
1612 @end deffn
1614 @c -----------------------------------------------------------------------------
1615 @anchor{zn_primroot}
1616 @deffn {Function} zn_primroot @
1617 @fname{zn_primroot} (@var{n})  @
1618 @fname{zn_primroot} (@var{n}, [[@var{p1}, @var{e1}], @dots{}, [@var{pk}, @var{ek}]])
1620 If the multiplicative group (Z/@var{n}Z)* is cyclic, @code{zn_primroot} computes the 
1621 smallest primitive root modulo @var{n}.  (Z/@var{n}Z)* is cyclic if @var{n} is equal to 
1622 @code{2}, @code{4}, @code{p^k} or @code{2*p^k}, where @code{p} is prime and 
1623 greater than @code{2} and @code{k} is a natural number.  @code{zn_primroot} 
1624 performs an according pretest if the option variable @mref{zn_primroot_pretest}
1625 (default: @code{false}) is set to @code{true}.  In any case the computation is limited 
1626 by the upper bound @mrefdot{zn_primroot_limit}
1628 If (Z/@var{n}Z)* is not cyclic or if there is no primitive root up to 
1629 @code{zn_primroot_limit}, @code{zn_primroot} returns @code{false}.
1631 The applied algorithm needs a prime factorization of @code{totient(n)}. This factorization 
1632 might be time consuming in some cases and it can be useful to factor first 
1633 and then to pass the list of factors to @code{zn_log} as an additional argument. 
1634 The list must be of the same form as the list returned by @code{ifactors(totient(n))} 
1635 using the default option @code{factors_only : false}.
1637 See also @mrefcomma{zn_primroot_p}  @mrefcomma{zn_order}  @mrefcomma{ifactors}  @mrefdot{totient}
1639 Examples:
1641 @code{zn_primroot} computes the smallest primitive root modulo @var{n} or returns 
1642 @code{false}.
1644 @example
1645 (%i1) n : 14$
1646 (%i2) g : zn_primroot(n);
1647 (%o2)                               3
1648 (%i3) zn_order(g, n) = totient(n);
1649 (%o3)                             6 = 6
1650 (%i4) n : 15$
1651 (%i5) zn_primroot(n);
1652 (%o5)                             false
1653 @end example
1655 The optional second argument must be of the same form as the list returned by 
1656 @code{ifactors(totient(n))}.
1658 @example
1659 (%i1) (p : 2^142 + 217, primep(p));
1660 (%o1)                             true
1661 (%i2) ifs : ifactors( totient(p) )$
1662 (%i3) g : zn_primroot(p, ifs);
1663 (%o3)                               3
1664 (%i4) [time(%o2), time(%o3)];
1665 (%o4)                    [[15.556972], [0.004]]
1666 (%i5) is(zn_order(g, p, ifs) = p - 1);
1667 (%o5)                             true
1668 (%i6) n : 2^142 + 216$
1669 (%i7) ifs : ifactors(totient(n))$
1670 (%i8) zn_primroot(n, ifs), 
1671       zn_primroot_limit : 200, zn_primroot_verbose : true;
1672 `zn_primroot' stopped at zn_primroot_limit = 200
1673 (%o8)                             false
1674 @end example
1676 @opencatbox{Categories:}
1677 @category{Number theory}
1678 @closecatbox
1679 @end deffn
1681 @c -----------------------------------------------------------------------------
1682 @anchor{zn_primroot_limit}
1683 @defvr {Option variable} zn_primroot_limit
1684 Default value: @code{1000} 
1686 If @mref{zn_primroot} cannot find a primitive root, it stops at this upper bound.
1687 If the option variable @mref{zn_primroot_verbose} (default: @code{false}) is 
1688 set to @code{true}, a message will be printed when @code{zn_primroot_limit} is reached. 
1690 @opencatbox{Categories:}
1691 @category{Number theory}
1692 @closecatbox
1693 @end defvr
1695 @c -----------------------------------------------------------------------------
1696 @anchor{zn_primroot_p}
1697 @deffn {Function} zn_primroot_p @
1698 @fname{zn_primroot_p} (@var{x}, @var{n})  @
1699 @fname{zn_primroot_p} (@var{x}, @var{n}, [[@var{p1}, @var{e1}], @dots{}, [@var{pk}, @var{ek}]])
1701 Checks whether @var{x} is a primitive root in the multiplicative group (Z/@var{n}Z)*.
1703 The applied algorithm needs a prime factorization of @code{totient(n)}. This factorization 
1704 might be time consuming and in case @code{zn_primroot_p} will be consecutively 
1705 applied to a list of candidates it can be useful to factor first and then to 
1706 pass the list of factors to @code{zn_log} as a third argument. 
1707 The list must be of the same form as the list returned by @code{ifactors(totient(n))} 
1708 using the default option @code{factors_only : false}.
1710 See also @mrefcomma{zn_primroot}  @mrefcomma{zn_order}  @mrefcomma{ifactors}  @mrefdot{totient}
1712 Examples:
1714 @code{zn_primroot_p} as a predicate function.
1716 @example
1717 (%i1) n : 14$
1718 (%i2) units_14 : sublist(makelist(i,i,1,13), lambda([i], gcd(i, n) = 1));
1719 (%o2)                     [1, 3, 5, 9, 11, 13]
1720 (%i3) zn_primroot_p(13, n);
1721 (%o3)                            false
1722 (%i4) sublist(units_14, lambda([x], zn_primroot_p(x, n)));
1723 (%o4)                            [3, 5]
1724 (%i5) map(lambda([x], zn_order(x, n)), units_14);
1725 (%o5)                      [1, 6, 6, 3, 3, 2]
1726 @end example
1728 The optional third argument must be of the same form as the list returned by 
1729 @code{ifactors(totient(n))}.
1731 @example
1732 (%i1) (p: 2^142 + 217, primep(p));
1733 (%o1)                         true
1734 (%i2) ifs: ifactors( totient(p) )$
1735 (%i3) sublist(makelist(i,i,1,50), lambda([x], zn_primroot_p(x,p,ifs)));
1736 (%o3)  [3, 12, 13, 15, 21, 24, 26, 27, 29, 33, 38, 42, 48]
1737 (%i4) [time(%o2), time(%o3)];
1738 (%o4)                  [[7.748484], [0.036002]]
1739 @end example
1741 @opencatbox{Categories:}
1742 @category{Predicate functions}
1743 @category{Number theory}
1744 @closecatbox
1745 @end deffn
1747 @c -----------------------------------------------------------------------------
1748 @anchor{zn_primroot_pretest}
1749 @defvr {Option variable} zn_primroot_pretest
1750 Default value: @code{false} 
1752 The multiplicative group (Z/@var{n}Z)* is cyclic if @var{n} is equal to 
1753 @code{2}, @code{4}, @code{p^k} or @code{2*p^k}, where @code{p} is prime and 
1754 greater than @code{2} and @code{k} is a natural number.  
1756 @code{zn_primroot_pretest} controls whether @mref{zn_primroot} will check 
1757 if one of these cases occur before it computes the smallest primitive root. 
1758 Only if @code{zn_primroot_pretest} is set to @code{true} this pretest will be 
1759 performed.
1761 @opencatbox{Categories:}
1762 @category{Number theory}
1763 @closecatbox
1764 @end defvr
1766 @c -----------------------------------------------------------------------------
1767 @anchor{zn_primroot_verbose}
1768 @defvr {Option variable} zn_primroot_verbose
1769 Default value: @code{false} 
1771 Controls whether @mref{zn_primroot} prints a message when reaching 
1772 @mrefdot{zn_primroot_limit}
1774 @opencatbox{Categories:}
1775 @category{Number theory}
1776 @closecatbox
1777 @end defvr