lunes, 18 de mayo de 2015

La función de Smarandache y los números de Kempner (1)

La función de Smarandache se define, para un número natural n, como el menor entero tal que su factorial es divisible entre n. La designaremos como S(n). Por ejemplo, para n=12, el menor valor de k tal que k! sea divisible entre 12 es el 4, ya que 4!=24 es el menor factorial divisible entre 12. Lo expresaremos como S(12)=4. Es fácil que entiendas que S(6)=3 o que S(7)=7. Plantéate otros ejemplos.

Esta función fue estudiada por Lucas y Kempner antes de que Smarandache le asignara su propio nombre. Por eso, la sucesión de sus valores recibe el nombre de “números de Kempner”, y es esta:

1, 2, 3, 4, 5, 3, 7, 4, 6, 5, 11, 4, 13, 7, 5, 6, 17, 6, 19, 5, 7, 11, 23, 4, 10, 13, 9, 7, 29, 5, 31, 8, 11,… (http://oeis.org/A002034)

Aprovecha estos valores para comprobar la definición de la función en cada uno de ellos. Pronto descubrirás casos particulares, que podrás ampliar en la próxima entrada de este blog. Por ejemplo, adelantamos que el valor de S(p) para un número primo p es el mismo número: S(p)=p para p primo, o que S(n!)=n. Lo veremos más adelante.

En las dos primeras entradas de esta serie nos dedicaremos sólo a intentar construir algoritmos que reproduzcan los valores de la función. Comenzaremos por el más ingenuo y seguiremos con otros que contienen más artificio. Ante todo hay que notar que S(N)<=N, ya que todo número es divisor de su propio factorial. Esto nos beneficia, porque las búsquedas terminan en N.

Algoritmo “ingenuo”

Para encontrar el valor de S(n) bastará recorrer todos los factoriales desde 1 hasta n! y detenernos en el primero que sea múltiplo de n. El gran inconveniente de este procedimiento es que pronto se sobrepasará la capacidad de cálculo de la herramienta que usemos, especialmente si es una hoja de cálculo. Lo intentamos para Excel:

Public Function smar1(x)
Dim n, f
Dim seguir As Boolean

If x < 3 Then smar1 = x: Exit Function ‘Para x=1,2 S(x)=x
n = 1: f = 1  ‘Recorremos naturales desde 2 hasta x y f es su factorial
seguir = True  ‘variable para controlar el WHILE-WEND
While n <= x And seguir  ‘mientras no se llegue a n (cota natural) ni a la solución
n = n + 1: f = f * n  ‘se incrementa n y su factorial
If f = Int(f / x) * x Then seguir = False  ‘se ha llegado a un factorial divisible entre n
Wend
smar1 = n  ‘el valor de la función es el entero cuyo factorial es divisible
End Function

El algoritmo es sencillo, pero como se usan factoriales, aunque sea de forma indirecta, comete errores muy pronto (en Excel): De hecho, los valores para n=23 y n=29 ya son erróneos (destacados en rojo en la imagen):



Así no llegaremos muy lejos. Podíamos intentarlo en PARI a ver si funciona mejor (hemos eliminado los casos x=1,2):

smar1(x)={local(n=1,f=1,s=1);while(n<=x&&s,n+=1;f*=n;if(f(x==f\x,s=0));return(n)}
{for(i=3, 90, print1(smar1(i), ", "))}

Es el mismo algoritmo traducido a PARI. Para los 90 primeros casos funciona bien:

3, 4, 5, 3, 7, 4, 6, 5, 11, 4, 13, 7, 5, 6, 17, 6, 19, 5, 7, 11, 23, 4, 10, 13, 9, 7, 29, 5, 31, 8, 11, 17, 7, 6, 37, 19, 13, 5, 41, 7, 43, 11, 6, 23, 47, 6, 14, 10, 17, 13, 53, 9, 11, 7, 19, 29, 59, 5, 61, 31, 7, 8, 13, 11, 67, 17, 23, 7, 71, 6, 73, 37, 10, 19, 11, 13, 79, 6, 9, 41, 83, 7,…

Hemos probado el algoritmo para valores mayores y parece funcionar, pero nosotros deseamos mejorar el proceso para hoja de cálculo.

Algoritmo con el MCD

Este algoritmo lo hemos creado para el blog, pero es posible que ya esté publicado con anterioridad. Así que no se reclamará ninguna autoría.

La idea base es la de que el número dado va tomando factores de los elementos del conjunto 1, 2, 3, 4, 5,… hasta agotarlos todos. Por ejemplo, para encontrar S(18) necesitamos contar con los factores 2, 3, 3. Si tomamos los primeros números, el 18 podrá tomar de cada uno de ellos el MCD. La idea que lo simplifica todo es que una vez encontrado un factor, dividimos entre él para ir disminuyendo el valor primitivo (en este caso el 18).

Lo vemos en esta tabla:



Repetimos la idea: Elegimos un número, lo comparamos con 1, 2, 3, 4, 5,… dividiendo entre el MCD de ambos. Cuando el número llegue a 1, hemos terminado, y la solución será el último término de 1, 2, 3, 4, … probado. Lo explicamos de nuevo con n=250. Si el MCD es 1, lo saltamos:

MCD(250,2)=2, luego dividimos entre 2 y nos queda N=125
El MCD con 3 y 4 es 1, luego los saltamos
MCD(125,5)=5, dividimos N=25
Saltamos a 10: MCD(25,10)=5 y dividimos N=5
Saltamos al 15: MCD(5,15)=5, y al dividir obtenemos N=1. Ya hemos terminado: la solución es 15: S(250)=15

La ventaja de este método estriba en que no se multiplica, sino que se divide, con lo que los valores disminuyen hasta 1, evitando el desbordamiento en hoja de cálculo. Aunque se puede usar el cálculo manual con la misma hoja (sería muy pedagógico intentarlo en la Enseñanza), hemos implementado la función SMAR2, mucho más rápida, al disminuir los datos y poder eliminar una sentencia IF:

Public Function smar2(x)
Dim n, x1, m

If x < 3 Then smar2 = x: Exit Function
n = 1: x1 = x ‘la variable x1 recoge los cocientes entre x y el MCD con 1, 2, 3, 4,…
While x1 > 1 ‘se sigue mientras el cociente no llegue a 1
n = n + 1 ‘nuevo valor para los primeros números
m = mcd(n, x1) ‘encontramos el MCD
x1 = x1 / m ‘no hay que usar un IF porque es divisible con seguridad
Wend
smar2 = n
End Function

Con esta nueva implementación podemos calcular S(x) para valores mayores. Lo hemos intentado para comprobar que S(200001)=409 y se ha conseguido de forma prácticamente instantánea. Coincide con el resultado obtenido con PARI.

El problema está en que necesitas la función MCD. Aquí tienes una:

Public Function mcd(a1, b1)
Dim a, b, r

'Halla el MCD de a1 y b1

r = 1
a = a1
b = b1
If b = 0 Then b = 1
If a = 0 Then a = 1
While r <> 0
r = a - b * Int(a / b)
If r <> 0 Then a = b: b = r
Wend
mcd = b
End Function

Puedes estudiar esta versión muy sintética en PARI:

a(n)={local(m=1,x=n);while(x>1,m++;x=x/gcd(x,m));m}

En la siguiente entrada experimentaremos con algoritmos basados en la descomposición factorial, siguiendo las ideas de Kempner y basados en las propiedades de la función que estudiamos.

No hay comentarios: