martes, 23 de abril de 2019

Simulaciones - La distribución binomial


Continuamos la serie de entradas sobre simulaciones de distribuciones estadísticas. Hoy llegamos a la distribución binomial. Como en toda la serie, y ese es uno de sus objetivos, usaremos la hoja de cálculo Simulador. Aunque la hayas descargado en otra ocasión, es conveniente que lo vuelvas a hacer, pues se le han añadido nuevas prestaciones, como la imitación de la máquina de Galton con probabilidad prefijada.

Se encuentra en las direcciones


Versión LibreOffice: 


y la puedes descargar para tu uso.

Distribución binomial

 Esta importante distribución se aplica a pruebas repetidas de la ley de Bernouilli


con las siguientes condiciones:

a) Se realizan experimentos repetidos del tipo Bernouilli, n en total.
b) La probabilidad p permanece constante en todos ellos
c) Cada experimento es independiente del resultado anterior.

Llamamos a n el número de intentos. Estamos interesados en estudiar el número de veces que aparece el suceso A (éxito). A su número de ocurrencias le llamaremos número de éxitos.

Por tanto la ley binomial se aplicará cuando repetimos un experimento cumpliendo las condiciones a), b) y c) establecidas y deseamos estudiar el número de éxitos que obtendremos. Son de este tipo las tiradas múltiples de monedas, de dados, de ruleta, ...

La probabilidad de obtener r éxitos en n intentos se demuestra que equivale a


En  ella el paréntesis es el número combinatorio n sobre r. Del hecho de que esta fórmula sea muy similar a la del Binomio de Newton proviene el nombre de binomial.

La media (esperanza matemática) de esta distribución viene dada por

y su varianza por


Consecuencia de esta es una fórmula que nos será muy útil, y es la de su desviación típica, que viene dada por

La distribución binomial de probabilidad p y número de intentos n se representa generalmente por B(n,p)

Puedes completar su estudio en


Uso de la hoja Simulador

 Comenzamos con un ejemplo, para estudiar la forma de plantear una simulación de este tipo. Habrá que concretar algunos parámetros, al igual que se hizo en simulaciones anteriores.

Tiramos 100 veces tres monedas. ¿En cuántas de ellas esperamos obtener tres caras?

La distribución binomial contiene decisiones automáticas en el Simulador, por lo que sólo hay que fijarle los siguientes parámetros:

·        Número filas y columnas: 100 filas y una columna (en este caso)
·        Tipo de simulación: Binomial (usa el desplegable)
·        Número de intentos: Lo escribes como parámetro A. En este caso son 3.
·        Probabilidad: Se escribe como parámetro B. Si deseas usar fracciones, escribe delante el signo =. Así, en el ejemplo escribiríamos =1/2.

El resto de parámetro lo rellena la hoja.

En la imagen observamos que ha fijado el número de intervalos en 4, para contar con el 0.


Los parámetros que no cambian se ignoran, como por ejemplo, 2 y 77.

Pulsamos sobre el botón “Simulación” y obtenemos los resultados de la simulación:

Número de veces en el que resultan tres caras:


Nos resultan 11 veces. Repetimos simulación y obtenemos 14, 16, 13, 14, 9, 17,…Esta variabilidad confirma lo peligroso de obtener conclusiones con solo 100 repeticiones. La simulación siempre es orientativa, pero se debe efectuar con más ensayos.

En este caso binomial también se pueden consultar los intervalos en la segunda hoja. Escribimos como extremos a y b el mismo número esperado de caras, 3, en la primera columna y nos devuelve el número de casos obtenido. En la imagen no aparece el 11, sino 16, porque corresponde a otra simulación:



Estudio con funciones de hoja de cálculo

Otra forma de responder a la cuestión es mediante las funciones estadísticas de Excel y Calc. Aquí estaría indicada DISTR.BINOM.N(X;N;P;Tipo) En ella escribimos los parámetros siguientes:

·        X es el punto en el que deseamos consultar la distribución, el resultado que esperamos. En este caso sería x=3, porque esperamos tres caras.
·        N es el número de intentos, que aquí también es 3.
·        P representa la probabilidad, que en monedas es 0,5
·        Tipo indica si la distribución es acumulada o no. Si su valor es 1, la distribución es acumulada y con 0 sin acumular. En el ejemplo deseamos no acumular. Sólo nos interesa el caso de 3 caras.

La escribimos en una celda y obtenemos: DISTR.BINOM.N(3;3;0,5;0)=0,125
Luego para 100 tiradas, el valor esperado sería 12,5. En la simulación obtuvimos 11, 14, 16, 13, 14, 9, 17,…Esto da idea de la variabilidad que presenta nuestra simulación.

Póquer con dados

En el juego familiar de póker con cinco dados, obtendremos póquer cuando cuatro de ellos marquen un mismo valor. Estudiaremos el caso en el que aparezcan en la primera tirada, sin comodín y sin acudir a otras tiradas, un póquer de reyes, por ejemplo. Se considera un suceso difícil. Lo simulamos:

·        Tipo: Binomial
·        Parámetro A: 5 dados
·        Parámetro B: Probabilidad 1/6 (lo escribimos como =1/6 y obtendremos 0,1666…)
·        Repeticiones: 500



Con ese número de repeticiones, no resultan en la simulación ni póquer de reyes ni repóquer:



Las frecuencias del 4 y el 5 son nulas, luego nuestra intuición de que son sucesos improbables no iba descaminada. Sobre el 500 escribimos un 2, para obligar a repetir la simulación dos veces.



Esta vez, con 1000 intentos sí se ha conseguido un póquer de reyes en una de las tiradas:


Aquí hemos trabajado con un solo valor en el resultado (reyes), pero el póquer puede salir con cualquier valor, lo que, al ser sucesos disjuntos, multiplicaría por 6 la probabilidad, pero es tan pequeña, que la simulación vale para darnos una idea de su dificultad.

Un ejemplo con intervalos

En un bombo de lotería se han introducido 100 bolas, numeradas del 00 al 99, con lo que todas las decenas figuran con 10 elementos. La probabilidad de obtener a ciegas una bola cuyo número comience por 3, será de 10/100=0,1
Imaginemos que 200 personas van sacando 10 bolas con reposición y contando las veces en las que obtienen un 3 en las decenas. ¿Cuántas de ellas esperaríamos que obtengan entre 3 y 10 éxitos?


El planteo sería:

·        Una repetición con 200 filas (200 personas)
·        Número de intentos: 10
·        Probabilidad: 1/10


Preparamos los intervalos de la segunda hoja para contar entre 3 y 10. Iniciamos la simulación y obtenemos:


Nos dan 9 casos, con una frecuencia relativa de 0,045. Fijamos ahora el experimento con 10 repeticiones (parte alta del cuadro de parámetros), pero en este caso el botón de intervalos no nos sirve, porque funciona sobre las filas de la simulación, pero no acumula. Mejor es leer las frecuencias de la tabla y sumar:



Entre 3 y 10 se han obtenido 124+17+3=144, que comparado con 2000 repeticiones nos da una frecuencia relativa de 144/2000=0,072, que parece más ajustada a la realidad que el 0,045 que obtuvimos con una sola tirada.

Podemos acudir a la función DISTR.BINOM.N(). En este caso restaremos la función acumulada en 10 de la acumulada en 2:

DISTR.BINOM.N(10;10;0,1;1)-DISTR.BINOM.N(2;10;0,1;1)=0,07019
Se observa, como era de esperar, que la frecuencia en la simulación repetida (0,072) se acercaba más al valor teórico (0,07019).

Aproximación a la normal

Se sabe que la distribución binomial se acerca a la normal bajo ciertas condiciones. Puedes repasar esta cuestión en 


Normalmente se piensa en la distribución normal cuando el valor de p es cercano a 1/2 y N tiende a infinito, aunque se suele obtener una buena aproximación práctica para N>30. Todo esto es un poco empírico, y se basa en el Teorema de Moivre, que puedes consultar en 


Ahí también se incluyen otros consejos para poder usar esta aproximación. Sólo queda indicar que la distribución normal límite poseería la misma media y desviación típica que la binomial.

Para ilustrar este ajuste, elegimos una binomial de 10 intentos y probabilidad 0,5, la simularemos en 100 filas con 10 repeticiones:



Con estas condiciones la media es 5 y la desviación típica 1,5811. Procedemos a la simulación y obtenemos, efectivamente, un resultado que se aproxima a la distribución normal:



Podemos estudiar el ajuste mediante la función =DISTR.NORM.N(x;M;D;0), que da la frecuencia para un valor dado en la distribución normal de media M y desviación típica D. Hemos incrementado el número de elementos a 1000, dejando en 10 el de repeticiones. A la tabla de frecuencias le hemos añadido esta función de Excel, con el siguiente resultado:



El buen ajuste se puede observar entre la columna de simulación y la de normal. La columna intermedia proviene de la función =DISTR.NORM.N(x;5;1.5811;0) y la última se ha creado multiplicando por N, que aquí es 10000.

Gráficamente se percibe mejor el buen ajuste:



La máquina de Galton como curiosidad

En la nueva versión del Simulador se ha añadido una hoja nueva con el experimento de Galton. Este modelo ya se ofrecía en nuestro curso de Estadística y en otros materiales, pero la nueva versión admite fijar una probabilidad que no tiene que ser necesariamente 1/2. De esta forma se puede visualizar una distribución binomial no centrada en la máquina. Lo puedes ver en esta imagen, correspondiente a p=0,7:



Observamos su sesgo hacia la derecha y su ajuste razonable a la distribución binomial teórica. En el gráfico se observa mejor el ajuste:



Con esto dejamos el tema de la simulación binomial. Se podía extender algo más, pero alargaría el texto.




jueves, 11 de abril de 2019

Suma de cuadrados de cifras (6) - Casos particulares


Llegamos hoy a la sexta entrada dedicada a la suma de los cuadrados de las cifras de un número. A esta seguirán otras, hasta completar una publicación más extensa.


Por tratar de casos particulares, es inevitable que esta entrada contenga mucha casuística, por lo que se desarrollarán más extensamente los primeros casos, reduciendo a meras referencias los últimos.

Suma cuadrada

En la tercera entrada de esta serie vimos los casos en los que la suma de cuadrados de las cifras es también cuadrada.


A continuación desarrollaremos algunos casos particulares de esta propiedad. 

Todas las sucesiones estudiadas serán partes de la general http://oeis.org/A175396

Cuadrados

Comenzamos con el caso más natural, y es que el número con suma cuadrada sea a su vez cuadrado. Estos números están publicados en http://oeis.org/A061090

