Introduction:
This is a continuing blog post for any Project Euler problems I've solved that I found interesting (or at least difficult!). Most of my solutions are written in Python, but that will start changing as I start learning more languages.
Problem 357:
Problem Statement:
Consider the divisors of 30: 1,2,3,5,6,10,15,30.
It can be seen that for every divisor d of 30, d+30/d is prime.
Find the sum of all positive integers n not exceeding 100 000 000
such that for every divisor d of n, d+n/d is prime.
A couple of inital observations are that any n that fits the description is not going to be: odd or contain squares in its prime factorization. For any odd n with one of its divisors d, the experession d+n/d evaluates to an even (not prime) value since d and n/d are both odd. If an n has a square p in its prime factorization then we end up with p+n/p which is divisible by p and therefore not prime.
This however doesn't really narrow it down much considering the size of the sum, plus the mobius function to identify non square functions isn't fast for the amount of even values we have to consider. Instead we make the identification that every number will generate the term 1+n/1 for divisor 1 and therefore only numbers one less than a prime should be considered. This does reduce our potential candidates by a lot but producing primes up to 100 000 000 is still expensive.
My last observation is that any candidate cannot be equivalent to 4 or 6 modulo 10. 4mod10 implies a 2^2 in the prime facotization and 6mod10 implies that for some d, d+n/d is equivalent to 5mod10 which is not prime. There are some 8mod10 numbers that satisfy the constraints since a non-square number equivalent to 8mod10 doesn't immedietly produce a contradiction.
#This code is written by Rohit Banerjee (2020)
#In reference to Project Euler problem 357
import sympy
def pgen(n,show=False):
divisor_list = sympy.divisors(n)
for i in divisor_list[:int(len(divisor_list)/2)+1]:
if show:
print("%d+%d=%d"%(i,int(n/i),i+int(n/i)))
if not sympy.isprime(i+int(n/i)):
return(False)
return(True)
def main(n):
"""Cases below 10 are added manually, after that we first make sure the number
isn't =4(mod10)(implies a 2^2 in prime factorization)
or =6(mod10)(implies there exists a d+(n/d) that is divisible
by 5). Then we screen for numbers for which the mobius function
returns 1 or -1 so we don't check numbers with non-distinct
prime factors. Any numbers left are manually checked by pgen(i-1)"""
output = 1 + 2 + 6
for i in sympy.primerange(8,n):
if (i-1)%10 in [2,8,0]:
if abs(sympy.mobius(i-1)):
if pgen(i-1):
output += (i-1)
return(output)
def brute_force(n):
"""This function was used for testing purposes and is
the most naive but surefire way to get the solution sum"""
output = 0
for i in range(1,n+1):
if pgen(i):
output+=i
return(output)
def problem():
"""Run this to get the solution sum as well as computation time
in seconds"""
start = time.time()
print(main(100000000))
end = time.time()
print(end-start)
This algorithm took about 7 minutes (which is by all means horrible) but it was the first >100 problem I'd solved so I was happy.