@c end concepts Number Theory
@menu
* Definitions for Number Theory::  
@end menu

@node Definitions for Number Theory,  , Number Theory, Number Theory
@section Definitions for Number Theory

@deffn {Function} bern (@var{n})
Returns the @var{n}'th Bernoulli number for integer @var{n}.
@c WELL, ACTUALLY bern SIMPLIFIES, LIKE FACTORIAL -- DO WE WANT TO GET INTO THAT ???
@c OR JUST PRETEND IT'S "RETURNED" ???
Bernoulli numbers equal to zero are suppressed if @code{zerobern} is @code{false}.

See also @code{burn}.

@example
(%i1) zerobern: true$
(%i2) map (bern, [0, 1, 2, 3, 4, 5, 6, 7, 8]);
                  1  1       1      1        1
(%o2)       [1, - -, -, 0, - --, 0, --, 0, - --]
                  2  6       30     42       30
(%i3) zerobern: false$
(%i4) map (bern, [0, 1, 2, 3, 4, 5, 6, 7, 8]);
            1  1    1   5     691   7    3617  43867
(%o4) [1, - -, -, - --, --, - ----, -, - ----, -----]
            2  6    30  66    2730  6    510    798
@end example

@end deffn

@deffn {Function} bernpoly (@var{x}, @var{n})
Returns the @var{n}'th Bernoulli polynomial in the
variable @var{x}.

@end deffn

@deffn {Function} bfzeta (@var{s}, @var{n})
Returns the Riemann zeta function for the argument @var{s}.
The return value is a big float (bfloat);
@var{n} is the number of digits in the return value.

@code{load ("bffac")} loads this function.

@end deffn

@deffn {Function} bfhzeta (@var{s}, @var{h}, @var{n})
Returns the Hurwitz zeta function for the arguments @var{s} and @var{h}.
The return value is a big float (bfloat);
@var{n} is the number of digits in the return value.

The Hurwitz zeta function is defined as

@example
sum ((k+h)^-s, k, 0, inf)
@end example

@code{load ("bffac")} loads this function.

@end deffn

@deffn {Function} binomial (@var{x}, @var{y})
The binomial coefficient @code{@var{x}!/(@var{y}! (@var{x} - @var{y})!)}.
If @var{x} and @var{y} are integers, then the numerical value of the binomial
coefficient is computed.
If @var{y}, or @var{x - y}, is an integer,
the binomial coefficient is expressed as a polynomial.

Examples:

@c ===beg===
@c binomial (11, 7);
@c 11! / 7! / (11 - 7)!;
@c binomial (x, 7);
@c binomial (x + 7, x);
@c binomial (11, y);
@c ===end===
@example
(%i1) binomial (11, 7);
(%o1)                          330
(%i2) 11! / 7! / (11 - 7)!;
(%o2)                          330
(%i3) binomial (x, 7);
        (x - 6) (x - 5) (x - 4) (x - 3) (x - 2) (x - 1) x
(%o3)   -------------------------------------------------
                              5040
(%i4) binomial (x + 7, x);
      (x + 1) (x + 2) (x + 3) (x + 4) (x + 5) (x + 6) (x + 7)
(%o4) -------------------------------------------------------
                               5040
(%i5) binomial (11, y);
(%o5)                    binomial(11, y)
@end example

@end deffn

@deffn {Function} burn (@var{n})
Returns the @var{n}'th Bernoulli number for integer @var{n}.
@code{burn} may be more efficient than @code{bern} for large, isolated @var{n}
(perhaps @var{n} greater than 105 or so), @c CLAIM MADE IN bffac.usg !!!
as @code{bern} computes all the Bernoulli numbers up to index @var{n} before returning.

@c STATEMENTS ABOUT TIMING NEED VERIFICATION !!!
@c CAN'T VERIFY NOW AS burn IS BROKEN IN 5.9.1 AND CVS BUILD AT PRESENT !!!
@c (BERN(402) takes about 645 secs vs 13.5 secs for BURN(402).
@c The time to compute @code{bern} is approximately exponential,
@c while the time to compute @code{burn} is approximately cubic.
@c But if next you do BERN(404), it only takes 12 secs,
@c since BERN remembers all in an array, whereas BURN(404) will take
@c maybe 14 secs or maybe 25, depending on whether Maxima needs to
@c BFLOAT a better value of %PI.)

