84 lines
1.7 KiB
Python
84 lines
1.7 KiB
Python
import os
|
|
import math
|
|
import numpy as np
|
|
|
|
|
|
# Find the sum of all prime numbers below two million
|
|
|
|
def isDivisible(n, ms):
|
|
res = False
|
|
for m in ms:
|
|
if n % m == 0:
|
|
res = True
|
|
break
|
|
return res
|
|
|
|
|
|
def getNextPrime(ps):
|
|
if len(ps) == 0:
|
|
return [2]
|
|
else:
|
|
p = ps[-1] + 1
|
|
while isDivisible(p, [q for q in ps if q <= np.sqrt(p)]):
|
|
p = p + 1
|
|
ps.append(p)
|
|
return ps
|
|
|
|
def getFirstNPrimesBad(n):
|
|
assert n >= 0
|
|
|
|
l = 0
|
|
ps = []
|
|
|
|
while l < n:
|
|
ps = getNextPrime(ps)
|
|
l = l + 1
|
|
|
|
return ps
|
|
|
|
def sieve(n):
|
|
assert n > 1
|
|
|
|
ns = [True] * n
|
|
|
|
for i in range(2, math.ceil(np.sqrt(n))):
|
|
if ns[i]:
|
|
j = pow(i, 2)
|
|
|
|
while j < n:
|
|
ns[j] = False
|
|
j = j + i
|
|
|
|
return [i for i,val in enumerate(ns) if val][2:]
|
|
|
|
|
|
def product(xs):
|
|
res = 1
|
|
|
|
for x in xs:
|
|
res = res * x
|
|
|
|
return res
|
|
|
|
|
|
# Found a formula for this exact question:
|
|
# https://oeis.org/A007504
|
|
# a(n) = (n^2/2)*(log n + log log n - 3/2 + (log log n - 3)/log n + (2 (log log n)^2 - 14 log log n + 27)/(4 log^2 n) + O((log log n/log n)^3))
|
|
|
|
# def getSumOfPrimes(n):
|
|
# return 0.5 * pow(n, 2) * (np.log(n) + np.log(np.log(n)) - 1.5 + (np.log(np.log(n)) - 3) / np.log(n) +
|
|
# (2 * pow(np.log(np.log(n)), 2) - 14 * np.log(np.log(n)) + 27) / (4 * np.log2(n)) + some O term)
|
|
|
|
# Unusable because of the o term, unfortunately. Also, I don't know exactly for which n prime numbers I want to calculate this sum.
|
|
# We will have to do it the hard way
|
|
|
|
|
|
def main():
|
|
print("Hello, this is Patrick")
|
|
|
|
ps = sieve(2000000)
|
|
print(sum(ps))
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main() |