lunes, 19 de diciembre de 2016

Numeradores de los números armónicos


Un número armónico Hn 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:




lunes, 5 de diciembre de 2016

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.