Simulaciones - Distribución uniforme (1/2)

Iniciamos hoy una serie, que nos tomaremos con calma, sobre simulaciones elementales de variables aleatorias. Nos basaremos en las prácticas de nuestro curso de Estadística, (http://www.hojamat.es/estadistica/iniestad.htm) adaptándolas al formato de un blog. Usaremos nuestro Simulador implementado para hojas de cálculo, el cual puede sufrir cambios a lo largo de la serie, por lo que se aconseja su recarga en caso de duda.

Comenzaremos con la distribución uniforme. Si no tienes claro el concepto puedes acudir a la Teoría correspondiente

(http://www.hojamat.es/estadistica/tema6/teoria/teoria6.pdf).

Por ahora te basta con la idea de que representa experimentos aleatorios en los que todos los elementos presentan una misma probabilidad de ocurrir, como tiradas de dados o las loterías. Se suele distinguir entre distribución uniforme discreta, cuando sólo existe un número finito de posibilidades (dados o monedas), o continua, cuando pueden aparecer infinitos sucesos, o al menos, tantos, que sea preferible tratarlos como infinitos.


Distribución uniforme discreta

En ella se trabaja sobre un conjunto finito de n elementos con la hipótesis de que todos ellos poseen la misma probabilidad de aparecer, que será, por tanto, 1/n. Un ejemplo es el de una tirada de dados. La distribución uniforme más útil es aquella en la que el conjunto está formado por los números comprendidos entre a y b ambos inclusive. En los dados a=1 y b=6
.
Puedes consultar en cualquier manual los valores de los principales estadísticos de esta distribución.

Media:

Varianza:

En el caso de las tiradas de dados m=3,5 y var=35/12=2,92 y su desviación típica 1,7321.

Podemos comprobar estos valores mediante nuestro simulador, alojado en estas direcciones:

http://www.hojamat.es/estadistica/tema1/open/simulador.ods (versión LibreOffice Calc)

http://www.hojamat.es/estadistica/tema1/open/simulador.xlsm (versión Excel)

Contiene dos hojas, la de criterios y simulación y la de los estadísticos. La mejor forma de aprender su funcionamiento es proceder a la primera simulación.

Deseamos saber si con 1000 tiradas de dados su media y varianza se acercan suficientemente a la teoría. Para ello, en la primera página del simulador concretamos lo siguiente:

Cinco repeticiones de 200 filas y una columna, para que se acumulen 1000 tiradas, mínimo=1 y máximo=6. Como criterios, “Uniforme”, “Entero”, para que la distribución sea discreta, y “Máximo-Mínimo”. También es conveniente fijar, unas celdas más abajo, el número de intervalos en 6. El resto de parámetros se puede ignorar. Con ello, al dar al botón Simulador obtendrás los resultados en la siguiente hoja de Estadísticos. En nuestro caso serían:

Se ha obtenido una aproximación apreciable. La herramienta también nos proporciona la asimetría y la kurtosis, pero prescindimos de ellas en esta simulación. A la derecha puedes observar la tabla de frecuencias y el diagrama de barras, que presenta una uniformidad de alturas bastante aceptable para el número de simulaciones que hemos fijado:

Como ya sabrás, es probable que, si aumentamos el número de repeticiones, el ajuste mejore, pero no lo des por seguro, que sólo existe una probabilidad. Si aumentamos a 50 repeticiones obtenemos:

Ha mejorado bastante el ajuste. En general, las simulaciones comienzan a ser útiles si las repites miles de veces. Con unas pocas no son útiles.


Distribución uniforme continua

En esta modalidad los datos se distribuyen de forma continua (en la práctica, con todos los decimales que deseemos) entre dos extremos a y b. Prácticamente no hay ejemplos en la vida diaria de distribuciones uniformes, ya que son más frecuentes otras, como la normal.  Suelen aparecer en experimentos diseñados o en instrumentos creados por nosotros, como puede ser el movimiento de las manecillas de un reloj, que recorre de manera uniforme toda la circunferencia.

Un ejemplo de distribución uniforme continua es el experimento de calcular π mediante simulación (método de Montecarlo). Consiste en simular dos coordenadas X e Y de manera uniforme entre 0 y 1 y contar aquellos pares en los que X^2+Y^2<1. De esa forma, su frecuencia relativa deberá ser ?/4=0,7854 aproximadamente. En la imagen sería como contar todos los puntos que caen dentro de la zona sombreada.

Lo organizaremos así:

Planteamos los criterios contenidos en la imagen:

Tomamos 500 filas y dos columnas (que representarán X e Y). No planteamos repeticiones porque añadiremos una columna nueva a la simulación. Concretamos un mínimo de 0 y un máximo de 1. En los restantes criterios elegimos “Uniforme”, “Decimal” (por ser continua) y “Máximo-Mínimo”.

Con ello obtenemos dos columnas de 500 valores de X e Y. Ahora, en una columna paralela, por ejemplo comenzando en J5, escribimos =SI(G5^2+H5^2<1;1;0). Esto significa que obtendremos un 1 si el punto (X,Y) pertenece al sector circular sombreado, y 0 si está fuera. Extendemos esa fórmula hacia abajo hasta abarcar los 500 valores:

Ahora basta sumar esa columna nueva y deberemos obtener un valor próximo a 500π/4≈393. Esto no se suele obtener en una simulación con tan pocos datos. En nuestro caso se ha obtenido 389. Podemos repetir el trabajo (de forma manual) varias veces y encontrar la media de resultados. Aquí tienes un ejemplo, con 7 repeticiones o 3500 casos:

389, 408, 383, 389, 392, 386, 384

Sumo y obtengo 2731, divido entre 3500 (para encontrar la frecuencia relativa) y multiplico por 4 para aproximar a ?: 2731/3500*4=3,1211. No es una extraordinaria aproximación a π, pero resulta aceptable si tenemos en cuenta las herramientas utilizadas.

Análisis de intervalos

En una de las actualizaciones del Simulador hemos añadido una simulación entre intervalos. En el caso de la distribución uniforme nos servirá para analizar la igualdad aproximada de las frecuencias en intervalos de igual longitud.

En la segunda hoja Estadísticos figura una tabla y un botón para obtener frecuencias f entre dos extremos a y b . Estos extremos se consideran alcanzables, por lo que en las distribuciones discretas estarán incluidos. En la última fila se calcula la frecuencia relativa h. Sólo se puede usar para tres columnas o menos. Basta fijar los extremos en cada columna, pulsar el botón Intervalos y comparar resultados.

La imagen corresponde a una simulación uniforme continua entre 10 y 20, con tres columnas. En cada una de ellas hemos fijado extremos con una diferencia de 4, con lo que las tres frecuencias relativas se acercan a 0,4.

Escaladas de Conway

El año pasado, 2017, se notificó que una conjetura de Conway, conocida como “escalada a un primo” había resultado ser falsa. Tienes la noticia, comentarios y contraejemplos en estos dos blogs:

http://francis.naukas.com/2017/06/14/contraejemplo-a-una-conjetura-de-conway-sobre-los-primos/

https://www.gaussianos.com/la-conjetura-de-la-escalada-hasta-un-primo/

Si los has leído (y si no, te bastará con mi explicación) entenderás que el proceso que propone Conway es el de descomponer el número en sus factores primos agrupados con exponentes y ordenados en orden creciente, como 144=2⁴*3², y después escribir seguidos y mezclados bases y exponentes de las potencias resultantes (2432).

Si un número primo está elevado a la unidad, esta se ignora y no se añade al nuevo número.

Si llamamos f(n) a esa función obtendremos:

f(144)=2432

La idea de Conway es la de ir reiterando esta función hasta llegar a un número primo p, en el que es evidente que f(p)=p, dando fin al proceso.
En el caso del 144:

F(144)=2432, f(2432)=2719, que es primo y con él termina el proceso.

Los dos blogs citados te darán más detalles. Nuestro objetivo en esta entrada es conseguir el algoritmo con el Visual Basic de las hojas de cálculo. Es, por tanto un interés operativo más que matemático.

Función fconway(n)

A continuación se inserta el listado para hoja de cálculo de la función de Conway, pero antes hay que acudir a la función ajusta(n), que elimina del número n los espacios en blanco que Excel o Calc añaden a los números naturales. Su código es el siguiente:

Function ajusta$(a)
Dim d$

d$ = Str$(a) ‘Convierte el número en un string
While Left$(d$, 1) = " " ‘Mientras esté precedido de un espacio, este se elimina
d$ = Right$(d$, Len(d$) - 1) ‘Se corta el string desde su segundo carácter
Wend
ajusta$ = d$
End Function 

Una vez contamos con esta función, se tratará ahora de descomponer el número en factores y concatenar factores primos y exponentes, eliminando aquellos iguales a la unidad. Puede ser así:

Public Function fconway(n)

Dim primo(20), expo(20), numomega ‘Reservamos 20 memorias para primos y exponentes
Dim s$
Dim f, a, e, i


a = n ‘Recibimos n en la variable a
f = 2: i = 0: numomega = 0 ‘numomega es el número de primos
While f * f <= a ‘Vamos extrayendo primos hasta la raíz cuadrada de a
e = 0
While a / f = Int(a / f)
e = e + 1 ‘Se toma nota del exponente, que va creciendo
a = a / f ‘Se elimina el factor primo encontrado
Wend
If e > 0 Then ‘Se incorpora el nuevo primo con su exponente a las memorias
numomega = numomega + 1
primo(numomega) = f
expo(numomega) = e
End If
If f = 2 Then f = 3 Else f = f + 2 ‘Se avanza en posibles primos
Wend
If a > 1 Then ‘Se recoge el último primo
numomega = numomega + 1
primo(numomega) = a
expo(numomega) = 1
End If
s$ = "" ‘Construimos la concatenación de primos y exponentes mayores que 1
For i = 1 To numomega
s$ = s$ + ajusta(primo(i)) ‘Se incorpora el primo
If expo(i) > 1 Then s$ = s$ + ajusta(expo(i)) ‘Se añade el exponente si es mayor que 1
Next i
fconway = Val(s$) ‘Convertimos el string en número
End Function

Si escribimos un número cualquiera (en los blogs citados no se recomienda usar el 20) en una celda de Excel y tenemos implementada la función anterior (debes entrar en Programador – VisualBasic. Puedes consultar http://hojamat.es/guias/descubrir/htm/macros.pdf), podemos calcularla en la celda inferior, y después extenderla hacia abajo hasta que el resultado se repita o alcancemos una magnitud para la que Excel pasa a notación científica. Puedes ver algún ejemplo en la imagen:

En el primero, se alcanza el primo inmediatamente. En el segundo, al segundo intento, como ya vimos más arriba. El siguiente sobrepasa la capacidad de Excel, y el cuarto llega al primo en cinco pasos.

Podemos convertir este proceso en una función, pero si llegamos a los límites de Excel habrá que devolver un valor que lo indique, como podría ser el cero. Hemos diseñado esta:

Public Function finconway(n)
Dim a, b

a = 1: b = n ‘Usamos dos variables para cada iteración
While a <> b And a < 100000000000# And a > 0 ‘Límites de la iteración
a = b ‘Guardamos b en la variable a
b = fconway(b) ‘Avanzamos un paso en la iteración
If a >= 100000000000# Then a = 0 ‘Si el número es muy grande, lo hacemos cero
Wend
finconway = a
End Function

Con esta función se resuelve rápidamente el ascenso a primo. Aquí tienes los resultados desde 5 hasta 15, por ejemplo:

Esta función, aplicada al 542, devolvería un 0, ya que sobrepasa el límite de Excel para la escritura de todas las cifras. Este inconveniente se salva usando otro lenguaje. Hemos adaptado y completado la función en PARI incluida en la sucesión http://oeis.org/A080670 para definir la función finconway en ese lenguaje. Su código es:

fconway(n)=if(n>1, my(f=factor(n), s=""); for(i=1, #f~, s=Str(s, f[i, 1], if(f[i, 2]>1, f[i, 2], ""))); eval(s), 1)
finconway(n)=my(a=1,b=n);while(a<>b&&a>>0,a=b;b=fconway(b));a

Con ella es fácil reproducir los resultados de la tabla anterior:

Con el lenguaje PARI sí podemos encontrar el primo al que llegan las iteraciones con inicio en 542. Sería 131811420855589. En la imagen se puede observar cómo lo presenta PARI:

La conjetura

Si has leído las entradas de blog recomendadas más arriba, sabrás que existe un contraejemplo en el que la iteración no asciende a un primo, sino a un compuesto. Se trata del valor 13532385396179=13*532*3853*96179

Con nuestra función en PARI se detecta su invariancia en la iteración:

fconway(n)=if(n>1, my(f=factor(n), s=""); for(i=1, #f~, s=Str(s, f[i, 1], if(f[i, 2]>1, f[i, 2], ""))); eval(s), 1)
finconway(n)=my(a=1,b=n);while(a<>b&&a>>0,a=b;b=fconway(b));a
print(finconway(13532385396179))

En la imagen observamos que el resultado es el mismo valor:

Con esto terminamos, ya que el único objetivo era usar medios elementales para reproducir los ascensos a primo de Conway y el contraejemplo descubierto recientemente.

Acercamiento entre potencias

Se sabe que sólo existen dos potencias que se diferencien en una unidad, que son 8=2³ y 9=3². Sí pueden existir otros casos con diferencias mayores entre cuadrado y cubo. Vemos los primeros ejemplos:

Dif=2: 5²+2=3³
Dif:=3: 2²-3=1³
Dif=4: 2²+4=2³; 11²+4=5³
Dif=5 y 6: No hay ejemplos elementales
Dif=7: 1²+7=2³; 181²+7=32³
Dif=8: 3²-8=1³; 4²-8=2³; 312²-8=46³
Dif=9: 6²-9=3³; 15²-9=6³; 253²-9=40³
Dif=10: No hay ejemplos elementales.

Todos son casos particulares de la ecuación x²=y³+k, que necesita conocimientos profundas de Teoría de Números para su resolución. Por eso, aquí se usarán técnicas de búsqueda de soluciones en casos particulares.

Diferencias entre cuadrados y cubos

Podemos usar una función que mida el acercamiento a un cubo. Así, si se construye una lista de cuadrados, se podrán elegir aquellos que se diferencien de un cubo en un número dado. Así se han encontrado los ejemplos del párrafo anterior.

La función adecuada puede ser esta, que está diseñada para cualquier valor del exponente. Actúa sobre un número (que puede ser otra potencia), una diferencia dada y el exponente de la potencia. Al número le sumamos y restamos la diferencia dada, y con otra función, espotentipo, se averigua si es potencia perfecta o no.


Public function difepot(n,dife,ex1) as boolean
'para un número n ver si dista dife unidades de una potencia dada
‘Variables: n es el número, dife, la diferencia dada y ex1 el exponente de la potencia.

dim a,b
dim es as boolean

es=false ‘Declaramos que la función es falsa
a=n+dife ‘Buscamos la potencia “por arriba”
if espotentipo(a,ex1) then ‘Más abajo se explica esta función
es=true ‘Si es potencia, la diferencia es válida.
else
b=n-dife ‘Buscamos la potencia “por arriba y procedemos del mismo modo”
if b>0 then
If espotentipo(b,ex1) then es=true
end if
end if
difepot=es
end function

La función espotentipo busca si un número es potencia con un exponente dado. Actúa sobre el número y el exponente y devuelve VERDADERO o FALSO. Su código es el siguiente:

Public Function espotentipo(n, k) As Boolean
Dim m, i
Dim e As Boolean

m = Log(n) / k
m = Int(Exp(m)) ‘Mediante logaritmos, encuentra un posible valor para la base de la potencia
e = False
For i = m - 1 To m + 1 ‘Por si existen errores de redondeo, prueba con m, m+1 y m-1
If i ^ k = n Then e = True ‘Si con uno de ellos se reconstruye la potencia, es válido
Next i
espotentipo = e
End Function

Con estas dos funciones, bastará crear una lista de cuadrados y ver si alguno de ellos se diferencia de un cubo (mayor o menor que él) en una cantidad dada.  Por ejemplo, si construimos la lista de los veinte primeros cuadrados, veremos que al llegar al 14, su cuadrado se diferencia en 20 unidades del cubo de 6: 14²+20=6³.

Potencia más cercana

Podemos encontrar, para un número cualquiera, la potencia de cierto exponente que esté más cercana, y evaluar su diferencia. Para ello podemos usar esta otra función:

Public function difepot2(n,ex1)
'busca la potencia ex1 más cercana y devuelve la diferencia
dim p,q

p=int(n^(1/ex1))
dife=abs(n-p^ex1)
q=abs(n-(p+1)^ex1) ‘Busca la potencia si es menor o mayor. Ambas valen.
if q<dife then dife=q ‘Se queda con la diferencia menor.
difepot2=dife
end function

En esta tabla puedes observar las diferencias con el cubo más cercano. Se descubre que esos cubos son el 8 o el 27. Hasta el 17, el cubo más cercano es el 8, y a partir del 18, el 27.

Versión en PARI

Estas funciones de Excel no tienen mucha potencia de cálculo. Por eso, en algún momento usaremos sus versiones en PARI, cuyo código es el siguiente:

espotentipo(n,k)=local(m,i,es);m=log(n)/k;m=truncate(exp(m));es=0;for(i=m-1,m+1,if(i^k==n,es=1));es
difepot(n,ex1)=local(p,q,dife);p=truncate(n^(1/ex1));dife=abs(n-p^ex1);q=abs(n-(p+1)^ex1);if(q<=dife,dife=q);dife

Con el uso de ambas, hemos reproducido la tabla anterior entre el 10 y el 20:

Como ejemplo del uso de esta función, se inserta a continuación la tabla de diferencias entre cuadrado y cubo menores que 10 (y mayores que 0) para potencias interiores a 1000000:

Otras potencias

Diferencias mínimas entre cubo y cuarta potencia, dife<30:

Entre cuartas y quintas potencias, dife<50:

Así podemos seguir con otros casos. Terminamos con cubos y séptimas potencias:

Para terminar, agrupamos en una misma tabla diferencias menores que 10 en algunos casos. Hemos eliminado en la segunda potencia los valores 1 y 2, ya conocidos y triviales:

Permutación de cifras al sumar su producto

Nuestras colaboraciones en Twitter (@connumeros) nos sugieren otros estudios similares con más profundidad. En el día 9/5/18 publicamos que el número 9158, al sumarle el producto de sus cifras, se convertía en 9518, que está formado por las mismas cifras en distinto orden.

Investigando un poco vimos que Derek Orr ya había publicado casos similares en OEIS en 2014 (Ver http://oeis.org/A247888 y http://oeis.org/A243102)

Los primeros casos son:

239+2*3*9=293
326+3*2*6=362
364+3*6*4=436

Al estar publicados, aquí sólo estudiaremos las técnicas de hoja de cálculo que nos pueden permitir reproducir resultados.

Función PRODUCIFRAS

Es evidente que la primera herramienta que necesitamos es la función PRODUCIFRAS, que devuelve el producto de las cifras de un número. Ese producto no debe ser cero, pues en ese caso el resultado sería el número original y no se produciría la permutación “no trivial” de las cifras. Esto supone que sólo nos resultarán números sin la cifra cero en las soluciones.

Sugerimos este listado para PRODUCIFRAS en VBA de Excel:

Public Function producifras(n)
Dim h, i, s, m

h = n  ‘la variable h recoge el valor de n
s = 1 ,‘El producto comienza con un 1 en la variable s
While h > 0  ‘Mientras queden cifras…
i = Int(h / 10) ‘Se queda con todas las cifras menos la última
m = h - i * 10 ‘La variable m recoge la última cifra
h = i ‘La variable h tiene una cifra menos
s = s * m ‘La última cifra se incorpora al producto
Wend
producifras = s
End Function

Necesitamos además una función que nos indique si dos números presentan el mismo conjunto de cifras y con la misma frecuencia, como 2344 y 4342. El Basic de Excel no incluye la función VECSORT, propia de otros lenguajes, que nos permite ordenar un vector de números. Con ella bastaría ordenar las cifras de uno y otro número y comparar. Así se soluciona el problema en PARI, como puedes ver en la página enlazada, http://oeis.org/A243102.

Lo solucionaremos para VBA de Excel creando las matrices ca(10) y cb(10), que acumulen las frecuencias de las cifras de cada número a o b. Para gestionar bien el 0, anotaremos la cifra k en el elemento k+1. Así, un 7 se contaría en el elemento ca(8) o cb(8), y el conjunto {0,1,2…8,9) en los elementos ca(1), ca(2),… o bien cb(1), cb(2),…Finalmente, comparamos la matriz ca con la cb y si son idénticas es que los números contienen las mismas cifras con la misma frecuencia.

El listado sería este:

Public Function cifras_identicas(m, n) As Boolean

Dim i, j, s
Dim ci As Boolean
Dim nn$, mm$, c$
Dim ca(10), cb(10)

For i = 1 To 10: ca(i) = 0: cb(i) = 0: Next i ‘Vaciamos las matrices ca y cb

h = m  ‘Extraemos las cifras de m y las volcamos en la matriz ca
While h > 0
i = Int(h / 10)
s = h - i * 10
h = i
ca(s + 1) = ca(s + 1) + 1
Wend

h = n  ’Igualmente, las de n las volcamos en cb
While h > 0
i = Int(h / 10)
s = h - i * 10
h = i
cb(s + 1) = cb(s + 1) + 1
Wend
 ‘Ahora comparamos las dos matrices, y si existe una sola discrepancia, los números no tienen iguales cifras e iguales frecuencias

ci = True ‘Comenzamos suponiendo que las cifras son idénticas
For i = 1 To 9
If ca(i) <> cb(i) Then ci = False ‘Basta una discrepancia para que sea falso
Next i
cifras_identicas = ci
End Function

Ahora sólo nos queda recorrer un rango de números, por ejemplo del 1 al 3000 y quedarnos con aquellos en los que al sumarles el producto de las cifras (si es distinto de cero) se convierten en otros de las mismas cifras. Podría ser así:

For i=1 to 3000
a=producifras(i)
if a<>0 and cifras_identicas(i,i+a) then msgbox(i)
Next i

Aquí tienes el listado de los primeros números con esa propiedad, que coinciden con los publicados por Derek Orr

Como comprobación de la potencia del lenguaje PARI, se incluye el código usado por dicho autor, en el que usa la función vecsort:

for(n=1, 10^5, d=digits(n); p=prod(i=1, #d, d[i]); v=digits(n+p); if(v!=d, v=vecsort(v); d=vecsort(d); if(v==d, print1(n, ", "))))

Nuestro objetivo ha sido el uso de una hoja de cálculo, que es más lenta, pero encaja en los objetivos de este blog.

Suma de números oblongos consecutivos

En los cálculos sobre fechas que publicamos en Twitter (@connumeros) se dio la casualidad de que dos fechas muy cercanas se podían expresar ambas como extensas sumas de números oblongos consecutivos. Así:

Día 14/5/18

14518 es suma de siete productos de números consecutivos que a su vez son consecutivos (es decir, siete oblongos consecutivos):
14518=42×43+43×44+44×45+45×46+46×47+47×48+48×49

Día 19/5/18

Se nos vuelve a presentar la suma de productos de números consecutivos que a su vez son consecutivos (es decir, oblongos consecutivos), pero esta vez son nueve, nada menos:
19518=42×43+43×44+44×45+45×46+46×47+47×48+48×49+49×50+50×51

Esta cercanía nos animó a estudiar la propiedad con cierta extensión, para ver, como es costumbre en este blog, hasta dónde nos llevarían las exploraciones sobre un tema. Es fácil ver que, al ser los números oblongos doble de los triangulares, este problema esté muy relacionado con el de sumas de números triangulares consecutivos.

Sabemos que la suma de los n primeros números triangulares es n(n+1)(n+2)/6.
(Ver https://es.wikipedia.org/wiki/N%C3%BAmero_triangular#Suma_de_los_primeros_n%C3%BAmeros_triangulares)

De aquí se deduce que en el caso de los oblongos será n(n+1)(n+2)/3. Si en la suma no se comienza desde el primer oblongo (1*2) podremos restar esta expresión aplicada al último sumando con la correspondiente al anterior al primer sumando.

Esto nos lleva a una generación de soluciones en forma de tabla de doble entrada y otra mediante una función:

Soluciones obtenidas mediante una tabla de hoja de cálculo

En la siguiente tabla hemos situado valores de n(n+1)(n+2)/3 tanto en fila como en columna. Si restamos después unos de otros nos resultarán los números que se pueden formar mediante sumas de oblongos consecutivos. Aparecerán desordenados y repetidos. Nos valen también los de la segunda fila:

Para obtener una lista ordenada es preferible una búsqueda algorítmica mediante una función. Lo vemos:

Búsqueda mediante una función

La suma entre n(n+1) y (n+k-1)(n+k), k sumandos, sería así:

 n(n+1)+(n+1)(n+2)+…..(n+k-1)(n+k)=(n+k-1)(n+k)(n+k+1)/3-(n-1)n(n+1)/3

Hemos restado la suma de los n+k-1 primeros con la correspondiente a n-1.

 Desarrollamos y queda:

Por ejemplo, 23*24+24*25+25*26+26*27+27*28=S(23,27)=3260 equivale a
5*23^2+5^2*23+(5-1)*5*(5+1)/3=3260

Este desarrollo nos da una primera condición para que un número P pueda desarrollarse como suma de oblongos consecutivos. Proseguimos:

Despejamos n en función de P:

El discriminante

ha de ser cuadrado perfecto. Esta es la primera condición para que P sea suma de oblongos. También el resultado para n ha de ser entero positivo.

Lo vemos para P=14518 y k=7, según lo publicado en Twitter:

D=12*7^2-3*7^4+36*7*14518=3651921=1911^2, luego se cumple que D es cuadrado perfecto. Sustituimos en la expresión de n y queda:

 n=(-3*7^2+1911)/(6*7)=42, que coincide con el desarrollo publicado para 14518.

Por una casualidad, 19518 también inicia su suma de oblongos en 42:

P=18518, k=9

D=12*9^2-3*9^4+36*9*19518=6305121=2511^2
  n=(-3*9^2+2511)/(6*9)=42

Este criterio nos puede servir para ver qué números se pueden descomponer en suma de oblongos consecutivos.

Tomamos el número P y probamos los criterios anteriores para los valores de k entre 1 y el máximo valor que cumpla k*(k^2-1)<3P, que es una cota fácil de deducir de las igualdades anteriores. Para cada k exigiremos que el discriminante sea cuadrado perfecto y calcularemos n para ver si es entero positivo.

Así encontraremos todos los números que sean suma de oblongos consecutivos:

2, 6, 8, 12, 18, 20, 30, 32, 38, 40, 42, 50, 56, 62, 68, 70, 72, 90, 92, 98, 104, 110, 112, 128, 132, 148, 156, 160, 162, 166, 168, 170, 182, 200, 210, 218, 220, 232, 238, 240, 242, 260, 272, 288, 290, 306, 310, 322, 328, 330,…

Para encontrarlos se puede usar esta función, que devuelve una cadena vacía si el número no se puede descomponer, o una cadena que contiene todos los valores de n y k para los que es suma de oblongos. Podría ser esta:

Function essumaob(n) As String
Dim e$
Dim k, p, q

e$ = "" ‘Comenzamos con una cadena vacía
k = 1
While k * (k ^ 2 - 1) <= 3 * n ‘Recorremos los valores posibles de k
p = 12 * k ^ 2 - 3 * k ^ 4 + 36 * k * n ‘Se calcula el discriminante
If escuad(p) Then ‘Si es cuadrado perfecto, se sigue el proceso
q = (Sqr(p) - 3 * k ^ 2) / 6 / k ‘Encontramos el valor de n inicial en la suma
‘Si n es entero positivo, se toma nota en la cadena
If q = Abs(Int(q)) and q>0 Then e$ = e$ + Str$(q) + ", " + Str$(k) + " "
End If
k = k + 1
Wend
essumaob = e
End Function

Esta función sólo se debe aplicar a números pares, que son los únicos que pueden coincidir con una suma de oblongos. Aquí tienes algunos:

Vemos que 12 presenta n=3, k=1, ya que 12=3*4. Existen dos soluciones para 20: n=4, k=1, pues 20=4*5, y también n=1, k=3, es decir, 20=1*2+2*3+3*4

Con esta función podemos analizar cualquier otro número. Aquí tenemos la solución para los ejemplo 14518, 19518:

Se confirman las soluciones encontradas más arriba.

Si todas las soluciones las dividimos entre 2, resultarán los números que equivalen a suma de triangulares consecutivos. Los tienes publicados en http://oeis.org/A034706 con el añadido de un cero:

Relación con los números combinatorios

Es fácil ver que n(n+1)(n+2)/3 es el doble del número combinatorio C(n;3). Por tanto, las fórmulas que hemos usado en párrafos anteriores se pueden resumir en

Tomamos la diagonal del 3 en el triángulo de pascal y formamos todas las diferencias mutuas multiplicadas por 2:

Deberemos tomar los valores 2, 8, 20, 40, 70, 112, 168,…dobles de la cuarta diagonal y restar “todos con todos”. Llegaríamos a los mismos valores 2, 6, 8, 12, 18, 20, 30, 32, 38, 40, 42, 50,…

Relación con los números piramidales

La suma de triangulares consecutivos daba lugar a los números piramidales triangulares, u ortoedros (Ver nuestra entrada http://hojaynumeros.blogspot.com/2017/04/numeros-piramidales-2-tetraedros.html)

Como en este caso no los sumamos todos, sino sólo a partir de un índice, en lugar de piramidales triangulares serían “troncopiramidales”. Como también nuestros sumandos son el doble de un triangular, lo que obtenemos son “troncopiramidales de base oblonga”. En la imagen tienes representado así el número 38=2*3+3*4+4*5

¿Qué números están repetidos en el listado?

Algunos números suma de oblongos aparecen repetidos en los listados. Por ejemplo, 128=7*8+8*9=56+72 y 128=5*6+6*7+7*8=30+42+56

Para descubrir el número de soluciones que puede presentar un número P como suma de oblongos consecutivos basta modificar ESSUMAOB para que nos devuelva el número de soluciones en lugar de los valores de n y k. Sólo hay que cambiar alguna línea. Puede quedar así:

Function numsumaob(n)
Dim k, p, q, v


v = 0 ‘En lugar de una cadena devuelve el número de soluciones
k = 1
While k * (k ^ 2 - 1) <= 3 * n ‘Esta parte coincide con ESSUMAOB
p = 12 * k ^ 2 - 3 * k ^ 4 + 36 * k * n
If escuad(p) Then
q = (Sqr(p) - 3 * k ^ 2) / 6 / k
If q = Abs(Int(q)) And q > 0 Then v = v + 1 ‘Incrementa el contador de soluciones
End If
k = k + 1
Wend
numsumaob = v ‘Devuelve el número de soluciones
End Function

Con esta función podemos identificar rápidamente los primeros números que coinciden al menos dos veces con una suma de oblongos consecutivos:

Por ejemplo, 328 presenta los valores n=7, k=4 y n=2, k=8, es decir:

328=7*8+8*9+9*10+10*11
328=2*3+3*4+4*5+5*6+6*7+7*8+8*9+9*10

El primer número que admite tres desarrollos es el 4360, con los valores n=31, k=4; n=27, k=5;  n=9,  k=15

Con más soluciones no existen ejemplos menores que 500000. Se podrían buscar con instrumentos más rápidos que las hojas de cálculo.

Números que contienen las cifras de sus divisores

No es infrecuente que un número expresado en un sistema de numeración (nos limitaremos a la base decimal) contenga las cifras bien ordenadas de alguno o de todos sus divisores. Por ejemplo, 1734 contiene a sus divisores 17 y 34 como  subconjuntos ordenados de cifras. Existen publicadas muchas listas, de las que veremos algunas, especialmente con números primos. Para efectuar búsquedas sería conveniente disponer de alguna función que nos devolviera si un número contiene las cifras (siempre consecutivas) de algún divisor o, mejor aún, que indicara esos divisores.

Funciones necesarias

El criterio de si las cifras de un número están contenidas dentro de otro es fundamental en esta cuestión, pero existen dos inconvenientes:

• Los números deberán convertirse en cadenas de texto
• El formato de hoja de cálculo puede añadir espacios en blanco delante de los números positivos.

Estos dos problemas los resuelve la función AJUSTA$, como podemos observar en su listado:

Function ajusta$(a)
Dim d$

d$ = Str$(a) ‘Convierte el número en String
While Left$(d$, 1) = " "
d$ = Right$(d$, Len(d$) - 1) ‘Le suprime los espacios en blanco
Wend
ajusta$ = d$ ‘El resultado es una cadena de caracteres
End Function

Con esta función es fácil construir otra, DENTROCIFRAS, que indique si las cifras de un número contienen a las de otro  de forma consecutiva, es decir que el segundo sea una subcadena del primero.

Public Function dentrocifras(a, b) As Boolean 'Ve si las cifras de b están en a
Dim aa$, bb$
aa$ = ajusta(a) ‘Ajusta ambos números
bb$ = ajusta(b) ‘después, con función InStr averigua si está contenido
If InStr(aa$, bb$) > 0 Then dentrocifras = True Else dentrocifras = False
End Function

Es una función tipo boolean, que devuelve VERDADERO (True) o FALSO (False). Puedes probarla: DENTROCIFRAS(2431;43) debe dar VERDADERO y DENTROCIFRAS(818;88), FALSO, porque las cifras han de ser consecutivas.

Función CONTIENEDIV

Con esta función ya estamos preparados para saber si las cifras de un número contienen las de un divisor determinado, como ocurre con 725 y 25. Bastará hacernos las preguntas de si es divisor y después de si sus cifras están contenidas. Podemos organizarlo así:

Public Function contienediv(a, b) As Boolean
contienediv = a / b = a \ b And dentrocifras(a, b)
End Function

En primer lugar determina si es divisor, mediante a/b=a\b, que significa que al dividirlos resulta el mismo cociente que en la división entera \ (no aparece resto, por lo que se puede usar también la expresión a MOD b = 0) y después, con dentrocifras, si está contenido. Así, contienediv(912;12)=VERDADERO, porque 12 es divisor y sus cifras están contenidas en 912. Por el contrario, contienediv(324;24) da FALSO, porque 24 está contenido pero no es divisor.

Están publicadas las listas de aquellos números múltiplos de alguno de los veinte primeros que contienen sus cifras. Por ejemplo, la de múltiplos de 17 es

17, 170, 1173, 1700, 1717, 1734, 1751, 1768, 1785, 2176, 3179, 3417, 5117, 6171, 6817, 7174, 8177, 8517, 10217, 11713,…

Todos contienen un 17 y son múltiplos de él. Los tienes publicados en http://oeis.org/A121037

Con las funciones que hemos presentado puedes emprender otras búsquedas. Aquí tienes los primeros que contienen 23 como divisor entre sus cifras:

23, 230, 2231, 2300, 2323, 2346, 2369, 2392, 4232, 4623,…

Aunque es algo complicado, se adjunta a continuación el código PARI correspondiente. En la última línea el resultado sería cero, ya que hemos visto que es falso que 324 contenga a 24 como divisor. Hay que recordar que en PARI el valor VERDADERO se codifica como 1 y el FALSO como 0.

indigit(a, b)={ local(u,v,indi=0,la,lb,i,d,x);u=Vec(Str(a)); v=Vec(Str(b)); la=#u; lb=#v; i=1; while(i<=la-lb+1&&indi==0, d=0; for(x=1, lb, if(v[x]==u[i+x-1], d+=1)); indi=(d==lb) ; i+=1); indi}
indigitdiv(a,b)=(a%b==0)&&indigit(a,b)
print(indigitdiv(324,24))

Números que contienen todos sus factores primos

En la propuesta anterior fijábamos el valor del divisor que debería estar contenido en las cifras del número, pero ese dato no tenemos por qué saberlo.

Podíamos introducir una función que recorriera los factores primos del número y fuera probando uno a uno para ver si están contenidos o no como subcadena.
El problema está en esa extracción de factores primos. Si deseas profundizar lee la entrada de este blog

http://hojaynumeros.blogspot.com.es/2013/09/tus-funciones-disponibles-en-todas-las.html

En ella se explica la rutina sacaprimos y la función factores. Si las usamos, los factores primos  se guardarán en la matriz primo(n), y es esta la que vamos a usar ahora. Crearemos la función CONTIENEDIVPRIM(A,TIPO). Ella recorrerá los factores primos de A, comprobará con CONTIENEDIV si las cifras de estos están contenidas en el total y las irá acumulando en un string, al que añadirá la expresión del número de primos contenidos localizado. Así,
CONTIENEDIVPRIM (432;0)=2  3 y k= 2

Significa que contiene los factores 2 y 3 y que se han encontrado 2.

El parámetro TIPO nos sirve para, si vale 0, exigir que aparezcan todos los factores primos del número, como en el ejemplo anterior, y si es distinto de  0, fijará el número de factores que deseemos.

Public Function contienedivprim(a, tipo) As String
Dim n, i, k
Dim c$

c$ = ""
n = sacaprimos(a)
k = 0
For i = 1 To n
If contienediv(a, primo(i)) Then c$ = c$ + Str$(primo(i)) + " ": k = k + 1
Next i
c$ = c$ + "y k=" + Str$(k)
If tipo = 0 Then
If k = n Then contienedivprim = c Else contienedivprim = ""
Else
If k >= tipo Then contienedivprim = c Else contienedivprim = ""
End If
End Function

Con el TIPO=0 podemos confeccionar un listado de aquellos números compuestos que contienen a todos sus factores (el caso primo carece de interés):
25, 32, 125, 128, 135, 175, 243, 250, 256, 324, 375, 432, 512, 625, 735, 875, 1024,...

http://oeis.org/A050694

Números que contienen algunos factores primos

Con la función CONTIENEDIVPRIM, modificando el TIPO desde 1 hasta el tope que deseemos, se obtendrán los números que al menos coinciden con tantos factores primos como indique el valor de TIPO. Así, para dos factores obtenemos:

Figuran en la tabla los números con al menos dos coincidencias, los factores con los que coinciden y el número k de ellos. Si prolongamos la lista, aparecerá algún número que coincida con tres factores o más. El primero, que figurará en la siguiente lista, será el 735, que contiene los factores 3, 5 y 7.

Con tres coincidencias o más, haciendo TIPO=3, obtenemos

735, 1326, 1365, 1785, 2346, 2510, 2570, 2730, 3162, 3192, 3276, 3570, 3675, 3792


Estos son los primeros que contienen cuatro primos:

21372, que contiene a 2, 3, 13 y 137
37296, respecto a los divisores 2, 3, 7 y 37

Y con cinco

271362, con divisores 2, 3, 7, 13 y 71
527310, con 2, 3, 5, 7 y 31

Por último, si hacemos TIPO=0, aparecerán los compuestos que contienen a todos sus factores primos (está publicada en http://oeis.org/A050694)

25, 32, 125, 128, 135, 175, 243, 250, 256, 324, 375, 432, 512, 625, 735, 875, 1024, 1250, 1352, 1372, 1593, 1675, 1715, 1792, 2048, 2176, 2304, 2500, 2510, 2560, 2570, 2744, 3072, 3087, 3125, 3375, 3645, 3675, 3792, 4232, 4375, 5120, 5210, 5230, 5832…

Con esta entrada despedimos el curso 2017-18. Volveremos en septiembre si el verano se da bien. Os deseo un buen descanso.

Diferencias mínimas entre divisores

Si ordenamos en orden creciente (o decreciente) los divisores de un número compuesto, las diferencias entre dos consecutivos son siempre menores que las existentes entre dos más alejados, como es fácil de razonar, ya que todas las del segundo tipo son sumas de las consecutivas. La cuestión que nos planteamos es saber qué diferencia entre divisores es la mínima y cuántas veces se repite.

Luego particularizaremos a la diferencia 1, que daría lugar a divisores que son números consecutivos. Para la primera parte de este estudio no consideraremos el divisor 1, porque así obtendremos propiedades más interesantes.

Por ejemplo, si escribimos ordenados los divisores de 135 mayores que 1, obtendremos:

135 45 27 15 9 5 3

Vemos que la diferencia mínima entre los consecutivos es 2 (entre 3 y 5), y que es menor que todas las demás diferencias, sean o no entre consecutivos. Esta diferencia mínima se puede repetir en la sucesión de divisores. Tenemos, por ejemplo el caso del 540:

540 270 180 135 108 90 60 54 45 36 30 27 20 18 15 12 10 9 6 5 4 3 2

En él se repite cinco veces la diferencia mínima 1:

1=3-2=4-3=5-4=6-5=10-9

Es fácil ver que ocurre esto por ser múltiplo de un factorial.

Esta cuestión da lugar a varios cálculos distintos. Los vemos uno a uno.

Búsqueda de la diferencia mínima

La primera cuestión, dado un número concreto, es averiguar cuál es la diferencia mínima entre divisores.

La siguiente función en VBasic resuelve la cuestión utilizando tan sólo las funciones predefinidas de Excel o Calc. Recorre los divisores del número tomando nota de las diferencias entre consecutivos y almacenando la menor que se encuentre:

Public Function mindifdivi(n)
Dim i, p, p1, a1, a2, es

p = n - 1: p1 = p: a1 = 1: a2 = n: es = 0 ‘Variables iniciales
For i = 2 To n / 2 ‘Recorre desde 2 hasta la mitad del número
If n / i = n \ i Then  ‘Comprueba si es divisor
a1 = i: es = 1 ‘Guarda memoria del divisor
p1 = Abs(a2 - a1) ‘Calcula la diferencia con el divisor anterior
If p1 < p Then p = p1 ‘Si es más pequeña que la almacenada, la sustituye
a2 = a1
End If
Next i
If es = 0 Then mindifdivi = 0 Else mindifdivi = p ‘Si no es primo, devuelve la diferencia mínima
End function

Como la gran mayoría de números es múltiplo de los primeros primos, 2, 3, 5 o 7, abundan los resultados 0 (los primos), 1 y 2, como podemos ver en los valores de 1 a 20:

Las mayores diferencias suelen presentarse en los números semiprimos, como el 14, o potencias de primos, que sería el caso del 9. Es lógico que sea así. Por ejemplo el año 2018 es semiprimo con los factores muy alejados (2018=2*1009), por lo que su mínima diferencia es 1007.

Las mínimas diferencias se dan entre los múltiplos de 2 y 3 o de 3 y 5 (o 5 y 7, por ejemplo), que tienen el valor de 1 o 2. Esto puede ser engañoso. Las diferencias entre factores primos pueden no ser las mínimas. En el caso de 2450 sus factores son  7, 5 y 2, lo que haría esperar una diferencia mínima de 2, pero en realidad es 1, diferencia entre 49 y 50, consecutivos en la lista de divisores:

2450 1225 490 350 245 175 98 70 50 49 35 25 14 10 7 5 2 1

¿Cuántas veces se repite una diferencia mínima?

Una vez que hemos encontrado la diferencia mínima, con una función similar a la anterior podemos contar las veces en las que aparece. Recordamos de nuevo que se excluye el divisor 1.

Su código puede ser el siguiente:

Public Function nummindifdivi(n)
Dim i, p, p1, aa, a1

p = mindifdivi(n): p1 = 0 ‘Se calcula la diferencia mínima p y se pone el contador p1 a cero.
If p = 0 Then nummindifdivi = 0: Exit Function ‘Caso en el que n sea primo
For i = 2 To n / 2 ‘Recorremos los divisores propios
aa = n / i
If aa = Int(aa) Then ‘Si es un divisor, seguimos
a1 = aa + p ‘Incrementamos el divisor en la diferencia mínima
If n / a1 = Int(n / a1) Then p1 = p1 + 1 ‘Si resulta otro divisor, incrementamos el contador
End If
Next i
nummindifdivi = p1 ‘El resultado es el contador
End Function

A continuación reproducimos la tabla anterior adjuntando una columna con el número de ocurrencias de la diferencia mínima y la lista de divisores propios:

Así vemos, por ejemplo, que la diferencia mínima en el 12 es 1, y que se presenta dos veces, entre el 2 y el 3 y entre el 3 y el 4. En el caso del 14, la mínima es 5, y sólo aparece una vez, entre 2 y 7.

Con esta función podemos determinar los números con mayor número de ocurrencias de la diferencia mínima.

Ejemplos:

Números con más de cuatro diferencias mínimas

Los primeros son

Vemos que, por ejemplo, 420 presenta 7 diferencias mínimas iguales a 1:
2-3, 3-4, 4-5, 5-6, 6-7, 14-15 y 20-21.

Números con diferencia 2 y varias repeticiones

Los números 3465 y 4095 son los primeros que presentan cinco diferencias mínimas con valor 2. La razón es que todos sus factores primos son impares y esto facilita la repetición del 2. Así, 3465=3*3*5*7*11, y presenta diferencia 2 en
3-5, 5-7, 7-9, 9-11, 33-35

Las primeras diferencias son previsibles, restando los factores primos, pero otras es difícil encontrarlas por razonamiento, como 33-35.

Caso en el que diferencia mínima es 1 (consecutivos)

Este caso es el más interesante, por lo que merece una función especial para él. Aquí sí vamos a incluir el divisor 1, por la siguiente razón:

El número de pares de divisores consecutivos coincide con el de divisores oblongos, del tipo n(n+1). Es claro que todo número divisible por n y n+1 también lo es entre n(n+1), e igual ocurre a la inversa. Si no tenemos en cuenta el 1, perderíamos el oblongo 1*2=2.

Otro detalle que nos facilita la búsqueda es que los oblongos comienzan en 2 y luego sus diferencias van siendo 4, 6, 8, 12,…los pares consecutivos, lo que facilita el recorrido entre oblongos.

Según estas consideraciones, se puede diseñar la siguiente función numconsec que cuente los divisores consecutivos.

Public Function numconsec(n)
Dim i, m, p

m = 0: i = 2: p = 2 ‘Inicios de números oblongos
While i <= n
If n / i = n \ i Then m = m + 1 ‘Si es un divisor oblongo, se incrementa el contador
p = p + 2: i = i + p ‘Se genera el siguiente oblongo
Wend
numconsec = m
End Function


En http://oeis.org/A088723 están publicados los números que al menos poseen dos divisores consecutivos (sin contar el 1):

6, 12, 18, 20, 24, 30, 36, 40, 42, 48, 54, 56, 60, 66, 72, 78, 80, 84, 90, 96, 100, 102, 108, 110, 112, 114, 120, 126, 132, 138, 140, 144, 150, 156, 160, 162, 168, 174, 180, 182, 186, 192, 198, 200, 204, 210, 216, 220, 222, 224, 228, 234, 240, 246, 252, 258, 260…

Tomamos uno de ellos, por ejemplo el 222. Le buscamos todos los divisores y resultan ser:

222 111 74 37 6 3 2 1, con consecutivos 2-3

Con nuestra función resultarían dos, ya que contamos 1-2 y 2-3, para que así coincida con oblongos, que en este caso serían 2 y 6.

Es claro que todos los múltiplos de 6 pertenecen al listado, pero hay otros, como el 110, no múltiplo de 6 pero es oblongo y sus divisores consecutivos son 10 y 11.

Estos son los que presentan más de dos pares de consecutivos:

12, 24, 30, 36, 42, 48, 60, 72, 84, 90, 96, 105, 108, 120, 126, 130, 132, 144,…

Casi todos son múltiplos de 6 o incluso de 12.

Con más de tres consecutivos aparecen:

60, 72, 84, 90, 120, 132, 144, 156, 168, 180, 210, 216, 240,…

Es fácil ver que al aumentar el número de consecutivos iremos obteniendo entre ellos múltiplos de 60, pues nos garantizamos los consecutivos 1, 2, 3, 4,  5 y 6. Un caso especial vemos que es el 132=2*2*3*11, cuyos consecutivos son 1-2, 2-3, 3-4 y 11-12.

Al llegar aquí podemos pensar en los factoriales, que tienen garantizados n-1 consecutivos. Veamos si aparecen muchos más pares:

Hasta el 4 resultan los previstos. Después van apareciendo más casos, como en el 6!, que podemos recorrer 8: 1-2, 2-3, 3-4, 4-5, 5-6, 8-9, 9-10 y 15-16.

Números “record”

En el párrafo anterior podíamos sospechar que el incremento del número de pares de divisores consecutivos no seguía una sucesión creciente. Hemos programado la búsqueda de los números en los que el número de pares se incrementa en una unidad, y hemos obtenido estos:

Vemos que el orden natural se rompe entre 72 y 60, ya que este segundo presenta más pares que el 72. Como era de esperar, no figuran apenas factoriales, que serían los mejores candidatos para poseer más pares de números consecutivos.

El resto de soluciones, salvo el ya conocido desfase en la primera,  lo tienes en http://oeis.org/A088726

A088726: Smallest numbers having exactly n divisors d>1 such that also d+1 is a divisor.

1, 6, 12, 72, 60, 180, 360, 420, 840, 1260, 3780, 2520, 5040, 13860, 36960, 41580, 27720, 55440, 83160, 166320, 277200, 491400, 471240, 360360, 1113840, 720720, 1081080, 3341520, 2162160, 2827440, 5405400, 4324320, 12972960, 6126120,…

Se vuelve a romper el orden en otros casos, como 3780 y 2520.