Funciones y variables para teoría de números
============================================

 -- Función: bern (N)
     Devuelve el N-ésimo número de Bernoulli del entero N. Los números de Bernoulli iguales a cero son suprimidos si zerobern vale false.

     Véase también burn.

          (%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

 -- Función: bernpoly (X, N)
     Devuelve el N-ésimo polinomio de Bernoulli de variable X.

 -- Función: bfzeta (S, N)
     Devuelve la función zeta de Riemann para el argumento S. El valor que devuelve es del tipo "big float" (bfloat) y N  es su número de dígitos.

     Es necesario cargar en memoria esta función haciendo load ("bffac").

 -- Función: bfhzeta (S, H, N)
     Devuelve la función zeta de Hurwitz para los argumentos S y H. El valor que devuelve es del tipo "big float" (bfloat) y N  es su número de dígitos.

     La función zeta de Hurwitz se define como

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

     Es necesario cargar en memoria esta función haciendo load ("bffac").

 -- Función: binomial (X, Y)
     Es el coeficiente binomial X!/(Y! (X - Y)!). Si X y Y son enteros, entonces se calcula el valor numérico del coeficiente binomial. Si Y o X - Y son enteros, el coeficiente binomial se expresa como un polinomio.

     Ejemplos:

          (%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)

 -- Función: burn (N)
     Devuelve el N-ésimo número de Bernoulli del entero N. La función burn puede ser más eficiente que bern para N grande (mayor que 105, por ejemplo), pues bern calcula todos los números de Bernoulli hasta N antes de devolver el resultado.

     La función burn se beneficia del hecho de que los números racionales de Bernoulli pueden aproximarse por funciones zeta con una eficiencia aceptable.

     Es necesario cargar en memoria esta función haciendo load ("bffac").

 -- Función: cf (EXPR)
     Transforma EXPR a fracciones continuas. La expresión EXPR debe contener fracciones continuas y raíces cuadradas de números enteros. Los operandos de la expresión pueden combinarse con operadores aritméticos. Además de fracciones continuas y raíces cuadradas, los factores de la expresión deben ser enteros o números racionales. Maxima no tiene más conocimiento sobre operaciones con fracciones continuas que el que aporta la función cf.

     La función cf evalúa sus argumentos después de asignar a la variable listarith el valor false, retornando una fracción continua en forma de lista.

     Una fracción continua a + 1/(b + 1/(c + ...)) se representa como la lista [a, b, c, ...], donde los elementos a, b, c, ... se evalúan como enteros. La expresión EXPR puede contener también  sqrt (n) donde n es un entero; en tal caso, cf devolverá tantos términos de la fracción continua como indique el valor de la variable cflength multiplicado por el período.

     Una fracción continua puede reducirse a un número evaluando la representación aritmética que devuelve cfdisrep. Véase también  cfexpand, que es otra alternativa para evaluar fracciones continuas.

     Véanse asimismo cfdisrep, cfexpand y cflength.

     Ejemplos:

       * La expresión EXPR contiene fracciones continuas y raíces cuadradas de enteros.

               (%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]

       * La variable cflength controla cuantos períodos de la fracción continua se calculan para números irracionales algebraicos.

               (%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]

       * Una fracción continua puede calcularse evaluando la representación aritmética que devuelve cfdisrep.

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

       * Maxima no sabe sobre operaciones con fracciones continuas más de lo que aporta la función cf.

               (%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]

 -- Función: cfdisrep (LISTA)
     Construye y devuelve una expresión aritmética ordinaria de la forma  a + 1/(b + 1/(c + ...)) a partir de la representación en formato lista de la fracción continua  [a, b, c, ...].

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

 -- Función: cfexpand (X)
     Devuelve la matriz con los numeradores y denominadores de la última (columna 1) y penúltima (columna 2) convergentes de la fracción continua X.

          (%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

 -- Variable opcional: cflength
     Valor por defecto: 1

     La variable cflength controla el número de términos de la fracción continua que devuelve la función cf, que será cflength multiplicado por el período. Así, el valor por defecto será el de un período.

          (%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]

 -- Función: divsum (N, K)
 -- Función: divsum (N)

     La llamada divsum (N, K) devuelve la suma de los divisores de N elevados a la K-ésima potencia.

     La llamada divsum (N) devuelve la suma de los divisores de N.

          (%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

 -- Función: euler (N)
     Devuelve el N-ésimo número de Euler del entero no negativo N.

     Para la constante de Euler-Mascheroni, véase %gamma.

          (%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]

 -- Constante: %gamma
     Es la constante de Euler-Mascheroni, 0.5772156649015329 ....

 -- Función: factorial (X)
     Representa la función factorial. Maxima considera factorial (X) lo mismo que X!. Véase !.

 -- Función: fib (N)
     Devuelve el N-ésimo número de Fibonacci. La llamada fib(0) es igual a 0, fib(1) devuelve 1 y fib (-N) es igual a (-1)^(N + 1) * fib(N).

     Después de llamar a fib, la variable prevfib toma el valor fib (X - 1), que es el número de Fibonacci que precede al último calculado.

          (%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]

 -- Función: fibtophi (EXPR)
     Expresa los números de Fibonacci en EXPR en términos de la razón áurea %phi, que es (1 + sqrt(5))/2, aproximadamente 1.61803399.

     Ejemplos:

          (%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) fibtophi (%);
                      n + 1             n + 1       n             n
                  %phi      - (1 - %phi)        %phi  - (1 - %phi)
          (%o3) - --------------------------- + -------------------
                          2 %phi - 1                2 %phi - 1
                                                    n - 1             n - 1
                                                %phi      - (1 - %phi)
                                              + ---------------------------
                                                        2 %phi - 1
          (%i4) ratsimp (%);
          (%o4)                           0

 -- Función: ifactors (N)
     Devuelve la factorización del argumento N, siendo éste un número entero positivo. Si n=p1^e1..pk^nk es la descomposición de N en números primos, ifactors devuelve [[p1, e1], ... , [pk, ek]].

     Los métodos de factorización son las divisones tentativas hasta el 9973, así como los  métodos rho de Pollard y el de la curva elíptica.

          (%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

 -- Función: inrt (X, N)
     Devuelve la raíz entera N-ésima del valor absoluto de  X.

          (%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]

 -- Función: inv_mod (N, M)
     Calcula el inverso de N módulo M. La llamada inv_mod (n,m) devuelve false si N es un divisor nulo módulo M.

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

 -- Función: jacobi (P, Q)
     Devuelve el símbolo de Jacobi para P y Q.

          (%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]

 -- Función: lcm (EXPR_1, ..., EXPR_N)
     Devuelve el mínimo común múltiplo de sus argumentos. Los argumentos pueden ser tanto expresiones en general como enteros.

     Es necesario cargar en memoria esta función haciendo load ("functs").

 -- Función: next_prime (N)
     Devuelve el menor de los primos mayores que N.

          (%i1) next_prime(27);
          (%o1)                       29

 -- Función: partfrac (EXPR, VAR)
     Expande la expresión EXPR en fracciones parciales respecto de la variable principal VAR. La función partfrac hace una descomposición completa en fracciones parciales. El algoritmo que se utiliza se basa en el hecho de que los denominadores de la expansión en fracciones parciales (los factores del denominador original) son primos relativos. Los numeradores se pueden escribir como combinaciones lineales de los denominadores.

          (%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)

 -- Función: power_mod (A, N, M)
     Utiliza un algoritmo modular para calcular a^n mod m, siendo A y N enteros cualesquiera y M un entero positivo. Si N es negativo, se utilizará inv_mod para encontrar el inverso modular.

          (%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

 -- Función: primep (N)
     Comprueba si el número entero N es o no primo, devolviendo true o false según el caso.

     Cuando el resultado de primep (n) es false, N es un número compuesto, y si es true, N es primo con alta probabilidad.

     Si N es menor que 341550071728321, se utiliza una versión determinística de la prueba de Miller-Rabin. En tal caso, si primep (n) devuelve true, entonces N es un número primo.

     Para N mayor que 34155071728321 primep realiza un número de pruebas de seudo-primalidad de Miller-Rabin igual a primep_number_of_tests y una prueba de seudo-primalidad de Lucas. La probabilidad de que N pase una prueba de Miller-Rabin es menor que 1/4. Con el valor por defecto de primep_number_of_tests, que es 25, la probabilidad de que N sea compuesto es menor que 10^-15.

 -- Variable opcional: primep_number_of_tests
     Valor por defecto: 25

     Número de pruebas de Miller-Rabin a realizar por primep.

 -- Función: prev_prime (N)
     Devuelve el mayor de los primos menores que N.

          (%i1) prev_prime(27);
          (%o1)                       23

 -- Función: qunit (N)
     Devuelve la unidad principal de sqrt (N), siendo N un entero; consiste en la resolución de la ecuación de Pell  a^2 - N b^2 = 1.

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

 -- Función: totient (N)
     Devuelve el número de enteros menores o iguales a N que son primos relativos con N.

 -- Variable opcional: zerobern
     Valor por defecto: true

     Si zerobern vale false, bern excluye los números de Bernoulli iguales a cero. See bern.

 -- Función: zeta (N)
     Devuelve la función zeta de Riemann para X entero negativo, 0, 1 o número par positivo. No se evalúa zeta (N) para cualesquiera otros argumentos, incluyendo racionales no enteros, números en coma flotante o argumentos complejos.

     Véanse también bfzeta y zeta%pi.

          (%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

 -- Variable opcional: zeta%pi
     Valor por defecto: true

     Si zeta%pi vale true, zeta devuelve una expresión proporcional a %pi^n si n es un número par positivo. En caso contrario, zeta no se evalúa y devuelve la forma nominal zeta (n).

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