Con nuestra función sumacifras(n;k) , que puedes consultar en la primera entrada de la serie, (http://hojaynumeros.blogspot.com/2018/11/suma-de-cuadrados-de-cifras-1-un.html) es fácil crear una búsqueda en la que tanto el número buscado como la suma de los cuadrados de sus cifras sean también cuadrados. Podría ser algo así:

For i = 1 to 100
a=i*i
b = sumacifras(a; 2)
If escuad(b) Then msgbox(a)
Next i

Con una rutina de búsqueda más elaborada hemos obtenido fácilmente el listado, que coincide con el publicado en A061090.



Listado publicado

1, 4, 9, 100, 400, 676, 841, 900, 1444, 4225, 10000, 24025, 40000, 42025, 42436, 43264, 66049, 67600, 84100, 90000, 109561, 119716, 144400, 155236, 239121, 244036, 248004, 252004, 335241, 355216, 362404, 373321, 422500, 643204, 664225

Como bien señala Charles R Greathouse IV, si un número pertenece a esta sucesión, también estará si lo multiplicamos por 100, ya que las cifras significativas no cambian.

El código PARI incluido en esta publicación resulta oscuro y no aprovecha mejoras habidas en este lenguaje, Proponemos otro:

for(p=1, 1000, b=p*p; a=norml2(digits(b)); if(issquare(a), print1(b,", ")))

Se comprende su funcionamiento con solo leerlo atentamente. Su resultado es



Primos

También están publicados los números primos con suma de cuadrados de cifras también cuadrada. Son estos:

2, 3, 5, 7, 43, 263, 269, 1153, 1531, 1933, 2063, 2069, 2287, 2609, 3319, 3391, 3511, 3931, 4003, 4441, 4801, 4889, 5113, 5399, 5939, 6029, 6067, 6203, 6469, 6607, 8849, 9133, 9539, 10111, 10177, 10513, 10531, 10771, 11149, 11213, 11273, 11321, 11491, 11503,…( http://oeis.org/A223035)

Los de una cifra pertenecen todos, como es evidente, y el siguiente, 43, corresponde a la terna pitagórica 42 + 32 = 52 = 25

Puedes reproducirlos con la función sumacifras(n;k):


También nos valdría este programa en PARI:

forprime(p=2, 10000, a=norml2(digits(p)); if(issquare(a), print1(p,", ")))

El uso de forprime simplifica mucho el código.

Triangulares

Para descubrir si un número N es triangular, basta ver que 8*N+1 sea cuadrado. Es una forma muy cómoda de descubrirlos. También podemos seguir la serie natural 1, 2, 3, 4,…y aplicarles la fórmula N(N+1)/2.

Aquí usaremos otro método, que es partir de 1 e ir añadiendo de forma recursiva, 2, 3, 4,…

1=1, 1+2=3, 3+3=6, 5+4=10, 10+5=15,…

Esta forma de generar números triangulares suele ser bastante rápida, y es la que usaremos en PARI:

p=1; q=2; while(p<=500000,  a=norml2(digits(p)); if(issquare(a), print1(p,", ")); p+=q; q+=1)

A la variable p se le va añadiendo q, que a su vez crece de 1 en 1. Así se generan los triangulares y después se suman los cuadrados de sus cifras y se investiga si la suma es cuadrada. El resultado es

0, 1, 3, 6, 10, 300, 2278, 6670, 10153, 27028, 28203, 33411, 35778, 40470, 44850, 58311, 62481, 66066, 105570, 113050, 143916, 144991, 151525, 155961, 181503, 198765, 199396, 204480, 213531, 220780, 239778, 241860, 248160, 271953, 279378, 333336, 348195, 381501, 392055, 431985, 461280, 479710,…(en la versión publicada se ha incorporado el cero, pues en OEIS se considera triangular, ya que 0=0*1/2)

Aquí también figuran, como ocurrió con los cuadrados, los triangulares de una cifra. Cualquier número de la sucesión será triangular y cuadrada la suma de sus cifras al cuadrado. Por ejemplo:

33411 es triangular, porque 33411=258*259/2, y la suma de los cuadrados de sus cifras es 9+9+16+1+1=36, un cuadrado.

Suma prima

 El caso general está publicado en http://oeis.org/A108662

11, 12, 14, 16, 21, 23, 25, 27, 32, 38, 41, 45, 49, 52, 54, 56, 58, 61, 65, 72, 78, 83, 85, 87, 94, 101, 102, 104, 106, 110, 111, 113, 119, 120, 126, 131, 133, 137, 140, 146, 159, 160, 162, 164, 166, 168, 173, 179, 186, 191, 195, 197, 199, 201, 203, 205, 207, 210,...

Como suele ocurrir en las sucesiones que estamos estudiando, si N está en la sucesión, todos los números cuyas cifras formen anagrama de las de N, también pertenecerán a la sucesión, Igual se puede afirmar de N*10^k. Así, si 38 pertenece a la sucesión, también estarán en ella 83, 380, 830, 3800, 8300,…Por tanto, esta sucesión es infinita.

Vemos algunos casos particulares.

Primos con suma prima

Están ya publicados en  http://oeis.org/A052034

11, 23, 41, 61, 83, 101, 113, 131, 137, 173, 179, 191, 197, 199, 223, 229, 311, 313, 317, 331, 337, 353, 373, 379, 397, 401, 409, 443, 449, 461, 463, 467, 601, 641, 643, 647, 661, 683, 719, 733, 739, 773, 797, 829, 863, 883, 911, 919, 937, 971, 977, 991, 997, 1013,...

Observa que los anagramáticos de un término, también pertenecen si son primos, como 113, 131 y 311.

Un ejemplo: 397 es primo, y 3^2+9^2+7^2=139, que es primo.

Podemos reproducir esta sucesión mediante el lenguaje PARI:

forprime(p=2, 10000, a=norml2(digits(p)); if(isprime(a), print1(p,", ")))


Triangulares con suma prima

Ya destacamos en la anterior cuestión que tenemos predilección por los números triangulares. Estos son los que presentan suma de cuadrados de cifras que es prima:

21, 45, 78, 120, 210, 276, 595, 780, 861, 1275, 1653, 2016, 2485, 2775, 3240, 3321, 3570, 3916, 4005, 4656, 6328, 7260, 7503, 9316, 9453, 9730, 10011, 11476, 12246, 12561, 13203, 15225, 15753, 16471, 16653, 17205, 17955, 18145, 19306, 19900, 20100, 20706, 21321, 21945, 22155, 22366, 22791, 23653, 23871, 24090, 24753,...

Vemos un ejemplo:

17205 es triangular, porque 17205=185*186/2, y 1^2+7^2+2^2+0^2+5^2=79 es primo.

El código PARI aprovecha la generación de triangulares como suma de los primeros números naturales:

(PARI) p=1; q=2; while(p<=25000,  a=norml2(digits(p)); if(isprime(a), print1(p,", ")); p+=q; q+=1)

Cuadrados con suma prima

Sorprendentemente, esta sucesión estaba inédita, a pesar de la sencillez de su definición. La inversa, primos con suma cuadrada sí está publicada en http://oeis.org/A223035

Los primeros números que cumplen esta condición son

16, 25, 49, 289, 324, 1369, 1521, 1600, 1936, 2209, 2304, 2500, 2601, 2809, 4900, 5929, 7569, 8649, 8836, 11025, 11881, 12321, 13225, 13689, 19881, 22201, 22801, 27556, 28900, 29929, 32400, 33489, 34596, 34969, 37636, 39601, 42849, 45369, 46656, 49284, 54756, 55225, 56169, 62001, 65536, 66564, 71289,...

for(p=1, 270, b=p*p; a=norml2(digits(b)); if(isprime(a), print1(b,", ")))

Suma triangular

Hay poco publicado, y tampoco parece que interese mucho, por lo que nos limitaremos a presentar los casos con algún comentario o código:

Los primeros números cuya suma de cuadrados es triangular son los siguientes

0,1, 6, 10, 13, 31, 36, 60, 63, 100, 103, 111, 112, 118, 121, 124, 130, 139, 142, 147, 174, 181, 193, 211, 214, 241, 244, 245, 254, 257, 275, 301, 306, 310, 319, 360, 391, 399, 412, 417, 421, 424, 425, 442, 452, 455, 458, 471, 485, 524, 527, 542, 545, 548, 554, 572, 584, 588, 600, 603, 630, 668, 669, 686, 696, 714, 725, 741,...

Se pueden generar mediante PARI:

for(p=0, 750, a=norml2(digits(p)); if(issquare(8*a+1), print1(p,", ")))

Dentro de ellos se pueden distinguir varios casos.

Triangulares con suma triangular

Los primeros son estos:

0, 1, 6, 10, 36, 630, 741, 1081, 2211, 7140, 10011, 10153, 13366, 15576, 17766, 23220, 24531, 25651, 28920, 33411, 42486, 47586, 52326, 59685, 61776, 69006, 112575, 113050, 121771, 125751, 128778, 129286, 146611, 156520, 163306, 165025, 167331, 185136, 202566, 207046, 211575, 219453, 221445, 222111, 231540, 235641, 237705, 243951,...

Ya están publicados en https://oeis.org/A094890

Por ejemplo, 741 es triangular (741=38*39/2) y la suma de cuadrados de cifras es 7^2+4^2+1^2=66, que también es triangular.

Puedes usar este código PARI

 p=0; q=1; while(p<=250000,  a=norml2(digits(p)); if(issquare(8*a+1), print1(p,", ")); p+=q; q+=1)

Primos con suma triangular

Los primeros son:

13, 31, 103, 139, 181, 193, 211, 241, 257, 421, 811, 1021, 1039, 1093, 1123, 1151, 1153, 1201, 1213, 1217, 1231, 1283, 1321, 1511, 1531, 1721, 1801, 1823, 1889, 2011, 2113, 2131, 2237, 2273, 2311, 2347, 2381, 2437, 2467, 2473, 2551, 2621, 2647, 2687, 2711, 2861, 3001, 3019, 3109, 3121, 3163, 3313, 3331, 3361, 3511, 3613, 3631, 3821,...

Su código:

forprime(p=2, 4000, a=norml2(digits(p)); if(issquare(8*a+1), print1(p,", ")))

Cuadrados con suma triangular

1, 36, 100, 121, 1024, 2401, 3136, 3600, 3844, 6724, 10000, 10201, 12100, 12769, 14161, 15376, 15625, 17689, 18769, 24649, 24964, 29584, 47089, 49729, 51529, 54289, 61504, 73441, 76176, 78961, 81796, 93636, 96721, 97344, 102400, 105625, 113569, 116964, 118336, 156025,...

Intenta construir un programa en PARI para ellos.

Aquí dejamos el recorrido de casos particulares en la suma de los cuadrados de las cifras. Como siempre, se ha preferido no abordar todos los casos a cansar a los seguidores con cálculos excesivamente parecidos.



martes, 2 de abril de 2019

Suma de potencias con bases en progresión aritmética



En el año 2017 publiqué en http://oeis.org/A292313 la sucesión de números equivalentes a una suma de tres cuadrados en progresión aritmética. Ahora he querido cambiar la cuestión, exigiendo a las bases que estén en progresión, pero no el resultado. Para un exponente general m, lo que buscaremos será aquellos números N que presenten la equivalencia con números enteros positivos:


En ella k representa la diferencia de la progresión y suponemos que es mayor que 0. Escrita así la igualdad podemos beneficiarnos de la simetría en los cálculos, como veremos más adelante.

Un ejemplo con cuadrados es 24318 = 872+902+932 = (90-3)2+902+(90+3)2

Comenzamos por ellos.

Suma de cuadrados en progresión aritmética

Encontrar los números que se puedan representar como suma de tres cuadrados con bases enteras positivas en progresión aritmética equivale a buscar aquellos que admitan la representación

N=(a-k)2 + a2 + (a+k)2

Esto parece fácil, pues basta crear una tabla de doble entrada para valores de a y k, siendo k<a.


El inconveniente radica en que están desordenados y que no se percibe a simple vista si existen duplicados (de hecho, 371 está repetido). Por eso, los valores 14, 29, 35, 50, 56, 66,… los encontraremos también con otras técnicas.

Como es costumbre en este blog, se caracterizarán estos números mediante una función. Para ello hay que considerar la expresión simplificada de la que los define:


En el problema propuesto, conocemos N, y el valor de a lo podemos ir cambiando entre 2 y la raíz cuadrada de N/3, que es una cota fácil de razonar. El valor de k depende de ellos, lo que nos evita un doble bucle de búsqueda, ya que es la raíz cuadrada de (N-3a2)/2. En la función que se usará nos preguntaremos si esa expresión es cuadrada, y en caso afirmativo, de ella obtendremos k y después N. Un ejemplo de función sería el siguiente:


Public Function basesenprog$(n)
Dim a, b, k, p, q, d
Dim s$

s$ = "" ‘Se usa un string para recoger las soluciones
k = Sqr(n / 3) ‘Cota para los valores de a
For a = 2 To k
b = (n - 3 * a ^ 2) / 2 ‘Se estudia el posible cuadrado de la diferencia de la p.a.
If escuad(b) Then d = Sqr(b) Else d = 0 ‘Si es cuadrado, se halla la diferencia d, y si no d=0
If d > 0 And d < a Then p = a - d: q = a + d: s$ = s$ + str$(p) + " " + str$(a) + " " + str$(q) + "&"
Next a
If s$ = "" Then s$ = "NO" ‘Si no hay solución, devuelve “NO”
basesenprog = s$
End Function

Con esta función y un bucle de búsqueda, obtenemos la lista ordenada de los números que obtuvimos con la tabla de doble entrada:



Cada uno viene acompañado de las tres bases en progresión aritmética. Por ejemplo, 140=22+62+102=4+36+100, con 6-2=10-6=4.

Con esta función también se detecta si un número presenta más de una solución. El primer número con esta propiedad es 371=12+92+172 y también 371=92+112+132.

Modificando la salida de la función, se pueden descubrir más números que admitan dos o más soluciones. Los primeros son estos:


El primero que presenta tres soluciones es 

2387=32+232+432=92+252+412=172+272+372

Las diferencias son 20, 16 y 10 respectivamente.

Si deseas las soluciones de los párrafos anteriores en forma de lista, puedes acudir al lenguaje PARI.

for(n=3,600,k=sqrt(n/3);a=2;v=0;while(a<=k&&v==0,b=(n-3*a^2)/2;if(b==truncate(b)&&issquare(b),d=sqrt(b);if(d>=1&&d<=a-1,v=1;print1(n,", ")));a+=1))

Obtendrás así el listado:

14, 29, 35, 50, 56, 66, 77, 83, 93, 107, 110, 116, 126, 140, 149, 155, 158, 165, 179, 194, 197, 200, 210, 219, 224, 242, 245, 251, 261, 264, 275, 290, 293, 302, 308, 315, 318, 332, 341, 350, 365, 371, 372, 381, 395, 398, 413, 428, 434, 435, 440, 450, 461, 462, 464, 482, 491, 504, 509, 515, 525, 530, 539, 557, 560, 563, 579, 590, 594, 596,…

Otra variante se inspira en la tabla de doble entrada y después elimina duplicados:

w=List();for(n=3,600,k=sqrt(n/3);for(a=2,k,for(c=1,a-1,v=(a-c)^2+a^2+(a+c)^2;if(v==n,listput(w,n)))));print(vecsort(Vec(w),,8))



Esta sucesión permanecía inédita y la hemos publicado en http://oeis.org/A306212


Suma de cubos en progresión aritmética

Acudiremos en este caso a las mismas técnicas que usamos con los cuadrados, pero de forma más breve. Buscamos números que presenten la equivalencia


En primer lugar, acudimos a una tabla de doble entrada para a y k:


Como se presenta el mismo inconveniente de los cuadrados, de presentación sin ordenar y sin depuración de repetidos, usaremos mejor el método de la función.

Para ello simplificaremos (a-k)3 + a3 + (a+k)3 = 3a3 + 6ak2 = N

Despejando k observamos que (N-3a3)/6a ha de ser cuadrado. Así que probaremos valores de a entre 2 y la raíz cúbica de N/3, tomando nota de cuando esa expresión sea cuadrada. Puede ser así:

Public Function basesenprog3$(n)
Dim a, b, k, p, q, d
Dim s$

s$ = ""
k = (n / 3) ^ (1 / 3) ‘Tope para la búsqueda
For a = 2 To k
b = (n - 3 * a ^ 3) / (6 * a) ‘Expresión que ha de ser cuadrada
If escuad(b) Then d = Sqr(b) Else d = 0
If d > 0 And d < a Then p = a - d: q = a + d: s$ = s$ + str$(p) + " " + str$(a) + " " + atr$(q) + "&"
‘Existe una solución
Next a
If s$ = "" Then s$ = "NO"
basesenprog3 = s$
End Function

Filtrando los primeros números naturales mediante esta función, obtenemos el listado:

36, 99, 153, 216, 288, 405, 408, 495, 645, 684, 792, 855, 972, 1071, 1197, 1224, 1407, 1548, 1584, 1701, 1728, 1968, 2079, 2241, 2304, 2403, 2541, 2673, 2736, 3051, 3060, 3240, 3264, 3537, 3540, 3888, 3960, 4059, 4131, 4257, 4500, 4587, 4833, 5049, 5160, 5256, 5472, 5643, 5832, 5940, 6336, 6369, 6669, 6840, 6903, 6984, 7227, 7293, 7776, 7839, 7860, 8217, 8316, 8541, 8568, 8712, 8988, 9339, 9399, 9576, 9792,...

Hemos publicado esta sucesión en http://oeis.org/A306213

El primer número que presenta dos soluciones es el 5643, ya que

5643=(9-8)3 + 93 + (9+8)3 = (11-5)3 + 113 + (11+5)3

Esta versión en PARI se basa en una tabla de doble entrada, eliminando repetidos:

w=List();for(n=3,10000,k=(n/3)^(1/3);for(a=2,k,for(c=1,a-1,v=(a-c)^3+a^3+(a+c)^3;if(v==n,listput(w,n)))));print(vecsort(Vec(w),,8))

Esta otra se basa en nuestra función Excel:

for(n=3,10000,k=(n/3)^(1/3);a=2;v=0;while(a<=k&&v==0,b=(n-3*a^3)/(6*a);if(b==truncate(b)&&issquare(b),d=sqrt(b),d=0);if(d>=1&&d<=a-1,v=1;print1(n,", "));a+=1))

Entre estas soluciones aparecen cuadrados destacables, como

62 = 13 + 23 + 33 y 482 = 43 + 83 + 123

Otra relación atractiva es 66 = (24-6)3 + 243 + (24+6)3 o 154 = (25-5)3 + 253 + (25+5)3

Podemos encontrar también muchos cubos. El primero es 63 = 3^3 + 4^3 + 5^3

Potencias cuartas

Actuamos de igual forma, iniciando una tabla de doble entrada y después una función. La tabla es la siguiente:


De ella se descubren los primeros elementos: 98, 353, 707, 962, 1568,…

Para construir la función necesitamos algo de Álgebra. Despejamos k, tal como procedimos en los casos anteriores:

(a-k)4 + a4 + (a+k)4 = 3a4 + 12a2k2 + 2k4 = N

De esta igualdad se deduce indirectamente que N>3a4 y que, por tanto, (N/3)1/4 es cota para a.
(Este resultado es general: Si (a-k)n + an + (a+k)n = N,
 la cota es (N/3)1/n. Se puede ver derivando la expresión para demostrar que es creciente respecto a k)

Seguimos:

2k4 + 12a2k2 - (N-3a4) = 0

Es una ecuación bicuadrada con solución positiva para k2:



Esto nos da una condición que debe cumplir esa expresión final, y es que sea cuadrado de un entero. Por eso, esta sería la función adecuada para descubrir los números que admiten estas sumas de potencias cuartas:

Public Function basesenprog4$(n)
Dim a, b, k, p, q, d
Dim s$

s$ = ""
k = Sqr(Sqr(n / 3)) ‘Cota superior N^(1/4)
For a = 2 To k
b = Sqr(2 * n + 30 * a ^ 4) / 2 - 3 * a ^ 2 ‘Valor de k^2
If escuad(b) Then d = Sqr(b) Else d = 0 ‘La variable d representa a k si existe solución
If d > 0 And d < a Then p = a - d: q = a + d: s$ = s$ + ajusta(p) + " " + ajusta(a) + " " + ajusta(q) + "&" ‘Formación de la solución en modo texto
Next a
If s$ = "" Then s$ = "NO"
basesenprog4 = s$
End Function

Con esta función se puede establecer la búsqueda, resultando, para las primeras soluciones:

98          1 2 3&
353        2 3 4&
707        1 3 5&
962        3 4 5&
1568      2 4 6&
2177      4 5 6&
2658      1 4 7&
3107      3 5 7&
4322      5 6 7&
4737      2 5 8&
5648      4 6 8&
7187      1 5 9&
7793      6 7 8&
7938      3 6 9&
9587      5 7 9&
11312    2 6 10&
12657    4 7 10&
13058    7 8 9&
15392    6 8 10&
15938    1 6 11&
17123    3 7 11&
19362    5 8 11&
20657    8 9 10&
23153    2 7 12&
23603    7 9 11&
25088    4 8 12&
28593    6 9 12&
30963    1 7 13&

La segunda columna presenta las tres bases en progresión aritmética.

Si ordenamos resulta el listado ordenado y sin repeticiones:

98, 353, 707, 962, 1568, 2177, 2658, 3107, 4322, 4737, 5648, 7187, 7793, 7938, 9587, 11312, 12657, 13058, 15392, 15938, 17123, 19362, 20657, 23153, 23603, 25088, 28593, 30963, 31202, 32738, 34832, 35747, 40962, 42528, 45233, 45377, 49712, 49763, 54722, 57153, 57267, 61250, 63938, 67667, 69152, 72113, 75792, 77922, 81473, 87713, 90083, 90368, 93602, 93827, 98787,...

Hemos publicado esta sucesión en http://oeis.org/A306214

Como en los casos anteriores,disponemos de dos versiones en PARI:

w=List();for(n=3,100000,k=(n/3)^(1/4);for(a=2,k,for(c=1,a-1,v=(a-c)^4+a^4+(a+c)^4;if(v==n,listput(w,n)))));print1(vecsort(Vec(w),,8))

y

for(n=3,100000,k=(n/3)^(1/4);a=2;v=0;while(a<=k&&v==0,d=sqrt(sqrt(2*n+30*a^4)/2-3*a^2);if(d==truncate(d)&&d>=1&&d<=a-1,v=1;print1(n,", "));a+=1))

Como curiosidad podemos destacar:

3923 = (56-28)4 + 564 + (56+28)4

Te invitamos a resolver el problema para potencias superiores. Ya tienes encauzados los cálculos y algoritmos.