In [1]:
numero = 78277485
In [2]:
a = numero % 1000; a
Out[2]:
In [3]:
b = numero % 100; b
Out[3]:
In [4]:
from mpmath import mp #Librería usada para modiuficar la precision decimal de PI
from math import sqrt
In [5]:
mp.dps = 1000 #Establecer la precisión que queremos
In [6]:
pi_1000 = mp.pi
In [7]:
print(pi_1000)
Convenimos que la posición decimal de las decimas está indicada con el 0.
In [8]:
besima = str(pi_1000)[b+2]
print(besima)
In [9]:
aesima = str(pi_1000)[a+2]
print(aesima)
In [10]:
#Generador de python que devuelve los divisores de un número n
def divisores(n):
if n < 1:
raise ValueError("El número ha de ser natural y mayor que 0")
div = [] #Complementos de numeros
for numero in range(1, int(sqrt(n) +1 )):
if n%numero == 0:
yield numero #Si es un divisor lo devolvemos
#En la siguiente iteración comprobamos si su cuadrado
#no es n. Simpre va a ser <= n ya
#que buscamos hasta la raiz del numero
if numero*numero != n:
#Si no es lo añadimos como divisor
div = div + [int(n/numero)]
for divisor in div:
yield divisor
In [11]:
list(divisores(a))
Out[11]:
In [12]:
# La suma de los divisores es
sum( list(divisores(a)) )
Out[12]:
In [13]:
# Un ejemplo para ver como sería con otro número
list( divisores( 24 ) )
Out[13]:
In [14]:
sum( list( divisores( 24 ) ) )
Out[14]:
In [15]:
primos = [2]
In [16]:
#Una lista para almacenar los primos que encontremos
def es_primo(n):
if n >= 2:
for y in range(2,int(sqrt(n))+1):
if not ( n % y ):
return False
else:
return False
return True
In [17]:
for i in range(1,920):
if es_primo(i):
primos.append(i)
In [18]:
print(primos)
Luego el mayor primo < a es
In [19]:
primos [-1]
Out[19]:
Vamos a calcular ahora el menor primo > que a. Para ello vamos a recurrir a itertools que es una librería del nucleo de python que nos proporciona funciones para hacer algo próximo a la programación funcional
In [20]:
from itertools import takewhile, count
In [21]:
[x for x in takewhile(lambda x: not es_primo(x), count(a+1))][-1]+1
Out[21]:
Luego el siguiente primo mayor que a es 937
Una manera de calcular el factorial es usar la librería que trae python.
In [22]:
from math import factorial
In [23]:
#factorial(1000+a)
Tambíen podemos hacernos nuestro propio factorial
In [24]:
from operator import mul
from functools import reduce
In [25]:
def factorial_mio(n):
return reduce(mul, range(1,n+1))
In [26]:
factorial(1000+a) == factorial_mio(1000+a)
Out[26]:
Vemos que da el mismo resultado
El numero de ceros en los que acaba el factorial es:
In [27]:
"""
Para buscar el número de 0 se ha procedido a convertir el número en una cadena de caracteres
invertir la posicion de los digitos, el último pasa a ser el primero y así con el resto de dígitos,
por último se recorre la cadena hasta que se encuentra un dígito distinto de 0
"""
total_ceros = 0
distinto_encontrado = False
numero = str( factorial_mio(1000+a) )[::-1]
while not distinto_encontrado:
if numero[total_ceros] == '0':
total_ceros += 1
else:
distinto_encontrado = True
In [28]:
total_ceros
Out[28]:
Ya tenemos una función para calcular los divisores de un número n
In [29]:
def divisores(n):
if n < 1:
raise ValueError("El número ha de ser natural y mayor que 0")
div = [] #Complementos de numeros
for numero in range(1, int(sqrt(n) +1 )):
if n%numero == 0:
yield numero #Si es un divisor lo devolvemos
#En la siguiente iteración comprobamos si su cuadrado
#no es n. Simpre va a ser <= n ya
#que buscamos hasta la raiz del numero
if numero*numero != n:
#Si no es lo añadimos como divisor
div = div + [int(n/numero)]
for divisor in div:
yield divisor
La suma de los divisores es (incluyendo al 1 y omitiendo el propio número):
In [30]:
def suma_divisores(n):
divis = list(divisores(n))
#Ahora quitamos el propio número
divis.remove(n)
return sum(divis)
In [31]:
suma_divisores(24)
Out[31]:
Los números perfectos entre los 1000 primeos naturales son:
In [32]:
[x for x in range(1,1001) if x == suma_divisores(x)]
Out[32]:
En total son
In [33]:
len([x for x in range(1,1001) if x == suma_divisores(x)])
Out[33]:
Demostremos que si $2^{p} -1 es primo necesariamente p es primo. Vamos a demostrarlo usando el contrarecíproco.
Supongamos que $p$ no es primo. Por tanto lo podemos expresar como producto de dos números. Sean $m,n \in \mathbb{N}$ tal que $2 \le m,n < p$ y $mn=p$. Tenemos que demostar que "Si $p$ no es primo entonces $2^{p}-1$ tampoco es primo".
Sea $p = mn$ para $m,n \in \mathbb{N}$ tales que $2 \le m,n < p$
$2^{mn}-1 = (2^m)^n -1 = (2^m -1 )((2^m)^{n-1} + (2^m)^{n-1} + ... + 2^m + 1)$
por nuestra hipótesis tenemos
$2 \le 2^m -1 < 2^p -1$
Así $2^m -1$ divide a $2^p -1$ y esto nos dice que $2^p -1$ es compuesto.
Hemos demostrado que si p no es primo entonces $2^p -1$ tampoco lo es. Podemos afirmar pues:
Si $2^p -1$ es primo, entonces $p$ es primo.
Sabiendo esto, primero buscamos los primos entre los primeros 257 numeros.
In [34]:
primos_en_258 = []
In [35]:
for i in range(1,258):
if es_primo(i):
primos_en_258.append(i)
In [36]:
print(primos_en_258)
Como comprobar si un primo de la forma $2^{p} -1$ es una tarea muy pesada cuando p crece vamos a usar el algoritmo de Lucas-Lehmer.
In [37]:
def LucasLehmer(p):
"""
Este algoritmo determina si 2**p-1 es primo para p >= 2
"""
if p == 2:
return True
s = 4
M = 2**p -1
for i in range(3, p+1):
s = (s**2 -2) % M
return s == 0
In [38]:
lista = []
for i in primos_en_258:
lista.append( (i, LucasLehmer(i) ) )
In [39]:
print(lista)
En esta lista tenemos tuplas de los primeros 258 números y si son primos de Mersenne. Los que lo cumplen son
In [40]:
print([x[0] for x in lista if x[1] == True ])
Luego Mersenne se había equivocado y la lista correcta es la anterior.
Vemos que no es cierto que para todo $p$ primo $2^{p} -1$ sea primo.
In [41]:
numero = 1000+2*a
print(numero)
creamos una lista de primos hasta 1000+2a
In [42]:
def lista_primos(n):
lista = []
for i in range(1,n+1):
if es_primo(i):
lista.append(i)
return lista
In [43]:
print(lista_primos(numero))
In [44]:
def suma_primos(numero):
primos = lista_primos(numero)
indice_i = 0
indice_m = len(primos) - 1
while indice_i <= indice_m:
suma = primos[indice_i] + primos[indice_m]
if suma < numero:
indice_i += 1
else:
if suma == numero:
yield primos[indice_i], primos[indice_m]
indice_m -= 1
In [45]:
list(suma_primos(numero))[0]
Out[45]:
In [46]:
print(list(suma_primos(numero)))
In [47]:
len(list(suma_primos(numero)))
Out[47]:
In [48]:
def golbach(n):
primos = lista_primos(n)
sumas = list(suma_primos(n))
return sumas[0]
In [49]:
golbach(12)
Out[49]:
In [50]:
def dgoldbach(n):
primos = lista_primos(n)
sumas = list(suma_primos(n))
return len(sumas)
In [51]:
dgoldbach(24)
Out[51]:
In [52]:
def ackermann(m, n):
if m == 0:
return n+1
if n == 0:
return ackermann(m-1, 1)
return ackermann(m-1, ackermann(m, n-1))
In [53]:
ackermann(3,2)
Out[53]:
Al intentar calcular el valor para ackermann(4,2) se desborda la pila por lo que hay que buscar otra manera de resolver el problema.
Propiedades de la función de Ackermann:
Así para el caso de ackerman(4,2) el resultado sería
In [54]:
numero = 2
n = 2
for i in range(n+2):
numero = numero**2
print(numero)
Así Akermann(4,2) = $2^{65536} - 3 $