@code{burn} exploits the observation that (rational) Bernoulli numbers can be
approximated by (transcendental) zetas with tolerable efficiency.

@code{load ("bffac")} loads this function.

@end deffn

@deffn {Function} cf (@var{expr})
Converts @var{expr} into a continued fraction.
@var{expr} is an expression
comprising continued fractions and square roots of integers.
Operands in the expression may be combined with arithmetic operators.
Aside from continued fractions and square roots,
factors in the expression must be integer or rational numbers.
Maxima does not know about operations on continued fractions outside of @code{cf}.

@code{cf} evaluates its arguments after binding @code{listarith} to @code{false}.
@code{cf} returns a continued fraction, represented as a list.

A continued fraction @code{a + 1/(b + 1/(c + ...))}
is represented by the list @code{[a, b, c, ...]}.
The list elements @code{a}, @code{b}, @code{c}, ... must evaluate to integers.
@var{expr} may also contain @code{sqrt (n)} where @code{n} is an integer.
In this case @code{cf} will give as many
terms of the continued fraction as the value of the variable
@code{cflength} times the period.

A continued fraction can be evaluated to a number
by evaluating the arithmetic representation
returned by @code{cfdisrep}.
See also @code{cfexpand} for another way to evaluate a continued fraction.

See also @code{cfdisrep}, @code{cfexpand}, and @code{cflength}.

Examples:

@itemize @bullet
@item
@var{expr} is an expression comprising continued fractions and square roots of integers.

@example
(%i1) cf ([5, 3, 1]*[11, 9, 7] + [3, 7]/[4, 3, 2]);
(%o1)               [59, 17, 2, 1, 1, 1, 27]
(%i2) cf ((3/17)*[1, -2, 5]/sqrt(11) + (8/13));
(%o2)        [0, 1, 1, 1, 3, 2, 1, 4, 1, 9, 1, 9, 2]
@end example

@item
@code{cflength} controls how many periods of the continued fraction
are computed for algebraic, irrational numbers.

@example
(%i1) cflength: 1$
(%i2) cf ((1 + sqrt(5))/2);
(%o2)                    [1, 1, 1, 1, 2]
(%i3) cflength: 2$
(%i4) cf ((1 + sqrt(5))/2);
(%o4)               [1, 1, 1, 1, 1, 1, 1, 2]
(%i5) cflength: 3$
(%i6) cf ((1 + sqrt(5))/2);
(%o6)           [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2]
@end example

@item
A continued fraction can be evaluated by evaluating the arithmetic representation
returned by @code{cfdisrep}.

@example
(%i1) cflength: 3$
(%i2) cfdisrep (cf (sqrt (3)))$
(%i3) ev (%, numer);
(%o3)                   1.731707317073171
@end example

@item
Maxima does not know about operations on continued fractions outside of @code{cf}.

@example
(%i1) cf ([1,1,1,1,1,2] * 3);
(%o1)                     [4, 1, 5, 2]
(%i2) cf ([1,1,1,1,1,2]) * 3;
(%o2)                  [3, 3, 3, 3, 3, 6]
@end example

@end itemize
@end deffn

@c NEEDS CLARIFICATION -- MAKE EXPLICIT HOW list IS RELATED TO a, b, c, ...
@c ALSO, CAN list CONTAIN ANYTHING OTHER THAN LITERAL INTEGERS ??
@deffn {Function} cfdisrep (@var{list})
Constructs and returns an ordinary arithmetic expression
of the form @code{a + 1/(b + 1/(c + ...))}
from the list representation of a continued fraction @code{[a, b, c, ...]}.

@example
(%i1) cf ([1, 2, -3] + [1, -2, 1]);
(%o1)                     [1, 1, 1, 2]
(%i2) cfdisrep (%);
                                  1
