viernes, 17 de abril de 2015

Sucesión de Padovan


En una entrada anterior estudiamos la sucesión de Perrin

(http://hojaynumeros.blogspot.com.es/2014/11/sucesion-de-perrin.html).

La de hoy, de Padovan, es muy parecida, por lo que se recomienda leer antes la entrada enlazada. Recordamos:

La sucesión de Perrin 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

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

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

Condiciones iniciales: X0=3   X1=0  X2=2 Ecuación de recurrencia: Xn= Xn-2+Xn-3

Pues bien, la sucesión de Padovan es similar, pero con distintos valores iniciales:

X0=1   X1=1  X2=1

Como con la anterior, podemos construirla con nuestra herramienta de hoja de cálculo adaptada a las sucesiones recurrentes de tercer orden.

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

Escribimos los coeficientes 0, 1,1 y los valores iniciales 1, 1, 1:



Y obtenemos:



Son los números espirales de Padovan contenidos en http://oeis.org/A134816. Existen otras variantes de esta sucesión, pero nos dedicaremos en esta entrada a la que comienza con 1, 1, 1. Por el carácter de este blog, omitiremos propiedades gráficas, como la espiral de triángulos, que puedes consultar en otras páginas.

Relaciones recurrentes

Para abreviar a los términos de esta sucesión los identificaremos como P(n).

En muchas páginas web podrás encontrar otras relaciones recurrentes además de la de la definición, P(n)=P(n-2)+P(n-3). Aquí sólo comentaremos alguna dejando como ejercicio el análisis de las demás.

(1)  P(n)=P(n-1)+P(n-5)

Se puede verificar por inducción: Se cumple en los primeros términos, como puedes comprobar con la misma hoja de cálculo:



Extensión a P(n+1)

P(n+1)=P(n-1)+P(n-2)=P(n-2)+P(n-6)+P(n-3)+P(n-7)=P(n)+P(n-4), luego se cumple la inducción completa.

(2) P(n)= P(n-2)+P(n-4)+P(n-8)

Sólo veremos los primeros términos con hoja de cálculo y dejaremos la demostración por inducción como ejercicio.



Hay más relaciones de este tipo. Las tienes en

http://es.wikipedia.org/wiki/Sucesi%C3%B3n_de_Padovan

Una interesante es la que relaciona la sucesión de Perrin con la de Padovan:

Perrin(n)=P(n+1)+P(n-10)

Con nuestra hoja hemos construido este esquema para que compruebes que se cumple para los primeros términos. El justificarlo por inducción es fácil por compartir ambas sucesiones la misma fórmula de recurrencia.



Ecuación característica

La ecuación característica correspondiente será x3-x-1=0, es decir, la misma que para la sucesión de Perrin. Con el botón Resolver de esa hoja obtienes las tres soluciones de la ecuación, una real y dos complejas



La solución real 1,32471… es el número plástico y,  que ya presentamos en el estudio de la sucesión de Perrin. También la sucesión de Padovan se acerca progresivamente a las potencias de este número, como puedes ver en este cálculo realizado con nuestra hoja:



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))

Crea un archivo de texto “sucesión.txt” en la misma carpeta de PARI y verás cómo te reproduce la sucesión:

x + x^2 + x^3 + 2*x^4 + 2*x^5 + 3*x^6 + 4*x^7 + 5*x^8 + 7*x^9 + 9*x^10 + 12*x^11 + 16*x^12 + 21*x^13 + 28*x^14 + 37*x^15 + 49*x^16 + 65*x^17 + 86*x^18 + 114*x^19 + O(x^20)

Los coeficientes del polinomio reproducen la sucesión de Padovan, con el índice desfasado en 1 porque hemos comenzado con el valor 0.

Relación con cuestiones combinatorias

Todas las sucesiones recurrentes suelen tener relación con particiones y composiciones (particiones con orden), porque su generación a partir de elementos anteriores puede coincidir. En el caso de la sucesión de Padovan también existen esas relaciones. Veamos:

P(n) coincide con las composiciones de n+2 en sumandos 2 y 3

En efecto, P(0)=P(1)=P(2) valen 1, que son las formas de descomponer 2, 3 y 4 en sumandos ordenados 2 y 2. P(3)=2 porque 5=2+3=3+2. P(4)=2, ya que 6=3+3=2+2+2.

Con nuestra hoja Cartesius (aún no publicada) se pueden comprobar estos desarrollos. Por ejemplo, para el caso de 8, plantearíamos:

Aunque no conozcas su sintaxis, basta explicarte que hemos pedido que desde 1 hasta 8, usando el conjunto {2,3} busque todas las sumas iguales a 8 con repetición.

Efectivamente, resultan 4=P(6)


