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, …


lunes, 23 de febrero de 2015

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


Esta entrada y las siguientes tienen el doble objetivo de presentar unas curiosidades numéricas (algo intrascendentes) y analizar cómo organizar búsquedas de cierto tipo intentando dar con el algoritmo más rápido posible, ya que llegan fácilmente al orden de 10^7. Comenzamos con un ejemplo:

Números triangulares que son producto de triangulares

Muchos números triangulares son  producto de otros dos también triangulares. Por ejemplo, 45=15*3, 210=21*10, todos triangulares. Los tienes publicados en https://oeis.org/A188630

36, 45, 210, 630, 780, 990, 1540, 2850, 3570, 4095, 4851, 8778, 11781, 15400, 17955, 19110, 21528, 25200,…

Esta búsqueda está resuelta, pero imagina que la deseamos reproducir. No es fácil, porque para cada número natural deberíamos buscar lo siguiente:


  • Ver si ese número N es triangular
  • En caso afirmativo, recorrer todos sus divisores d.
  • Para cada uno de ellos, investigar si tanto d como N/d son ambos triangulares, y en caso afirmativo, parar la búsqueda para ese valor N y proseguir con N+1.

Es fácil darse cuenta de que se perderá mucho tiempo recorriendo números de uno en uno, que no son triangulares o bien que no poseen muchos divisores triangulares (o ninguno). Con  búsquedas de ese tipo, llamemos “ingenuas”, nuestro ordenador se pasaba minutos y minutos cuando llegaba a números grandes. Una solución es encaminar la búsqueda para que, hasta donde sea posible, sólo se recorran los números de cierto tipo. En el caso de triangulares, cuadrados, oblongos o primos, es posible realizar ese filtro. Lo concretamos:

Generación de triangulares

Los números que usaremos, salvo los primos, se pueden engendrar a partir de los precedentes.

Comenzaremos en esta entrada por explicar distintas formas de generar algunos tipos de números, y así ya las tenemos preparadas para cuestiones posteriores.

En el caso de los triangulares manejamos dos variables: Número N e incremento D. Comenzamos haciendo N=1 (primer triangular) e incremento D=2, para que 1+2 genere el siguiente triangular, el 3. Luego, en cada paso, N se convierte en N+D y D en D+1. Así funciona, ya que los números triangulares se forman añadiendo una fila con un elemento más.


Observa estos valores, generados con hoja de cálculo:



Recordamos: Los triangulares se generan tomando incrementos con una unidad más en cada paso. Ya utilizaremos esto más adelante.

Generación de cuadrados

Aquí tenemos dos soluciones, ambas prácticas, según el contexto, para recorrer sólo números cuadrados: la primera es trivial, declarar una variable k y luego usar k^2 en los cálculos. Está bien, pero a veces es lenta y no admite con facilidad ciertas acotaciones. La segunda es similar a la de los triangulares, pero incrementando D en D+2, en dos unidades en lugar de en una.


Observa el esquema:



Para quienes conozcáis estas propiedades, esto parecerá trivial, pero no está mal recordarlo, porque más adelante dará velocidad a nuestras búsquedas.

Generación de oblongos

Un número es oblongo cuando tiene la forma N=k(k+1), es decir, doble de un triangular. Es fácil ver que el siguiente oblongo será (k+1)(k+2). Su diferencia (k+1)(k+2)-k(k+1) = 2(k+1), es decir, el doble del mayor en el producto, luego el incremento adecuado será par y crecerá de 2 en 2. Esto nos permite generar oblongos: comenzamos con N=2 y D=4, Así generamos los siguientes: N=2+2*2=6, N=6+2*3=12, N=12+2*4=20,…También podemos declarar una variable y después trabajar con k*(k+1). Así hemos procedido en nuestra primera búsqueda con hojas de cálculo.



Así que, mientras los cuadrados se generan sumando impares, los oblongos sumando pares (y los triangulares sumando todos)

Caracterización de estos números

Necesitaremos también en las búsquedas que emprenderemos una forma de caracterizar estos números, cómo saber si un resultado es cuadrado u oblongo, por ejemplo. Aunque es sencillo y conocido, lo recordamos aquí:

Para saber si un número natural es cuadrado, la mejor prueba es que la parte entera de su raíz cuadrada, elevada a su vez al cuadrado, nos dé como resultado el número primitivo. En hoja de cálculo:

Public Function escuad(n) As Boolean
If n < 0 Then
escuad = False
Else
If n = Int(Sqr(n)) ^ 2 Then escuad = True Else escuad = False
End If
End Function

Si trabajas con las celdas, sin macros, el procedimiento es el mismo, pero con distinto lenguaje. Si tienes, por ejemplo, en la celda C12, un número del que deseas saber si es cuadrado, escribe en otra celda esta fórmula: =SI((ENTERO(RAIZ(C12)))^2=C12;1;0), y te devolverá un 1 si es cuadrado y un 0 si no lo es.

La caracterización de un número como triangular se basa en lo anterior, ya que ser triangular N es equivalente a que 8*N+1 sea cuadrado. Por tanto, podemos definir esta función:

Function estriangular(n) As Boolean
If escuad(8 * n + 1) Then estriangular = True Else estriangular = False
End Function

En lenguaje de celdas sería algo más complejo:
=SI((ENTERO(RAIZ(C12*8+1)))^2=C12*8+1;1;0),
Por último, los oblongos, al ser doble de triangulares, se descubren fácilmente:

Public Function esoblongo(n) As Boolean
If escuad(4 * n + 1) Then esoblongo = True Else esoblongo = False
End Function

Sin macros, =SI((ENTERO(RAIZ(C12*4+1)))^2=C12*4+1;1;0)

Algoritmos de búsqueda

Ya podemos pasar a nuestras búsquedas. Comenzaremos generando la sucesión https://oeis.org/A188630, con la que comenzamos esta entrada: Números triangulares que son producto de otros dos triangulares.

Nos organizaremos así:

Iniciamos triangulares: N=3, D=3 (Comenzamos por el 3)
1 Mientras no lleguemos al tope que nos hayamos marcado
    Iniciamos otros triangulares para ver si son divisores: K=3, P=3 
      2 Mientras no lleguemos a N ni encontremos un producto de triangulares
         Para cada K vemos si:
         1.-Es divisor de N
         2.- Si el cociente N/k es triangular
         Si cumple ambas condiciones, cerramos la búsqueda para N e imprimimos.
         Si no, generamos otro triangular convirtiendo K en K+P y P en P+1
     Fin de 2 Mientras
  Generamos otro triangular convirtiendo N en N+D y D en D+1
Fin del 1 Mientras

 Hemos comenzado los triangulares en el 3 para evitar trivialidades.

Para quienes no manejen mucho los algoritmos puede resultar complicado, pero hay que repasar hasta entenderlo. Se puede traducir al Visual Basic de Excel:

Sub productriang()
Dim i, j, k, p, c
i = 3: j = 3 Iniciamos la búsqueda en 3, para eliminar trivialidades
While i <= 10 ^ 4
k = 3: p = 3: c = 0 También iniciamos los divisores en 3
While c = 0 And p < i
If i / k = i \ k And estriangular(i / k) And i / k > 1 Then c = k
En la línea anterior buscamos que sea divisor, cociente triangular y no trivial
If c <> 0 Then MsgBox (i) Si cumple todo, presentamos el resultado
k = k + p: p = p + 1 Generamos el siguiente divisor
Wend
i = i + j: j = j + 1 Generamos el siguiente triangular
Wend
End Sub

Elimina comentarios, copia el resto como rutina para Visual Basic y al ejecutar verás aparecer los valores 36, 45, 210,…

Si te apetece explorar, aquí tienes la versión para PARI

{i=3;j=3; Iniciamos la búsqueda en 3, para eliminar trivialidades
while(i<=10^4,k=3;p=3;c=0; También iniciamos los divisores en 3
while(k<i&&c==0,if(i/k==i\k&&ispolygonal(i/k,3)&&i/k>1,c=k);
En la línea anterior buscamos que sea divisor, cociente triangular y no trivial
if(c>0,print1(i,", ")); Si cumple todo, imprimimos
k+=p;p+=1); Generamos el siguiente divisor
i+=j;j+=1) Generamos el siguiente triangular
}

Igualmente, si eliminas los comentarios y ejecutas este código PARI verás reproducida la sucesión https://oeis.org/A188630



Nos hemos detenido mucho en la generación de algunos tipos de números, su caracterización  y en la estructura general del algoritmo de búsqueda, pero en las siguientes entradas nos servirá todo esto para entender mejor los procesos.

lunes, 16 de febrero de 2015

Números “Tribonacci”


Los números “tribonacci” son análogos a los de Fibonacci, pero generados mediante recurrencias de tercer orden homogéneas.  Existen muchas sucesiones con este nombre, según sean sus condiciones iniciales. Aquí comenzaremos con la contenida en

http://mathworld.wolfram.com/TribonacciNumber.html

pero podemos cambiar más tarde si surgen propiedades interesantes para su estudio con hoja de cálculo.

En estos números la fórmula de recurrencia posee todos sus coeficientes iguales a la unidad

xn= A*xn-1+B*xn-2+C*xn-3 se convertiría en

xn= xn-1+xn-2+xn-3

Al igual que en el caso de Fibonacci, los dos valores iniciales también valen 1, y el tercero, 2, pero ya hemos explicado que existen otras variantes. Dejamos los enlaces de algunas de ellas:

http://oeis.org/A000073  comienza con a(0)=a(1)=0, a(2)=1
http://oeis.org/A000213  con a(0)=a(1)=a(2)=1
http://oeis.org/A001590 con a(0)=0, a(1)=1, a(2)=0
http://oeis.org/A081172 comienza con 1,1,0.

Y hay más.

Como ya hemos indicado, nosotros comenzaremos con:

Condiciones iniciales: x0=1, x1=1x2=2 Ecuación de recurrencia: xn= xn-1+xn-2+xn-3

