Algoritmo de Euclides binario

Esta entrada participa en la Edición 4.12 del  Carnaval de Matemáticas cuyo anfitrión es el blog High Ability Dimension.

En este blog nos motiva mucho el construir un esquema en hoja de cálculo que explique de la mejor forma posible el funcionamiento de un proceso o un algoritmo. Lo haremos hoy con la variante binaria del Algoritmo de Euclides.

Restar en lugar de dividir

Ya se sabe que el algoritmo de Euclides por divisiones se puede sustituir por otro efectuado a base de restar los dos números. Parece más lento, pero al ahorrar divisiones puede resultar más eficiente. Hace sesenta años lo usábamos los alumnos de Bachillerato (con once añitos) para comprobar si dos segmentos eran inconmensurables o no. Llevábamos uno sobre otro y nos quedábamos con la diferencia. Si al reiterar llegábamos al segmento nulo, es que tenían una medida común. Aquello era Geometría de Euclides pura.

Si el algoritmo clásico se organizaba a partir de la división euclídea, en esta variante se usa la resta. Así, para encontrar MCD(96,36) por divisiones sería: 96=2*36+24; 36=1*24+12; 24=2*12+0, luego el MCD(96,36)=12. Por restas formaríamos estas parejas: (96,36), (60,36), (36,24), (24,12), (12,12), (12,0), llegando al mismo resultado.

Eliminar el factor 2 siempre que se pueda.

En computación con sistema de numeración binario la división y la multiplicación por 2 se reducen a trasladar un lugar a la derecha o izquierda del dígito que se multiplica. Por eso, es interesante dar protagonismo al número 2 en los cálculos. Una primera idea es que si expresamos los dos números A y B de la forma A=2^m*p, B=2ⁿ*q, con p y q los mayores divisores impares, el M.C.D se puede encontrar con dos cálculos por separado. Por una parte quedándonos con el menor exponente del 2 y por otra hallando MCD(p,q).

En el algoritmo que estamos presentando, se elimina en primer lugar la potencia de 2 común a ambos números, y se toma nota de ella. Una vez eliminada, el factor 2 no va a influir en el resultado, y cuando aparezca en alguno de los dos números se podrá igualmente suprimir. En términos binarios suprimir un 2 es trasladar los dígitos una posición hacia la derecha.

En Wikipedia y otras páginas que hemos consultado expresan lo anterior en forma de tres reglas. http://en.wikipedia.org/wiki/Binary_GCD_algorithm. No son difíciles de razonar.

(1) Si A y B son ambos pares se cumple MCD(A,B)=2*MCD(A/2,B/2)

Esto justifica que el primer paso que demos sea el de separar las potencias de 2.

(2) Si A es par y B impar se tiene MCD(A;B)=MCD(A/2,B)

Igualmente se aplicaría si B es par y A impar. En virtud de esta regla eliminaremos todos los factores 2 una vez que se ha separado la potencia común.

(3) Si ambos A y B son impares y A<B, MCD(A,B)=MCD(B-A,A). Si es B<A, MCD(A,B)=MCD(A-B,B)

Esta es la esencia del algoritmo de Euclides, restar ambos números.

Eliminar la potencia de 2 común (Regla 1)

Hemos creado una demo en hoja de cálculo para seguir visualmente los pasos del algoritmo binario. Lo puedes descargar desde la dirección

http://hojamat.es/sindecimales/divisibilidad/herramientas/herrdiv.htm#euclibin

Como nuestro objetivo es visual, desarrollaremos ahora los pasos sugeridos pero en forma de esquema. En una hoja de cálculo se escribirán los dos números y con un botón iterativo se irán avanzando pasos hasta llegar al MCD. Esta primera fase de eliminar el factor común lo puedes ver en las imágenes siguientes:

En primer lugar hemos escrito los números 192 y 84. Al pulsar sobre el botón Aceptar números se han convertido ambos en binario y se han reconstruido a la derecha. La potencia de 2 común figura al principio con el valor 1.

Observamos que ambos números terminan en dos ceros (en binario), luego compartirán el factor 2²=4.
Según estamos explicando, el primer paso del algoritmo binario es eliminar el factor 2 común. Si ahora usamos el botón Paso a paso dos veces veremos que los dígitos de ambos números se mueven dos posiciones a la derecha y que la potencia de 2 se convierte en 4.

A la derecha figuran los números ya simplificados, 48 y 21, y la potencia de 2, que nos servirá al final.

Resto de pasos

Ahora la estrategia es triple:

(1) Si el primer número contiene aún potencias de 2, se eliminan (Regla 2)

Observa que en nuestro ejemplo el número 48 termina en cuatro ceros, luego al cabo de cuatro toques del botón Paso a paso desaparecerán.

(2) Se opera igual con el segundo: se le elimina el factor 2 (Regla 2)

En nuestro ejemplo ya no quedan factores 2 (por ahora)

(3) Si ambos son impares, se sustituye el mayor de ellos por la diferencia entre ambos.

Este es el núcleo del algoritmo de Euclides por diferencias. En la práctica quizás haya que intercambiar los valores de A y B, pero no entraremos en esos detalles.

Al restar dos impares se producirán nuevos factores 2, por lo que en la siguiente pasada del algoritmo los eliminará. Lo ves en las siguientes dos imágenes:

Reiteramos estas tres reglas hasta que lleguemos al valor 0. En ese momento el algoritmo construye el M.C.D. multiplicando el resultado por la potencia de 2 que tenía almacenada.
Aquí lo tienes:

Visto así, en hoja de cálculo, no parece ser nada del otro mundo, pero todas las operaciones que realiza son altamente eficientes en el sistema de numeración binario. Por algo lo introdujo el programador israelí Stein en 1967. Aquí sólo se nos queda como un tema de cultura matemática, pero es divertido implementarlo.

Versión recursiva

Lo que hemos desarrollado, con un botón que en cada paso reacciona según lo que se encuentra, nos permite sospechar que todo esto se resuelve también mediante una función recursiva. Es cierto, y está publicado. Lo que haremos aquí es adaptarlo al Basic de Excel y Calc:

Public Function mcdbin(m, n)
Dim mb

'Si son iguales, hemos llegado al MCD
If m = n Then
mb = m

'Si ambos son divisibles entre 2, se saca ese factor
ElseIf m / 2 = m \ 2 And n / 2 = n \ 2 Then
mb = 2 * mcdbin(m / 2, n / 2)

'El primero contiene un 2. Se elimina
ElseIf m / 2 = m \ 2 Then
mb = mcdbin(m / 2, n)

'Operamos de igual forma con el segundo
ElseIf n / 2 = n \ 2 Then
mb = mcdbin(m, n / 2)

'Ambos son impares. Se restan
Else
If m > n Then mb = mcdbin(m - n, n)
If n > m Then mb = mcdbin(m, n - m)
End If

'La función recoge el valor de mb
mcdbin = mb
End Function

Tiene toda la elegancia de las funciones recursivas
(ver http://hojaynumeros.blogspot.com.es/2012/03/funciones-recursivas-en-las-hojas-de.html)

En este caso resulta un poco complicada, pero funciona muy bien. Si te apetece, estudia la versión clásica y recursiva:

Public Function mcdeuclid(m, n)
Dim mb

'Si uno es múltiplo de otro, obtenemos el MCD
If m / n = m \ n Then
mb = n
ElseIf n / m = n \ m Then
mb = m
Else
'En caso contrario, se restan
If m > n Then mb = mcdeuclid(m - n, n)
If n > m Then mb = mcdeuclid(m, n - m)
End If
'La función recoge el valor de mb
mcdeuclid = mb
End Function

Ambas están implementadas en la herramienta que ofrecemos.

Funciones generatrices en Combinatoria– La teoría

En la entrada anterior presentábamos como un artificio sustituir los elementos de un producto cartesiano condicionado de conjuntos numéricos por polinomios en una indeterminada con exponentes idénticos a los números a combinar.

Hoy lo convertiremos en un desarrollo teórico.

Función generatriz de un conjunto numérico

Dado un conjunto ordenado de números reales o complejos (aquí usaremos casi exclusivamente los enteros positivos) {a₀, a₁, a₂, a₃,…a_n} llamaremos función generatriz (ordinaria) o generadora del mismo al polinomio de la forma

Si el conjunto tiene infinitos términos sustituiremos el polinomio por una serie de potencias, pero en este caso la igualdad

sólo tendrá sentido si dicha serie posee un radio de convergencia no nulo y la función está definida dentro de ese radio. No obstante, aquí  no trataremos la convergencia, sino las relaciones entre coeficientes. Si no converge, P(x) no será una función, pero las técnicas siguen valiendo.

Casos sencillos

El caso más sencillo de función generatriz es la que corresponde al conjunto {1,1,1,1,…1}
Si este es finito con n elementos, sabemos que su función generatriz se puede obtener mediante la fórmula de la suma de una progresión geométrica

Ya usamos esta fórmula en la entrada anterior

Si el conjunto es infinito {1,1,1,1,…}también es sencillo verificar que

Aquí nos encontramos con la potencia que tiene este método. Si derivamos miembro a miembro (omitimos detalles y también la cuestión de la convergencia) nos encontraremos con que

Por tanto esta es la función generatriz del conjunto {1,2,3,4,….}

La fórmula del binomio de Newton nos proporciona otro ejemplo de función generatriz. Los números combinatorios para un n dado tienen por función generatriz (1+x)ⁿ ya que

Podíamos seguir dando ejemplos hasta llenar páginas enteras, pero destacaremos especialmente dos técnicas

Manipulación algebraica

Con las técnicas algebraicas y sin plantearnos ahora el problema de la convergencia podemos encontrar la función generatriz de muchos conjuntos de coeficientes.

Por ejemplo, es fácil encontrarla para los números de Fibonacci F0, F1, F2,  F3,  F4…, de los que suponemos se conocen sus valores y propiedades. Observa estas manipulaciones:

F(x)=F₀+F₁x+F₂x²+ F₃x³+ F₄x⁴+…=
F₀+F₁x+(F₀+F₁)x²+(F₁+F₂)x³+(F₂+F₃)x⁴+… =
F₀+F₁x+(F₀ x²+ F₁ x³+ F₂ x⁴+…) + (F₁ x²+ F₂ x³+ F₃ x⁴+…)

Pero en los paréntesis se está reconstruyendo F(x) de alguna forma, por lo que podemos escribir (pon tú los detalles)

F(x)= F₀+F₁x+F(x) x²+(F(x)- F₀)x

Despejamos F(x), sustituimos F₀ por su valor 0 (a veces se toma 1) y F₁ por 1 y ya la tenemos:

Sólo hemos usado técnicas algebraicas sencillas. Más adelante comprobaremos este resultado.

Manipulación analítica

Si consideramos la derivación e integración formales podemos encontrar fácilmente funciones generatrices. Ya hemos considerado un ejemplo de derivación.

De la misma forma podemos usar la integración. Por ejemplo en la geométrica podemos integrar

Nos resultaría entonces

En cualquier manual puedes encontrar muchos ejemplos similares de este tipo de manipulación. No olvides que podemos mezclar las dos técnicas, analítica y algebraica, así como sumar, multiplicar y otras.

Problema inverso

Si la serie que define la función generatriz converge y conocemos esta, encontrar los coeficientes de la serie siempre es posible por la fórmula de McLaurin

Es un camino muy pesado, pero seguro. Sin embargo el problema contrario de dar los coeficientes y encontrar la expresión de P(x) quizás no lo puedas resolver. Es el clásico problema de la suma de series.

Con la ayuda de un ordenador se puede simplificar el proceso. Damos un ejemplo:

Antes vimos que los números de Fibonacci tenían como función generatriz (si se toma F0=0. A veces se toma F0=1 y entonces tiene una expresión ligeramente distinta)

Si tuviéramos posibilidad de desarrollar por McLaurin en algún lenguaje o programa, nos ahorraríamos mucho trabajo. Nosotros lo hemos hecho con el lenguaje PARI. Se entiende fácilmente aunque no se haya usado nunca:

{write("final.txt",taylor(x/(1-x-x^2), x,12))}

Ordenamos que escriba en el archivo “final.txt” 12 términos del desarrollo de Taylor (es en x=0) de la función dada. El resultado es:

0+x + x^2 + 2*x^3 + 3*x^4 + 5*x^5 + 8*x^6 + 13*x^7 + 21*x^8 + 34*x^9 + 55*x^10 + 89*x^11 + O(x^12)

Como vemos, los coeficientes son los números de Fibonacci. Si quisiéramos hacer F0=1 nos daría otro resultado, pues la función generatriz seria 1/(1-x-x2) (intenta demostrarlo)

Programaríamos en PARI esta variante:

{write("final.txt",taylor(1/(1-x-x^2), x,12))}

Obtendríamos la sucesión comenzando en 1:

1 + x + 2*x^2 + 3*x^3 + 5*x^4 + 8*x^5 + 13*x^6 + 21*x^7 + 34*x^8 + 55*x^9 + 89*x^10 + 144*x^11 + O(x^12)

Lo podemos intentar con el ejemplo de la entrada anterior, el de las bolas de colores, cuya función generatriz sin desarrollar era

Deberíamos escribir en PARI

{write("final.txt",taylor(x^5*(x^4-1)^2*(x^5-1)/(x-1)^3, x,20))}

Obtendríamos

x^5 + 3*x^6 + 6*x^7 + 10*x^8 + 13*x^9 + 14*x^10 + 13*x^11 + 10*x^12 + 6*x^13 + 3*x^14 + x^15 + O(x^20)

Identificamos los resultados de la entrada anterior: Con grado 9 el coeficiente es el esperado, 13. Para el grado 14 sólo 3 y para el grado 7 los 6 que ya se obtuvieron.

En otra entrada posterior veremos las funciones generatrices de combinaciones, variaciones y demás. Hoy seguiremos con ejemplos:

¿De cuántas formas se puede descomponer el número 28 como suma de primos distintos?

Los primos inferiores a 28 son 2, 3, 5, 7, 11, 13, 17, 19, 23. Cada uno de ellos o está en la suma una vez o no está. Entonces podemos usar términos del tipo (1+x^7) o (1+x^11) para combinarlos en la función generatriz:

F(x)= (1+x²) (1+x³) (1+x⁵) (1+x⁷) (1+x¹¹) (1+x¹³) (1+x¹⁷) (1+x¹⁹) (1+x²³)

Desarrollamos con wxMmaxima y nos da (hemos recortado la imagen):

Vemos que para el exponente 28 el coeficiente es 6, luego existen 6 formas distintas de expresar 28 como suma de primos distintos

Lo comprobamos con nuestra herramienta PARTLISTA (http://hojamat.es/sindecimales/aritmetica/herramientas/herrarit.htm#reprenum)

Con ella vemos las seis sumas: 28=11+17=5+23=3+5+7+13=2+7+19=2+3+23=2+3+5+7+11

Con la función generatriz hemos conseguido también otros muchos coeficientes, pero cuidado con la lista de primos 2, 3, 5, 7, 11, 13, 17, 19, 23. Este desarrollo sólo valdría hasta el 28. Para números mayores deberíamos añadir (1+x^29)(1+x^31)…

Luego con la función generatriz podemos resolver varios problemas a la vez.

Estos desarrollos algebraicos son pesados incluso con ayuda de los CAS. Sería bueno elegir un solo coeficiente en ellos. El lenguaje PARI que estamos usando últimamente (es gratuito aunque de gestión poco amigable) posee la función POLCOEFF(P(x),E) en la que P(x) es un polinomio y E un exponente dado, y nos devuelve el coeficiente de ese polinomio correspondiente al grado E. Nuestro problema del 28 se calcularía así:

print(polcoeff((x^2+1)*(x^3+1)*(x^5+1)*(x^7+1)*(x^11+1)*(x^13+1)*(x^17+1)*(x^19+1)*(x^23+1),28))

Daría un resultado de 6. Si cambiamos 28 por un número inferior obtendremos más resultados. Para números superiores deberíamos incrementar los primos.

Funciones generatrices en Combinatoria

Nota: Esta entrada es la primera de una serie que dedicaremos a las funciones generatrices y a los productos cartesianos condicionados. Aunque estarán enlazadas dentro de una secuencia lógica, para no convertir este blog en un tratado primaveral sobre Combinatoria, las iremos alternando con varias referentes a otros temas, en los que no abandonaremos nuestros queridos números ni a las hojas de cálculo.


Antes de iniciar cualquier planteamiento teórico sobre las funciones generadoras (o generatrices), muy usadas en Combinatoria y en el estudio de las sucesiones de números, las introduciremos mediante un ejemplo. Después, en sucesivas entradas, estudiaremos el concepto con más profundidad. Al principio no se ve bien la utilidad de estas funciones, pero si lees la serie entera que vamos a ir publicando descubrirás que constituyen un buen instrumento de cálculo. Paciencia, pues.

Problema

Deseamos elegir nueve cuentas de colores. Disponemos de cantidad suficiente de cuentas rojas, amarillas y verdes, pero queremos elegir entre 2 y 5 rojas, menos de 4 amarillas y al menos 3 verdes. ¿De cuántas formas podemos efectuar la elección?

Este problema nos plantea el desarrollo de particiones condicionadas de un número.

Ya hemos tocado este tema en el blog (http://hojaynumeros.blogspot.com.es/2011/02/particiones-de-un-numero.html y en
http://hojaynumeros.blogspot.com.es/2011/02/funciones-de-particion-de-un-numero.html)

Independientemente de que se trate de particiones, intentaremos resolver el problema con varias técnicas, y entre ellas la del uso de una función generatriz.

Con un producto cartesiano

Como de las rojas hemos de elegir entre 2 y 5, de las amarillas de 0 a 3 y de las verdes un mínimo de 3 y un máximo de 7 (¿por qué 7?), bastará formar el producto cartesiano  {2,3,4,5}{0,1,2,3}{3,4,5,6,7} e ir eligiendo las ternas que sumen 9. Un problema totalmente elemental que se puede resolver en enseñanzas medias.

A poco que nos pongamos obtendremos: 9= 2+0+7 = 2+1+6 = 2+2+5 = 2+3+4 = 3+0+6 = 3+1+5 = 3+2+4 = 3+3+3 =4+0+5 = 4+1+4 = 4+2+3 = 5+0+4 = 5+1+3.

En total 13 particiones. Para este problema no se necesitarían más técnicas, pero lo estamos tomando como modelo sencillo de introducción.

Con la hoja de cálculo “Cartesius”, aún en fase beta y que presentaremos más adelante, se pueden construir productos cartesianos condicionados. En el problema que nos ocupa plantearíamos esto, que se comprende sin más explicación:

Y obtendríamos

Como era de esperar, son los mismos resultados. Veamos ahora el método que deseamos explicar hoy.

Función generatriz

La idea revolucionaria de la función generatriz consiste en sustituir los distintos elementos numéricos por potencias de una indeterminada, y los conjuntos convertirlos en polinomios. Así el producto cartesiano {2,3,4,5}{0,1,2,3}{3,4,5,6,7} se puede escribir en forma de producto de polinomios:

Es fácil  darse cuenta de que si multiplicamos algebraicamente estos polinomios, el término de grado 9 tendrá como coeficiente el número de particiones pedido, en este caso 13.

Hemos sustituido una técnica de conteo por otra de tipo algebraico o analítico (esto último lo veremos más adelante)

En este caso tan sencillo parece que esto es una complicación, pero en casos generales veremos que puede resultar muy útil. Por dos motivos:

•  Las técnicas algebraicas y analíticas permiten simplificar estos productos y  encontrar directamente el coeficiente deseado.
•  Al desarrollar estos polinomios no sólo resolvemos el problema para 9 cuentas, sino para cualquier otro total posible.

Disponemos en este caso de una ayuda, y es que sabemos sumar muy bien las progresiones geométricas. Así, el producto de polinomios dado se puede expresar como

Este a su vez se puede simplificar como

Con tres pasadas del algoritmo de Ruffini encontraremos todos los coeficientes (omitimos los que no influyen en el problema)

Hemos destacado el que nos interesa: con grado 9 aparece el coeficiente 13, que es la solución, pero este procedimiento nos devuelve mucho más. Por ejemplo, con suma 7 sólo son posibles estas elecciones: 7=2+0+5=2+1+4=2+2+3=3+0+4=3+1+3=4+0+3, seis en total, como marca el esquema, y con suma 14 sólo deberán aparecer 3. Son estas: 4+3+7 = 5+3+6 = 5+2+7

Hemos resuelto varios problemas en uno

Si dispones de un CAS puedes ahorrarte bastante trabajo. Aquí tienes el resultado con la calculadora Wiris (hemos recortado la imagen) Compara los coeficientes con los que resultan con Ruffini.

Es normal que pienses que es mucha complicación, pero se trataba de un problema elemental. En sucesivas entradas daremos un enfoque teórico a las funciones generatrices y nos complicaremos un poco, descubriendo así su utilidad.

Carnaval de triangulares

Esta entrada se ha desbordado, como una serpentina que al arrojarla ya no puede volver a ser rollo. Comenzamos estudiando variantes de la entrada anterior

(http://hojaynumeros.blogspot.com.es/2013/02/de-los-triangulares-alojados-los-primos.html)

y a la multiplicidad de divisores triangulares le siguió su suma, las coincidencias en esa suma y la reconstrucción de otro triangular. Por ello, aunque sea infrecuente en un blog, comenzamos con un esquema:

* Número de divisores triangulares (http://oeis.org/A007862)

• Cálculo general
• Números con un solo divisor triangular propio (http://oeis.org/A203468)

* Suma de divisores triangulares (http://oeis.org/A185027)

* Curiosidades

• Números con suma de triangulares también triangular (http://oeis.org/A209309)
• Números triangulares con la misma propiedad (http://oeis.org/A209310)
• Otras curiosidades menores
• Caso en el que la suma de divisores triangulares es otro divisor (http://oeis.org/A209311)

Número de divisores triangulares

Vimos en la entrada anterior que el número 40 posee una parte triangular igual a 10, que le permite ser representado como un prisma triangular.

Esta representación 40=4*T₄ es única (en toda la entrada no consideramos el triangular 1, por lo que no volveremos a citarlo). Ningún otro número triangular menor o igual que 40 (3, 6, 10, 15, 21, 28 o  36) lo divide salvo el 10. Esto por lo que se refiere al 40, pero existen otros números que admiten varias representaciones. El 30 admite cuatro: 30=10*T₂ = 5*T₃ = 3*T₄ = 2*T₅

No es difícil contar los divisores triangulares que posee un número N (al menos, el 1). Basta cambiar el algoritmo que publicamos en la anterior entrada para que cuente en lugar de quedarse con el mayor

Public Function numdivtriang(n)
Dim p, i, t, tr

p = Int((Sqr(n * 8 + 1) - 1) / 2) ‘Calcula el máximo orden
t = 1
For i = 2 To p
tr = i * (i + 1) / 2 ‘forma todos los triangulares menores o iguales a n
If n / tr = n \ tr Then t = t+1  ‘si es divisor, incrementa el contador
Next i
numdivtriang = t ‘se queda con el mayor
End Function

Esta función cuenta el 1, por lo que para 30 dará 5 posibilidades y para 40 sólo 2. En la siguiente tabla parcial lograda con hoja de cálculo lo puedes comprobar

30    5
31    1
32    1
33    2
34    1
35    1
36    4
37    1
38    1
39    2
40    2

Este resultado lo tienes en http://oeis.org/A007862 y es interesante leer los comentarios que se incluyen.

Números con un solo divisor triangular propio mayor que 1

El caso del 40 no es único. Hay muchos números que sólo pueden representarse de una sola forma como un prisma triangular con base y altura mayores que uno (para evitar trivialidades). Son estos:

6, 9, 15, 20, 21, 27, 33, 39, 40, 50, 51, 56, 57, 69, 70, 80, 81, 87, 93, 99, 100, 111, 112, 117, 123, 129, 130, 141, 153, 159, 160, 170, 171, 177, 182, 183, 190, 196, 200, 201, 207, 213, 219, 224, 230, 237, 243…

Los hemos publicado en http://oeis.org/A203468

Prueba con cualquiera de ellos, el 182=2*7*13. Puedes usar la propiedad que vimos de que su doble ha de tener dos divisores consecutivos. 364=2*2*7*13 y su conjunto de divisores es {364, 182, 91, 52, 28, 26, 14, 13, 7, 4, 2, 1}. Los únicos divisores consecutivos son 13 y 14, que dan lugar a un único divisor triangular de 182, el 91.

Por cierto, su consecutivo 183 presenta la misma situación: su único divisor triangular es el 3. No es el único par de consecutivos contenido en la sucesión. Por ejemplo, tenemos 170 y 171.

Dentro de esta sucesión figuran números triangulares. Todos ellos presentarán tres divisores triangulares: ellos, un divisor propio y la unidad. Así, 351 tiene como únicos divisores triangulares 1, 3 y el propio 351.

Suma de divisores triangulares

Además de considerar la suma de todos los divisores de un número, puede resultar curioso sumar sólo los de un tipo. Por ejemplo, el número 720 tiene como suma de divisores 2418, pero si sólo consideramos los que son cuadrados, sumarían 210=144+36+ 16+9+4+1 y con los triangulares 236=120+45+36+15+10+6+3+1. Se pueden considerar otros tipos de divisores: los pares, los oblongos…

Un algoritmo un poco burdo, pero que funciona, es el de recorrer todos los posibles divisores y someter a cada uno a una condición antes de incorporarlo a la suma. Aquí tienes el que hemos usado para cuadrados y triangulares:

Public Function sumadiv(nume, tipo)
'tipos
'0 da todos los divisores
'1 los cuadrados
'2 los triangulares

Dim i, s

s = 0
For i = 1 To nume
If esmultiplo(nume, i) Then
If tipo = 0 Then s = s + i
If tipo = 1 And escuad(i) Then s = s + i
If tipo = 2 And estriangular(i) Then s = s + i
End If
Next i
sumadiv = s
End Function

Con un algoritmo similar hemos publicado en OEIS la función que recoge la suma de los divisores triangulares de los primeros números naturales:

1, 1, 4, 1, 1, 10, 1, 1, 4, 11, 1, 10, 1, 1, 19, 1, 1, 10, 1, 11, 25, 1, 1, 10, 1, 1, 4, 29, 1, 35, 1, 1, 4, 1, 1, 46, 1, 1, 4, 11, 1, 31, 1, 1, 64, 1, 1, 10, 1, 11, 4, 1, 1, 10, 56, 29, 4, 1, 1, 35, 1, 1, 25, 1, 1, 76…
(http://oeis.org/A185027)

Curiosidades

Esta sucesión da lugar a varias curiosidades:

La suma de triangulares puede ser triangular.

Excluimos el caso en que sea igual a 1 por trivial. Estos son los números que lo cumplen:

6, 12, 18, 24, 48, 54, 96, 102, 110, 114, 138, 162, 174, 186, 192, 204, 220, 222, 228, 246, 258, 282, 315, 318, 348, 354, 364, 366, 372, 384, 402, 414, 426, 438, 440, 444, 456, 474, 486, 492, 498, 516, 522, 534, 550…
También la acabamos de publicar (http://oeis.org/A209309)

Por ejemplo, 444 tiene como divisores triangulares 6, 3 y 1, y su suma es 10 que es triangular. Más complejo sería el caso de 1320, cuyos divisores triangulares, 120, 66, 55, 15, 10, 6, 3 y 1 suman 276, que es triangular igual a 23*24/2.

Similares a esta, pero menos exigentes, son estas condiciones:

(1) La suma de los divisores ordinarios es triangular
1, 2, 5, 8, 12, 22, 36, 45, 54, 56, 87, 95, 98, 104, 116, 152, 160, 200,… (A045746)

(2) La que es triangular es la suma de las partes alícuotas, y mayor que 1
2,4,6,14,16,18,24,25,28,33,36,51,54,66,91,112,...(no publicada)

(3) Números triangulares en los que la suma de sus divisores propios es también triangular
1, 3, 6, 28, 36, 66, 91, 231, 496, 8128, 14196, 15225, 129795, 491536,… (A083675)

 (4) Números triangulares cuya suma de divisores es también triangular
1, 36, 45, 23220, 105111, 135460, 2492028, 5286126, 6604795, 14308575, 45025305, 50516326, 54742416, 99017628, 108125865, 152486916,... (A083674)

Ahora viene la nuestra, la más exigente:

Números triangulares cuya suma de divisores triangulares es mayor que 1 y triangular

Esto es ya mucho exigir, por lo que las soluciones crecen rápidamente:

6, 4186, 32131, 52975, 78210, 111628, 237016, 247456, 584821, 750925, 1464616, 3649051, 5791906, 11297881, 16082956, 24650731, 27243271, 38618866, 46585378, 51546781, 56026405, 76923406, 89880528, 96070591…(http://oeis.org/A209310)

El estudio del código PARI de esta sucesión te enseñará técnicas útiles:

istriangular(n)=issquare(8*n+1)
{t=0; for(n=1, 10^8, if(istriangular(n), k=sumdiv(n, d, istriangular(d)*d) ; if(istriangular(k)&&k>>1, t+=1; write("b209310.txt", t, " ", n))))}

Y por último, para no cansar (si es que has llegado hasta aquí), la última curiosidad

Números en los que la suma de divisores triangulares es mayor que 1 y divisor del número

285, 1302, 1425, 1820, 2508, 3640, 3720, 4845, 4956, 5016, 5415, 7125, 7280, 9100, 9114, 9912, 11685, 12255, 12740, 14508, 15105, 16815, 17385, 18200, 19095, 19824, 20235, 20805, 22134, 22515, 23655, 23660, 24021, 24738…http://oeis.org/A209311

Aquí tienes dos ejemplos:

285.-Divisores triangulares:1, 3 y 15 y su suma, 19,  es divisor de 285

1302.- Divisores triangulares: 21 + 6 + 3 + 1 = 31 que es divisor de 1302.

Código PARI

 istriangular(n)=issquare(8*n+1)
{t=0; for(n=1, 10^7, k=sumdiv(n, d, istriangular(d)*d); if(n/k==n\k&&k>>1, t+=1; write("b209311.txt", t, " ", n)))}

Nos queda algo en el tintero, porque en esta última el cociente puede ser también triangular, pero esto queda para otro día.

De los triangulares alojados a los primos de Sophie Germain

Esta entrada participa en la Edición 4.1 del Carnaval de Matemáticas cuyo anfitrión es Tito Eliatron Dixit.

Si tomamos 40 cubos, los podemos apilar en forma de prisma con base un triángulo isósceles y rectángulo, o en términos aritméticos, un número triangular mayor que 1. Excluimos la unidad porque en ese caso se pierde la forma triangular.

Este ejemplo es válido porque 40=4*10, y 10 es el cuarto número triangular.

No todos los números enteros se pueden representar así, pues han de ser múltiplos de un número triangular y eso no siempre ocurre. Por ejemplo, el 14, ya que entre sus divisores no figuran 3, 6 ó 10, que son los triangulares menores que él (recuerda que excluimos el 1). Esto ya nos divide el conjunto  de los números naturales entre los que tienen divisores triangulares mayores que 1 y los que no.

Los segundos, que no admiten la representación propuesta, son 1, 2, 4, 5, 7, 8, 11, 13, 14, 16, 17, 19, 22, 23, 25, 26, 29, 31, 32, 34, 35, 37, 38… (http://oeis.org/A112886) y les llamaremos libres de triangulares. Verás que están los primos, algunos semiprimos, potencias de primos y otros a los que volveremos más adelante.

Parte triangular y parte libre de triángulos

Sabemos que los primeros admiten un divisor triangular pero, como pueden ser varios, nos quedaremos con el mayor: llamaremos parte triangular (PTR) de un número al mayor divisor triangular que posea. Si has leído sobre estos temas, te recordará esto a la parte cuadrada y la parte libre de un número (http://hojaynumeros.blogspot.com.es/2011/05/parte-cuadrada-y-parte-libre.html)

El mayor divisor triangular puede ser 1 o el mismo número, como se comprueba en la lista de todos ellos (http://oeis.org/A115017):

N 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11,  12,  13, 14, 15, 16,  17,  18,  19,  20, 21, 22…
PTR 1, 1, 3, 1, 1, 6, 1, 1, 3, 10,   1,   6,    1,   1,  15,   1,   1,   6,    1,  10, 21,   1…

En ella están los libres de triangulares, que son los que se corresponden con un 1, como el 4 y el 5, los triangulares, cuya PTR son ellos mismos, como 6 y 10, y el resto, en el que se tiene una parte triangular y otra libre ambas mayores que la unidad. Es el caso de 12 o 40. La parte libre de estos últimos está recogida en http://oeis.org/A121289

Una idea: dos números con la misma parte libre y partes triangulares consecutivas formarán un prisma cuadrado. Imagina el prisma de la primera imagen y su complementario.

Búsqueda de la parte triangular

Un algoritmo simple es el de ir recorriendo los números naturales k, formar con ellos los triangulares mediante k(k+1)/2 e ir verificando si el número dado N es múltiplo de alguno. El mayor de todos ellos será la PTR(N).

Previamente es bueno calcular el orden del máximo triangular que es menor o igual que N, para acortar el ciclo de búsqueda. Se deja a los lectores la demostración de que ese orden k se calcula mediante

En hoja de cálculo sería =ENTERO((RAIZ(8*N+1)-1)/2)

Por ejemplo, para N=14534, k=169 y el mayor triangular menor que N, 169*170/2 = 14365. A nosotros nos interesaría el 169, porque entre 2 y 169 estaría el orden del triangular buscado. Todo esto se puede plasmar en una función:

Public Function partetriang(n)
Dim p, i, t, tr

p = Int((Sqr(n * 8 + 1) - 1) / 2) ‘Calcula el máximo orden
t = 1
For i = 2 To p
tr = i * (i + 1) / 2 ‘forma todos los triangulares menores o iguales a n
If n / tr = n \ tr Then t = tr ‘si es divisor, toma nota
Next i
partetriang = t ‘se queda con el mayor
End Function

El algoritmo busca los triangulares entre el menor 3 y el mayor k(k+1)/2 y se va quedando con los divisores. El último encontrado será PTR(A).

Si en lugar de recoger el valor de i*(i+1)/2 hubiéramos ido recogiendo i, nos hubiera resultado el orden de PTR. Los tienes en http://oeis.org/A083312

¿Qué números dan alojamiento a un triangular?

Para que N tenga un divisor triangular mayor que 1 se ha de poder escribir de la forma N=k(k+1)*M/2 con k>1. Esto da lugar a varias interpretaciones:

(a) N tiene un divisor triangular mayor que 1, si y sólo si 2N posee dos divisores consecutivos mayores que 1.

Es condición necesaria, pues la expresión de 2N sería 2N=k(k+1)*M con k>1, con lo que k y k+1 son los divisores pedidos.

Por ejemplo, el triangular 21 divide a N=8883, con lo que el doble 17766=6*7*423 contiene a los consecutivos 6 y 7.

La condición es suficiente: Si 2N posee dos divisores consecutivos h y h+1 con h>1, estos serán primos entre sí, luego su MCM(h,h+1) será su producto h(h+1). Como 2N es múltiplo de h y h+1, lo será de su MCM, es decir de su producto. Por tanto 2N=h(h+1)P y será múltiplo del triangular h(h+1)/2, ya que uno de los dos h o h+1 es par.

Los números 2N de este tipo los tienes en http://oeis.org/A132895. Son el doble de un número libre de triángulos.

Sería interesante que pensaras en un algoritmo que descubriera esos números.

(c) Los semiprimos N=p*q son números libres de triángulos salvo que uno de sus factores  sea 3, o bien q=2p±1 

En efecto, si N=p*q con p y q primos, 2N=2pq ha de contener dos divisores consecutivos. Si p o q fueran iguales a 3, ya se cumpliría, porque 2N=2*3*k, pero entonces N sería múltiplo de 3.

Si ni p ni q son iguales a 3, lo serán a 2 o a un primo mayor que 3. Si por ejemplo p=2 entonces 2N=2*2*q y q se ve obligado a ser 3, con lo que pasamos al primer caso. Seria múltiplo de 3.

Así que sólo nos queda que N=p*q con p y q primos mayores que 3  y no son números consecutivos (porque son impares). En ese caso es claro que 2N=2pq no podría tener divisores consecutivos salvo que q=2p+1 o bien q=2p-1 (o simétricamente, p=2q+1 o 2q-1). En el primer caso p sería un primo de Sophie Germain.

Recuerda que los primos de Sophie Germain son aquellos en los que 2*p+1 también es primo: 2, 3, 5, 11,…

(d) Los números primarios (potencias de primos) están libres de triángulos salvo el caso N=3k

Esta es trivial: Si N=pk con k>1, entonces 2N=2pppp…sólo contendría divisores consecutivos en el caso 2N=2*3*3*3...

¿Se te ocurren más propiedades? A nosotros por ahora no.

Convertir esquemas de cálculo en tablas

Desde este blog y nuestra página hojamat.es hemos promovido siempre el uso de esquemas de cálculo y apuntes interactivos (http://hojamat.es/contenidos/apuntes.htm)

La limitación que presentan es que al provenir de cálculos de cierta complejidad es difícil recogerlos en forma de función o tabla. En esta entrada presentaremos una herramienta sencilla para recoger resultados de esquemas en forma de tabla.

La herramienta

Al principio intentamos recogerlos como funciones, pero ni Excel ni OpenOffice ni LibreOffice lo permiten. Hay algo en la programación de funciones que hace que si se altera el valor de una celda cualquiera dentro del proceso de cálculo de la función, esta no recoja el valor que ha de devolver. Continuamente da mensajes de error.

Lo que sí podemos es construir una tabla que altere los parámetros del esquema y recoja el valor final del cálculo. Como estamos abusando de generalidades, lo explicaremos mejor con un ejemplo.
Partiremos del cálculo del fósil de un número (http://hojaynumeros.blogspot.com.es/2008/10/dndole-vueltas-2.html)

Este valor se halla multiplicando las cifras de un número, volviendo a realizar esta operación en el resultado y en los siguientes hasta llegar a un número de una cifra al que llamamos fósil.

Por ejemplo, partimos de 876, multiplicamos 8*7*6=336. Volvemos a multiplicar y reiteramos. 3*3*6=54, 5*4=20, 2*0=0, luego el fósil de 876 es 0.

Este cálculo lo podemos tener implementado en una hoja mediante un esquema (en este momento no nos va a interesar qué fórmulas se han usado)

Lo que deseamos es poder extraer de este esquema una tabla de valores y un gráfico si vamos cambiando el 876 por ejemplo en el intervalo 860…880. Esta tarea la realiza la hoja de cálculo esquefun, que está alojada en http://hojamat.es/sindecimales/aritmetica/herramientas/herrarit.htm#esquefun


Tabla simple

Si la abres verás que la primera hoja estará en blanco o contendrá un esquema de cálculo cualquiera. En el segundo caso puedes borrarlo todo y construir tu propio proceso. No sobrepases el tamaño de una pantalla o algo más (unas treinta filas por treinta columnas).

En la segunda hoja se te ofrece la posibilidad de construir la tabla

¿Cómo se consigue esto? Lo explicaremos paso a paso para una tabla simple X,FX:

(1) En la primera hoja, en la celda que está a la izquierda del valor de la variable independiente escribes una X. Debes procurar tener esa celda siempre libre. También es conveniente que uses las primeras filas y columnas.

 (2) También a la izquierda del resultado que te interese como valor de la función (en nuestro caso el fósil) escribes una F.

Si tu resultado no está en ellas siempre puedes copiarlas dinámicamente. Por ejemplo, si tienes el resultado en la celda AH44, basta que en una celda más arriba y a la izquierda escribas la fórmula =AH44 y así todos los resultados se copiarán a esa celda.

Con eso ya has terminado la preparación de la hoja 1: Escribir “X” a la izquierda de la variable independiente y una “F” al lado de la variable dependiente.

(3) En la segunda hoja define el mínimo, máximo y salto de la tabla que deseas para concretar los valores de X en la misma (puedes hacerlo de forma manual, pero sería cosa tuya borrar en la columna de la X todos los datos sobrantes)

(4) Pulsa el botón FX y obtendrás tabla y gráfico de los resultados. En el volcado de pantalla de más arriba puedes observar que la tabla va de 860 a 880 y que el comportamiento del fósil es totalmente irregular.

Ya está. No hay que trabajar más. Después tabla y gráfico los puedes exportar a cualquier otro documento.

Tabla simple con dos funciones

Lo explicamos también con un ejemplo:

Deseamos hacer entender a nuestro alumnado que la raíz de una suma no da el mismo resultado que la suma de raíces. Para ello hemos pensado en usar

Preparamos un esquema de cálculo en el que se manifiesten las diferencias. Podía ser este:

Si ahora deseamos ver en un gráfico cómo evolucionan las diferencias necesitaremos definir dos funciones. Así que rotulamos con X el número, con F el primer cálculo y con G el segundo (estos nombres son obligatorios)

A partir de este esquema ya rotulado podemos crear tabla y gráfico

Hemos construido una tabla del 10 al 2000 con saltos de 10, con la ¿sorpresa? de que las diferencias tienden a estabilizarse a valores cercanos al 2

En la imagen aparece una tercera columna de diferencias que se han creado manualmente. El objetivo de construir la tabla a partir del esquema se ha conseguido.

En cursos algo más avanzados puedes intentar demostrar que efectivamente el límite de la diferencia entre ambos resultados es 2.

Tabla doble

Sería también útil estudiar una función que dependiera de dos variables. Para eso dispones de la tercera hoja de esta herramienta. No se ha incluido el gráfico para no tener que insertar otra cuarta hoja, pero nuestros lectores sabrán cómo construirlo.

En este caso deberemos rotular con X e Y las dos variable independientes y con F la función. Lo explicaremos con un ejemplo que no tiene más interés que la mera curiosidad:

Sabemos que los pasos necesarios en el algoritmo de Euclides para obtener el MCD de dos números varía mucho según los datos usados. Intentemos formar una tabla de doble entrada con ellos.
Imagina que hemos trasladado a la primera hoja el algoritmo de nuestra herramienta Euclides (http://hojamat.es/sindecimales/congruencias/herramientas/herrcong.htm):

Debemos ahora, después de comprobar que funciona bien, borrar  los rótulos “Primer número” y “Segundo número” y sustituirlos por X e Y respectivamente. Abajo también sustituiremos “Número de cocientes” por F, para recoger su valor como una función. En este ejemplo tenemos un problema, y es que esas celdas están combinadas. Debes primero anular la combinación y después escribir X,Y,F de forma contigua a su valor.

Pasamos a la tercera hoja y definimos intervalos y saltos para X e Y, por ejemplo, de 20 a 30 con saltos de 1 (el carácter optativo se incluye porque se puede efectuar un relleno manual, aunque no es muy aconsejable).

Pulsamos el botón Fxy y se formará la tabla de doble entrada:

Llama la atención que no es simétrica para intercambios entre X e Y, pero es que si el primer número es menor que el segundo nos cuesta un paso más en el algoritmo.

Con los procedimientos habituales podemos traducirla a un gráfico 3D:

Estas son las tres modalidades de creación de tablas que hemos incluido en esquefun. Con ellas basta para encontrar usos en la enseñanza y como herramienta de búsqueda. Que os sea útil.

Pandigitales, cromos y un poco de Benford (2)

Esta es la segunda parte de nuestra participación en el  Carnaval de Matemáticas, Edición 3.1415926535, cuyo anfitrión es La Aventura de la Ciencia.

Estudiamos en la anterior entrada cuándo unos resultados de potencias se convierten todos en pandigitales en sentido amplio (con repetición) Llegamos a la conclusión de que esto ocurre aproximadamente cuando la potencia alcanza unas 50 cifras. A partir de este resultado sospechamos que la distribución de cifras no presentan uniformidad.

Estadísticas de cifras

Según lo anterior, hemos de tener cuidado en considerar como casi uniforme la aparición de cifras en las potencias de un número. Si así aparecieran, deberíamos encontrar más similitud entre lo que se espera de un fenómeno aleatorio y este que nos ocupa.

La uniformidad

Para estudiar de forma empírica la distribución de cifras en las potencias hemos elegido las bases entre 2 y 9, a cada una la hemos elevado a todos los exponentes comprendidos entre 1 y 50, pues son los cálculos antecedentes de la región en la que desaparecen los resultados que no son pandigitales. Para ello nos ha sido útil nuestra calculadora STCALCU para hoja de cálculo (que presentaremos próximamente). Hemos obtenido este resultado:

 (1) Las desviaciones típicas mayores se corresponden con los divisores del 10, 2 y 5. Después baja algo en los números no coprimos con 10: 4, 6 y 8. Por último, son más homogéneas en los coprimos, 3, 7 y 9. Esto tiene cierto sentido, pero no seguiremos por ahí.

(2) La distribución por cifras presentan un máximo en la cifra 1 (como en la Ley de Benford) de 11,2%, muy alejado del 8,7% de la cifra 0. Además, las cifras impares aparecen más que las pares. Esto lo afirmamos descriptivamente, pues una prueba chi-cuadrado no da significación.

(3) Casi todas las cifras presentan un máximo en la base igual a ellas. Las hemos destacado en rojo. Llama la atención el 14% del 5 y del 6. A eso no es ajeno el que en esas cifras coincidan las terminaciones de sus potencias sucesivas: 5, 25, 125, 625, … y 6, 36, 216,…Te puedes divertir intentando analizar otros casos.

Resumiendo, no aparece una uniformidad clara en los resultados, que más bien parecen sesgados hacia el 1 y los impares. ¿Se mantendrá esta tendencia para potencias mayores?

Hemos acumulado los resultados desde exponente 1 al 200, para ver cómo evoluciona la distribución de cifras, llegando a esto:

Aquí el panorama cambia algo: se percibe más uniformidad, aunque el 1 es la cifra que presenta mayor frecuencia. Por tanto debemos pensar que en las primeras potencias las cifras aparecen con frecuencias más alejadas del 10% y que eso es lo que produce que se tenga que llegar a unas 50 cifras para llegar a completar el carácter pandigital

La herencia

Otra pregunta sería pertinente: estas desviaciones de la uniformidad ¿se mantienen de cierta forma entre unas potencias y las posteriores dentro de una misma base? Si una cifra presenta una frecuencia en 2^N sería interesante saber cómo se comporta en 2^N+1. Pues bien, aquí tampoco se ve relación clara y significativa entre las frecuencias de un exponente  con el siguiente. Tomamos como ejemplo la base 5 haciendo trampa, porque podía esperase que la cifra 5 y la 0 se mantuvieran en sus frecuencias al crecer N.

Basta ver los máximos y mínimos para darnos cuenta de lo alejada de la uniformidad que está la distribución de cifras. Respecto a la herencia, si recorres los porcentajes correspondientes a cada cifra sí se percibe una cierta constancia en la tendencia. No es importante. Le hemos aplicado la prueba Chi-cuadrado y no nos da una diferencia significativa respecto a la homogeneidad máxima.

Así que, por si acaso, no uses potencias para extraer números psudoaleatorios, que te puedes llevar sorpresas.

Nuestra Ley de Benford

Y ya puestos, ¿cómo se comportan las primeras cifras de cada potencia? Recuerda que según la Ley de Benford (en la Red tienes muchas referencias a ella, por ejemplo en http://www.estadisticaparatodos.es/taller/benford/benford.html),

se podría esperar un  30% para el 1, un 17% para el 2, 12% para el 3 y así disminuyendo para el resto, como se ve en la gráfica incluida en la página recomendada.

Lo intentamos: elevaremos las distintas bases de 2 a 9 (podían ser otras) a todos los exponentes comprendidos entre 1 y 250 y recogeremos las estadísticas de la primera cifra.
Son estas:

Esto quiere decir que respecto a la Ley de Benford las potencias se comportan admirablemente. Hemos comparado nuestras frecuencias con la fórmula de Benford LOG((d+1)/d) y nos ha resultado:

No necesita comentario. El comportamiento de las estadísticas globales viene dado más por las cifras intermedias que por la primera, que sigue la distribución esperada. A partir de aquí puedes emprender un estudio del que sólo hemos esbozado el principio.