Damos vueltas a los triangulares cuadrados (1) Generación.

Generación de la sucesión

Dos entradas del blog de John D. Cook (http://www.johndcook.com/blog/2015/08/20/when-is-a-triangle-a-square/ y siguiente) me han animado a volver a tomar el tipo de entrada al que llamé “dar vueltas” a un tema o concepto. Lo haré sobre estos números:

0, 1, 36, 1225, 41616, 1413721, 48024900, 1631432881,… (http://oeis.org/A001110)

Evidentemente,  tienen en común el ser triangulares y cuadrados a la vez. Puedes leer un desarrollo sencillo y claro en este documento:

http://www2.caminos.upm.es/Departamentos/matematicas/revistapm/revista_impresa/vol_IV_num_1/jue_mat_num_triang.pdf

Me he dado cuenta de que es un concepto sencillo pero que da lugar a bastantes reflexiones, algoritmos y repasos de teoría.

Nosotros seguiremos en parte este documento para iniciar el tema.

Búsqueda de números triangulares cuadrados

Un número triangular tiene por fórmula n(n+1)/2 y un cuadrado m². Aquellos números que participen de las dos características tendrán que cumplir la igualdad

De esta igualdad deducimos esta otra mucho más práctica:

O bien

Con cambio de variable se convierte en una ecuación de Pell:

Esta ecuación la tenemos muy estudiada

http://hojamat.es/parra/pell.pdf (documento de Rafael Parra)

http://hojaynumeros.blogspot.com.es/2010/02/ecuacion-de-pell.html

Disponemos además de una hoja de cálculo para ayudar a resolverla:

http://www.hojamat.es/sindecimales/aritmetica/herramientas/herrarit.htm#pell

Usamos esta herramienta para el coeficiente 8 y el segundo miembro 1 y nos dan las primeras soluciones:

(1,0) (3,1) (17,6) (99,35) (577,204) (3363,1189) (19601,6930) (114243, 40391),…

Por recurrencia:

Según la teoría de la ecuación de Pell, las soluciones aparecen con las recurrencias (en este caso) x_n=3*x_n-1+8*y_n-1  y_n=3*y_n-1+1*x_n-1. Por ejemplo, 99=3*17+8*6, 35=3*6+1*17.

Ahora sólo nos queda elevar al cuadrado las soluciones de y (que equivalen a la variable m de la primera igualdad que planteamos) y nos resultarán los triangulares cuadrados:

0²0, 1²=1, 6²=36, 35²=1225, 204²=41616, 1189²=1413721,…

Primer algoritmo

El estudio que acabamos de desarrollar nos da una pista para la generación de términos triangulares cuadrados: Iniciamos dos variables X=1, Y=0, y en  cada paso del algoritmo convertimos X en 3X+8Y y la Y en 3Y+X. Terminado el cálculo presentamos el valor de Y² como siguiente triangular cuadrado. En el Basic de las hojas de cálculo quedaría así:

Sub triangcuad()
Dim x, y, x1, y1, i, t, fila

x = 1: y = 0 ‘Valores de inicio
fila = 3 ‘Fila inicial
Cells(fila, 4).Value = 0 ‘El primer valor es un cero
For i = 1 To 8 ‘Calculamos sólo ocho
x1 = 3 * x + 8 * y ‘Iteración para x
y1 = 3 * y + x ‘Iteración para y
x = x1
y = y1
t = y * y ‘Número triangular cuadrado
fila = fila + 1
Cells(fila, 4).Value = t ‘Se presenta el resultado
Next i
End Sub

Obtendríamos:

El que el último se nos ofrezca en coma flotante nos da idea de las limitaciones de la hoja para cálculos con enteros de muchas cifras. Si acudimos a PARI no nos encontraremos con esas limitaciones. Prueba este código:

{x=1;y=0;print(0);while(x<10^10,x1=3*x+8*y;y1=3*y+x;x=x1;y=y1;t=y^2;print(t))}

En pocos segundos te presenta los triangulares cuadrados menores que 10^10.

Relación de recurrencia con una sola variable

Por la naturaleza de su definición podemos esperar que estos números sigan una relación de recurrencia de segundo orden. Para encontrar su expresión, que será del tipo

A_n=aA_n-1+bA_n-2+c usaremos los valores iniciales 0, 1, 36, 1225, 41616 para plantear:

36=a*1+b*0+c
1225=a*36+b*1+c
41616=a*1225+b*36+c

Resolvemos

1189=35a+b
40391=1189a+35b

1224=36a y a=34, b=-1 y c=2

Según los cálculos anteriores, la relación de recurrencia será

A_n=34A_n-1-A_n-2+2

Es la misma que propone John D. Cook en su blog.

En este blog no olvidamos la hoja de cálculo. Intenta una resolución como la de la imagen usando cálculo matricial:

A partir de esta escritura matricial del sistema de ecuaciones, creamos debajo la matriz inversa de los coeficientes con MINVERSA, y a su derecha su producto por los términos independientes con MMULT:

Conseguimos así la misma solución 34, -1, 2

Segundo algoritmo

La relación de recurrencia nos permite un segundo algoritmo para encontrar los triangulares cuadrados. El que describimos a continuación presenta los nueve primeros (después existen problemas de coma flotante)

Sub triancuad1()
Dim m, n, p, k, fila
m = 0: n = 1 ‘Valores iniciales
fila = 3
Cells(1, 3).Value = m ‘presenta los dos primeros términos
Cells(2, 3).Value = n
For k = 1 To 7
p = 34 * n - m + 2 ‘relación de recurrencia
Cells(fila, 3).Value = p: fila = fila + 1 ‘presenta los siguientes términos
m = n: n = p ‘cada término se convierte en el anterior
Next k
End Sub

Los términos rellenarán una columna de hoja de cálculo:

Siguiendo nuestra costumbre, lo traducimos a PARI para conseguir más términos:

{x=0;y=1;print(0);print(1);for (k=1, 20, z=34*y-x+2;print(z);x=y;y=z)}

No resistimos la tentación, al igual que propone J.C. Cook, de intentar una versión recursiva en forma de función. Funciona muy bien en hoja de cálculo:

Public Function ftriangcuad(n)
If n < 2 Then
ftriangcuad = n
Else
ftriangcuad = 34 * ftriangcuad(n - 1) - ftriangcuad(n - 2) + 2
End If
End Function

No necesita explicación. La tabla siguiente se forma con gran rapidez de cálculo:

Fórmula directa

Si lees el capítulo sobre sucesiones recurrentes en nuestra publicación Sucesiones (http://www.hojamat.es/publicaciones/sucesiones.pdf) entenderás que a partir de la fórmula de recurrencia es posible encontrar la expresión directa de cada término (fórmula del término general).  Sólo insertamos la captura de pantalla de nuestra hoja de cálculo Recurrencias
(http://www.hojamat.es/sindecimales/aritmetica/herramientas/herrarit.htm#recurre2) en la parte homogénea de la recurrencia:

Con un ligero retoque y la interpretación de los decimales llegamos a la propuesta por John D. Cook:

El uso de la raíz cuadrada de 2 le quita utilidad en nuestro trabajo, por lo que intentaremos prescindir de ella. Para ver la influencia del formato de coma flotante, la implementamos en hoja de cálculo, y el resultado es similar al de los algoritmos anteriores:

Función generatriz

Para quien no lo sepa, diremos que la función generatriz de una sucesión, si se desarrolla como una serie de potencias, poseerá como coeficientes de esas potencias de x los términos de la sucesión.

En el caso de los números triangulares cuadrados la función generatriz es

 (ver http://oeis.org/A001110)

 Con  esta sencilla orden de PARI podemos comprobar su desarrollo.

print(taylor(x*(1+x)/((1-x)*(1-34*x+x^2)),x,20))

Tercer algoritmo

Finalizamos la entrada con la presentación de un algoritmo de los que llamamos “ingenuos”, que no usan la teoría para simplificar los cálculos, pero sí la fuerza bruta de la velocidad de proceso. En este caso obligaremos a los números naturales a ir creciendo hasta alcanzar un cuadrado, y después a la inversa, que los cuadrados avancen hasta alcanzar un triangular. Cuando se llegue a una igualdad se imprime el resultado. A pesar de su simplicidad, no resulta lento. Es éste:

Sub triancuad2()
Dim i, j, m, n, k, fila
m = 3: n = 4: i = 2: j = 3 ‘Se inician las variables
k = 10 ^ 7
fila = 3
Cells(fila, 3).Value = 1: fila = fila + 1
While m < k ‘Busca soluciones menores que k
While m <= n ‘Los triangulares crecen
If m = n Then Cells(fila, 3).Value = m: fila = fila + 1 ‘Hay igualdad
i = i + 1: m = m + i
Wend
While n <= m ‘Los cuadrados crecen
If m = n Then Cells(fila, 3).Value = m: fila = fila + 1 ‘Hay igualdad
j = j + 2: n = n + j
Wend
Wend
End Sub

Dejamos a los lectores el estudio de  por qué funciona este algoritmo para descubrir los triangulares cuadrados. Como los anteriores, llega a los mismos resultados, en este caso  hasta 10^7.