Los primeros términos son:

 1, 1, 2, 4, 7, 13, 24, 44, 81, 149, 274, 504, 927, 1705, 3136, 5768, 10609, 19513, 35890, 66012, 121415, 223317, 410744,… http://oeis.org/A000073

Como en otras entradas sobre el mismo tema, podemos acudir a nuestra herramienta de hoja de cálculo para sucesiones recurrentes

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

En la imagen puedes identificar los coeficientes y valores iniciales



Con el botón “Ver sucesión” podremos obtener el listado de estos números:


Ecuación característica

Al igual que en otras sucesiones recurrentes, su ecuación característica se formará a partir de sus coeficientes, en este caso todos iguales a 1, luego será x3-x2-x-1=0

Con nuestra herramienta podemos encontrar sus raíces:



La misma solución obtenemos con WxMaxima



Recordemos que los elementos de las sucesiones recurrentes se pueden expresar como suma de potencias de las tres soluciones, pero con estos números ocurre como con algunos similares (los de Fibonacci, Perrin o Narayana), y es que las raíces complejas, al tener módulo inferior a la unidad, tienden a cero si prolongamos la sucesión. Por ello, las potencias de la raíz real, 1,839286…generan con bastante aproximación los números Tribonacci, y, lo que es lo mismo, esta constante coincidirá aproximadamente con el cociente entre dos de estos números consecutivos. Lo vemos con hoja de cálculo:



Por ello, al número 1,839286…se le llama Constante Tribonacci.

Función generatriz

Todas las variantes de las sucesiones Tribonacci comparten los mismos coeficientes de recurrencia, y por tanto también el denominador de su función generatriz. La que estamos estudiando en esta entrada, de inicio 1, 1, 2, se genera con la siguiente:


Al igual que con otras sucesiones, la comprobaremos con PARI:

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

Te escribirá en un archivo sucesión.txt su desarrollo (este archivo lo deberás tener vacío en la misma carpeta que PARI), y aparecerán como coeficientes los términos de la sucesión Tribonacci:

x + x^2 + 2*x^3 + 4*x^4 + 7*x^5 + 13*x^6 + 24*x^7 + 44*x^8 + 81*x^9 + 149*x^10 + 274*x^11 + 504*x^12 + 927*x^13 + 1705*x^14 + 3136*x^15 + 5768*x^16 + 10609*x^17 + 19513*x^18 + 35890*x^19 + O(x^20)

Una excursión por la hoja de cálculo

Podemos usar la versión matricial de la generación de estos números para recordar algunos detalles sobre hojas de cálculo.

Es elemental comprobar que las ternas de números consecutivos de Tribonacci. T(n), T(n+1), T(n+2) pueden engendrar matricialmente la terna siguiente T(n+1), T(n+2), T(n+3), mediante la siguiente fórmula matricial:


Esta fórmula es adecuada para repasar las fórmulas matriciales de las hojas de cálculo. Comenzamos construyendo un esquema como el de la imagen:



Para efectuar el producto matrical deberemos usar la función MMULTI, con parámetros la primera matriz y la columna de la primera terna:

{=MMULT(D4:F6;H4:H6)}

Observa que como multiplicamos rangos de celdas, usamos el separador :

Para que la hoja entienda que se trata de una multiplicación matricial, cuando termines de escribir la fórmula, en lugar de terminar con INTRO, usaremos Ctrl+Mayúscula+INTRO. La aparición de las llaves es la señal de que la fórmula ha sido introducida correctamente.

Una vez efectuado el cálculo sobre una terna, basta que copies el resultado como dato, usando Copiar y Pegado especial como valores, y proseguirán apareciendo ternas nuevas.

Uno de los autovalores de la matriz que hemos usado es la constante de Tribonacci, 1,839286…La razón es que el polinomio característico de la matriz es el mismo que el de la ecuación característica de la recurrencia, x3-x2-x-1=0.

Curiosidades

En esta serie sobre sucesiones recurrentes solemos presentar en cada una de ellas propiedades curiosas, no todas las conocidas, que llenarían libros, sino las que más nos llamen la atención o se adapten mejor a las herramientas que usamos. Para la de Tribonacci presentaremos una propiedad combinatoria.

Particiones de un número en sumandos no mayores que 3

Los números de Tribonacci (salvo los iniciales) cumplen que T(N) coincide con las particiones de N-1 en sumandos que se pueden repetir, en cualquier orden y con los sumandos menores o iguales a 3. Por ejemplo, T(5)=7, que coincide con las particiones del número 4 en partes no superiores a 3:
Lo comprobamos con el listado obtenido con nuestra hoja no publicada “Cartesius”:



Observamos que resultan 7 particiones distintas.

Para T(4)=4 obtenemos el mismo resultado con particiones del número 3:


La razón de que esto funcione así es que cualquier partición de este tipo con N elementos ha resultado a adjuntar un 1 a las particiones de N-1, un 2 a las de N-2 y un 3 a las de N-3, con los que se cumple xn= xn-1+xn-2+xn-3 . Para que lo entiendas mejor hemos coloreado estos tres sumandos para el caso de T(6)=13: