Factores primos de la parte libre de cuadrados (y 4)

Función P(n), la omega de G(n) 

En los tuits citados en las anteriores tres entradas, de @republicofmath y @jamestanton en los que hemos basado los desarrollos de las mismas se introduce también la función P(n), que es el número de factores primos de G(n). Para entender mejor lo que sigue es conveniente releer esas  entradas. Son las tres anteriores a esta.

Las funciones P y Q aplicadas a G(n)

Si en lugar de multiplicar los factores primos de la parte libre de cuadrados del factorial, los contamos (función OMEGA), obtendremos la función P(n), que ya estudiamos en anterior entrada en su versión general.

Definiremos, pues, P(n)=omega(partelibre(n!). En código PARI se escribiría

P(n) = omega(core(n!))

Así se han encontrado los primeros valores de P(n): 0, 1, 2, 2, 3, 1, 2, 3, 3, 1, 2, 3, 4, 4, 4, 4, 5, 4, 5, 4, 6, 6, 7, 5, 5, 5, 6, 5, 6, 5, 6, 7, 9,… recogidos en http://oeis.org/A055460

Por ejemplo P(5)=3, porque 5!=120= 2³*3*5 contiene tres factores primos con exponente impar. Sin embargo P(7)=2 porque su factorial contiene primos elevados a un exponente par salvo el  5 y el 7.

En el Basic de las hojas de cálculo se evalúa esta función de forma idéntica a la de g(n), usando la fórmula de Polignac, solo que se cuentan factores en lugar de multiplicarlos:

Public Function p(n)
Dim i, s
s = 0
For i = 1 To n
If Not esnumpar(polignac(n, i)) Then s = s + 1
Next i
p = s
End Function

Así hemos reproducido sin dificultad los primeros valores:

Al igual que g(n), la función p(n) crecerá en los números primos y se mantendrá constante en los cuadrados. En los demás podrá aumentar o disminuir. Recorre la tabla para verificarlo.

Su crecimiento queda claro en el gráfico

Ajustes de P(n)

Esta función presenta una clara tendencia lineal. Si aumentamos el número de términos y añadimos una línea de tendencia obtenemos un ajuste bastante bueno a una recta de pendiente 0,1236 con R²=0,9853

¿Podríamos afinar más?

En los tuits citados se sugiere un crecimiento potencial suave. Proponen la fórmula potencial P(n) »  0.307*n^0.854. Hemos creado dos columnas paralelas, una con P(n) y otra con la formula 0.307*n^0.854.

La hemos prolongado a más de 1000 filas y hemos pedido una función que no se suele usar mucho en las hojas de cálculo: COEFICIENTE.R2. Esta función te devuelve el coeficiente de determinación, que evalúa la parte explicada de la función P(n) mediante esa aproximación. Resulta, tal como afirman los autores, R²=0,996998973, impresionante en su ajuste.

Volvemos al Solver

Como uno de los objetivos de este blog es el aprendizaje del uso de las hojas de cálculo, acudimos a la herramienta Solver para ver si Excel (en este caso) puede aproximar los valores 0.307 y 0.854 de la formula.  Al igual que operamos en la entrada anterior, asignamos dos celdas a estos parámetros, y los iniciamos, por ejemplo, en 0,3 y 0,8, a ver qué ocurre. A su derecha construimos la columna de las diferencias al cuadrado

Sumamos la columna de diferencias en la celda J1146. Con todo ello acudimos a Solver para ponerlo a prueba:

La solución que nos da no es la óptima

Para una herramienta no dedicada a usos científicos no está mal, pero vemos que no es fiable si se le exige mucho. Para comprobaciones serviría, pero sólo para eso.

Factores primos de la parte libre de cuadrados (3)

Ajustes de la función g(n) 

La entrada anterior la dedicamos a la parte libre de cuadrados de los factoriales. La llamamos g(n)=core(n!) e indicábamos que sus valores estaban contenidos en http://oeis.org/A055204. En dicha página señala Charles R Greathouse IV que log g(n) ~ n log 2. Comencemos por ahí:

Como en la entrada anterior se ofrecía una forma de evaluar g(n), podemos crear dos columnas paralelas, una con log(g(n)) y otra con n*log(2). El gráfico correspondiente a los primeros números nos indica que esta aproximación es siempre por exceso, y con un ajuste bastante alto: R²=0,99

La función log(g(n)) tiende a infinito con n de forma sensiblemente lineal

No he encontrado desarrollo teórico sobre esta aproximación, pero es algo que llama la atención. También se puede expresar como g(n) » 2ⁿ. También es sorprendente que g(n) se ajuste bastante bien al número de subconjuntos de un conjunto de n elementos.

James Tanton propone como aproximación inferior en media g(n) » 1,85n.  ¿Qué podríamos afirmar nosotros con una hoja de cálculo? No mucho, pero lo intentaremos:

Ajuste por mínimos cuadrados y Solver

Preparamos cuatro columnas de datos, en la imagen desde la I hasta la L

En la columna I escribimos los primeros números naturales, en la siguiente el logaritmo de G(N), y su aproximación mediante N*LOG(2) en la columna K. Observa que el 2 está escrito en la celda K1. A continuación calculamos en la última columna la diferencia de ambas expresiones elevada al cuadrado. Esta columna la sumamos al final, en la imagen en la celda L1057.

Ahora interviene Solver: le pedimos que elija el valor mínimo en la celda K1 (para sustituir el 2) que consiga minimizar la suma de diferencias al cuadrado contenida en L1057 con lo que habremos realizado un ajuste por mínimos cuadrados:

Nos da como mejor valor 1,94 aproximadamente, muy cercano al 2 de partida.

Este pequeño cambio hace que el ajuste en el gráfico se aprecie mejor:

El ajuste no está sesgado como en el caso del 2.

Esta técnica que acabamos de usar es sencilla, pero no muy usada. La ventaja que tiene es que tú puedes elegir la función de ajuste, que en las líneas de tendencia está obligada a ser lineal, exponencial o potencial. Recuerda los pasos:

* Situamos en dos columnas paralelas la función a estudiar y la que deseamos sirva de ajuste
* Si la función de ajuste depende de unos parámetros, tomamos nota de en qué celdas están situados.
* Creamos una tercera columna con las diferencias al cuadrado entre las dos primeras. La sumamos en una celda cuya referencia recordaremos.
* Acudimos a Solver y solicitamos minimizar la celda de la suma de diferencias al cuadrado a partir de las celdas que contienen los parámetros. Se usará un Solver no lineal.
* Si el ajuste es posible, aparecerán los nuevos valores de los parámetros.

Podemos, por pura curiosidad, intentar un ajuste lineal al LOG(G(N)) (neperiano). Resulta una coincidencia bastante fuerte, porque, además del sumando -7,2383, descubrimos que el coeficiente que da para la X es 0,6738, que es el logaritmo de 1,96, luego la expresión  log g(n) ~ n log 2 da un ajuste ligeramente superior.

G(n) y el primorial

Podemos conseguir otra aproximación comparando G(n) con los primoriales:

Recordamos que G(n) es la parte libre de cuadrados del factorial de n. Es un divisor del primorial n#, que es el producto de todos los números primos menores o iguales que n

 (ver http://hojaynumeros.blogspot.com.es/2012/02/el-primorial.html)

G(n) elige del primorial sólo los factores primos que presentan exponente impar en n. Basta recordar los esquemas que usamos cuando presentamos la función:

En el esquema, si multiplicamos los elementos de la primera columna nos resultará un primorial, y como en la segunda se marcan los que entran en G(n), si sólo multiplicamos los que figuran con 1, resultará, como hemos afirmado, que G(n) es un divisor de n#, y es claro que este, a su vez, es un divisor de n!. Esto nos lleva a unas acotaciones claras:

G(n) divide a n# y este a n!

Los cocientes tienen valores altos en el caso de los factoriales, como vemos en esta tabla.

Sin embargo, los correspondientes a N#/G(N) parecen más asequibles a nuestro estudio. Sabemos que los logaritmos de los primoriales se ajustan bien al valor de N. Veamos el ajuste del logaritmo del cociente N#/G(N)

Así que log(N#/G(N)) se acerca a 0,2928N y log(N#) a N. Se tendrá entonces:

Log(G(N)) » log(N#)-0,2928N » N-0,29N » 0,7072N >» Nlog(2)

Hemos llegado a un ajuste muy cercano al que obtuvimos anteriormente, pero por exceso. Lo más llamativo es que los distintos logaritmos presentan una tendencia lineal.

Factores primos de la parte libre de cuadrados (2)

La función G(n), cuadrados y factoriales

Hace unas semanas, navegando por Twitter encontré unos comentarios de Republic of Math (@republicofmath)  sobre resultados relativos a esta función. Me interesaron bastante y decidí estudiarla mediante hojas de cálculo, que es donde nos movemos en este blog. En la anterior entrada se incluyó un estudio sobre los factores primos de las partes cuadrada y libre como introducción al que se inicia hoy.

En dichos textos de Twiter se define g(n) como el mínimo número que multiplicado por el factorial de n lo convierte en un cuadrado. Ahora bien, según razonamos en la entrada

http://hojaynumeros.blogspot.com.es/2011/12/emparedado-de-cuadrados-2.html

esa función g(n) es, simplemente la parte libre de cuadrados del factorial de n. Si la parte libre la representamos como PL, la fórmula adecuada sería g(n)=PL(n!).

En lenguaje PARI esta función se representaría por core(n!), y así es como se ha engendrado la sucesión de valores de g(n) en http://oeis.org/A055204:

1, 2, 6, 6, 30, 5, 35, 70, 70, 7, 77, 231, 3003, 858, 1430, 1430, 24310, 12155, 230945, 46189, 969969, 176358, 4056234, 676039, 676039, 104006…

Desafortunadamente, en hoja de cálculo, si usamos la expresión equivalente con funciones nuestras: PARTELIBRE(FACT(N)), el cálculo se ralentiza hasta llegar a hacerse inútil. Para conseguir la tabla que sigue, hemos tenido que esperar varios minutos.

Para resolver esto, y entrando ya en un tema de algoritmos, podemos contar con una ayuda:

Fórmula de Polignac

Esta útil fórmula la estudiamos en la entrada http://hojaynumeros.blogspot.com.es/2009/02/formula-de-polignac.html, a la que remitimos para su definición y estudio.

La fórmula recorre todas las potencias de los factores primos menores que n y para cada una de ellas evalúa la parte entera del cociente de n entre cada una de las potencias.

El resultado equivale al exponente del factor primo dentro del factorial. Esto nos da una oportunidad para encontrar la parte libre de dicho factorial:

* Recorremos todos los números primos menores que n
* A cada uno le aplicamos la fórmula de Polignac
* Si su exponente es par, pertenece a la parte cuadrada del factorial, y no nos interesa.
* Si es impar, pertenecerá la parte libre, es decir, a g(n), tomándolo con exponente la unidad.

No es difícil programar como función estos cálculos. Este listado lo entenderás bien. Devuelve un cero si el número no es primo y su exponente dentro del factorial si lo es:

Public Function polignac(n, p)  n es el número y p el primo
Dim pol, pote

pol = 0 El valor se inicia en cero. Si no es primo p, se queda así
If esprimo(p) Then
pote = p Recorrerá las potencias de p menores que n
While pote <= n
pol = pol + Int(n / pote)  Sumando de la fórmula de Polignac
pote = pote * p  Se pasa a otra potencia del primo.
Wend
End If
polignac = pol
End Function

Puedes comprobar con esta fórmula la descomposición de 22! que publicamos en la entrada sobre Polignac:

Con esta función podemos encontrar los valores de la parte libre de cuadrados del factorial. En el ejemplo obtendríamos g(22)=2*3*7*13*17*19=176358.

Seguimos las operaciones sugeridas más arriba: recorrer los primos y tomar tan sólo aquellos que presenten un valor impar en la fórmula de Polignac:

Este segundo listado es más simple, y nos da el valor de g(n):

Public Function g(n)
Dim i, s
s = 1
For i = 1 To n Recorre los menores que n, sean o no primos
If Not esnumpar(polignac(n, i)) Then s = s * i Multiplica sólo los de exponente impar (si no es primo suma cero)
Next i
g = s
End Function

Ahora el proceso es mucho más rápido. Este listado se ha conseguido en segundos:

A primera vista hay algo que llama la atención, y es que la función no es creciente, aunque sí tenga esa tendencia a la larga, y que el valor para un cuadrado es idéntico al de su número anterior. Esto último se comprende porque un cuadrado no aporta nada a la parte libre de cuadrados del factorial. El que no sea creciente se explica porque la aportación del nuevo número puede ser de exponente impar que se acumule a otro impar ya existente y que entre ambos formen uno par, que por ser cuadrado se elimina. Pensemos en esto con más detenimiento.

Proceso recursivo

Si disponemos de la descomposición en factores primos de g(n) y la de n+1 entenderemos mejor por qué la función g(n) a veces crece otras decrece y en algunos casos queda igual. Usaremos el siguiente esquema:

En él hemos representado los exponentes (todos iguales a 1) de g(15) que es el producto 2*5*11*13. En la siguiente columna se han situado los exponentes de 16, que en este caso sólo figura el 4 correspondiente a 2⁴. Al pasar de 15! a 16!, el factor nuevo tiene exponente  par, luego el exponente del 2 no cambia, con lo que g(15)=g(16)=1430. Esto ocurriría en todos los cuadrados.

El valor de g(n) es igual al de g(n+1) cuando n+1 sea un cuadrado.

Si n+1 es un número primo, la situación es la opuesta, Observa el paso de 18 a 19:

g(18)=5*11*13*17. Como 19 es primo, no se combinará con los anteriores, y aparecerá como factor nuevo en g(19)= 5*11*13*17*19. Así ocurrirá con todos los números primos:

Si n+1 es primo, se cumplirá g(n+1)=(n+1)*g(n)

Recorre la tabla, detente en un número primo N  y observarás que g(N)=g(N-1)*N

En los demás casos, crece cuando el producto de los nuevos factores es superior al de los que se eliminan. Vemos dos ejemplos:

Paso del 19 al 20

Aquí los factores nuevos que aporta el 20 son 2 y 5. El 2 no cuenta porque está elevado al cuadrado, y se elimina. El 5 tampoco cuenta, porque con el 5 que ya está presente en g(19) forma un cuadrado y también se elimina. El resultado es que se pierde un 5 y la función disminuye.

Paso del 14 al 15

Aumenta, según el esquema. Estúdialo bien:

En la siguiente entrada estudiaremos los ajustes de esta función y veremos que existe una tendencia lineal para ella bastante clara.

Factores primos de la parte libre de cuadrados (1)

Esta entrada es la primera de cuatro consecutivas que dedicaremos al mismo tema. Introduciremos las funciones P(n), Q(n) y G(n), que nos servirán como excusa para profundizar en el conteo y producto de los factores primos de la parte libre de un número y de su factorial. Un apartado interesante es el de sus ajustes, pues revisaremos técnicas de hoja de cálculo,


Factores de la parte libre y de la parte cuadrada

Ya vimos en otra entrada

http://hojaynumeros.blogspot.com.es/2011/05/parte-cuadrada-y-parte-libre.html

que todos los números naturales poseen una parte cuadrada PC(N) y otra libre de cuadrados PL(N). La primera contiene como divisores todos los de N que son cuadrados. Si un factor primo está elevado a un exponente par pertenecerá a la parte cuadrada, pero si es impar, el par mayor contenido en él pasará a la parte cuadrada, y quedará en la parte libre el mismo factor elevado a la unidad.

Todos los factores primos de la parte libre de cuadrados están elevados a la unidad.

Puedes seguir la teoría en la citada entrada y también en nuestra publicación sobre funciones multiplicativas.

http://www.hojamat.es/publicaciones/multifun.pdf

También puede ser interesante contar los factores primos de la parte cuadrada, sin repetición.

Llamaremos  función Q(N) al resultado de contar esos primos. Así, por ejemplo, en el número 2520=2³×3²×5×7 tendríamos:

Parte cuadrada 2²×3²=36, Parte libre de cuadrados: 2×5×7=70, Q(2520)= 2, porque la parte cuadrada contiene dos primos distintos.

Los valores de esta función Q(N) los tienes en http://oeis.org/A056170

0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0,…

Puedes leer ahí algunos comentarios y desarrollos. El valor 0 aparece en los números libres de cuadrados. Verifícalo en la sucesión. Es sencillo de entender.

Presentarán valor 1 aquellos números cuya parte cuadrada posee un solo factor primo, como 4, 8, 9, 12, 16, 18, 20, 24, 25, 27, 28,…( http://oeis.org/A190641). El primer valor Q(n)=2 ocurre en el 36, y, en general, esta función cuenta los factores no unitarios de N.

Aprenderás bastante si ejecutas y analizas este código PARI que engendra esos valores. Ahí te lo dejamos. Recuerda que OMEGA cuenta los factores primos sin repetirlos y que CORE es la parte libre.

{for(i=2,36,print1(omega(i/core(i)),", "))}

Podíamos efectuar idéntica operación con la parte libre, contar sus factores primos. Llamaremos al resultado P(N). Sus valores son:

0, 1, 1, 0, 1, 2, 1, 1, 0, 2, 1, 1, 1, 2, 2, 0, 1, 1, 1, 1, 2, 2,…

y están contenidos en la sucesión OEIS http://oeis.org/A162642. En ella los valores 0 se corresponden con los cuadrados, porque en ellos la parte libre es 1 y no tiene factores primos.
Como en la anterior recomendamos la lectura del desarrollo de este enlace de OEIS y el que generes la sucesión mediante el código PARI

{for(i=1,36,print1(omega(core(i)),", "))}

Recuerda que core es la parte libre de cuadrados

Las funciones P(N) y Q(N) no actúan sobre conjuntos disjuntos de factores y pueden contar ambas el mismo factor, como ocurría con el 2 en el ejemplo de más arriba, el del 2520, que pertenecía a la parte cuadrada y también a la libre. Por tanto, la suma  P(N)+Q(N) es igual o mayor que OMEGA(N).

En la tabla siguiente podemos observar que en los números que contienen cubos, como 8, 24 y 27, presentan esa desigualdad P(N)+Q(N) > OMEGA(N).

Puedes reflexionar sobre qué números presentan esa desigualdad además de los cubos.

P y Q como funciones aditivas

En Teoría de Números una función f(n) se llama aditiva cuando se cumple

F(ab) = f(a) + f(b) siempre que a y b sean coprimos

En efecto, si a y b son coprimos, tanto su parte cuadrada como su parte libre poseerán factores primos diferentes en ambos números. Por tanto,  P y Q aportarán al producto factores que no pertenecerán a la otra función. En ese producto figurarán los que aporta cada uno sin coincidencias, por lo que sus cuentas se sumarán. Lo puedes verificar en la tabla de más arriba, por ejemplo:

P(2)=1, P(9)=0 y P(2*9)=P(18)=1=P(2)+P(9)

Prueba también con otros pares (coprimos) y con Q(n), y comprobarás la aditividad.

Al igual que las funciones multiplicativas, las aditivas se definen sólo para potencias de primos. En este caso la definición adecuada de Q(p^m) sería

Q(p^m)=0 si m=1, y Q(p^m)=1 en los demás casos. 

Lo puedes expresar también como p^sg(m-1), donde sg es la función signo, que vale 1 en los positivos y 0 en el cero.

Para la función P tendríamos la situación opuesta:

P(p^m)=1 si m es impar, y P(p^m)=0 si m es par. 

También se puede resumir como P(p^m)=(m mod 2)

La falta de simetría en las definiciones viene dada por el hecho de que si un primo está elevado a exponente 2 o mayor, se cuenta en Q y no en P, tanto si es par o impar.

Sobre estas funciones volveremos en las siguientes entradas, cuando estudiemos los factores primos de las parte libre de cuadrados en el caso del factorial.

Conjeturas con hoja de cálculo: Primos tipo n^2+1

Es uno de los problemas de Landau, y en el momento de redactar este texto sigue sin conocerse si es verdadera o no la siguiente conjetura:

Existen infinitos primos de la forma n²+1

Hardy y Littlewood supusieron que la conjetura era verdadera, y aproximaron el número de tales números primos menores que n, P (n), asintóticamente a

Con C una constante adecuada.

El listado de los primeros primos de este tipo lo puedes consultar en http://oeis.org/A002496

2, 5, 17, 37, 101, 197, 257, 401, 577, 677, 1297, 1601, 2917, 3137, 4357, 5477, 7057, 8101, 8837, 12101, 13457, 14401,…

Si la conjetura es cierta, esta sucesión deberá poseer infinitos términos.

¿Qué estudios podríamos abordar sobre este tema con una hoja de cálculo?

El primer objetivo razonable es el de comprobar que, dado un número cualquiera, existe un número primo del tipo n²+1 que es mayor que él.

 Usaremos la herramienta de hoja de cálculo conjeturas, alojada en

http://www.hojamat.es/sindecimales/divisibilidad/herramientas/herrdiv.htm#global

Para encontrar ese primo mayor que el dado, reiteraremos el uso de la función PRIMPROX hasta que encontremos un número primo p tal que p-1 sea un cuadrado.

(A) Planteamiento manual

Basta estudiar este esquema brevemente para descubrir su funcionamiento:

El primer número primo de la lista es el PRIMPROX(N), en la imagen 1511. Los siguientes se obtienen como los próximos primos del de la fila superior. Esta lista se puede extender hacia abajo todo lo que se desee.

En la segunda columna hemos usado una fórmula del tipo

=SI(ESCUAD(C9-1);C9;""),

es decir, si C9 u otro primo de la lista cumple que al restarle la unidad se convierte en un cuadrado, lo escribimos, y , si no, dejamos la celda en blanco. Así descubrimos que el primer primo de este tipo es 1601. Si la conjetura es cierta, siempre llegaremos a un número de ese tipo.

Este método puede necesitar muchas filas hasta dar con el primo esperado. Por eso, se puede plantear como una función:

(B) Estudio mediante una función

Si suponemos cierta la conjetura, para cada número existirá un primo mayor que él con la forma n²+1. Entonces lo podemos plantear como una función. Su listado lo entenderás fácilmente:

Public Function proxn2mas1(n)
Dim p
p = primprox(n)
While Not escuad(p - 1)
p = primprox(p)
Wend
proxn2mas1 = p
End Function

De esta forma, la búsqueda manual que emprendimos en el caso anterior la podemos reducir al planteamiento de esta función:

Se comprende que para números grandes esta función tardará algo en calcularse. Lo hemos intentado con 10^7:

Tarda unos segundos, aunque no es un retraso desesperante. Hemos añadido la raíz cuadrada del primo menos uno, 3174.

Lista de primos de este tipo

Con esa función proxn2mas1 podemos reproducir toda la lista de OEIS. Basta escribir un 2, debajo de él proxn2mas1(2) y nos resultará un 5. Le aplicamos de nuevo proxn2mas1 y obtendremos el 17, y así seguimos hasta donde deseemos.

Si en lugar de comenzar con el 2 inicias con un número cualquiera, se escribirá la continuación de la lista, salvo quizás el primero, que no tiene que ser primo de ese tipo. Aquí tienes los siguientes a 10000:

Aproximación asintótica

Para comprobar la aproximación de Hardy y Littlewood necesitamos contar los primos de este tipo anteriores a N. Algo parecido a la función PRIME(N), pero quedándonos sólo con los primos de forma n²+1

Entenderás a la primera esta definición:

Public Function ppn2mas1(n)
Dim pp, i
' para valores de n superiores s 2

i = 2: pp = 0
While i <= n
pp = pp + 1
i = proxn2mas1(i)
Wend
ppn2mas1 = pp
End Function

Esta función cuenta los primos del tipo n²+1 inferiores o iguales a N. Como nos interesan valores grandes por cuestiones asintóticas, suponemos, para simplificar la programación, que N es mayor que 2. Observa esta tabla en la que se percibe que tratamos con una función escalonada, y que los cambios ocurren en 5 y 17, primos del tipo estudiado.

Para comprobar la aproximación asintótica y evaluar la constante C crearemos una tabla a partir del 10 en progresión geométrica hasta llegar a 10000000:

Hemos evaluado la constante C como cociente entre la función de distribución de los primos de tipo  n²+1 inferiores a N y la aproximación RAIZ(N)/LN(N). No usamos una herramienta adecuada, pero se ve que los valores de C presentan una cierta convergencia.

Variante de la conjetura

La más sencilla es la que busca primos de la forma  n²+a. Podemos crear una función similar a la que hemos usado, pero añadiendo un parámetro A

Public Function proxn2masa(n,a)
Dim p

p = primprox(n)
While Not escuad(p - a)
p = primprox(p)
Wend
proxn2masa = p
End Function

Con esta función se puede comprobar que,  dado cualquier valor de N (primo o no) y elegida una constante A, existe un número primo del tipo  N²+A superior a N

Observa un posible esquema de búsqueda:

En él elegimos N y A y se nos devuelve el primo adecuado y la raíz cuadrada de N-A

Sucesión de Perrin

En la pasada temporada dedicamos varias entradas a las sucesiones definidas mediante una recurrencia de segundo orden. Ahora comenzaremos una serie sobre las de tercer orden. Entre ellas son muy populares las de Perrin y Padovan. Como en las anteriores, nuestro planteamiento no será teórico, pues ya existe mucho publicado sobre ellas. El objetivo será crear esquemas y cálculos que faciliten  la comprensión de sus propiedades.

Sucesión de Perrin

La teoría fundamental sobre esta serie la puedes consultar en

http://mathworld.wolfram.com/PerrinSequence.html
http://en.wikipedia.org/wiki/Perrin_number

Aquí la describiremos con la ayuda de la herramienta que hemos ofrecido en entradas anteriores, alojada en

http://hojamat.es/sindecimales/aritmetica/herramientas/herrarit.htm#recurre2

Definición

Esta sucesión es recursiva de tercer orden homogénea, por lo que necesita tres valores iniciales y que X(n) dependa de los tres valores anteriores X(n-1), X(n-2) y X(n-3) mediante la relación

X_n= A*X_n-1+B*X_n-2+C*X_n-3

En este caso particular sólo depende de los dos últimos, y no de X(n-1). Concretando:

Condiciones iniciales: X₀=3   X₁=0  X₂=2 Ecuación de recurrencia: X_n= X_n-2+X_n-3

Es como una sucesión del tipo Fibonacci pero “con retraso”, pues los que se suman no son los dos anteriores, sino los que están un paso más atrás.

En nuestra hoja de cálculo se define así (segunda hoja del libro):

El primer coeficiente es nulo, que es lo que produce el “retraso”, y debajo tienes los tres valores iniciales.

La sucesión resultante la vemos pulsando el botón correspondiente:

Esta popular sucesión la tienes disponible en http://oeis.org/A001608, donde les llaman números skiponacci, quizás por los saltos o retardos que presentan: 3, 0, 2, 3, 2, 5, 5, 7, 10, 12, 17, 22, 29, 39, 51, 68, 90, 119, 158, 209, 277, 367, 486, 644, 853,…

Ecuación característica

La ecuación característica correspondiente será X³-x-1=0. Con el botón Resolver de esa hoja obtienes las tres soluciones de la ecuación, una real y dos complejas

Coinciden con las soluciones que da WxMaxima

La solución real 1,32471…(aquí sólo aproximada) es el número plástico y, cuyo nombre se eligió como afín al número de oro o el de plata. En estas páginas puedes estudiarlo más a fondo:

http://es.wikipedia.org/wiki/N%C3%BAmero_pl%C3%A1stico
http://revistasuma.es/IMG/pdf/57/055-064.pdf http://cscmates.blogspot.com.es/2010/11/el-numero-de-plastico.html

Recordemos que, como en sucesiones anteriores, todo número de Perrin es combinación lineal de las potencias de las tres soluciones de la ecuación característica, pero las dos complejas tienen módulo menor que la unidad, por lo que sus potencias tenderán a cero en valor absoluto. Por tanto, X(n) se acercará asintóticamente a yⁿ

Se puede construir una tabla doble en la que se observe este acercamiento:

A partir de un cierto orden basta redondear la potencia para obtener el número de Perrin  correspondiente. Lo puedes comprobar en las últimas filas de la tabla.

Función generatriz

Usando procedimientos similares a los que explicamos para las recurrentes de segundo orden, se puede demostrar que la función generatriz es

Puedes comprobar que esta es la F.G. adecuada efectuando este desarrollo en PARI

write("sucesion.txt",taylor((3-x^2)/(1-x^2-x^3),x,20))

Te escribirá en un archivo sucesión.txt su desarrollo, y aparecerán como coeficientes los términos de la sucesión de Perrin:

3 + 2*x^2 + 3*x^3 + 2*x^4 + 5*x^5 + 5*x^6 + 7*x^7 + 10*x^8 + 12*x^9 + 17*x^10 + 22*x^11 + 29*x^12 + 39*x^13 + 51*x^14 + 68*x^15 + 90*x^16 + 119*x^17 + 158*x^18 + 209*x^19 + O(x^20)

Sucesión de Perrin y números primos

La propiedad más conocida de estos números es que si p es primo, p divide a X(p). Por ejemplo, X(11)=22, que es múltiplo de 11. Podemos construir una tabla en la que dividamos X(n) entre n y los cocientes enteros se corresponderán con los números primos:

A pesar de su carácter algo extraño, la propiedad ha sido demostrada para todos los números primos. La contraria no es cierta. X(n) puede ser múltiplo de n sin que este sea primo. A estos términos se les suele llamar pseudoprimos de Perrin (http://oeis.org/A013998):

271441, 904631, 16532714, 24658561, 27422714, 27664033, 46672291,…

Otras propiedades

La paridad de X(n) recorre el ciclo {1, 0, 0, 1, 0, 1, 1} Es fácil de ver: las tres primeras vienen determinadas por la definición (en color rojo en la imagen). Las siguientes dependen de dos anteriores. Por tanto, existirá ciclo si se vuelve a repetir el par 1 0, y esto ocurre siete lugares más adelante (color verde):

Para ampliar el tema puedes visitar http://www.mathpages.com/home/kmath345/kmath345.htm en el que se incluye la espiral triangular creada con estos números.

Propiedades matriciales

Estas entradas sobre sucesiones recurrentes también se plantean el objetivo de un mayor conocimiento de las hojas de cálculo. Por eso vamos a aprovechar las propiedades matriciales de la sucesión de Perrin para  repasar este tipo de funciones.

La primera propiedad matricial se resume en la siguiente fórmula para n>2:

Recuerda que la traza es la suma de los elementos de la diagonal principal de una matriz cuadrada.
Para comprobarlo con una hoja de cálculo organizaremos este esquema:

Comenzamos escribiendo a la izquierda la matriz M dos veces, y a la derecha las multiplicamos.
Para ello usaremos la función matricial MMULT, pero como es de tipo matricial deberás seleccionar la matriz de la derecha (debajo del rótulo “Potencia n de M), después escribir una fórmula similar a esta: =MMULT(C3:E5;G3:I5), tomando como rangos los de las matrices de la izquierda. Cuando escribas la fórmula no termines con Intro, sino con la combinación Ctrl+Mayúsc+Intro, para indicar que la fórmula es de tipo matricial. Notarás que lo has escrito bien porque la fórmula se verá entre corchetes.

A la derecha de las matrices puedes incluir la traza de la tercera, que en la imagen te da 2. Después copia la tercera sobre la primera matriz con copia sólo de valores, y te resultará el siguiente número de Perrin, en este caso 3, porque esta propiedad genera la sucesión a partir del tercer término.
Seguirían 2, 5, 5, 7, 10,…

Variante de la anterior expresión

Si en lugar de usar la traza empleamos un producto por la matriz (en vertical) (3, 0, 2), obtenemos tres términos en lugar de uno. La expresión sería ahora:

Bastaría borrar la traza en el anterior esquema y sustituirla por otro nuevo producto matricial con la (3, 0, 2). Lo dejamos como ejercicio. Aquí tienes la generación de los términos 5, 7 y 10

Comprobar conjeturas con hoja de cálculo. Goldbach.

La formulación más simple de la Conjetura de Goldbach es:

Todo número par mayor que 2 es suma de dos primos

Fue propuesta por Goldbach el 7 de Junio de 1742, en una carta dirigida a Euler.  En realidad, su propuesta se refería a la conjetura ternaria: " Todo número impar es la suma de tres primos" y Euler le respondió con la propuesta binaria que todos conocemos.

Ha sido comprobada hasta números muy grandes, pero no se ha podido demostrar. No obstante, se han logrado resultados provisionales:

Cualquier número par es suma de 6 o menos números primos.(Ramaré 1995)

Todo número par suficientemente grande es suma de un primo y del producto de dos primos.(Chen 1966)

Todo número impar N mayor que 5 es suma de tres primos. (Demostración de la conjetura ternaria a cargo de Vinogradov en 1937).

En esta entrada sólo nos plantearemos, como en toda la serie que vamos desarrollando sobre conjeturas, la comprobación de algunos aspectos de la misma mediante el uso de la hoja de cálculo.

Primer nivel

Comprobaremos la conjetura en tres niveles distintos, según el uso que se haga del lenguaje de macros. En primer lugar lo efectuaremos con las técnicas usuales de las hojas de cálculo. Usaremos la hoja Conjeturas, alojada en la página

http://www.hojamat.es/sindecimales/divisibilidad/herramientas/herrdiv.htm#global

(Búscala en la relación de herramientas)

Organizaremos la comprobación según un esquema similar a este:

Escribimos un número par mayor que 2 en una celda. En la imagen es el 612. Después ordenamos los números primos en columna, hasta el límite que queramos. Para ello escribimos el 2, debajo primprox(2) (función implementada en esta hoja, y que encuentra el primo siguiente a uno dado).

Rellenamos hacia abajo y nos resultará la lista de primos.

En una segunda columna escribiremos una fórmula similar a a siguiente, que copiaremos de arriba a abajo:

=SI(Y(F15<=H$11/2;esprimo(H$11-F15));H$11-F15;"")

En ella H11 es la celda donde hemos escrito el 612. En tu caso podrá ser otra. La F15 en nuestro esquema apunta al número primo que tiene a su izquierda. De esa forma, la podemos interpretar así: “Si el primo no llega a la mitad del número par probado (aquí el 612) y su diferencia con él es otro primo, escribo esa diferencia, pero en caso contrario dejo la celda en blanco”.

Es sencillo de entender y funciona escribiendo los pares de primos en los que se descompone el 612. En la imagen 5+607, 11+601, 19+593,…hasta un total de 26 pares. Si no logras ese número, deberás rellenar hacia abajo las dos columnas hasta llegar a la mitad de 612

Este esquema puede aclarar, probando con varios pares, el sentido de la conjetura. También te da confianza en ella, pues no sólo existe un par de primos para cada número par, sino muchos. ¡Pero no se ha probado aún!

Nivel 2

Ya que con el esquema anterior nos han resultado varias descomposiciones en primos para cada número par, podíamos simplificar mucho si lo plasmáramos en una función. En la hoja de cálculo que estamos usando hemos implementado NUMGOLDBACH(N), que devuelve un cero si N no es par y el número de descomposiciones si es par. En el caso del 612 devuelve correctamente 26.

Aquí tienes los primeros resultados. Si la conjetura es cierta, deberán ser todos mayores que 0. Están recogidos en http://oeis.org/A045917

Merece la pena recorrer la codificación de esta función y así entenderás mejor las cuestiones.

Public Function numgoldbach(n)
Dim ng, i

If n <> 2 * Int(n / 2) Then  ‘si es impar devuelve un cero (valor de ng
ng = 0
Else
i = 2: ng = 0
While i <= n / 2  ‘si es par recorre todas las posibles sumas de primos
If esprimo(n - i) Then ng = ng + 1 ‘si el segundo sumando es primo, incrementa el contador ng
i = primprox(i) ‘esta línea asegura que el primer sumando sea primo
Wend
End If
numgoldbach = ng
End Function

Nivel 3

Podemos dejar que sea la hoja de cálculo la que recorra automáticamente los primeros números hasta un tope o hasta que numgoldbach dé un cero. Como lo segundo es imposible para números pequeños (ya está comprobada la conjetura), el resultado final será siempre un cero.

Podíamos usar un esquema similar al siguiente:

Escribimos un tope, pulsamos el botón e irán apareciendo valores de Numgoldbach, ninguno nulo, hasta finalizar la búsqueda. Si uno fuera cero, se interrumpiría el proceso con un solemne mensaje. La programación del botón podría ser similar a esta:

Sub buscagoldbach()
Dim i, g, p

p = ActiveWorkbook.Sheets(1).Cells(5, 3).Value ‘lee el tope
i = 4 ‘inicio búsqueda
g = 1’inicio valor de numgolbach
While g <> 0 And i <= p
i = i + 2 ‘busca de 2 en 2
g = numgoldbach(i)
ActiveWorkbook.Sheets(1).Cells(6, 3).Value = g ‘escribe el valor de g
If g = 0 Then MsgBox ("¡Contraejemplo!") ‘Esto no va a ocurrir
Wend
End Sub

Variantes

Variante ternaria

“Todo número impar mayor que 5 es la suma de tres primos"

No vamos a repetir con ella los tres niveles anteriores. El primer nivel necesitaría un estructura de datos tridimensional, poco intuitiva en una hoja de dos dimensiones. El tercero sería semejante al del primer caso. Así que sólo desarrollaremos un esquema con todas las posibles descomposiciones en tres sumandos primos:

Como en el caso anterior, no vamos a analizar si el número es impar o no. Simplemente hemos programado un botón que lo descompone en esos sumandos de todas las formas posibles (lo haremos con sumandos decrecientes)

Para quien le guste la programación, ahí tiene explicado el algoritmo que hemos usado:

Sub goldbach3()
Dim fila, a, b, c, n

fila = 7
n = ActiveWorkbook.Sheets(1).Cells(fila, 3).Value ‘lee el número, que se encuentra en C7
a = 2 ‘primer sumando primo
While a < n ‘ el primer sumando llega hasta n en lo posible
b = 2
While b < a ‘ el segundo es inferior al primero
c = n - a – b ‘tercer sumando
If esprimo(c) And c <= b Then ‘si el tercer sumando también es primo, se presenta el resultado
fila = fila + 1
ActiveWorkbook.Sheets(1).Cells(fila, 3).Value = a
ActiveWorkbook.Sheets(1).Cells(fila, 4).Value = b
ActiveWorkbook.Sheets(1).Cells(fila, 5).Value = c
End If
b = primprox(b) ‘se incrementa el segundo primo
Wend
a = primprox(a) 'se incrementa el primero
Wend
End Sub 

Como era de esperar, siempre aparecen los tres sumandos primos. Se deja a los lectores el definir una función que cuente las soluciones. Siempre existirá al menos una.

Expresión mediante equidistancia

Un comentario a la entrada http://culturacientifica.com/2013/06/26/la-conjetura-de-goldbach/ me ha dado la idea de organizar una comprobación distinta.

Si la conjetura es cierta, para todo número par 2N, si es la suma de dos primos p y q, con p>q, cumplirán que p+q=2N, o bien que p-N=N-q:

Todo número entero positivo mayor que 4 es equidistante de dos primos

Es fácil ver que es otra formulación distinta de la conjetura de Goldbach. En los párrafos anteriores hemos visto la consecuencia directa. A la inversa, si es cierto que todo N equidista de dos primos, dado un par 2N aplicamos p-N=N-q para cierto par de primos, con lo que 2N=p+q. El exigir que sea mayor que 4 es porque no habría primos inferiores para números menores.

Es muy fácil organizar la comprobación con esta variante. Lo efectuaremos en el Nivel 1, de cálculo manual:

Escribimos la lista de números consecutivos 1, 2, 3,…y los sumamos y restamos con el número dado. Después, en la tercera columna, escribimos una fórmula similar a

=SI(Y(ESPRIMO(D9);ESPRIMO(E9));"SI";"") 

que nos devuelve un “SI” si los equidistantes son primos.

En la imagen, 17=29-12 y 41=29+12.