5. Primes as sums of squares

VCMNA221: level 6: Design algorithms involving branching and iteration to solve specific classes of mathematical problems

  • implementing algorithms such as the Euclidean division algorithm

5.1. Primes 1 mod 4

If a prime number is divided by 4 and its remainder is 1, then it can be written as the sum of 2 squares.
Steps:

  • Define the upper limit of the prime numbers to check

  • Generate the primes using the sieve_of_eratosthenes function

  • Iterate over the primes and check if they are 1 mod 4 (have a remainder of 1 when divided by 4)

  • Find the sum of squares for each prime:

  • Print the list of primes with the sums of squares.

 1# Import math library
 2import math
 3
 4
 5def is_one_mod_four(n):
 6    # Define a function to check if a number is 1 mod 4
 7    return n % 4 == 1
 8
 9
10def is_sum_of_squares(n):
11    # Define a function to check if a number is a sum of two squares
12    # Iterate over the possible square roots from 1 to sqrt(n)
13    for i in range(1, int(math.sqrt(n)) + 1):
14        # Check if the difference between n and i^2 is also a perfect square
15        j = math.sqrt(n - i**2)
16        if j == int(j):
17            # Return the two square roots as a tuple
18            return (i, int(j))
19    # Return None if no such pair exists
20    return None
21
22
23def sieve_of_eratosthenes(n):
24    # Define a function to generate primes using the Sieve of Eratosthenes algorithm
25    # Create a boolean array of size n+1, initialized to True
26    is_prime = [True] * (n + 1)
27    # Mark 0 and 1 as not prime
28    is_prime[0] = is_prime[1] = False
29    # Iterate over the numbers from 2 to the square root of n
30    for i in range(2, int(math.sqrt(n)) + 1):
31        # If i is prime, mark all its multiples as not prime
32        if is_prime[i]:
33            for j in range(i * i, n + 1, i):
34                is_prime[j] = False
35    # Return a list of primes by filtering the array
36    return [i for i in range(2, n + 1) if is_prime[i]]
37
38
39# Define the upper limit of the range
40n = 100
41
42# Generate the primes using the sieve function
43primes = sieve_of_eratosthenes(n)
44
45# Iterate over the primes and check if they are 1 mod 4 and sum of squares
46for p in primes:
47    if is_one_mod_four(p):
48        result = is_sum_of_squares(p)
49        if result:
50            # Print the prime and the two squares
51            print(f"{p} = {result[0]}^2 + {result[1]}^2 = {result[0]**2} + {result[1]**2}")
Sample output is:
5 = 1^2 + 2^2 = 1 + 4
13 = 2^2 + 3^2 = 4 + 9
17 = 1^2 + 4^2 = 1 + 16
29 = 2^2 + 5^2 = 4 + 25
37 = 1^2 + 6^2 = 1 + 36
41 = 4^2 + 5^2 = 16 + 25
53 = 2^2 + 7^2 = 4 + 49
61 = 5^2 + 6^2 = 25 + 36
73 = 3^2 + 8^2 = 9 + 64
89 = 5^2 + 8^2 = 25 + 64
97 = 4^2 + 9^2 = 16 + 81