Find the first 10-digit prime in the decimal expansion of 17π.
The first 5 digits in the decimal expansion of π are 14159. The first 4-digit prime in the decimal expansion of π is 4159. The first 10-digit prime in the expansion of e is 7427466391. You are asked to find the first 10-digit prime in the decimal expansion of 17π. Write unit tests for each of these three functions.
Solution :
According to the prime number theorem , the proportion of prime numbers in all 10-digit numbers are approximately 4.5%. This is because:
\[\frac {1} {\ln(10^9)} \approx \frac {1} {\ln(10^{10})} \approx 4.5\%\]
As the decimal expansion of π comprises random numbers, we would expect a probability of observing at least one prime number in the first 100 digits of the decimal part of π of \(Pr = 1- (1 - 4.5\%)^{100} \approx 99\%\). Therefore, the first prime number would likely occur in the first 100 digits of the decimal part of π.
As indicated above, we first need to write a function to generate π with high precision. This goal can be achieved with the library sympy
:
#! python3 -m pip install --quiet sympy
def value_generator ( exp , digits ):
'''This function is used to generate n-digits pi or Euler's number e by evaluating
argument "exp".
PARAMETER:
exp: string, a math expression to be evaluated
digits: integer, how many digits of pi that you wanted.
RETURN:
A pi number of class decimal.
'''
import sympy as sp
value = sp . N ( exp , digits )
return ( value )
With sympy
, we can evaluate the value of a mathematical expression to any digit we want:
def unit_tests_1 ():
'''This function is to test if other functions go well
RETURN:
If all functions are on the right track, return 0.
'''
test1 = str ( value_generator ( "pi" , 1 + 2 )) #"3.14", 1+2 means 1 digit for integer part and 2 digits for the decimal part.
test2 = str ( value_generator ( "E" , 1 + 2 )) #"2.72"
test3 = str ( value_generator ( "17 * pi" , 2 + 2 )) #"53.41"
assert test1 == "3.14"
assert test2 == "2.72"
assert test3 == "53.41"
return 0
unit_tests_1 () #0, no problem
Thereafter, we will need some functions to examine if a number is prime or not. According to the definition of prime numbers, a prime number greater than 1 needs to be co-prime with all the prime numbers less than its square root. This is because if \(x= a \times b\), either \(a\) or \(b\) would be smaller than \(\sqrt{x}\). If \(x\) can not be divided by a number \(a \leq \sqrt{x}\), it would neither be divided by \(b > \sqrt{x}\). Thus, I write one function to generate the list of these prime numbers less than a given limit and one to check if the current number is prime.
def prime_list ( limit ):
'''This function returns a list of primes smaller than a given limit
PARAMETER:
number: an integer required to be checked if prime
RETURN:
prime_table: a list of prime numbers
'''
prime_table = [ 2 , 3 , 5 , 7 ] #a list to store the prime numbers smaller than the limit
for i in range ( 3 , - 1 , - 1 ):
if prime_table [ i ] > limit :
prime_table . pop ()
else :
break
if limit <= 10 :
return prime_table
#start from here, the limit should >10
for i in range ( 11 , limit + 1 ):
is_composite = False
for primes in prime_table : # note that the prime_table may be appended after each i loop
if i % primes == 0 : #i.e. i is composite
is_composite = True
break #break the "primes" loop
if not is_composite : prime_table . append ( i ) #if i is prime, add it to the prime_table
pass #the next "i" loop
return prime_table
def prime_checker ( number , prime_table = None ):
'''This function is to check if a given number is coprime with the number in the prime_table. If
prime_table is not given, it returns whether the number is prime or not.
PARAMETER:
number: an integer required to be checked if prime
prime_table: a list of primes given to check if the number is co-prime with them.
Using this argument with caution: all the elements should be prime numbers less
than "number"
RETURN:
flag: boolean, if the number is prime, return TRUE, otherwise FALSE
'''
import math
#exception handling
if not isinstance ( number , int ):
coprime_flag = False
return coprime_flag
if number <= 1 :
coprime_flag = False
return coprime_flag
if number == 2 :
#exclude 2, because ceil(sqrt(number)) will equal to itself
coprime_flag = True
return coprime_flag
if prime_table is None :
limit = math . ceil ( math . sqrt ( number )) # we just need to check if number can be divided by
#the prime numbers less than its square root
prime_table = prime_list ( limit )
coprime_flag = True
for primes in prime_table :
if number % primes == 0 : # that is "i is composite"
coprime_flag = False
break #break the "primes" loop
return coprime_flag
#the end of coprime_checker
The results of the two functions can be like:
def unit_tests_2 ():
'''This function is to test if other functions go well
RETURN:
If all functions are on the right track, return 0.
'''
test4 = prime_list ( 5 )
test5 = prime_list ( 10 )
test6 = prime_list ( 11 )
test7 = len ( prime_list ( 100 ))
assert test4 == [ 2 , 3 , 5 ]
assert test5 == [ 2 , 3 , 5 , 7 ]
assert test6 == [ 2 , 3 , 5 , 7 , 11 ]
assert test7 == 25 #should be 25 (means 25 primes less than 100)
pass
test8 = prime_checker ( 2 ) #to check if work well for small numbers
test8_1 = prime_checker ( 1 )
test9 = prime_checker ( 3 )
test10 = prime_checker ( 10 )
test11 = prime_checker ( 14159 ) #to check if work well for large numbers
test12 = prime_checker ( 14159 , prime_table = prime_list ( 1000 )) #when provided a prime list
test13 = prime_checker ( 14160 , prime_table = prime_list ( 1000 ))
assert test8 == True
assert test8_1 == False
assert test9 == True
assert test10 == False
assert test11 == True
assert test12 == True
assert test13 == False
pass
return 0
unit_tests_2 () #0, no problem
With the help of these functions, we finally need to write a wrapper function that integrates these parts and create a sliding window along the decimal expansion of π.
def find_first_prime ( exp , width , digit_limit ):
'''This function is a wrapper function which can find the first prime in the decimal
expansion of a math expression "exp"
PARAMETERs:
exp: string, a math expression to be evaluated, e.g. "17 * pi"
width: integer, the digits of prime that we are looking for, e.g. 10
digit_limit: integer, the function will be stopped if reaching a digit limit of
digit_limit, e.g., digit_limit = 100 means the function will stopped
at the 100th decimal digit if a prime number has not been found.
RETURN:
number: the first prime number in the decimal expansion of "exp". If None, meaning
no prime number has been found before reaching "digit_limit"
'''
import sympy
import math
value = value_generator ( exp , digit_limit + width + 10 )
#plus 10 to take the integer part into account and avoid round errors
dec_part = str ( value % 1 )
limit = math . ceil ( math . sqrt ( 10 ** width ))
prime_table = prime_list ( limit ) #to get all the prime numbers less than sqrt(10^width)
for i in range ( 2 , digit_limit + 2 ):
#start from the 3rd position of the value_dec string, right after "0."
#while Cliburn asks to create a function to generate this sliding windows,
#I feel it makes more sense to integrate it in the main function
#because it will make the program simpler and more logical.
number = int ( dec_part [ i :( i + width )])
flag = prime_checker ( number = number , prime_table = prime_table )
if flag : break
if ( flag ):
return ( number )
else :
return ( None )
Tested this function with some known examples:
#Here comes the unit tests
def unit_tests_3 ():
'''This function is to test if other functions go well
RETURN:
If all functions are on the right track, return 0.
'''
test14 = find_first_prime ( exp = "E" , width = 1 , digit_limit = 5 )
#check if the first digit is included
test15 = find_first_prime ( "pi" , 1 , 3 )
#check if stops at a desired digit_limit
test15_1 = find_first_prime ( "pi" , 1 , 4 )
test16 = find_first_prime ( "pi" , 5 , 10 ) #check the given examples
test17 = find_first_prime ( "pi" , 4 , 10 )
test18 = find_first_prime ( "E" , 10 , 100 )
assert test14 == 7
assert test15 == None
assert test15_1 == 5
assert test16 == 14159
assert test17 == 4159
assert test18 == 7427466391
pass
return 0
unit_tests_3 () #0, no problem
And finally solve the first prime in the deciaml expansion of 17π:
#Here comes the final result
print ( "The first prime number is" , find_first_prime ( exp = "17 * pi" , width = 10 , digit_limit = 100 ) )
#with digit_limit = 100, we have >98% probability to get a prime number
#The first prime number is 8649375157
Reference :
Prime number theorem. https://en.wikipedia.org/wiki/Prime_number_theorem.