Numeradores de los números armónicos

Un número armónico H_n es un número racional formado mediante la suma de los inversos de los números naturales hasta n, es decir

La sucesión de los números armónicos tiene límite infinito, por ser divergente la serie armónica a la que pertenecen los sumandos. Es fácil ver que los denominadores de los números armónicos son, salvo simplificaciones, los primeros factoriales, pero aquí nos van a interesar los numeradores.

Con hoja de cálculo no es difícil encontrarlos, pues situamos en columna los distintos valores hasta n y en otra columna vamos dividiendo n! entre 1, 2, 3,… para después sumar la columna de cocientes (que serán todos enteros). Aquí tienes el desarrollo para n=6

Se ha dividido el factorial de 6, 720, entre 1,2,…6, sumando los cocientes hasta conseguir 1764, luego el sexto número armónico tendrá la expresión 1764/720.

Nos basamos en este cálculo en la operación de sacar el denominador común.

Si nos atenemos a la definición primera, vemos que H(n+1)=H(n)+1/(n+1). Si  imaginamos mentalmente el resultado de esa suma entenderemos que el numerador de H(n) quedará multiplicado por n+1 con la suma posterior del primitivo denominador, es decir, n!. Esto nos lleva a una fórmula de recurrencia fácil para los numeradores: a(n+1)=a(n)*(n+1)+n!

Si aplicamos la fórmula de recurrencia al 1, obtenemos la sucesión 1, 3, 11, 50, 274, 1764, 13068, 109584,…que coincide con los números de Stirling de primera clase, publicados en http://oeis.org/A000254 y que puedes estudiar allí. Nosotros simplificaremos las fracciones de los números armónicos, con lo que los valores anteriores se verán sustituidos por

 1, 3, 11, 25, 137, 49, 363, 761, 7129, 7381, 83711, 86021, 1145993,...