En general, cualquier suma correspondiente a N resultará de añadir un 2 a las composiciones de N-2 y un 3 a las de N-3, por lo que su generación es idéntica a la de la sucesión de Padovan. Tal como nos ocurrió con la sucesión de Narayana (http://hojaynumeros.blogspot.com.es/2015/01/sucesion-de-las-vacas-de-narayana.html), esta descomposición da lugar a la expresión de los números de Padovan como suma de números combinatorios. En http://en.wikipedia.org/wiki/Padovan_sequence tienes uno de ellos:



Así, por ejemplo, en el desarrollo para k=11 con Cartesius vemos clara la descomposición en números combinatorios (recuerda que las permutaciones con repetición y dos elementos equivalen a esos números)



Para quienes apreciáis las técnicas de programación, insertamos esta función por si queréis implementarla en vuestra hoja de cálculo:

Public Function padovan(n)
Dim p, q, t, s, i, nn

 nn = n + 2: p = Int(nn / 2): q = nn - 2 * p: t = 0

While p >= q
s = 1: For i = 0 To q - 1: s = s * (p - i) / (q - i): Next i 'Calcula el número combinatorio
t = t + s 'Suma los números combinatorios
p = p - 1: q = q + 2
Wend
padovan = t
End Function






martes, 7 de abril de 2015

Algunas curiosidades sobre la antisigma

Se ha publicado algo, no mucho, sobre algunas propiedades y curiosidades acerca de la antisigma. Destacamos algunas y aportaremos otras.

Antisigmas cuadradas

La antisigma de un número puede ser un cuadrado. Esto ocurre en los siguientes números:

1, 2, 5, 6, 14, 149, 158, 384, 846, 5065, 8648, 181166, 196366, 947545, 5821349, 55867168, 491372910, 4273496001, 40534401950,… http://oeis.org/A076624

En el caso de números primos, como el 2 y el 5, deberemos resolver una ecuación diofántica de segundo grado, ya que (p+1)(p-2)/2=k2. Donovan Johnson ha encontrado la recurrencia que genera otros casos en la página de OEIS enlazada. El siguiente primo sería 8946229758127349.

Nosotros también nos aproximaremos al tema mediante una ecuación de Pell. Transformamos la igualdad multiplicando todo por 4 y desarrollando:

4p2-4p-8=8k2

 Completamos un cuadrado y queda: (2p-1)2-8k2=9 Cambiamos de variables haciendo X=(2p-1)/3  Y=k/3, obteniendo la ecuación de Pell

X2-8Y2=1

Tenemos una herramienta para resolver esta ecuación, en la dirección

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

Aplicándola para el caso D=8 obtenemos varias soluciones para X e Y



Si a ellas añadimos la solución trivial X=1, Y=0 y deshacemos el cambio X=(2p-1)/3 mediante p=(3X+1)/2, obtendremos todas las soluciones para p. En la imagen que sigue hemos añadido una columna para saber si son primos o no los valores obtenidos, y vemos (los de color rojo) que coinciden con los valores primos de la sucesión:



Con una herramienta más potente podemos seguir con las iteraciones y llegar a la siguiente solución prima dada por Donovan Johnson e incluso sobrepasarla. No damos muchos detalles por no alargar.

La iteración de Pell en este caso es
(ver la teoría en http://hojaynumeros.blogspot.com.es/2010/02/ecuacion-de-pell.html)


La aplicamos reiteradamente en PARI comenzando con X=1 Y=0, y tomamos nota de las soluciones de X que son primas. Obtendremos un listado en el que aparecerán los cuatro primos obtenidos aquí, más el propuesto por Donovan Johnson, 8946229758127349, y alguno más.

Programa en PARI

{x=1;y=0;for(n=1,60,m=(3*x+1)/2;if(isprime(m),print1(m,", "));p=3*x+8*y;q=x+3*y;x=p;y=q)}

Resultado obtenido:

2, 5, 149, 5821349, 8946229758127349, 13748537282247342677718149, 828287615476676026361062299923143963349, 32470531080787945457870876690417952137154149,

Aparecen tres nuevas y enormes soluciones. Este tipo de descubrimientos hacen que sigamos con estos temas a pesar del cansancio de los años.

Antisigmas triangulares

Sólo hemos encontrado el caso ya estudiado de las potencias de 2. No parece que haya ningún número que no sea potencia de 2 y tenga antisigma triangular.

Antisigmas primas

Estos son los números con antisigma prima:

3, 4, 10, 21, 34, 46, 58, 70, 85, 93, 118, 129, 130, 144, 178, 201, 226, 237, 262, 298, 310, 322, 324, 325, 333, 334, 346, 382, 406,... http://oeis.org/A200981

Sólo figura entre los primos el 3, porque si (p+1)(p-2)/2 ha de ser primo, si p es mayor que 3 aparecerían dos factores al menos en la antisigma.

Llama la atención la abundancia de semiprimos. Según la fórmula que vimos en la entrada anterior, deberá darse la casualidad de que si N=pq, se dé que pq(pq+1)/2-(p+1)(q+1) sea primo. Por ejemplo, 46=2*23 y su antisigma sería 46*47/2-3*24=1009, que es primo.

Antisigma y sigma coprimas

Terminamos por ahora con otra curiosidad: Números en los que sigma y antisigma son primos entre sí:

4, 9, 10, 16, 21, 22, 25, 34, 36, 46, 49, 57, 58, 64, 70, 81, 82, 85, 93, 94, 100, 106, 118, 121, 129, 130, 133, 142, 144,…

Al principio parece que pertenecerán a ella los cuadrados, pero a partir de 196 hay muchos que no cumplen esta propiedad: 441, 1521, 1764, 3249, 3600,...

En esta sucesión todos son compuestos, pues un primo p tiene como sigma p+1 y como antisigma (p+1)(p-2)/2, con lo que el factor (p+1)/2 divide a ambas para un primo mayor que 2. Y en el caso del 2, su sigma es 3 y su antisigma 0, que no tiene divisores


viernes, 27 de marzo de 2015

Antisigma de un número natural


Al igual que se ha definido la función SIGMA(N) como la suma de todos los divisores de N (incluido él mismo), podemos definir la ANTISIGMA(N), que es la suma de los números menores que N y que no lo dividen, Por ejemplo, la antisigma de 8 sería la suma de 3+5+6+7=21, y sigma(8) es igual a 1+2+4+8=15.

Los valores de esta función antisigma son los siguientes, que están incluidos en https://oeis.org/A024816

0, 0, 2, 3, 9, 9, 20, 21, 32, 37, 54, 50, 77, 81, 96, 105, 135, 132, 170, 168, 199, 217, 252, 240, 294, 309, 338, 350,…

Comprueba que para el 8 se da el valor de 21, que es el que hemos calculado.

No se debe confundir con la indicatriz de Euler, que cuenta los números primos con N. En la suma correspondiente al 8 figura el 6, que no divide al 8, pero tampoco es primo con él. Con estos números primos con N se puede formar también una suma. Nosotros, para entendernos, la llamaremos S_EULER. Puedes consultar la página https://oeis.org/A023896. En el caso del 8, la suma sería 1+3+7=11 (por convención, se considera el 1 primo con todos los demás números. Esto facilita los cálculos). Es evidente que S_EULER(N) será siempre menor o igual que ANTISIGMA(N), porque sus sumandos están incluidos en la otra suma.

La suma de SIGMA(N) y ANTISIGMA(N) es muy fácil de calcular, ya que se trata de sumar todos los números desde 1 hasta N, y esto sabemos que es igual a N(N+1)/2.

Relación fundamental: SIGMA(N)+ANTISIGMA(N)=N(N+1)/2

Si no se dispone de la función SIGMA, también se puede encontrar ANTISIGMA. Por ejemplo, puedes usar este código Basic de Excel:

Public Function antisigma(n) 'suma los no divisores
Dim i, a

a = 0
For i = 1 To n
If n / i <> n \ i Then a = a + i  ‘no es divisor de n, y se suma
Next i
End If
antisigma = a
End Function

Si ya tienes implementada SIGMA, el desarrollo es mucho más simple:

Public Function antisigma(n)
antisigma = n * (n + 1) / 2 - sigma(n)
End Function

Así lo haremos en PARI

antisigma(x)=x*(x+1)/2-sigma(x)

Antisigmas calculables mediante una fórmula 

Nos referimos a una fórmula sencilla, sin tener que proceder a una descomposición complicada en factores primos.

Antisigma de un número primo

Si p es primo, es posible encontrar una fórmula para la antisigma. En efecto, por la relación anterior, antisigma(p)=p(p+1)/2-sigma(p)=p(p+1)/2-(p+1)=(p2+p-2-2p)/2=( p2-p-2)/2=(p+1)(p-2)/2, Hemos usado el hecho de que la sigma de un primo p equivale a p+1, como es evidente.

Así que si p es primo, es válida la fórmula

Por ejemplo, antisigma(7)=8*5/2=20, antisigma(13)=14*11/2=77

Es curioso el hecho de que esta función sea evaluable directamente en este caso. Constituye una relación cuadrática, y su diagrama conjunto forma una parábola.


En los números compuestos hay que descomponer en factores primos previamente, y se pierde así una relación tan directa.

Antisigma de una potencia de 2

Si N=2k, entonces sigma(N)=1+2+4+8+…2k=(2k+1-1)/(2-1)=2k+1-1. Aplicamos la relación fundamental y nos queda:
Antisigma(2k)=2k(2k+1)/2-(2k+1-1)=(22k+2k-2*2k+1+2)/2=(22k-3*2k+2)/2
(2k-2)(2k -1)/2=(N-1)(N-2)/2  y también (2k-1-1)(2k -1) (ver http://oeis.org/A134057)

La antisigma de una potencia de 2 es un número triangular. Si la potencia es N, su antisigma hemos visto que es (N-1)(N-2)/2, el triangular de orden N-1. Lo puedes comprobar en este listado:


Antisigma de semiprimos con factores diferentes

Un caso también sencillo es el de semiprimos producto de dos primos diferentes. Si N=pq, con p y q primos, es posible encontrar una fórmula sencilla para la antisigma. Los valores son los siguientes:

N      Antisigma(N)
6 9
10 37
14 81
15 96
21 199
22 217
26 309
33 513r
34 541
35 582
38 681

Busquemos la fórmula que los genera: La sigma de un número primo p es p+1, luego de q será q+1. Como es una función multiplicativa, la sigma del producto equivaldrá a (p+1)(q+1) y por tanto la antisigma pq(pq+1)/2-(p+1)(q+1). Parece que queda mejor así sin intentar simplificarla. La comprobamos: 35=5*7, luego su antisigma será 35*36/2-6*8=35*18-48=582

Antisigma de la potencia de un primo

Es otro caso sencillo. La sigma de una potencia de primo pr viene dada por 1+p+p2+p3+…+ pr, es decir:



Por tanto, la antisigma vendrá dada por



Lo comprobamos: según el primer listado, la antisigma de 27 es 338, y según esta fórmula se obtendría

Asig(33)=27*28/2-(81-1)/2=338

jueves, 19 de marzo de 2015

Autonúmeros (2)


En la entrada anterior descubrimos que todos los números se pueden clasificar en generados por la digitadición de Kaprekar y autonúmeros, que no pueden ser generados. Nos dedicaremos en primer lugar a los generados, y con ellos descubriremos algo más sobre los autonúmeros.

Números generados

Señalamos que los números 2, 4, 6, 8, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 21, 22, 23, 24, 25, 26, 27, 28,… (http://oeis.org/A176995) son los complementarios de los autonúmeros. Todos ellos se pueden expresar como N+SUMACIFRAS(N), por lo que Kaprekar les llamó “generados”. A este número N le podemos llamar generador del mismo, pero desafortunadamente no es único, por lo que no podemos considerarlo una función dependiente del número. En efecto, existen números, como el 103, que poseen más de un desarrollo de este tipo, ya que 103=92+9+2 y también 103=101+1+1.

Podíamos definir una función GENERADOR si para cada N eligiéramos el mayor K que cumple que K+SUMACIFRAS(K)=N, y asignar el valor 0 como generador de los autonúmeros. Así sí sería una función. La función AUTO descrita más arriba y debidamente adaptada nos servirá para este propósito:

Public Function generador(n)
Dim k, g

g = 0
While k < n
If k + sumacifras(k) = n Then g = k
k = k + 1
Wend
generador = g
End Function

En esta tabla puedes comprobar que si a cada número de la derecha le sumas la suma de sus cifras, resulta el de la izquierda, salvo el caso del autonúmero 97 al que le hemos asignado un cero.



Si construimos un diagrama de dispersión entre N y su generador, obtenemos un resultado similar a este:



Abajo los puntos de valor 0 representan a los autonúmeros, y los de arriba parecen presentar pautas lineales escalonadas, pero es una ilusión, ya que esos tramos no se podrían solapar. Lo que ocurre es que en este proceso los números consecutivos aparecen muy cercanos, y dan la impresión de sucederse. Analizaremos estos números con más detalle.

Si representamos un número N de k+1 cifras en base 10, lo estamos generando mediante un polinomio del tipo


Si ahora le sumamos sus cifras (digitadición) se convertirá en


Todos los números generados se podrán expresar de esta forma, y los autogenerados, no.
Según esta fórmula, los dobles de las cifras 1…9 serán todos números generados, y si dos generadores se diferencian sólo en una unidad en las decenas, sus generados se diferenciarán en 11, como ya hemos observado. Si se diferencian en una unidad de las centenas, sus generados se diferenciarán en 101, y así. Estas correspondencias también se dan en casi todos los autonúmeros, pues van la par de estos. No se da en todos.

En este gráfico hemos descompuesto cada autonúmero en relación con los menores de 100, y vemos una repetición de tramos lineales a unas alturas de 1001, 2002, 3003 y 4004. Esto es consecuencia de la fórmula que hemos desarrollado.



Fundamento del algoritmo de Kaprekar


La fórmula

nos da una forma paralela a la de Kaprekar para decidir si N es autonúmero en pocos pasos (es el mismo proceso con otra orientación). Observa que todos los coeficientes son congruentes con 2 módulo 9, con lo que si llamamos G al posible generador de N y SC(G) a la suma de sus cifras, se cumplirá que Nº2SC(G) (mod 9. Para saber esto no hacía falta la fórmula, pues como G es congruente con SC(G) módulo 9 y se cumple que N=G+SC(G), es evidente que N será congruente con el doble de SC(G).

Como 2 es primo con 9, en la ecuación Nº2SC(G) (mod 9 se podrá encontrar una solución S, con lo que G=9*k+S para valores de k no superiores al número de cifras. Lo vemos con el 30: Si N=30, será 30º3º2*SC(G). En este caso SC(G)º6, porque 6+6º3 (mod 9. De esta forma G puede ser 6, 15 0 24. Recorremos de mayor a menor y descubrimos que el 24 vale, porque 24+2+4=30, luego 30 no es autonúmero, sino generado.

Probamos con un número mayor, como 4327, del que ya sabemos por otros medios que es autonúmero:

Hallamos el resto módulo 9 de 4327 (puedes dividir mentalmente o usar la función RESIDUO de Excel) y resulta ser 7. Planteamos 7º2SC(G) (mod 9. También, mentalmente, vemos que SC(G) º8 (mod 9. Por tanto, vamos restando de 4327 números del tipo 9*k+8 hasta ver si la suma de las cifras de la diferencia es ese número:

K=0: 4327-8=4321, y su suma de cifras no es 8.
K=1: 4327-17=4310 y no suman 17
K=2: 4327-26=4301. No suman 26,
K=3: 4327-35=4292. No vale
Acabamos con k=4, porque la suma de cifras ya no puede ser mayor:
K=4: 4327-44=4283, de suma 17, luego tampoco es válido.

Hemos comprobado, con la misma técnica de Kaprekar y distinta presentación, que 4327 es autonúmero. No puede ser generado.

Para los aficionados a la programación, esta es la función que podemos usar:

Public Function esauto3(n)
Dim nc, a, r, h, k
Dim es As Boolean

'número de cifras
a = 1: nc = 0
While a <= n
a = a * 10: nc = nc + 1
Wend
'módulo 9
r = n Mod 9
For k = 0 To 8: If (r - 2 * k) Mod 9 = 0 Then h = k
k = 0: es = True
While k <= nc And es
If sumacifras(n - 9 * k - h) = 9*k+h Then es = False
k = k + 1
Wend
esauto3 = es
End Function

Hemos construido un esquema de hoja de cálculo en el que vemos la equivalencia de las tres funciones que hemos programado en estas entradas. En la imagen tienes un ejemplo:



Sucesiones recurrentes de números generados

Kaprekar estudió también las sucesiones que se forman si a un número cualquiera le vamos aplicando la digitadición tanto a él como a los resultados sucesivos. Por ejemplo, si comenzamos con el 12 tendríamos la sucesión 12, 15, 21, 24, 30,…Evidentemente es infinita, creciente y todos sus elementos, salvo quizás el primero, son generados. Kaprekar destacó que si restamos el último del primero y sumamos las cifras del último, resulta la suma de cifras de todos ellos. Por ejemplo, aquí, 30-12+3=21=1+2+1+5+2+1+2+4+3+0

Aunque el autor le dio mucha importancia, basta recorrer la generación para darse cuenta de su validez: 12+3+6+3+6=30, luego al restar resultarán las sumas de cifras menos la del último. Es una obviedad.


viernes, 13 de marzo de 2015

Autonúmeros (1)


Autonúmeros o números colombianos

Estos números, llamados también de Devlali, fueron descritos por el matemático indio Kaprekar (el de la constante 6174, http://es.wikipedia.org/wiki/Constante_de_Kaprekar). Son muy conocidos por haberlos presentado Martin Gadner en uno de sus libros. En esta dirección puedes leer su artículo.

http://librosdemates.blogspot.com.es/2013/01/viajes-por-el-tiempo-y-otras.html

Su definición es algo extraña, porque son aquellos números enteros positivos que no pueden ser expresados como la suma de otro entero con la suma de sus cifras (Kaprekar llamó a esta operación digitadición). Por ejemplo, 20, no puede generarse con números más pequeños a los que les sumamos sus cifras. Con los de una cifra se ve que es imposible, y con los de dos: 11+1+1=13, 12+1+2=15, 13+1+3=17, 14+1+4=19, 15+1+5=21, y el resto tampoco daría como resultado 20.

Con esta definición se comprende que existan infinitos tipos de autonúmeros, dependiendo de la base elegida, y así se ha demostrado. Nosotros nos limitaremos a la base 10. En este caso los tienes en http://en.wikipedia.org/wiki/Self_number y son estos:

1, 3, 5, 7, 9, 20, 31, 42, 53, 64, 75, 86, 97, 108, 110, 121, 132, 143, 154, 165, 176, 187, 198, …

También los puedes estudiar en http://oeis.org/A003052

Estos números proceden de una criba. Puedes ir calculando todos los resultados posibles si se suma cada número con la suma de sus dígitos en la base dada. Resultarían estos:

2, 4, 6, 8, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 21, 22, 23, 24, 25, 26, 27, 28,…, que están contenidos en http://oeis.org/A176995 y los que nos ocupan serían su complemento.

Algoritmo de fuerza bruta

Una cuestión interesante es cómo saber si un número es autonúmero o no. El primer acercamiento a este tema puede consistir en recorrer todos los números inferiores a N, añadirles la suma de sus cifras y comprobar si el resultado es N. Desembocamos entonces en la programación con hojas de cálculo. Necesitaríamos:


  •  La función SUMACIFRAS
  •  Un bucle que recorra los números k de 1 a N y pare si el resultado de k+sumacifras(k) es N.
  •  Si la prueba anterior es positiva, declaramos que el número N no es autonúmero, y si es negativa para todos los números inferiores a N, diremos que sí lo es.


La función SUMACIFRAS debemos tenerla publicada en algún documento. Por si no fuera así, la reproducimos en Basic de hoja de cálculo:

Public Function sumacifras(n)

'No analiza si el número es entero positivo

Dim h, i, s, m

h = n ‘De la variable h se irán extrayendo las cifras
s = 0 ‘Esta variable recogerá la suma de cifras
While h > 9 ‘Bucle para extraer las cifras una a una
i = Int(h / 10)
m = h - i * 10
h = i
s = s + m ‘La nueva cifra se suma a la variable
Wend
s = s + h ‘La cifra residual se suma a la variable
sumacifras = s
End Function

Una vez tenemos la función sumacifras podemos organizar el test en forma de función booleana:

Public Function esauto(n)
Dim es As Boolean
Dim k

es = True ‘Se declara de entrada que el número es “colombiano”
k = 0  ‘Comenzamos a recorrer los números menores que n
While k < n And es  ‘No paramos hasta llegar a n o hasta que es sea falso
If k + sumacifras(k) = n Then es = False ‘Si hay igualdad, paramos y declaramos “False”
k = k + 1
Wend
esauto = es ‘La función recoge el valor de es
End Function

Aplicamos esta función a los primeros números y obtendremos el valor VERDADERO en los
autonúmeros:



Por ejemplo, busquemos el primer autonúmero que sigue a 1000:


Resulta ser el 1006.

Test de Kaprekar

El anterior algoritmo recorre demasiados números y es ineficiente para enteros grandes. El mismo Kaprekar demostró un test en el que sólo hay que analizar enteros no superiores al número de dígitos. Lo copiamos de la página http://www.numbersaplenty.com/set/self_number/



Lo hemos implementado como función para Excel. No detallamos el código. Sólo lo copiamos con algún comentario:

Public Function esauto2(n)
Dim nc, a, r, h, k
Dim es As Boolean

'número de cifras
a = 1: nc = 0
While a <= n
a = a * 10: nc = nc + 1
Wend
'Se preparan las variables del test de Kaprekar
If nc = 0 Then esauto2 = False: Exit Function
r = 1 + ((n - 1) Mod 9)
If r / 2 = r \ 2 Then h = r / 2 Else h = (r + 9) / 2
'Bucle del test
es = True
k = 0
While k <= nc And es
If sumacifras(Abs(n - h - 9 * k)) = h + 9 * k Then es = False
k = k + 1
Wend
esauto2 = es
End Function

Hemos preparado un esquema en el que se puede elegir el inicio y compara las dos funciones, ESAUTO y ESAUTO2. Por si se hubiera deslizado algún error se han efectuado varias comprobaciones y coinciden ambas funciones, pero con gran lentitud en la primera.


En la siguiente entrada justificaremos este test ofreciendo simultáneamente otro similar.

Distribución de los autonúmeros

Es fácil ver que el incremento entre dos autonúmeros consecutivos es frecuentemente 11. Por ello esperamos tramos lineales en su distribución. Los veremos creando un gráfico con los primeros, pero si elegimos un gráfico con muchos puntos, se nos presenta una distribución prácticamente lineal, con pendiente 10,2, cercana al 11, que como hemos visto, aparece en muchos tramos.



Para ver mejor los tramos lineales elegimos un rango de números más pequeño:



Normalmente el salto de 11 se produce porque equivale a rebajar en una unidad las decenas cuando es posible. Así la suma total se rebaja en 11. Si el primero es autonúmero puede obligar a que el segundo también lo sea, pues si éste admitiera una suma, al rebajar esa unidad podría no ser autonúmero el primero. Esto vale sólo para la mayoría de los casos. No es una regla general. Ocurre algo parecido con 101, 1001, 10001, … Estos números nos aparecerán en la siguiente entrada.
Hay una forma sencilla de engendrar autonúmeros (no todos), y es comenzar por 9 y después usar la fórmula recurrente

Ck=8*10k-1+Ck-1+8

El siguiente número sería 8*10+9+8=97, el siguiente 800+97+8=905.

Con hoja de cálculo es fácil construir esta generación recurrente. Inténtalo si quieres. En la imagen hemos añadido a la derecha la función esauto.




miércoles, 4 de marzo de 2015

Números especiales que son un producto especial (3)


Con números cuadrados

Proseguimos el tema de buscar números especiales que son producto de dos factores también especiales. Hemos estudiado algunas variantes con triangulares, cuadrados, oblongos y primos. En esta entrada descompondremos cuadrados como producto de otros dos de algún tipo especial.

Un caso que no necesita estudio es el de cuadrados producto de cuadrados, pues esta condición la cumplen todos salvo los cuadrados de primos. Así que pasamos a otros productos. En todos ellos exigiremos que los factores del producto sean distintos, para evitar trivialidades, ya que si son iguales su producto sería cuadrado sin necesidad de buscar más.

Los algoritmos de esta entrada se basarán todos en la generación de cuadrados mayores que 1 que ya estudiamos, es decir, iniciar P=4 y K=5 y después, en cada pasada del bucle, hacer P=P+K y K=K+2, ya que sabemos que los incrementos de los cuadrados crecen de 2 en 2.

Producto de triangulares distintos

Hemos encontrado esta sucesión, que contiene muchos cuadrados de términos de https://oeis.org/A175497:

900, 7056, 32400, 44100, 88209, 108900, 298116, 705600, 1368900, 1498176, 2924100, 5336100, 8643600, 8820900, 9217296, 10432900, 15210000, 24147396, 37088100, 50893956, 50979600, 52490025, 55353600, 80568576, 114704100, 160123716, 200930625, 219632400, 265559616, 268304400, 296528400, 394657956,…

Pronto se llega a números grandes, por la fuerte condición que les exigimos. Por ejemplo, 88209 es el cuadrado de 297, y se puede descomponer en el producto de dos triangulares distintos: 88209=29403*3 y 29403 es el triangular número 242 (29403 =242*243/2) y 3 es el segundo (3=2*3/2).

Hemos comprobado los cálculos con dos métodos. Aquí tienes el de PARI

{i=4;j=5; ‘Inicia cuadrados
while(i<=5*10^8,k=3;p=3;c=0; ‘Inicia factores triangulares
while(k<sqrt(i)&&c==0,if(i/k==i\k&&ispolygonal(i/k,3)&&i/k>1,c=k);
‘Sólo llegamos en la búsqueda hasta la raíz cuadrada, para que ambos factores sean distintos
if(c>0,print(i);write1("final.txt",i,", ")); ‘Si c>0 se ha encontrado una solución
k+=p;p+=1);
i+=j;j+=2)
}

Ambos factores triangulares han de tener algún factor en común para que se forme un cuadrado.



Sus partes libres de cuadrados serán iguales, ya que es la única forma de que al final resulte un cuadrado. Lo puedes comprobar con algunos ejemplos de la tabla. Por ejemplo, 705600 es el triangular número 1187, y se descompone en el triangular 28, con parte libre 7, y el triangular 25200=60^2*7, que tiene también un 7 como parte libre.

Cuadrados producto de dos oblongos distintos

Este caso es equivalente al que exige que un cuadrado sea igual a cuatro veces el producto de dos triangulares distintos. Así, si el cuadrado 3600 es el producto de los dos oblongos 6=2*3 y 600=24*25, también cumple, como es evidente, que su cuarta parte es el producto de dos triangulares 3600/4=900=3*300.

Los primeros cuadrados que cumplen esto son:

144, 3600, 4900, 28224, 129600, 166464, 176400, 352836, 435600, 1192464, 2822400, 5475600, 5654884, 5992704, 11696400, 21344400, 34574400, 35283600, 36869184, 41731600, 60840000, 96589584, 148352400, 192099600, 203575824, 203918400, 209960100, 221414400, 322274304, 458816400,…

Aquí tienes la descomposición en producto de oblongos de los primeros:


Como todos son pares, ya tienen un factor común, y como en el caso anterior, sus partes libres de cuadrados han de ser iguales.

Los puedes conseguir en PARI con

{i=4;j=5;while(i<=5*10^8,k=2;p=4;c=0;while(k<sqrt(i)&&c==0,if(i/k==i\k&&issquare(4*i/k+1)&&i/k>1,c=k);if(c>0,print(i));k+=p;p+=2);i+=j;j+=2)}

Si usamos nuestra imaginación podemos recorrer muchas variantes, pero nos quedaremos con esta última:

Cuadrados producto de triangular y primo

9, 196, 225, 900, 2601, 3249, 4225, 15376, 53361, 88209, 136161, 176400, 181476, 191844, 324900, 450241, 461041, 1032256, 2152089, 2873025, 3960100, 7027801, 8643600, 11826721, 12744900, 17791524, 19193161, 28515600, 43956900, 45360225, 61230625, 63282025, 96216481, 108680625,…

Aquí es como si el número primo aportara el factor que se necesita para convertir un triangular en un cuadrado. Por tanto, el factor primo es igual a la parte libre de cuadrados del otro factor triangular.

Por tener esta propiedad disponemos de dos códigos distintos para encontrar esos números. Uno es el “largo”, que sigue lo explicado en estas entradas:

{i=4;j=5;
while(i<=5*10^8,k=3;p=3;c=0;
while(k<i&&c==0,if(i/k==i\k&&isprime(i/k)&&i/k>1,c=k);
if(c>0,print(i));
k+=p;p+=1);
i+=j;j+=2)
}

Podemos usar este otro, mucho más corto, pero que no nos ordena los números en orden creciente:

{i=1;j=2;while(i<=10^5,m=core(i);if(isprime(m),print1(i*m,", "));i+=j;j+=1)}

Como hemos indicado, se podría seguir con otros casos. Intenta reproducir este:

Cuadrados producto de un triangular y un oblongo
36, 900, 3600, 7056, 15876, 41616, 44100, 54756, 69696, 108900,…

viernes, 27 de febrero de 2015

Números especiales que son un producto especial (2)


Números triangulares como producto de otros 

En la entrada anterior planteamos el problema de buscar los números triangulares que son a su vez producto de otros dos triangulares. Presentamos una forma de generar cuadrados, oblongos y triangulares de la forma más rápida posible, ya que estas búsquedas se hacen lentas para números grandes, y construimos un algoritmo para generar estos números, contenidos en la sucesión de OEIS https://oeis.org/A188630

En esta entrada generaremos triangulares con otros tipos de números, y llegaremos a sucesiones que permanecían inéditas hasta ahora.

Triangulares producto de dos cuadrados

Aquí no hay que buscar mucho, ya que basta considerar que el número que cumple esto es aquel que es cuadrado (producto de cuadrados) y triangular a la vez. Esto está muy estudiado. Son estos:

1, 36, 1225, 41616, 1413721, 48024900, 1631432881,… y está publicados en https://oeis.org/A001110

Vemos, pues, otros casos.

Triangulares que son producto de un triangular y un cuadrado

Ahora tenemos ocasión de aplicar lo expuesto en la entrada anterior: cómo generar de forma rápida cuadrados y triangulares y cómo saber si son o no de ese tipo. Si esto te quedó claro, entenderás el algoritmo que sigue. Primero en Visual Basic de Excel:

Sub productriang()
Dim i, j, k, p, c
i = 3: j = 3 Generamos el primer triangular
While i <= 10 ^ 4
k = 3: p = 3: c = 0 Aquí generamos posibles divisores triangulares
While c = 0 And p < i  Cuando c<>0 se para la búsqueda y se comunica el resultado
If i / k = i \ k And escuad(i / k) And i / k > 1 Then c = k
La línea de arriba investiga si k es divisor y si i/k es cuadrado
If c <> 0 Then MsgBox (i)
k = k + p: p = p + 1  Genera el siguiente posible divisor triangular
Wend
i = i + j: j = j + 1  Genera el siguiente triangular para seguir buscando
Wend
End Sub

Si repasas la entrada anterior, es el mismo algoritmo propuesto en ella, pero cambiando estriangular(i/k) por escuad(i/k). Puedes repasar los comentarios que se incluían. Con un ligero cambio en el código, hemos situado los resultados en columna:


A partir de este valor se produce desbordamiento, por lo que acudimos a PARI y al código

{i=3;j=3;
while(i<=10^6,k=3;p=3;c=0;
while(k<i&&c==0,if(i/k==i\k&&issquare(i/k)&&i/k>1,c=k);
if(c>0,print1(i,", "));
k+=p;p+=1);
i+=j;j+=1)
}
No insertamos comentarios porque es el mismo que presentamos en la entrada anterior, salvo el uso de la función issquare (“es cuadrado”)

Con él puedes obtener los primeros números triangulares que son producto de un triangular y un cuadrado:


300, 1176, 3240, 7260, 14196, 25200, 29403, 41616, 64980, 97020, 139656, 195000, 228150, 265356, 353220, 461280, 592416, 749700, 936396, 1043290, 1155960, 1412040, 1708476, 2049300, 2438736, 2881200, 3381300, 3499335, 3943836, 4573800,…

Esta sucesión no se había publicado, y lo hicimos en http://oeis.org/A253650

Si los escribimos en columna podemos desarrollarlos como producto del tipo deseado:



Antes de seguir adelante, hay que advertir que estas descomposiciones en producto no tienen que ser únicas. Por ejemplo, 2881200 admite estos dos productos: 3*960400 y 300*9604, ambos formados por un triangular y un cuadrado. Por eso en la tabla se puede tener la falsa idea de que el factor triangular es más pequeño que el cuadrado. No es así, sino que hemos detenido la búsqueda al encontrar un ejemplo.

Una curiosidad

Los números triangulares de esta sucesión, tendrán, como todos la forma T(P)=P(P+1)/2. Pues bien, ni P ni P+1 pueden ser primos, porque si se descompone en un producto de un triangular y un cuadrado, tendríamos P(P+1)=k(k+1)m2 y si P o P+1 fueran primos, tendrían que dividir a k o a k+1 o a m, y los tres son menores que P. Estúdialo. Lo vemos en una tabla con los primeros casos:



Subconjunto interesante

Los números triangulares de orden k2-1, siendo k impar y mayor que 1, pertenecen a esta sucesión. En efecto, si k es impar tendrá la forma 2m+1, luego el orden del triangular será

 (2m+1)2-1=4m2+4m.

El triangular formado sobre él tendrá la forma (((2m+1)2*(4m2+4m)/2 y se puede descomponer en un cuadrado y un triangular: 4((2m+1)2 * m(m+1)/2.

Otra forma de expresarlo es que son triangulares construidos sobre 8 veces otro número triangular, ya que 4m2+4m=8*(m(m+1)/2). Estos números los tienes en http://oeis.org/A185096

Casi todos los elementos de la sucesión tienen esta forma: 300, 1176, 3240, 7260, 14196, 25200,… y pertenecen a la sucesión http://oeis.org/A083374 pero otros no, como 29403, 1043290 y 3499335.

Triangulares que son producto de un triangular y un primo

En esta búsqueda nos organizaremos de forma similar a la precedente pero al llegar a investigar si (i/k) es cuadrado o triangular, deberemos sustituirlo por la pregunta de si es primo. Esta es más difícil de responder, pero disponemos de la función isprime en PARI y de esprimo en nuestra hoja Conjeturas (http://www.hojamat.es/sindecimales/divisibilidad/herramientas/herrdiv.htm#global)

Puedes adaptar el algoritmo usado en la búsqueda anterior cambiando esas funciones: Observa que ahora, con lo explicado antes, las operaciones se entienden mejor, por lo que omitimos los comentarios y usamos color rojo en las novedades. Para hoja de cálculo podría servir esta rutina de Visual Basic:

Sub productriang()
Dim i, j, k, p, c


i = 3: j = 3
While i <= 10 ^ 3
k = 3: p = 3: c = 0
While c = 0 And p < i
If i / k = i \ k And esprimo(i / k) And i / k > 1 Then c = k
If c <> 0 Then MsgBox (i)
k = k + p: p = p + 1
Wend
i = i + j: j = j + 1
Wend
End Sub

En PARI

{i=3;j=3;while(i<=10^6,k=3;p=3;c=0;while(k<i&&c==0,if(i/k==i\k&&isprime(i/k)&&i/k>1,c=k);if(c>0,print1(i,", "));k+=p;p+=1);i+=j;j+=1)}

Los triangulares obtenidos son, en este caso, los siguientes:

6, 15, 21, 45, 66, 78, 105, 190, 210, 231, 435, 465, 630, 861, 903, 1035, 1326, 2415, 2556, 2628, 3003, 3570, 4005, 4950, 5460, 5565, 5995, 7140, 8646, 8778, 9870, 12246, 16471, 16836, 17205, 17391, 17766, 20100, 22155, 26565, 26796, 28680, 28920, 30381, 32131, 33411, 33930, 36856, 40755,… (los publicamos en http://oeis.org/A253651)

Podíamos pensar que, al ser el segundo factor primo, el primero será el máximo divisor triangular propio que tiene el número que se descompone
 (ver http://hojaynumeros.blogspot.com.es/2013/02/de-los-triangulares-alojados-los-primos.html).

Esto es así en la gran mayoría de casos. Así, 4005 se descompone como 45*89, siendo 45 es el máximo divisor triangular propio de 4005 y 89 un factor primo. Sin embargo, en el número 3570, su divisor triangular máximo es 595, que daría lugar al producto 3570=595*6, y el 6 no es primo. El producto válido sería 3570=210*17, que sí es primo.

Triangulares producto de triangular y oblongo

Estos son los primeros:

6, 36, 120, 210, 300, 630, 1176, 2016, 3240, 3570, 4950, 7140, 7260, 10296, 14196, 19110, 23436, 25200, 32640, 39060, 41616, 52326, 61776, 64980, 79800, 97020, 116886, 139656, 145530, 165600, 195000, 228150, 242556, 265356, 304590, 306936, 349866, 353220, 404550, 426426, …

Si recuerdas la caracterización de los oblongos en la entrada anterior, entenderás este código en PARI:

{i=3;j=3;
while(i<=10^6,k=2;p=4;c=0;
while(k<i&&c==0,if(i/k==i\k&&issquare(4*(i/k)+1)&&i/k>1,c=k);
if(c>0,write1("final.txt",i,", "));
k+=p;p+=2);
i+=j;j+=1)
}

Destacamos que el doble de cualquiera de estos números tiene la forma N=P(P+1)Q(Q+1).  Así, 2*2016=4032=7*8*8*9. Al contrario, como uno de los factores es oblongo, el producto será par y se podrá dividir entre dos y su mitad será producto de dos triangulares.

Todos los términos no nulos de la sucesión A083374 (6, 36, 120, 300, 630, 1176, 2016, 3240,…) pertenecen a esta, pues si tienen la forma n2(n2-1)/2 se pueden descomponer en el oblongo n(n+1) y el triangular n(n-1)/2, o bien el oblongo n(n-1) y el triangular n(n+1)/2.

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

Triangulares producto de un cuadrado y un primo

Al llegar a este punto simplificamos la exposición, ya que tendrás una buena idea de cómo se generan. Los triangulares que se forman multiplicando un cuadrado y un primo son estos:

28, 45, 153, 171, 300, 325, 496, 2556, 2628, 3321, 4753, 4851, 7381, 8128, 13203, 19900, 25200, 25425, 29161, 29403, 56953, 64980, 65341, 101025, 166753, 195625, 209628, 320400, 354061, 388521, 389403, 468028, 662976, 664128, 749700, 750925, 780625, 781875, 936396,…

Los números perfectos lo cumplen

Entre los números encontrados figuran los perfectos 28, 496, 8128,… https://oeis.org/A000396
La razón es que estos números tienen la forma 2k-1(2k-1)= 2k(2k-1)/2, y son, por tanto, triangulares. Además, como (2k-1) es primo de Mersenne en esta expresión, k ha de ser impar, y k-1, par, con lo que 2k-1 es un cuadrado y (2k-1) un primo, cumpliéndose así la condición.

Para no cansar más con este tema, sólo incluimos el código PARI con las novedades en rojo:

PARI
{i=3;j=3;
while(i<=10^6,k=4;p=5;c=0;
while(k<i&&c==0,if(i/k==i\k&&isprime(i/k)&&i/k>1,c=k);
if(c>0,write1("final.txt",i,", "));
k+=p;p+=2);
i+=j;j+=1)
}

Pues ya está bien. Ahora te toca a ti. ¿Sabrías prolongar las siguientes sucesiones con las técnicas que hemos desarrollado en las dos entradas?:

Triangulares producto de un cuadrado y un oblongo: 120, 300, 378, 528, 990, 1176, 2016,…

Y por último, los triangulares producto de un oblongo y un primo: 6, 10, 36, 66, 78, 210, 276, …