(%o2)                     1 + ---------
                                    1
                              1 + -----
                                      1
                                  1 + -
                                      2
@end example

@end deffn

@deffn {Function} cfexpand (@var{x})
Returns a matrix of the numerators and denominators of the
last (column 1) and next-to-last (column 2) convergents of the continued fraction @var{x}.

@example
(%i1) cf (rat (ev (%pi, numer)));

`rat' replaced 3.141592653589793 by 103993//33102 = 3.141592653011902
(%o1)                  [3, 7, 15, 1, 292]
(%i2) cfexpand (%); 
                         [ 103993  355 ]
(%o2)                    [             ]
                         [ 33102   113 ]
(%i3) %[1,1]/%[2,1], numer;
(%o3)                   3.141592653011902
@end example

@end deffn

@defvr {Option variable} cflength
Default value: 1

@code{cflength} controls the number of terms of the continued
fraction the function @code{cf} will give, as the value @code{cflength} times the
period.  Thus the default is to give one period.

@example
(%i1) cflength: 1$
(%i2) cf ((1 + sqrt(5))/2);
(%o2)                    [1, 1, 1, 1, 2]
(%i3) cflength: 2$
(%i4) cf ((1 + sqrt(5))/2);
(%o4)               [1, 1, 1, 1, 1, 1, 1, 2]
(%i5) cflength: 3$
(%i6) cf ((1 + sqrt(5))/2);
(%o6)           [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2]
@end example

@end defvr

@deffn {Function} divsum (@var{n}, @var{k})
@deffnx {Function} divsum (@var{n})

@code{divsum (@var{n}, @var{k})} returns the sum of the divisors of @var{n}
raised to the @var{k}'th power.

@code{divsum (@var{n})} returns the sum of the divisors of @var{n}.

@example
(%i1) divsum (12);
(%o1)                          28
(%i2) 1 + 2 + 3 + 4 + 6 + 12;
(%o2)                          28
(%i3) divsum (12, 2);
(%o3)                          210
(%i4) 1^2 + 2^2 + 3^2 + 4^2 + 6^2 + 12^2;
(%o4)                          210
@end example

@end deffn

@deffn {Function} euler (@var{n})
Returns the @var{n}'th Euler number for nonnegative integer @var{n}.

For the Euler-Mascheroni constant, see @code{%gamma}.

@example
(%i1) map (euler, [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]);
(%o1)    [1, 0, - 1, 0, 5, 0, - 61, 0, 1385, 0, - 50521]
@end example

@end deffn

@defvr {Constant} %gamma
The Euler-Mascheroni constant, 0.5772156649015329 ....
@c DOUBTLESS THERE IS MORE TO SAY HERE.

@end defvr

@deffn {Function} factorial (@var{x})
Represents the factorial function. Maxima treats @code{factorial (@var{x})} the same as @code{@var{x}!}.
See @code{!}.

@end deffn

@deffn {Function} fib (@var{n})
Returns the @var{n}'th Fibonacci number.
@code{fib(0)} equal to 0 and @code{fib(1)} equal to 1,
and
@code{fib (-@var{n})} equal to @code{(-1)^(@var{n} + 1) * fib(@var{n})}.

After calling @code{fib},
@code{prevfib} is equal to @code{fib (@var{x} - 1)},
the Fibonacci number preceding the last one computed.

@example
(%i1) map (fib, [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]);
(%o1)         [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55]
@end example

@end deffn

@deffn {Function} fibtophi (@var{expr})
Expresses Fibonacci numbers in terms of the constant @code{%phi},
which is @code{(1 + sqrt(5))/2}, approximately 1.61803399.

@c SEEMS MISPLACED -- THIS WANTS TO BE UNDER %phi !!!
By default, Maxima does not know about @code{%phi}.
After executing @code{tellrat (%phi^2 - %phi - 1)} and @code{algebraic: true},
@code{ratsimp} can simplify some expressions containing @code{%phi}.

@example
(%i1) fibtophi (fib (n));
                           n             n
                       %phi  - (1 - %phi)
(%o1)                  -------------------
                           2 %phi - 1
(%i2) fib (n-1) + fib (n) - fib (n+1);
(%o2)          - fib(n + 1) + fib(n) + fib(n - 1)
(%i3) ratsimp (fibtophi (%));
(%o3)                           0
@end example

@end deffn

@deffn {Function} ifactors (@var{n})
For a positive integer @var{n} returns the factorization of @var{n}. If
@code{n=p1^e1..pk^nk} is the decomposition of @var{n} into prime
factors, ifactors returns @code{[[p1, e1], ... , [pk, ek]]}.

Factorization methods used are trial divisions by primes up to 9973,
Pollard's rho method and elliptic curve method.

@example
(%i1) ifactors(51575319651600);
(%o1)     [[2, 4], [3, 2], [5, 2], [1583, 1], [9050207, 1]]
(%i2) apply("*", map(lambda([u], u[1]^u[2]), %));
(%o2)                        51575319651600
@end example

@end deffn

@deffn {Function} inrt (@var{x}, @var{n})
Returns the integer @var{n}'th root of the absolute value of @var{x}.

@example
(%i1) l: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]$
(%i2) map (lambda ([a], inrt (10^a, 3)), l);
(%o2) [2, 4, 10, 21, 46, 100, 215, 464, 1000, 2154, 4641, 10000]
@end example

@end deffn

@deffn {Function} inv_mod (@var{n}, @var{m})
Computes the inverse of @var{n} modulo @var{m}. 
@code{inv_mod (n,m)} returns @code{false}, 
if @var{n} is a zero divisor modulo @var{m}.

@example
(%i1) inv_mod(3, 41);
(%o1)                           14
(%i2) ratsimp(3^-1), modulus=41;
(%o2)                           14
(%i3) inv_mod(3, 42);
(%o3)                          false
@end example

@end deffn

@deffn {Function} jacobi (@var{p}, @var{q})
Returns the Jacobi symbol of @var{p} and @var{q}.

@example
(%i1) l: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]$
(%i2) map (lambda ([a], jacobi (a, 9)), l);
(%o2)         [1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0]
@end example

@end deffn

@deffn {Function} lcm (@var{expr_1}, ..., @var{expr_n})
Returns the least common multiple of its arguments.
The arguments may be general expressions as well as integers.

@code{load ("functs")} loads this function.

@end deffn

@deffn {Function} minfactorial (@var{expr})
Examines @var{expr} for occurrences of two factorials
which differ by an integer.
@code{minfactorial} then turns one into a polynomial times the other.

@c I CAN'T TELL WHAT THIS IS SUPPOSED TO MEAN. !!!
@c minfactorial DOESN'T SEEM TO DO ANYTHING binomial DOESN'T DO BY ITSELF !!!
@c LOOKING AT THE minfactorial CODE DOESN'T HELP !!!
@c If exp involves binomial coefficients then they will be
@c converted into ratios of factorials.

@example
(%i1) n!/(n+2)!;
                               n!
(%o1)                       --------
                            (n + 2)!
(%i2) minfactorial (%);
                                1
(%o2)                    ---------------
                         (n + 1) (n + 2)
@end example

@end deffn

@deffn {Function} power_mod (@var{a}, @var{n}, @var{m})
Uses a modular algorithm to compute @code{a^n mod m} 
where @var{a} and @var{n} are integers and @var{m} is a positive integer. 
If @var{n} is negative, @code{inv_mod} is used to find the modular inverse.

@example
(%i1) power_mod(3, 15, 5);
(%o1)                          2
(%i2) mod(3^15,5);
(%o2)                          2
(%i3) power_mod(2, -1, 5);
(%o3)                          3
(%i4) inv_mod(2,5);
(%o4)                          3
@end example

@end deffn

@deffn {Function} next_prime (@var{n})
Returns the smallest prime bigger than @var{n}.

@example
(%i1) next_prime(27);
(%o1)                       29
@end example

@end deffn

@deffn {Function} partfrac (@var{expr}, @var{var})
Expands the expression @var{expr} in partial fractions
with respect to the main variable @var{var}.  @code{partfrac} does a complete
partial fraction decomposition.  The algorithm employed is based on
the fact that the denominators of the partial fraction expansion (the
factors of the original denominator) are relatively prime.  The
numerators can be written as linear combinations of denominators, and
the expansion falls out.

@example
(%i1) 1/(1+x)^2 - 2/(1+x) + 2/(2+x);
                      2       2        1
(%o1)               ----- - ----- + --------
                    x + 2   x + 1          2
                                    (x + 1)
(%i2) ratsimp (%);
                                 x
(%o2)                 - -------------------
                         3      2
                        x  + 4 x  + 5 x + 2
(%i3) partfrac (%, x);
                      2       2        1
(%o3)               ----- - ----- + --------
                    x + 2   x + 1          2
                                    (x + 1)
@end example
@end deffn

@c IS IT POSSIBLE TO MAKE A DECLARATION SUCH THAT primep RETURNS true ??
@deffn {Function} primep (@var{n})
Primality test. If @code{primep (n)} returns @code{false}, @var{n} is a
composite number and if it returns @code{true}, @var{n} is a prime number
with very high probability.

For @var{n} less than 341550071728321 a deterministic version of Miller-Rabin's
test is used. If @code{primep (n)} returns @code{true}, then @var{n} is a
prime number.

For @var{n} bigger than 34155071728321 @code{primep} uses
@code{primep_number_of_tests} Miller-Rabin's pseudo-primality tests
and one Lucas pseudo-primality test. The probability that @var{n} will
pass one Miller-Rabin test is less than 1/4. Using the default value 25 for
@code{primep_number_of_tests}, the probability of @var{n} beeing
composite is much smaller that 10^-15.

@end deffn

@defvr {Option variable} primep_number_of_tests
Default value: 25

Number of Miller-Rabin's tests used in @code{primep}.
@end defvr

@deffn {Function} prev_prime (@var{n})
Returns the greatest prime smaller than @var{n}.

@example
(%i1) prev_prime(27);
(%o1)                       23
@end example
@end deffn

@deffn {Function} qunit (@var{n})
Returns the principal unit of the real quadratic number field
@code{sqrt (@var{n})} where @var{n} is an integer,
i.e., the element whose norm is unity.
This amounts to solving Pell's equation @code{a^2 - @var{n} b^2 = 1}.

@example
(%i1) qunit (17);
(%o1)                     sqrt(17) + 4
(%i2) expand (% * (sqrt(17) - 4));
(%o2)                           1
@end example

@end deffn

@deffn {Function} totient (@var{n})
Returns the number of integers less than or equal to @var{n} which
are relatively prime to @var{n}.

@end deffn

@defvr {Option variable} zerobern
Default value: @code{true}

When @code{zerobern} is @code{false},
@code{bern} excludes the Bernoulli numbers which are equal to zero. 
See @code{bern}.

@end defvr

@deffn {Function} zeta (@var{n})
Returns the Riemann zeta function if @var{x} is a negative integer, 0, 1,
or a positive even number,
and returns a noun form @code{zeta (@var{n})} for all other arguments,
including rational noninteger, floating point, and complex arguments.

See also @code{bfzeta} and @code{zeta%pi}.

@example
(%i1) map (zeta, [-4, -3, -2, -1, 0, 1, 2, 3, 4, 5]);
                                     2              4
           1        1     1       %pi            %pi
(%o1) [0, ---, 0, - --, - -, inf, ----, zeta(3), ----, zeta(5)]
          120       12    2        6              90
@end example

@end deffn

@defvr {Option variable} zeta%pi
Default value: @code{true}

When @code{zeta%pi} is @code{true}, @code{zeta} returns an expression
proportional to @code{%pi^n} for even integer @code{n}.
Otherwise, @code{zeta} returns a noun form @code{zeta (n)}
for even integer @code{n}.

@example
(%i1) zeta%pi: true$
(%i2) zeta (4);
                                 4
                              %pi
(%o2)                         ----
                               90
(%i3) zeta%pi: false$
(%i4) zeta (4);
(%o4)                        zeta(4)
@end example

@end defvr