Estos son los numeradores que consideraremos en principio. Están también publicados (http://oeis.org/A001008)

Según las ideas contenidas en los párrafos anteriores, la primera suma que consideramos y la operación de simplificar, el algoritmo adecuado para generar estos numeradores simplificados puede ser el siguiente, que lo hemos orientado a una función de VBA:

Public Function numearmonico(n)
Dim i, v, f
v = 0
f = 1
For i = 1 To n: f = f * i: Next i ‘Se construye el factorial
For i = 1 To n
v = v + f / i ‘Se suman los cocientes del factorial con 1..n
Next i
v = v / mcd(v, f) ‘Simplificación de la fracción
numearmonico = v
End Function

Con esta función es fácil reproducir de nuevo los valores de los numeradores:

Todo lo anterior tenía por objeto que practicáramos un poco con cálculos y algoritmos, pero si hubiéramos acudido al lenguaje PARI, la generación de los numeradores es tan sencilla que se limita a su definición:

for(n=1,30,print1(numerator(sum(i=1, n, 1/i)),", "))

Obtenemos así rápidamente los primeros 30 términos:

Propiedades

Estos numeradores son también interesantes por sus propiedades. Algunas sorprenden, como que el primo p al cuadrado divide al término a(p-1) para p>3 (teorema de Wolstenholme)

Lo vemos en la tabla:

En ella se observa cómo los cuadrados de 5, 7, 11, 13 y 17 dividen a los términos de orden 4, 6, 10, 12 y 16 respectivamente. Sorprendente.

En http://oeis.org/A001008 puedes leer otras propiedades similares.

Desarrollo por Taylor

La función generatriz de los números armónicos es log(1-x)/(x-1). Podemos aprovecharla para generar de nuevo estos numeradores mediante PARI. Basta desarrollar la F.G usando el desarrollo de Taylor:

print(taylor(log(1-x)/(x-1),x,20))

Obtenemos:

Y basta leer los denominadores.

Números armónicos generalizados

Todo lo que hemos estudiado se puede extender al caso en el que la suma de fracciones se efectúe con potencias en los denominadores:

Si repetimos el algoritmo que hemos usado, pero con potencias, resultarían los numeradores generalizados, que se obtienen cambiando los sumandos 1/n por 1/n^m en el algoritmo presentado más arriba.

En Hoja de Cálculo:

Public Function numearmonico(n, m) ‘Valores de n y el exponente m
Dim i, v, f
v = 0
f = 1
For i = 1 To n: f = f * i ^ m: Next i
For i = 1 To n
v = v + f / i ^ m
Next i
v = v / mcd(v, f)
numearmonico = v
End Function

Y en PARI:

for(n=1,30,print1(numerator(sum(i=1, n, 1/i^m)),", "))

Resultan una serie de sucesiones, a las que OEIS da el nombre de “Wolstenholme numbers”, y que tienen propiedades similares a los numeradores estudiados. Las primeras, según su exponente son:

Con el 2 

http://oeis.org/A007406

1, 5, 49, 205, 5269, 5369, 266681,… y en ellas el primo p divide a a(p-1). Por ejemplo, el 5 divide a a(4)

Con el 3

http://oeis.org/A007408

1, 9, 251, 2035, 256103, 28567, 9822481,…

Con 4

http://oeis.org/A007410

1, 17, 1393, 22369, 14001361, 14011361,…

Con 5

http://oeis.org/A099828

1, 33, 8051, 257875, 806108207, 268736069,…

Todos tienen propiedades similares, según puedes consultar en los enlaces. Podemos resumir todos en una tabla de hoja de cálculo usando la función antes definida. Aquí tienes los primeros:

Trocear y desplazar

En esta entrada estudiaremos algunos números enteros que mantienen una propiedad (ser primos, cuadrados, triangulares,…) si rotamos algunas de sus cifras, es decir, que algunas situadas a la derecha las desplazamos a los primeros lugares de la izquierda, o la operación inversa, rotar de izquierda a derecha. Por brevedad, sólo consideraremos el desplazamiento de derecha a izquierda.

Por ejemplo, 622521 es cuadrado perfecto, porque 622521=789^2.  Si ahora rotamos las cifras de la derecha 21 y las situamos al principio, resulta otro cuadrado: 216225=465^2. Recorreremos algunos casos concretos, diferenciándolos por la propiedad que conservan y por el número de cifras rotadas

Cuadrados con rotación de una cifra

El cuadrado de 13516 es 182682256. Si desplazamos el 6 al primer lugar queda otro cuadrado: 618268225=24865^2. Todos los números que poseen una propiedad similar están recogidos en http://oeis.org/A045877 y los primeros son: 1, 2, 3, 16, 21, 31, 129, 221, 247, 258, 1062, 1593, 1964, 2221, 13516, 17287, 18516,… En esta lista figuran las bases de los cuadrados que presentan esa propiedad, no los cuadrados.

Este ejemplo ya publicado nos servirá de guía para presentar la función rotate(m,n), que rota n cifras de la derecha en el número m. La idea de la función es calcular en primer lugar el módulo de m respecto a 10^n, porque ese será el número formado por las cifras que rotan. Una vez obtenido, lo multiplicamos por 10^(número de cifras –n), lo que sitúa esas cifras en la parte izquierda. Para rellenar la parte derecha, tomaremos las otras cifras, que son el cociente entero del número entre 10^n. Esas no las multiplicamos porque quedan a la derecha. La expresión completa sería a = (m Mod 10 ^ n) * 10 ^ (k - n) + m \ 10 ^ n.

Adjuntamos la función completa en Basic VBA:

Public Function rotate(m, n)

Dim a, k

If m < 10 ^ n Then rotate = m: Exit Function

k = Int(Log(m) / Log(10) + 1)

a = (m Mod 10 ^ n) * 10 ^ (k - n) + m \ 10 ^ n

rotate = a

End Function

Hemos usado módulo y división entera. Otra técnica sería convertir el número en una cadena, trocearla, rotarla después y volverla a convertir en número. Esa técnica no nos viene bien en hoja de cálculo, porque introduce unos molestos espacios en blanco. Sí es adecuada en PARI. El listado siguiente lo podrás entender con facilidad:

rotate(a,p)=if(a<10^p, a, eval(Str(a%10^p, Str(a\10^p))))

La función Str convierte en cadena tanto el módulo a%10^p como el cociente entero a\10^p. El paréntesis que los une sirve para concatenar las cadenas en orden inverso, y eval las convierte de nuevo en número.

Con el siguiente código podemos reproducir la sucesión presentada arriba, de números cuyos cuadrados siguen siéndolo si rotamos una cifra:

rotate(a,p)=if(a<10^p, a, eval(Str(a%10^p, Str(a\10^p))))

for(i=1,10^5,s=i*i;r=rotate(s,1);if(issquare(i)&&issquare(r),print1(i,", ")))

Aquí tienes el resultado, que coincide con el de arriba.

Cuadrados con rotación de dos cifras

Si la función rotate la usamos con dos cifras podemos descubrir qué cuadrados siguen siéndolo si trasladamos sus dos últimas cifras de derecha a izquierda. Para simplificar el problema eliminaremos de la lista los que terminen en dos ceros, que introducirían demasiados elementos en el listado. Basta concretar el parámetro de cifras en dos y exigir que el módulo de cada candidato respecto a 100 no sea cero.

Con hoja de cálculo descubrimos sólo tres casos de tres cifras con esas condiciones: 144, 196 y 625, que se convierten en los cuadrados 441, 961 y 256 respectivamente. Con cuatro cifras no aparecen, y con cinco 21609, 38416, 60025 y 86436.

Para mayor número de cifras usaremos PARI. El código es similar al que usamos más arriba, con dos ligeros retoques:

rotate(a,p)=if(a<10^p, a, eval(Str(a%10^p, Str(a\10^p))))

for(i=10,10^5,if(i%10>>0,n=i*i;r=rotate(n,2);if(issquare(r),print1(n,", "))))

Con él aparecen todos los casos hasta 10^10:

144, 196, 625, 21609, 38416, 60025, 86436, 120409, 161604, 236196, 272484, 363609, 436921, 481636, 622521, 638401, 646416, 870489, 904401, 1560001, 6240004, 14010049, 20602521, 29822521, 77422401, 82410084, 1632321604, 3672723609, 6023622544, 6529286416, 9089524921,…

Es atractivo 60025=245^2, que se convierte en 25600=160^2

Puedes adaptar el código para encontrar cuadrados que con la rotación de tres cifras sigan siendo cuadrados. Los primeros son estos:

16384, 25600, 36864, 61009, 2896804, 7049025,…

Y aquí tienes los primeros que rotan cuatro cifras. Los incluimos como un reto para reproducirlos:

166464, 210681, 214369, 216225, 364816, 690561, 842724, 898704, 962361, 1560001, 2010724, 6240004, 8042896, 14010049, 20412324, 23242041, 32410249, 56040196, 81649296, 82410084, 92968164,…

Triangulares con rotación de cifras

Se puede repetir el estudio con los triangulares. Con rotación de una cifra ya están publicados en http://oeis.org/A068071:

1, 3, 6, 10, 55, 66, 210, 253, 666, 780, 4851, 8911, 15400, 17766, 125751, 303810, 407253, 906531, 1125750, 2648451, 5925403, 6706953, 15772536, 22207780, 50085036, 78406503, 438094800, 1623331710, 1764803755,…

Por ejemplo, 407253 es el triangular número 902, porque  407253=902*903/2, y al desplazar la cifra 3 queda 340725, triangular número 825.

Para que podamos reproducirlos necesitamos recordar que un número N es triangular si 8*N+1 es un cuadrado. Todo lo que hemos presentado en PARI hasta ahora nos vale, sustituyendo issquare(i) por issquare(8*i+1)

Así, la sucesión anterior la podemos reproducir así (en sus primeros términos):

rotate(a,p)=if(a<10^p, a, eval(Str(a%10^p, Str(a\10^p))))

for(i=1,10^4,n=i*(i+1)/2;r=rotate(n,1);if(issquare(8*r+1),print1(n,", ")))

Con rotación de dos cifras, a partir de 100 resultan:

300, 325, 666, 5050, 5151, 5565, 6105, 6555, 13530, 14196, 36046, 78606, 187578, 395605, 508536, 837865, 975106, 1400301, 7340196, 9135675, 14426506, 24580566, 36265386, 38434528, 209582101, 338715378, 495070311, 2936270028, 4201373611, 5664257830, 5735794065, 8997105153, 9652787040, …

Con el código en PARI

rotate(a,p)=if(a<10^p, a, eval(Str(a%10^p, Str(a\10^p))))

for(i=24,10^5,n=i*(i+1)/2;r=rotate(n,2);if(issquare(8*r+1),print1(n,", ")))

Prueba a reproducir esta sucesión, y comprueba algún caso, como el de 187578, que es el triangular número 612, mientras 781875 es el número 1250.

Y con tres a partir de 1000:

1035, 1485, 1891, 30135, 46360, 60031, 75078, 96141

Rotación de primos

Con una cifra están publicados en A234901:

2, 3, 5, 7, 11, 13, 17, 31, 37, 71, 73, 79, 97, 113, 131, 173, 197, 199, 271, 277, 311, 313, 337, 373, 379, 397, 419, 479, 491, 571, 577, 593, 617, 631, 673, 719, 733, 811, 839, 877, 911, 919, 971, 977…

Para no cansar, reproducimos los primos que rotan dos cifras:

Con dos: 101, 103, 107, 113, 127, 131, 149, 157, 163, 181, 191, 197, 199, 307, 311, 317, 331, 337, 359, 367, 373, 701, 709, 719, 727, 733, 739, 757, 761, 787, 797, 907, 919, 937, 941, 947, 971, 983, 991

Podemos usar la siguiente rutina en PARI para ir descubriendo rotaciones de primos con más cifras. Basta ir cambiando el valor del segundo parámetro de rotate. La siguiente sirve para obtener las rotaciones de dos cifras:

rotate(a,p)=if(a<10^p, a, eval(Str(a%10^p, Str(a\10^p))))

forprime(i=101,10^5,r=rotate(i,2);if(isprime(r),print1(i,", ")))

Comprueba estos listados:

Primos con rotación de tres cifras:

1103, 1109, 1123, 1163, 1181, 1193, 1301, 1303, 1319, 1321, 1327, 1361, 1777, 1783, 1907, 1913, 1931, 1933, 1949, 1951, 1979, 1987, 1993, 1997, 2113, 2131, 2161, 2311, 2333, 2339, 2347, 2377, 2381, 2389, 2393, 2399, 2707, 2713, 2729, 2741, 2777, …

Primos con rotación de cuatro cifras:

10007, 10069, 10091, 10099, 10103, 10151, 10181, 10193, 10211, 10253, 10259, 10267, 10271, 10273, 10301, 10333, 10337, 10357, 10369, 10391, 10427, 10459, 10487, 10501, 10559, 10601, 10613, 10627, 10631, 10657, 10687, 10691, 10733,…

Y así podemos seguir.

Máximo producto en la partición de un número (2)

En la entrada anterior estudiamos la función que asigna a cada número entero positivo el máximo producto formado entre los sumandos de todas su particiones. Mediante recurrencia calculamos los valores de la función correspondiente, a la que llamamos MPP, pudiéndose formar tablas similares a la siguiente:

El objetivo de esta segunda entrada sobre el tema es demostrar que todos los valores de la función se calculan con una de estas tres expresiones: 3^n, 2*3^n, 4*3^n, de forma cíclica salvo el valor inicial 1.

Factores del producto máximo

Razonaremos a continuación que los valores de mpp(n) (n>4) solo poseen los factores 2 y 3. Todo se basa en que si un número N se descompone en dos sumandos a y b, el máximo producto a*b se produce en el interior del intervalo (1,N), debido a que la función x(N-x) presenta el máximo en  N/2. Consecuencia de esto es que si una partición recibe un sumando nuevo para crear otra partición superior, ese sumando se puede ir descomponiendo en sumandos 1, 2 y 3, de forma que el producto aumente. Por ejemplo, si el nuevo sumando es 6, se puede sustituir por 3+3, cuyo producto es 9, superior al 6. Como el 6 es factor del producto, también lo será el 9, y el producto aumentará.

Con sumandos mayores se puede proceder de igual forma, Por ejemplo, el 10 se puede sustituir por 3+3+3+1, con lo que el producto final, en lugar de ser multiplicado por 10, lo hará en 27 (obsérvese que el 1 no interviene en el proceso).

Todos los valores del producto máximo presentarán la descomposición del tipo 2^p*3^q.

Concretemos algo más:

Caso 1: números del tipo 3k (múltiplos de 3)

Los primeros valores, según la tabla anterior, son: mpp(3)=3, mpp(6)=9, mpp(9)=27, es decir, las primeras potencias de 3. Si pasamos al valor 12, según el razonamiento de los párrafos anteriores, el producto máximo recibirá el factor 3, luego se dará que mpp(12)=3^4=81 y mpp(15)=3^5=243.

Si N=3*k, su función mpp(N)=3^k

Caso 2: Números del tipo 3k+2

Para ellos es inútil acudir al anterior 3k+1, pues el producto sólo perdería el factor 1 que no incrementa su valor, pero si acudimos al máximo de 3k=3^k, el nuevo factor máximo es el 2, y queda mpp(3k+2)=2*3^k.

Por ejemplo, mpp(8)=mpp(3*2+2)=2*3^2=18, como ya sabemos por la entrada anterior. Otro: mpp(14)=mpp(3*4+2)=2*3^4=162.

Si N=3*k+2, su función mpp(N)=2*3^k

Caso 3: Números del tipo 3k+1

Si acudimos al valor 3(k-1)+2, basta usar como nuevo sumando el 2 para ver que el máximo que buscamos es 4*3^(k-1).

Vemos algún ejemplo: mpp(13)=mpp(3*4+1)=4*3^(4-1)=4*27=108, como ya sabemos por casos anteriores. El valor de mpp(31) = mpp(3*10+1) = 4*3^9 = 4*19683 = 78732.

Lo hemos comprobado con la definición de función de la entrada anterior:

Si N=3*k+1, su función mpp(N)=4*3^(k-1)

Nueva definición de la función

Estas propiedades que acabamos de estudiar nos permiten simplificar mucho la definición de la función MPP. La presentamos en varias versiones para que las analices. La primera está construida con las funciones de Excel o Calc. Es esta, que la hemos copiado actuando sobre la celda D11:

=3^(ENTERO((D11)/3)-(RESIDUO(D11;3)=1))*2^RESIDUO(-D11;3)

Parece difícil. Intenta comprenderla. La primera parte decide a qué exponente elevaremos el 3, y la segunda el mismo problema para el 2. Para el primero hallamos el valor de k que acompaña al 3 y le restamos un 1 en el caso 3k+1. Para el segundo elegimos 1, 2 o 4 según el residuo respecto al 3, pero con signo menos. Ya, es complicado.

Con ella hemos construido esta tabla, para los números 40 a 45:

Código en Basic VBA:

Es quizás más fácil de entender la definición (segunda ya) de mpp en VBA. Como ya tenemos una versión, le llamaremos a esta MPP_2. Su desarrollo puede ser este:

Public Function mpp_2(n)
Dim k, m, p

k = Int(n / 3) ‘se determina el valor de k en n=3k+b
If n Mod 3 = 1 Then m = 1 Else m = 0 ‘si es múltiplo de 3, el exponente de 3 disminuye en 1
p = 3 ^ (k - m) ‘parte correspondiente al factor 3
m = 3 - n Mod 3: If m = 3 Then m = 0 ‘se prepara el exponent del 2
mpp_2 = p * 2 ^ m ‘se ensambla la función
End Function

Este código es más simple y rápido que el anterior. Hemos comprobado su equivalencia, y se demuestra aquí el poder simplificador del razonamiento matemático.

Para quienes entiendan el lenguaje PARI, copio aquí la solución de M. Somos:

 mpp(n) = floor( 3^(n - 4 - (n - 4) \ 3 * 2) * 2^( -n%3))

Es muy parecida a la que proponemos.

Otra recurrencia

Si repasas la página http://oeis.org/A000792 que nos viene ayudando en este estudio, podrás leer más propiedades y fórmulas sobre estos productos máximos.

Una muy sencilla y curiosa es la de Ivan Neretin:

a(n) = a(n-1) + largest proper divisor of a(n-1), n > 2

En efecto, estudiamos los tres casos:

Si a(n-1)=3*k, sabemos que mpp(a(n-1))=3^k, y su mayor divisor propio md=3^(k-1). Sumamos y obtenemos mpp(a(n))=mpp(3*k+1)=3^k+3^(k-1)=4*3^(k-1), como ya sabemos.

Si a(n-1)=3*k+1, se tendrá: mpp(a(n-1))=4*3^(k-1), su mayor divisor propio md=2*3^(k-1). Sumamos: mpp(a(n))=mpp(3*k+2)=4*3^(k-1)+2*3^(k-1)=2*3^k, expresión ya conocida.

Si a(n-1)=3*k+2, mpp(a(n-1))=2*3^k, su mayor divisor propio md=3^k y al sumar obtenemos la esperada 3^(k+1).

Si dispones de la función en VBA de “mayor divisor”, en una columna de hoja de cálculo puedes engendrar la sucesión completa. En este blog disponemos de la función MAYORDIV (no la desarrollamos porque acude a otras más complicadas), lo que nos permite desarrollar la recurrencia:

Se ha adosado la recurrencia a la derecha de una tabla de MPP, para ver la equivalencia.

Si no dispones de la función MAYORDIV, puedes copiar la que sigue:

Public Function mayordiv(n)
Dim f, m
Dim es As Boolean

f = 2
m = 1
If n / 2 = n \ 2 Then es = True: m = 2 Else es = False
While f * f <= n And Not es
If n / f = n \ f Then es = True: m = f
f = f + 1
Wend
If m = 1 Then mayordiv = 1 Else mayordiv = n / m
End Function

Máximo producto en la partición de un número (1)

Ya sabemos que una partición de un número entero positivo N es una suma de números también enteros positivos cuyo resultado es ese número. Ya hemos tratado en este blog el tema de las particiones, y lo volveremos a desarrollar próximamente. Si no tienes claro el concepto puedes acudir a las direcciones

http://hojaynumeros.blogspot.com.es/2011/02/particiones-de-un-numero.html

https://es.wikipedia.org/wiki/Partici%C3%B3n_(teor%C3%ADa_de_n%C3%BAmeros)

Está muy estudiado el tema del desarrollo de las particiones y el del cálculo de su número. Aquí nos interesará el máximo producto que se puede lograr multiplicando los sumandos de cada partición. Lo introducimos con ejemplos:

Máximo producto logrado con particiones

Tomemos el número 6. Mentalmente se pueden escribir sus particiones (no se tiene en cuenta el orden): 6=5+1=4+2=3+3=4+1+1=… En cada partición calculamos el producto entre sumandos: 6, 5, 8, 9, 4,… y nos quedamos con el máximo. En el esquema lo verás mejor:

Figuran las once particiones del 6 y en la columna de la derecha los productos de sumandos. En la partición de sumando único lo elegimos como producto. Se observa que el máximo producto es 9. Como el resultado es único, constituye una función del número elegido, que podríamos escribir como MPP (Máximo Producto en Particiones) y se tendría que MPP(6)=9.

Otro ejemplo: En PARI las particiones se obtienen con la función partitions. Si pedimos las particiones del número 8 las obtenemos como vectores independientes:

Si multiplicas los sumandos dentro de cada vector descubrirás que el máximo producto es 18=2*3*3, luego MPP(8)=18

Estos resultados figuran en la página de OEIS http://oeis.org/A000792 con una definición recursiva que ya trataremos:

1, 1, 2, 3, 4, 6, 9, 12, 18, 27, 36, 54, 81, 108, 162, 243, 324, 486, 729, 972, 1458, 2187, 2916, 4374, 6561, 8748, 13122, 19683, 26244, 39366, 59049, 78732, 118098, 177147, 236196, 354294, 531441, 708588, 1062882, 1594323, 2125764, 3188646, 4782969,…

Como era de esperar, los valores son crecientes (cada uno es un máximo que se apoya en los anteriores) y pronto adquieren una buena tasa de crecimiento. La página citada contiene una fórmula recursiva que es fácil de entender.

a(n) = max{ (n-i)*a(i) : i<n}; a(0) = 1

Se define a(0) como 1, y también es fácil entender las siguientes: a(1)=1, a(2)=2, a(3)=3,… La recursividad tampoco es difícil de captar. Se trata de multiplicar cada valor menor que n por el máximo correspondiente a su diferencia con n. En efecto, las particiones de n se forman eligiendo esos valores i:1…n-1 para luego unirlos con las particiones de n-i. Por ejemplo, en las particiones del 7, si elegimos el valor 4, se deberá combinar con las particiones del 3 para formar las particiones del 7 que contengan un 4. Igual ocurre con todos los valores: 5 se añadirá a las particiones del 2 y 3 se unirá a las particiones de 4.

Si conservamos los valores máximos de cada partición de n-i, al multiplicarlos por i resultarán productos de la partición superior, candidatos a ser máximos. Al recorrer todos los valores menores que n dispondremos de n-1 posibles máximos, y uno de ellos será el MPP.

Algoritmo de construcción de la función MPP

Las ideas anteriores nos permitirán construir la función  MPP. En VBA de Excel se puede usar esta definición de función:

Public Function mpp(n)
Dim mx, i, m, j, mm
Dim a(50) ‘Está preparado para n<=50. Se puede ampliar a otro número

If n = 0 Or n = 1 Then mpp = 1: Exit Function
If n = 2 Then mpp = 2: Exit Function ‘Casos particulares
a(0) = 1: a(1) = 1: a(2) = 2: mx = 2
If n > 2 Then
For i = 1 To n ‘Se recorren los valores anteriores para la recursión
m = 1 ‘Valor provisional del máximo
For j = 1 To i
mm = j * a(i - j)
If m < mm Then m = mm ‘Se busca un máximo nuevo mediante los productos con los anteriores
Next j
mx = m: a(i) = m ‘Se incorpora el máximo a la lista
Next i
End If
mpp = mx ‘Máximo final
End Function

Con esta función podemos encontrar el máximo producto entre particiones de cualquier número. Si es mayor que 50 bastará cambiar la dimensión del vector de máximos. En la tabla siguiente hemos recogidos los valores de MPP para los números comprendidos entre 30 y 40:

Para quienes conozcan el lenguaje PARI (gratuito y muy recomendable para estos temas) se inserta un código para esta función, que también devuelve los valores entre 30 y 40:

mpp(n)=my(a=vector(50), m, mm, mx=2,mp=1);a[1]=1;a[2]=2;if(n<2,mp=1,if(n==2,mp=2,for(i=3,n,m=1;for(j=1,i, d=i-j;if(d>0,mm = j * a[d],mm=j);if(m<mm,m=mm));mx=m;a[i]=m));mp=mx;mp)
for(k=30,40,print1(mpp(k),", "))

Aquí tienes el resultado, que coincide con el de Excel:

Interpretación algebraica

Estos valores coinciden con los cardinales máximos de los subgrupos del grupo simétrico S(n). Usando la descomposición en ciclos se les puede dar un significado (ver http://hojaynumeros.blogspot.com.es/2013/10/ciclos-2-descomposicion-en-ciclos.html)

Por ejemplo, en el caso del 8 visto más arriba, mpp(8)=18, y puede dársele el significado de que es el cardinal máximo de un subgrupo propio del grupo de permutaciones de 8. En concreto, usando la descomposición en ciclos, podría ser G=(1,2)(3,4,5)(6,7,8) o cualquiera de sus isomorfos. En el caso del 6 es fácil ver que el subgrupo maximal es el GM=(1,2,3)(4,5,6), de cardinal 3*3=9, que es el valor de mpp(6).

Existe una forma directa y simple para calcular mpp(n), sin recurrencias ni algoritmos. Como es un cambio importante en el desarrollo que hemos llevado hasta ahora, lo dejamos para la próxima entrada.

Conjetura de Rassias

Esta conjetura recibe el nombre de su autor, M. Th. Rassias, que la enunció siendo muy joven, mientras preparaba una Olimpiada Matemática. Se puede formular de varias formas, pero la que preferimos es la siguiente:

Para cada número primo p>2 existen dos primos p₁ y p₂, con p₁<p₂ tales que

(p-1)p₁=p₂+1

Es decir, que si el primer primo lo multiplicamos por p-1, conseguimos un número al que precede otro número primo. Por ejemplo:

Para el número 17, el par de primos puede ser 2 y 31, porque (17-1)*2=32=31+1. Para el primo 47 los primos pueden ser 3 y 137, porque (47-1)*3=138=137+1

La conjetura afirma que siempre se pueden encontrar esos dos primos para uno dado.

Obtención con hoja de cálculo

En teoría se podría comprobar esta conjetura mediante una tabla de doble entrada con los primeros números primos, pero sería un procedimiento costoso en espacio y tiempo. Es preferible acudir al VBASIC o lenguaje similar. En Excel puedes intentar una función que nos devuelva el más pequeño q de los dos primos, suficiente para comprobar la conjetura, ya que el otro lo podemos calcular mediante (p - 1) * q - 1

Public Function rassias(p)
Dim a, q

If Not esprimo(p) Then rassias = 0: Exit Function  ‘Si no es primo nos devuelve un cero
q = 2 ‘Posible valor del primo más pequeño
a = 0 ‘Si a=0 significa que aún no se han encontrado los primos
While a = 0
If esprimo((p - 1) * q - 1) Then ‘Prueba para saber que se encontraron los primos
a = q
End If
q = primprox(q) ‘Se prueba con el siguiente primo
Wend
rassias = a ‘Se encontró el primo menor
End Function

Las funciones ESPRIMO y PRIMPROX las puedes copiar desde nuestra entrada

http://hojaynumeros.blogspot.com.es/2012/04/proposito-de-ormiston.html

Con esta función y el cálculo posterior podemos construir una tabla en la que para cada número primo contenga los dos primos más pequeños que verifican la conjetura. Sería similar a esta:

Programa en PARI

Para quienes conozcan el lenguaje PARI, con este programa comprobamos la conjetura para los primos inferiores a 200:

p=2;while(p<200,p=nextprime(p+1);q=2;a=0;while(a==0,b=(p-1)*q-1;if(isprime(b),a=q);q=nextprime(q+1));print(p,", ",a,", ",b))

Resultado:

Con los cambios oportunos se puede lograr la comprobación para otros conjuntos de primos.

Otros puntos de vista

En la tabla anterior destaca la frecuencia con la que aparecen los valores 2 y 3 para el primo más pequeño. Es una indicación de que la conjetura no es algo complicado, sino que se comprueba fácilmente para valores pequeños. Podemos plantear una búsqueda para saber cuándo aparecerán otros valores, si es que lo hacen. Aquí tienes los resultados para la primera aparición de otros primos:

Esta tabla sugiere que la conjetura también se cumple para todo p₁. El problema radica en que no hay tope en la búsqueda de p y de  p₂ , por lo que de no cumplirse para algún valor, entraríamos en un bucle sin fin. No obstante, lo intentamos con esta función:

Public Function rassias2(p)
Dim q, b, a

If Not esprimo(p) Then rassias2 = 0: Exit Function
q = 2
a = 0
While a = 0
b = (q - 1) * p - 1
If esprimo(b) Then
a = b
End If
q = primprox(q)
Wend
rassias2 = a
End Function

Con ella vemos que a todo valor de  p₁ le corresponde otro de  p₂. No tienen que resultar los mismos valores anteriores, porque al cambiar el punto de vista se encuentran otros mínimos, pero lo importante es que existe siempre una solución.

Aquí tienes un resultado:

Por ejemplo, para 2081 como  p₁, el valor de  p₂ es 20809, calculado mediante p=11, ya que 20809=2081*10-1, y 20809 es primo.

No resistimos la elaboración de una función para p₂:

Public Function rassias3(p)
Dim q, b, a

If Not esprimo(p) Then rassias3 = 0: Exit Function
q = 2
a = 0
While a = 0
b = (p + 1) / q + 1
If esprimo(b) Then
a = b
End If
q = primprox(q)
Wend
rassias3 = a
End Function

Con ella se puede construir una tabla que relaciones  p₂ con  p₁ y p:

Todas estas tablas se podrían prolongar hasta números mucho mayores, y siempre existe una solución de dos primos respecto al dado, luego se puede dar por comprobada la conjetura dentro de la herramienta que hemos usado.

Primos relacionados con uno fijo

Por último, nos podríamos plantear si para cada valor de p podemos encontrar infinitos pares  p₁ y  p₂ que cumplan la conjetura. Lo dejamos como ejercicio. En la tabla observamos pares de seis cifras que cumplen la conjetura para p=11:

Como los primos de la primera columna se multiplican por 10, los de la segunda terminan todos en 9. Este ejercicio lo podemos repetir para cualquier valor, y dentro del rango que deseemos. Terminamos con los primos de siete cifras que corresponden a p=137:

Conjetura de Collatz

Ya se trató esta conjetura en este blog, pero desde el punto de vista de su experimentación en un Taller de Matemáticas

http://hojaynumeros.blogspot.com.es/2011/05/la-conjetura-de-collatz-en-un-taller-de.html

En esta ocasión se buscarán rutinas y funciones que nos ayuden a comprobar la conjetura para muchos números, así como encontrar sus cúspides y órbitas. Se ha escrito mucho sobre esta conjetura, por lo que aquí se desarrollará sólo ese aspecto.

Planteamiento

Para quienes no conozcan esta conjetura recordaremos su planteamiento:

Se toma un número entero positivo N cualquiera, por ejemplo el 13, y se le aplica la siguiente operación, a la que llamaremos función COLL(N):

• Si el número es par, se divide entre 2.
• Si el número es impar, se multiplica por 3 y se le suma 1.

En el caso del 13, como es impar, se le aplicará la segunda, y quedará COLL(13)=13*3+1=40.

La idea de la conjetura es que sigamos aplicando esta operación a todos los resultados que obtengamos, En nuestro caso sería COLL(40)=20 (por ser par), COLL(20)=10, COLL(10)=5, COLL(5)=3*5+1=16, COLL(16)=8, COLL(8)=4, COLL(4)=2, COLL(2)=1, y a partir del 1 se entra en el ciclo {4, 2, 1}

La conjetura afirma que este final en el 1 y el ciclo posterior ocurre para cualquier otro entero positivo. Sea cual sea el comienzo, se llegará al número 1. Todas las sucesiones construidas así terminarán en el ciclo 4, 2, 1.

Lo vemos con otro ejemplo:

COLL(6)=3, COLL(3)=10, COLL(10)=5, COLL(5)=16, COLL(16)=8, COLL(8)=4, COLL(4)=2, COLL(2)=1

Insistimos en que existen muchas publicaciones sobre esta conjetura y aquí sólo nos limitaremos a pequeñas comprobaciones. La más sencilla es mediante celdas en la hoja de cálculo.

Comprobación con celdas de hoja de cálculo

Escribimos el entero positivo inicial (lo podemos nombrar como semilla) en una celda, sea por ejemplo, la D4. En la celda inferior D5 escribimos SI(RESIDUO(D4;2)=0;D4/2;3*D4+1), que divide entre 2 el valor de la D4 si es par (porque entonces se verifica RESIDUO(D4;2)=0) y lo multiplica por 3 añadiendo 1 si es impar. En la imagen vemos el resultado para el número 132:

Al ser par el 132 se ha dividido entre 2. Ahora lo único que tenemos que hacer es rellenar esa fórmula hacia abajo y parar cuando aparezca un 1:

(Troceamos la imagen porque aparecen muchos números)

La conjetura ha sido comprobada hasta números muy grandes, por lo que puedes tener la seguridad de que llegarás siempre al valor 1. Al conjunto de números que se recorren hasta llegar a ese valor le podemos llamar órbita del número dado, que aquí son los 29 números que aparecen en la imagen. Al número mayor que hayamos alcanzado en la órbita le llamaremos cúspide. En este ejemplo la cúspide es el mismo 132.

De esta forma tan simple podemos comprobar la conjetura dentro del alcance de la herramienta que usamos. Si la órbita tiene una longitud grande este procedimiento puede alargarse. Por ello acudiremos ahora a la definición de funciones:

Funciones sobre la conjetura de Collatz

Ya hemos presentado COLL(N). Sería bueno introducir su versión en VBA para poder construir sobre ella otras funciones más complicadas. Su código es muy sencillo:

Public Function coll(n)
If n / 2 = n \ 2 Then coll = n / 2 Else coll = 3 * n + 1
End Function

No necesita grandes explicaciones. La condición n/2=n\2 equivale a indicar que n es par, ya que entonces el resultado de la división n/2 es idéntico al de la división entera n\2. El resto se entiende bien. Con esta función podemos reproducir las órbitas en columna de forma idéntica a como procedimos en el primer ejemplo.

En PARI el código es similar:

coll(n)=if(n/2==n\2,n/2,3*n+1)

En la imagen vemos el resultado de pedir coll(132)

Función orbicoll

A cualquier número entero le podemos asignar una cadena que contenga todos los números por los que “pasa” hasta llegar al 1, es decir, su órbita. En VBA podía ser esta:

Public Function orbicoll(n)
Dim b
Dim s$ ‘Cadena (string) para recoger los resultados
b = n: s$ = Str$(b) ‘La cadena comienza con el número inicial (semilla)
While b <> 1 ‘Se trabaja hasta llegar al 1
b = coll(b) ‘ En cada paso se aplica la función COLL
s$ = s$ + Str$(b) ‘Se incorpora el resultado al string
Wend
orbicoll = s$
End Function

En la imagen puedes comprobar la creación de la órbita del número 132:

Después puedes usar la prestación de convertir texto en columnas (antes debes copiar la órbita en otra celda mediante copiar valores) y crear un gráfico que abarque toda la órbita:

Puedes observar que arranca en 132 y termina en 1. Los tramos ascendentes representan números impares, y los descendentes a los pares.
Puedes seguir este proceso con cualquier otro número entero y seguir la evolución de su órbita.

En la imagen aparece la órbita del número 127, más compleja y con una cúspide cercana a 4500:

Función orbicoll en PARI

Es fácil la traducción de esta función a PARI:

coll(n)=if(n/2==n\2,n/2,3*n+1)
orbicoll(n)=my(b=n,s=Str(n));while(b<>1,b=coll(b);s=concat(concat(s," "),Str(b)));s

En la imagen se ha pedido la órbita de 127:

Estudia al código siguiente para la función lorbicoll, que devuelve el número de elementos de una órbita:

Public Function lorbicoll(n)
Dim a, b
a = 1: b = n
While b <> 1
b = coll(b)
a = a + 1
Wend
lorbicoll = a
End Function

Con ella podemos comprobar lo que ya sabemos, que 132 tiene una órbita de 29 elementos.

Con la versión PARI puedes abordar casos con números mayores:

coll(n)=if(n/2==n\2,n/2,3*n+1)
lorbicoll(n)=my(b=n,a=1);while(b<>1,b=coll(b);a+=1);a

Te proponemos comprobar que el número  871 es el que posee la órbita de más longitud entre los de tres cifras, y que contiene 179 elementos. De los de cuatro cifras el de órbita de más longitud es el número 6171, con 262 elementos.

Función cuspicoll

Del mismo modo que construimos la órbita de un número entero positivo, podemos encontrar su cúspide. El procedimiento será similar, pero, en lugar de añadir resultados a un string, tomaremos nota en cada paso del máximo valor que ha aparecido:

Public Function cuspicoll(n)
Dim a, b
a = n: b = n
While b <> 1
b = coll(b)
If b > a Then a = b
Wend
cuspicoll = a
End Function

La instrucción clave es If b > a Then a = b, que convierte a en el nuevo máximo si aparece un elemento b mayor que los precedentes.

Con las funciones definidas podemos construir un esquema en el que se analice el comportamiento de la conjetura de Collatz para una semilla dada:

Finales previsibles

Algunos tipos de números presentan una órbita bastante previsible. Por ejemplo:

Potencias de 2: a partir de ellos se entra en una ruta descendente y previsible que finaliza en el 1.

Números tipo (2^n-1)/3: Desembocan en una potencia de 2, por lo que también inician una ruta directa,  y esto ocurre una potencia sí y otra no, porque 2^n es del tipo 3k+1 o 3k+2 y al multiplicar por Si es 3k+1 es candidato a que (2^n-1)/3 sea entero, y si es del tipo 3k+2, la siguiente potencia será 6k+4, o sea del tipo 3k+1

Otros tienen recorrido corto, como el 6, el 10 o el 20.

Es normal que pensemos en que muchas órbitas pasarán por ellos, y existan pares de órbitas  que pasan ambas por el mismo punto de entrada, más o menos primario.

Podemos intentar ver si dos números presentan alguna coincidencia en sus órbitas, porque entonces compartirán final. No es difícil programar una función que nos devuelva un punto de coincidencia en las órbitas de dos números. En primer lugar necesitamos una función que nos indique si el número n pertenece a la órbita del número m. Puede ser esta:

Public Function enlacoll(m, n) As Boolean
Dim e As Boolean
Dim p
If m = n Then
e = True
Else
e = False
p = m  ‘p recorrerá la órbita de m
While Not e And p <> 1
p = coll(p)
If p = n Then e = True  ‘Si p es igual a n, sí pertenece
Wend
End If
enlacoll = e
End Function

Con esta función enlacoll(m,n) podemos saber si n pertenece a la órbita de m. Se puede organizar un esquema de cálculo:

En la imagen se ha verificado que 52 pertenece a la órbita de 57.

Coincidencia en las órbitas

Con la anterior función ya estamos preparados para encontrar la primera coincidencia entre dos órbitas. Si no existe una coincidencia anterior, se nos devolverá un cero. El código de la función es:

Public Function coincicoll(m, n)
Dim c, q
q = m
c = 0
While q > 1 And c = 0 ‘q recorre la órbita de m
If enlacoll(n, q) Then c = q ‘si q pertenece a la órbita de n, hay coincidencia
q = coll(q)
Wend
coincicoll = c
End Function

Con ella también podemos construir otro esquema de cálculo. En la imagen se comprueba que 125 y 126 comparten una subórbita de 95 elementos, que comienzan en 364.

Con estas ideas puedes construir otras muchas funciones, o emprender otras búsquedas. Sólo se ha pretendido en esta entrada dar ideas para comprobaciones y experimentaciones sobre la conjetura, ya de por sí bastante estudiada en otros aspectos.

Hipotenuseando

En mis exploraciones por la página OEIS (Enciclopedia On-line de sucesiones de números enteros, http://oeis.org/?language=spanish), me encontré con la sucesión http://oeis.org/A104863

10, 30, 31, 43, 53, 68, 86, 109, 138, 175, 222, 282, 358, 455, 578, 735, 935, 1189,…

En ella, a partir de los valores a(1)=10 y a(2)=30, se van formando los siguientes como la parte entera de la raíz cuadrada de la suma de cuadrados de los dos anteriores.

Así, el tercer elemento es igual a

ENTERO(RAIZ(10^2+30^2))=ENTERO(RAIZ(1000))=ENTERO(31,62)=31.

Prueba a justificar que el siguiente es 43.

Esta definición equivale a que cada término es la hipotenusa truncada correspondiente a los términos anteriores tomados como catetos. Por eso le hemos llamado a esta sucesión la de “hipotenusear”. Su expresión recurrente sería:

Esta sucesión presenta algunas características que le hacen merecer esta entrada:

(1) El uso de la parte entera obliga a renunciar a las fórmulas teóricas. De hecho, en términos avanzados de la sucesión, el cumplimiento del Teorema de Pitágoras es deplorable. Observa este ejemplo, en el que a(n-2)=8348, a(n-1)= 10618 y a(n)=13506. Si elevamos al cuadrado obtenemos:

13506^2=182412036
8348^2+10618^2=182431028

Restamos y nos queda un error de 18992 unidades. Así que las hipotenusas que obtengamos con esta recurrencia no son tales, aunque seguiremos llamándolas así.

(2) Como también ocurre en las recurrencias lineales, el cociente entre dos términos consecutivos se va estabilizando y tiende a un límite. Esto nos permite “razonar en el límite”, aunque sepamos que es una técnica aproximada, que sólo nos valdrá para explicar (y no demostrar) algunas conjeturas que se han afirmado para esta sucesión y otras similares.

(3) La sucesión presentada es sólo un caso particular de toda una familia en la que podemos fijar a(1) y a(2) como deseemos. La mayoría de las propiedades se mantendrán. Vemos la primera:

Conjetura: El límite del cociente a(n+1)/a(n) es la raíz cuadrada del número áureo.

La llamamos conjetura por causa de la parte entera, que nos impide mejores razonamientos. Esta cuestión, de manejarnos entre aproximaciones y conjeturas, es uno de los objetivos de esta entrada.

Podemos comprobar lo anterior con hoja de cálculo. Escribimos dos catetos uno debajo de otro, como 2 y 5, y después, en columna rellenamos la fórmula

=ENTERO(RAIZ(CATETO1^2+CATETO2^2)).

No es difícil de organizar. Después, en la columna de la derecha calculamos los cocientes entre dos términos consecutivos. Algo así:

Al llegar al término 14 ya se adivina el valor deseado. Si seguimos bajando, la aproximación mejora mucho

Podemos razonarlo en el límite. Llamamos k al cociente a(n+1)/a(n). Por tanto, en la expresión de a(n) podemos escribir:

O bien, pasando a(n-1) al primer miembro,

Elevando al cuadrado y agrupando, tenemos que k se debe aproximar a la solución de la ecuación k⁴ – k² – 1 = 0, una bicuadrada cuya solución es el límite sugerido, la raíz cuadrada del número áureo.

Incluimos las cuatro soluciones tal como las da WolframAlpha:

Elegimos la real positiva, y, efectivamente, resulta 1,27201964951407…

Manteniendo el razonamiento en el límite, si a(n-1) y a(n) se comportan como cateto e hipotenusa respectivamente con esa razón dada, el otro cateto, a(n-2), se podrá aproximar (también en el límite) de esta forma:

Plantéate como ejercicio demostrar el último paso. Recuerda que F-1=1/F
En esta sucesión a(n) tiende en el límite a a(n-2)*F

Lo hemos demostrado en el párrafo anterior. También lo podemos razonar mediante la idea de que si el cociente entre dos términos consecutivos se aproxima a la raíz del número áureo, el correspondiente a a(n) y a(n-2) será dicho número F.

Por tanto, en el límite, cada tres términos consecutivos forman un triángulo rectángulo cuyos lados son proporcionales a (1, 1,272019…, 1,618033…) y cuyo ángulo menor es de 38,17º.

Un ejercicio: ¿Cuál es, en el límite, el cociente entre el área del triángulo (a(n+1), a(n), a(n-1)) y el correspondiente a (a(n), a(n-1), a(n-2))?

Para quienes conozcáis el lenguaje PARI, con una línea de código similar a esta podéis estudiar la sucesión hasta términos más avanzados:

a=1;b=7;for(i=1,30,c=truncate(sqrt(a^2+b^2));a=b;b=c;print1(c,", "))

Podéis estudiar los cocientes añadiendo el código adecuado.

Conjetura: A partir de un término mínimo, a(n) se diferencia de a(n-2)+a(n-4) en a lo sumo una unidad.

Esta conjetura está publicada en la página OEIS citada para el caso a(1)=10 y a(2)=30, en el que la diferencia se estabiliza en 1. Su verificación no depende de los términos iniciales, salvo, quizás, el tope inferior de 1. Por ejemplo, lo comprobaremos con hoja de cálculo y los términos iniciales a(1)=4 y a(2)=7:

En este caso vemos que a(n) tiende a coincidir con a(n-2)+a(n-4)-1

En el límite se puede justificar usando todos los cocientes presentados más arriba:

a(n-2)+a(n-4) = a(n)/F+a(n)/F^2 = a(n)*(F+1)/F^2 = a(n)

Así que en el límite la coincidencia es exacta: a(n-2)+a(n-4) = a(n), y la unidad como error aparece por los truncamientos.

Puedes cambiar la función ENTERO por la de REDONDEAR. Así lo hacen las sucesiones A104803,  A104804, A104805 y A104806, con resultados similares.

¿Variaciones o combinaciones?

En el mes de julio pasado descubrí que el número 1716 equivale a un número de variaciones sin repetición y también de combinaciones sin repetición, ambas con el mismo índice superior. En efecto, 1716 = 13*12*11 = V(13,3), pero también equivale a C(13,6) o C(13,7), ya que

Se ha producido la feliz casualidad de que 13*12*11 = 6*5*4*3*2*1, y por eso se ha podido simplificar.

En esta propiedad el verdadero protagonista, a efectos de construcción de algoritmos, es el número 13, que es el que participa en ambas fórmulas, de combinaciones y variaciones. No existen muchos índices que cumplan esto. Los primeros son 8, 13, 27, 124, 725 y 5046, si imponemos la condición razonable de que el índice inferior sea mayor que 1, para evitar trivialidades.

Búsqueda “ingenua”

Para encontrar estos índices superiores podemos acudir a la definición que hemos insinuado: “Números n para los que existen dos índices k y h tales que V(n,k) = C(n,h)”. Si disponemos de las funciones C y V, basta recorrer índices para cada candidato y parar cuando se dé una coincidencia. Podemos acotar la búsqueda eligiendo para h el intervalo (2, n/2), por cuestión de simetría en los números combinatorios. Por otra parte, es claro que k ha de ser menor que h, para que tenga lugar la igualdad. En Basic de hojas de cálculo podía usarse algo así:

Public Function vari(n, k)
Dim v, i
v = 1
For i = 0 To k - 1: v = v * (n - i): Next i
vari = v
End Function

Public Function combi(n, k)
Dim v, w, i
v = 1: w = 1
For i = 0 To k - 1: v = v * (n - i): w = w * (i + 1): Next i
combi = v / w
End Function

Código para cada valor de n (aquí representado por la variable i):
a = 0 ‘variable para parar la búsqueda
k = 2
While k <= i / 2 And a = 0 ‘ se recorren los valores de k
b = combi(i, k) ‘ se encuentra el número de combinaciones b
h = 2
While h < k And a = 0 ‘se recorren los valores de h
c = vari(i, h) ‘ se encuentra el número de variaciones c
If b = c Then a = 1: m = k: n = h ‘en caso de igualdad, se para y toma nota
h = h + 1
Wend
k = k + 1
Wend

La salida será el valor m del índice k y el n del índice h. En forma de tabla estos serían los primeros valores:

Los comprobamos (salvo el 13 que ya se ha visto):

V(8,2)=8*7 = 56, C(8,3)=(8*7*6)/(3*2*1)=56

V(27,3)=27*26*25=17550, C(27,4)=(27*26*25*24)/(4*3*2*1) =27*26*25=17550

V(124,4)=124*123*122*121=225150024 = C(124,5)

Algoritmo con recursividad

En este primer intento estamos realizando más operaciones de lo debido. No es necesario calcular C(n,k) y V(n,h) en cada paso. Es mejor generar cada intento recursivamente a partir del anterior:

a = 0
k = 2
b = i
While k <= i / 2 And a = 0
b = b * (i - k + 1) / k ‘recursividad para combinaciones
h = 2
c = i
While h < k And a = 0
c = c * (i - h + 1) ‘recursividad para variaciones
If b = c Then a = 1: m = k: n = h
h = h + 1
Wend
k = k + 1
Wend

Con hoja de cálculo se llega pronto al desbordamiento de decimales. Deberemos cambiar a PARI:

for(i=3,1000,a=0;j=2;m=i;while(j<=i/2+1&&a==0,m=m*(i-j+1)/j;k=2;n=i;while(k<j&&a==0&&n<=m,n*=i-k+1;if(m==n,a=1);k+=1);j+=1);if(a==1,print(i,", ",j,", ",k,", ",m)))

Es poco legible. Contiene las mismas ideas desarrolladas con VBA, pero escritas de forma excesivamente compacta.

El resultado te será familiar, y aparece un nuevo índice, el 725:

Para seguir avanzando se requiere ya mucha paciencia, porque los cálculos se van haciendo lentos y complejos. El siguiente índice superior en aparecer es el 5046:

El valor de V(5046,7) = C(5046,8) da idea de cómo se va complicando esto. Sin embargo, nos conduce al hecho de que en

V(5046,7)=5046*5045*5044*5043*5042*5041*5040,

el último factor es el factorial de 7, lo que permite la simplificación que da lugar a la igualdad entre variaciones y combinaciones.

Casos particulares

El último ejemplo nos da una idea de la naturaleza de algunos de los índices superiores con la propiedad buscada. Es fácil entender que todo índice del tipo n!+(n-1) da lugar a un número m en el que el número de variaciones de n-1 elementos coincide con el de combinaciones de n elementos. Así ha ocurrido con muchas de las soluciones presentadas. En la siguiente tabla se han destacado en rojo las que ya conocíamos. Nos hemos detenido en el último factorial que Excel puede expresar de forma entera:

Por tanto, el número de índices adecuados es infinito, y crece a ritmo de factorial.
Los casos que faltan, como el 13, provienen de la casualidad de que un producto de números consecutivos equivalga a un factorial, que es lo que ocurre con 10*9*8 = 6! ¿Existirán más casos? Mediante una búsqueda manual descubrimos: 6*5*4=5!, lo que nos da de nuevo el candidato 8, pero esta vez con la expresión C(8,5) y la solución 56 ya vista. De ella podemos extraer la coincidencia 10*9*8*7=7!, que nos llevaría al índice 13.

Este estudio es un ejemplo más de una forma clásica de abordar problemas:

•  Usar un algoritmo sencillo, que no hace uso de propiedades especiales. Es un buen método para comenzar, pero suele ser largo y poco interesante. Así ha sido nuestro procedimiento “ingenuo”.
•  Perfeccionar el algoritmo a fin de conseguir mayor velocidad de búsqueda. En este caso se ha logrado con la recursividad.
•  Acudir a la teoría o el razonamiento. Hemos descubierto así que existen infinitos casos con la fórmula n!+n-1, más unos cuantos casos aislados. Así le hemos quitado el misterio a la cuestión planteada.

Expresión cuadrática X^2+kY^2 = N

Hace unos meses publiqué en Twiter, como una curiosidad, esta tabla de desarrollos para el número 4516

No existen muchos números que admitan esas diez expresiones cuadráticas con enteros. En concreto estos son los primeros:

1009, 1129, 1201, 1801, 2521, 2689, 3049, 3361, 3529, 3889, 4036, 4201, 4516, 4561, 4729, 4804, 5209, 5569, 5881, 6841, 7204, 7561, 7681, 8089, 8521, 8689, 8761, 8929, 9081, 9241, 9601, 9769,…

¿Qué hay detrás de esta lista?

Realmente, el problema radica en resolver la ecuación X²+kY²=N, con k>0 buscando soluciones enteras X>1 Y>1, para evitar trivialidades.

Despejando X en X²+kY²=N nos queda

Por tanto, para averiguar si un número se puede desarrollar de esta forma, bastará recorrer los valores de Y entre 1 y la raíz cuadrada de N/k. El valor que dé como resultado un cuadrado en el radicando será válido.

Por ejemplo, hemos afirmado arriba que 4516 se puede desarrollar como X²+9Y², con X>1 e Y>1. Según lo anterior, buscamos la raíz cuadrada de 4516/9, que resulta ser 22 (si tuviera decimales, truncaríamos). Por tanto habrá que ir probando desde Y=1 hasta Y=22 para ver qué valor da un cuadrado perfecto. Si se posee la función ESCUAD (puedes copiarla desde nuestra entrada http://hojaynumeros.blogspot.com.es/2015/02/numeros-especiales-que-son-un-producto.html ) para ver si un número es cuadrado perfecto, lo podemos organizar en una hoja de cálculo.

No hemos encontrado una solución hasta el valor 18, que junto a la raíz cuadrada de 1600 forma la expresión vista en el primer párrafo 40^2+9*18^2=4516

Función esforma(N;k)

Esta búsqueda que hemos efectuado se podría automatizar. Dado un número N y un coeficiente k podemos diseñar una función que devuelva TRUE si es posible la expresión N=X2+KY2 y FALSE en caso contrario. Existe una variante más útil, y es que si la expresión es posible, la devuelva en forma de String, y en caso contrario la frase “NO”.

En el Basic de VBA podría tener este código (añadimos la función ESCUAD por si no la has encontrado)

Public Function escuad(n) As Boolean
If n < 0 Then
escuad = False
Else
If n = Int(Sqr(n)) ^ 2 Then escuad = True Else escuad = False
End If
End Function

Public Function esforma(n, k) As String
Dim a, b, i
Dim es As Boolean
If k <= 0 Then esforma = "NO": Exit Function ‘Para k<=0 no hay solución
a = Int(Sqr(n / k)) ‘Tope de búsqueda
es = False
i = 1
While i <= a And Not es
b = n - k * i * i
If escuad(b) And b > 0 Then es = True: b = Sqr(b) ‘Se encuentra una solución
i = i + 1
Wend
If es Then
esforma = Str$(b) + "^2+" + Str$(k) + "*" + Str$(i - 1) + "^2" ‘Construcción del String
Else
esforma = "NO" ‘No es expresable
End If
End Function

Con esta función podemos construir un esquema para ver si un número admite la expresión N=X²+kY² con un valor de k dado. En esta imagen encontramos el desarrollo de 4516 con coeficiente 5:

Hay que advertir que la función ESFORMA sólo da la primera solución posible, sin descartar que existan otras.

En esta otra imagen se comprueba que el número 1298 no admite la expresión N=X²+7Y²,

Ordenando un poco los cálculos podemos reproducir la imagen con la que comenzamos la entrada (con formato de hoja de cálculo)

En el caso de k=3 nos devuelve una solución distinta, ya que hay más de una, y hemos prolongado la tabla para comprobar que para los valores k=11 y k=13 no existe solución.

Listado de números expresables

Con esta función y un bloque FOR-NEXT podemos encontrar la lista de los primeros números que se pueden expresar como N=X²+kY²  para un valor de k determinado, e incluso con un doble bucle, los que son expresables para varios valores. Así hemos construido la lista de los que son expresables para los valores k=1..10: 1009, 1129, 1201, 1801, 2521, 2689,…

Números que se puedan expresar con los valores de k=1..11 existen muchos menos. Los primeros son: 7561, 10756, 14116, 14281,…

El número 21961 satisface las expresiones para k=1..13 y los números 32356, 35044 y 35281 llegan al valor 14. Llegan hasta el 15 los números 32356, 35044 y 35281. Y así podríamos seguir, con valores cada vez más altos y escasos.

En estos listados están incluidos valores de k que pueden resultar redundantes. Por ejemplo, si un número es expresable como X^2+8Y^2, también lo es como X^2+2(2Y)^2. Así que si comprobamos la expresión X^2+8Y^2 lo estamos haciendo también con X^2+2Y^2. Como los cálculos no eran muy lentos, hemos preferido dejarlos. En otros trabajos similares se suelen estudiar tan solo los valores primos de k.

Aspecto modular

Si nos fijamos en los restos módulo k, es fácil ver que para que N=X²+kY²,  ha de ser N congruente con X²  módulo k, es decir, que N ha de ser un resto cuadrático módulo k. Si se dispone de un listado de esos restos, o un programa que los genere, podemos averiguar para qué polinomios del tipo dado es expresable un número. La hoja que hemos alojado en esta dirección

http://www.hojamat.es/sindecimales/congruencias/herramientas/hoja/congruencias2.xlsm

contiene restos cuadráticos para cada caso

Hemos adaptado provisionalmente esta herramienta para este caso particular. A cada número que probemos le calculamos el resto respecto a k y lo comparamos con la lista de cuadráticos. Esta prueba sólo la efectuaremos para valores de k que sean primos impares. Los demás casos se reducen a este. Vemos unos ejemplos mediante imágenes:

Aquí vemos que el resto 5 no figura en el listado de restos cuadráticos (columna de la izquierda), 1, 4, 6, 9, 10,…módulo 13, por lo que no es expresable para ese coeficiente.

El mismo resultado obtendríamos mediante ESFORMA(4516,13).

Por el contrario, el número 35281 sí es expresable mediante X²+13Y², , ya que su resto respecto al 13 es 12, y ese valor sí figura como resto cuadrático (el último de la primera columna). Lo dejamos aquí por si te apetece profundizar en la teoría de los restos cuadráticos.