# Pythagorean triples by factorizing

This simple C program finds primitive Pythagorean triples by factorizing the square of the lowest integer of the triple.

## C code

`/* Find primitive Pythagorean triples by factorizing`

`Version 1.1 (2014-01-25, revised 2019-03-04)`

`We want positive integers a <= b <= c such that a squared + b squared = c squared.`

`This is equivalent to a squared = c squared - b squared = (c-b) (c+b). That is, we want a squared to have factors (c-b) and (c+b). We want b and c to be positive integers with b < c, so we want (c+b) and (c-b) to be positive integers. So, for each positive integer a, we want pairs of factors of a squared.`

`We have (c-b) = (c+b) - 2b = 2c - (c+b). If (c+b) is even, (c-b) is the sum of two even numbers and the difference between two even numbers, so (c-b) is even.`

`We have (c+b) = (c-b) + 2b = 2c - (c-b). If (c-b) is even, (c+b) is the sum of two even numbers and the difference between two even numbers, so (c+b) is even.`

`So, if either of the factors (c+b) or (c-b) is even, both factors are even.`

`For a given a, a is either even or odd.`

`If a is even, a squared is even. So at least one of the factors is even. But we know that if one of the factors is even, both factors are even. So, if a is even, both factors are even. So, we can consider the factors ((c-b) / 2) and ((c+b) / 2) of (a / 2) squared.`

`If a is odd, a squared is odd. If a squared is odd, all the factors of a squared are odd. So, if a is odd, both factors are odd.`

`We have b = ((c+b) - (c-b)) / 2. We want be to be positive, so (c+b) > (c-b). So, the factors (c+b) and (c-b) must not be equal.`

`We have a < b = ((c+b) - (c-b)) / 2 = (((a squared) / (c-b)) - (c-b)) / 2.`

`So, 2a (c-b) < a squared - (c-b) squared.`

`So, 0 < a squared - 2 a (c-b) - (c-b) squared = ((a (sqrt (2) - 1) - (c-b)) (a (sqrt (2) + 1) + (c-b)).`

`(a (sqrt (2) + 1) + (c-b)) is positive, so, we need 0 < a (sqrt (2) - 1) - (c-b), or (c-b) < a (sqrt (2) - 1). This gives an upper bound for the factor (c-b).`

`The minimum factor is 1 so 1 <= (c-b) < a (sqrt (2) - 1. So, sqrt (2) + 1 <= a. So, a > 2. This gives a lower bound for a.`

`When we have found a value for (c-b), we can find the corresponding factor (c+b) = a squared / (c-b).`

`Similarly, when a is even and we have found a value for ((c-b) / 2), we can find the corresponding factor ((c+b) / 2) = (a / 2) squared / ((c-b) / 2).`

`When we have a suitable pair of factors (c-b) and (c+b), we can derive b and c to give the Pythagorean triple.`

`b = ((c+b) - (c-b)) / 2`

`c = ((c+b) + (c-b)) / 2`

`Similarly, when a is even and we have a suitable pair of factors ((c-b) / 2) and ((c+b) / 2), we can derive b and c to give the Pythagorean triple.`

`b = ((c+b) / 2) - ((c-b) / 2)`

`c = ((c+b) / 2) + ((c-b) / 2)`

`We only want primitive triples, so a, b, and c cannot have a common factor. If any two of a, b, and c have a common factor, the other one of them will have that common factor and so we need to make sure that no two of them share a common factor. No two of a, b, and c are equal and at most one of them is even. If at most one of a, b, and c is even, at least two of them are odd. If two of them are odd, the other must be even, because the square of an odd number is odd and the sum of two odd numbers is even and the difference between two odd numbers is even. So, only one of a, b, and c is even.`

`If a and b were odd and c was even, we could find integers i, j, and k such that a = 2i + 1, b = 2j + 1, and c = 2k. We could write a squared + b squared = c squared as (2i + 1) squared + (2j + 1) squared = (2k) squared, which we can write as 4 i squared + 4i + 1 + 4 j squared + 4 j + 1 = 4 k squared or 4 (i squared + j squared - k squared + i + J) + 2 = 0 or 2 (i squared + j squared - k squared + i + j) + 1 = 0. No integers i, j, and k can satisfy this equation, because the left-hand side is an odd number, but the right-hand side is even. So, c cannot be even.`

`So, in a primitive triple, one of a and b is even and c is always odd.`

`If the value a or b that was even was the product of 2 and an odd number, we could write it as 2(2i + 1) and we could write the other two values of the triple as 2j + 1 and 2k + 1. We could write a squared + b squared = c squared as 4 (4 i squared + 4 i + 1) + 4 j squared + 4 j + 1 = 4 k squared + 4 k + 1 or 4 (i squared + i) + 1 + j squared + j) - (k squared + k) = 0. Regardless of whether j or k are odd or even, j squared + j and k squared + k are always even. No integers i, j, and k can satisfy this equation, because the left-hand side is an odd number, but the right-hand side is even. So, the even integer a or b is not the product of 2 and an odd number. So, when a or b is even, it must be divisible by 4.`

`The factors (c-b) and (c+b) cannot share an odd prime factor, because if they did, b and c would share that common factor too.`

`If a is even, both the factors (c-b) and (c+b) are even. But, for c to be odd, one of the factors ((c-b) / 2) and ((c+b) / 2) is even and the other is odd. So, the factors ((c-b) / 2) and ((c+b) / 2) share no common prime factors. */`

`/* Standard header file for fprintf, fputs, printf, puts */`

`#include <stdio.h>`

`/* Standard header file for EXIT_FAILURE, EXIT_SUCCESS, free, malloc, NULL, realloc */`

`#include <stdlib.h>`

`/* Standard header file for strcmp */`

`#include <string.h>`

`/* Standard header file for sqrt */`

`#include <math.h>`

`/* Returns the position of a character in a character array */`

`unsigned char position_in_character_array (char const character)`

`{`

`char const character_array [] = "0123456789ABCDEFabcdef";`

`unsigned char position = 0;`

`while ((character_array [position] != character) && (character_array [position] != 0))`

`{`

`position ++;`

`}`

`return position;`

`}`

`/* Returns the maximum position of a character in a character array across all characters in a string */`

`unsigned char maximum_position_in_character_array (char const * const string)`

`{`

`unsigned char position = 0;`

`unsigned char maximum_position = 0;`

`while ((string [position] != 0) && (maximum_position < 22))`

`{`

`unsigned char character_array_position = position_in_character_array (string [position]);`

`if (character_array_position > maximum_position)`

`{`

`maximum_position = character_array_position;`

`}`

`position ++;`

`}`

`return maximum_position;`

`}`

`/* Sets the default base indicator */`

`unsigned char default_base_indicator (int const argument_count, char const * const argument_array [])`

`{`

`unsigned char base_indicator = 'u';`

`unsigned char const zero = 1 << 0;`

`unsigned char const base_8 = 1 << 1;`

`unsigned char const uppercase_base_16_with_prefix = 1 << 2;`

`unsigned char const lowercase_base_16_with_prefix = 1 <<3;`

`unsigned char const uppercase_base_16_without_prefix = 1 << 4;`

`unsigned char const lowercase_base_16_without_prefix = 1 << 5;`

`unsigned char const ambiguous_base = 1 << 6;`

`unsigned char const found_base = 1 << 7;`

`unsigned char status = 0;`

`int argument_index;`

`for (argument_index = 1; argument_index < argument_count; argument_index ++)`

`{`

`char const * const string = argument_array [argument_index];`

`unsigned long int position = 0;`

`unsigned char maximum_position;`

`if (string [position] == '-')`

`{`

`position ++;`

`}`

`if (string [position] == '0')`

`{`

`position ++;`

`if (string [position] == 0)`

`{`

`status |= (found_base | zero);`

`}`

`else if ((string [position] != 'O') && (string [position] != 'o') && (string [position] != 'X') && (string [position] != 'x'))`

`{`

`maximum_position = maximum_position_in_character_array (& (string [position]));`

`if (maximum_position <8)`

`{`

`status |= (found_base | base_8);`

`}`

`else if (maximum_position < 10)`

`{`

`status |= (found_base | ambiguous_base);`

`}`

`else if (maximum_position < 16)`

`{`

`status |= (found_base | uppercase_base_16_without_prefix);`

`}`

`else if (maximum_position < 22)`

`{`

`status |= (found_base | lowercase_base_16_without_prefix);`

`}`

`}`

`}`

`if ((status & found_base) == 0)`

`{`

`if ((string [position] == 'O') || (string [position] == 'o'))`

`{`

`position ++;`

`maximum_position = maximum_position_in_character_array (& (string [position]));`

`if (maximum_position <8)`

`{`

`status |= base_8;`

`}`

`}`

`else if ((string [position] == 'X') || (string [position] == 'x'))`

`{`

`position ++;`

`maximum_position = maximum_position_in_character_array (& (string [position]));`

`if (maximum_position < 16)`

`{`

`status |= uppercase_base_16_with_prefix;`

`}`

`else if (maximum_position < 22)`

`{`

`status |= lowercase_base_16_with_prefix;`

`}`

`}`

`else`

`{`

`maximum_position = maximum_position_in_character_array (& (string [position]));`

`if (maximum_position < 10)`

`{`

`status |= ambiguous_base;`

`}`

`else if (maximum_position < 16)`

`{`

`status |= uppercase_base_16_without_prefix;`

`}`

`else if (maximum_position < 22)`

`{`

`status |= lowercase_base_16_without_prefix;`

`}`

`}`

`}`

`else`

`{`

`status &= ~ found_base;`

`}`

`}`

`if ((status == base_8) || (status == (zero | base_8)))`

`{`

`base_indicator = 'o';`

`}`

`else if ((status & lowercase_base_16_without_prefix) != 0)`

`{`

`base_indicator = 'x';`

`}`

`else if ((status & uppercase_base_16_without_prefix) != 0)`

`{`

`base_indicator = 'X';`

`}`

`else if ((status == lowercase_base_16_with_prefix) || (status == (uppercase_base_16_with_prefix | lowercase_base_16_with_prefix)))`

`{`

`base_indicator = 'x';`

`}`

`else if (status == uppercase_base_16_with_prefix)`

`{`

`base_indicator = 'X';`

`}`

`else`

`{`

`base_indicator = 'u';`

`}`

`return base_indicator;`

`}`

`/* Returns the value corresponding to the position of the character in the character array */`

`unsigned char value (unsigned char position)`

`{`

`unsigned char return_value;`

`if (position > 15)`

`{`

`return_value = position - 6;`

`}`

`else`

`{`

`return_value = position;`

`}`

`return return_value;`

`}`

`/* Converts a string representing an integer into an integer */`

`unsigned long int string_as_integer (char const * const string, unsigned char const base_indicator, unsigned char * const error)`

`{`

`unsigned char const false = 0;`

`unsigned char const true = 1;`

`unsigned char negative = false;`

`unsigned long int position = 0;`

`unsigned char maximum_position;`

`unsigned char base = 0;`

`unsigned long int integer = 0;`

`* error = 0;`

`if (string [position] == '-')`

`{`

`negative = true;`

`position ++;`

`}`

`if (string [position] == '0')`

`{`

`position ++;`

`if (string [position] == 0)`

`{`

`base = 1;`

`}`

`else if ((string [position] != 'O') && (string [position] != 'o') && (string [position] != 'X') && (string [position] != 'x'))`

`{`

`maximum_position = maximum_position_in_character_array (& (string [position]));`

`if (maximum_position <8)`

`{`

`base = 8;`

`}`

`else if ((maximum_position < 10) && (base_indicator != 'X') && (base_indicator != 'x'))`

`{`

`base = 10;`

`}`

`else if (maximum_position < 22)`

`{`

`base = 16;`

`}`

`else`

`{`

`* error = 1;`

`}`

`}`

`}`

`if (base == 0)`

`{`

`if ((string [position] == 'O') || (string [position] == 'o'))`

`{`

`position ++;`

`maximum_position = maximum_position_in_character_array (& (string [position]));`

`if (maximum_position <8)`

`{`

`base = 8;`

`}`

`else`

`{`

`* error = 8;`

`}`

`}`

`else if ((string [position] == 'X') || (string [position] == 'x'))`

`{`

`position ++;`

`maximum_position = maximum_position_in_character_array (& (string [position]));`

`if (maximum_position < 22)`

`{`

`base = 16;`

`}`

`else`

`{`

`* error = 16;`

`}`

`}`

`else`

`{`

`maximum_position = maximum_position_in_character_array (& (string [position]));`

`if ((maximum_position < 22) && ((base_indicator == 'X') || (base_indicator == 'x')))`

`{`

`base = 16;`

`}`

`else if (maximum_position < 10)`

`{`

`base = 10;`

`}`

`else if (maximum_position < 22)`

`{`

`base = 16;`

`}`

`else`

`{`

`* error = 1;`

`}`

`}`

`}`

`if (base > 0)`

`{`

`while ((string [position] > 0) && (* error == 0))`

`{`

`unsigned long int previous_integer = integer;`

`if (integer > 0)`

`{`

`integer *= base;`

`}`

`integer += value (position_in_character_array (string [position]));`

`if ((previous_integer != 0) && (integer <= previous_integer))`

`{`

`* error = 255;`

`}`

`position ++;`

`}`

`if (negative == true)`

`{`

`integer = (- integer);`

`}`

`}`

`return integer;`

`}`

`int main (int const argument_count, char const * const argument_array [])`

`{`

`int return_value = EXIT_SUCCESS;`

`char const * const program_file_name = argument_array ;`

`unsigned char const option_help = 1 << 0;`

`unsigned char const option_version = 1 << 1;`

`char const * const option_string_array [] = { "?", "/?", "/H", "/h", "-?", "-H", "-h", "-help", "--help", "--version" };`

`unsigned char const option_value_array [] = { option_help, option_help, option_help, option_help, option_help, option_help, option_help, option_help, option_help, option_version };`

`unsigned char const option_string_count = sizeof (option_string_array) / sizeof (option_string_array );`

`unsigned char const option_value_count = sizeof (option_value_array) / sizeof (option_value_array );`

`unsigned char status = 0;`

`if (option_string_count != option_value_count)`

`{`

`unsigned char option_index;`

`fprintf (stderr, "%s: Internal error: The size of the option value array (%u) does not match the size of the option string array (%u)\n", program_file_name, option_value_count, option_string_count);`

`fputs ("index: value, string\n", stderr);`

`for (option_index = 0; (option_index < option_string_count) && (option_index < option_value_count); option_index ++)`

`{`

`fprintf (stderr, "%u: %u, %s\n", option_index, option_value_array [option_index], option_string_array [option_index]);`

`}`

`if (option_string_count > option_value_count)`

`{`

`for (option_index = option_value_count; option_index < option_string_count; option_index ++)`

`{`

`fprintf (stderr, "%u: [none], %s\n", option_index, option_string_array [option_index]);`

`}`

`}`

`else`

`{`

`for (option_index = option_string_count; option_index < option_value_count; option_index ++)`

`{`

`fprintf (stderr, "%u: %u, [none]\n", option_index, option_value_array [option_index]);`

`}`

`}`

`return_value = EXIT_FAILURE;`

`}`

`else`

`{`

`int argument_index;`

`/* Set options */`

`for (argument_index = 1; (argument_index < argument_count) && (status == 0); argument_index ++)`

`{`

`char const * const current_argument = argument_array [argument_index];`

`unsigned char const valid_option = 1 << 2;`

`unsigned char option_index;`

`status &= ~ valid_option;`

`for (option_index = 0; (option_index < option_string_count) && ((status & valid_option) == 0); option_index ++)`

`{`

`if (strcmp (option_string_array [option_index], current_argument) == 0)`

`{`

`status |= valid_option;`

`status |= option_value_array [option_index];`

`}`

`}`

`}`

`/* Print information */`

`if ((status & option_help) != 0)`

`{`

`printf ("Usage: %s [list of nonnegative integers separated by spaces]\n", program_file_name);`

`puts ("Print each primitive Pythagorean triple whose lowest integer is greater than or equal to the lowest integer in the list and less than or equal to the highest integer in the list.");`

`}`

`else if ((status & option_version) != 0)`

`{`

`puts ("Triples by factorizing 1.1");`

`}`

`else`

`{`

`/* Read the arguments that are integers */`

`unsigned char base_indicator = default_base_indicator (argument_count, argument_array);`

`unsigned int integer_count = 0;`

`unsigned long int minimum_a = 0;`

`unsigned long int maximum_a = 0;`

`for (argument_index = 1; argument_index < argument_count; argument_index ++)`

`{`

`unsigned char error = 0;`

`unsigned long int integer = string_as_integer (argument_array [argument_index], base_indicator, & error);`

`if (error == 0)`

`{`

`integer_count ++;`

`if (integer_count == 1)`

`{`

`minimum_a = integer;`

`maximum_a = integer;`

`}`

`else if (integer < minimum_a)`

`{`

`minimum_a = integer;`

`}`

`else if (integer > maximum_a)`

`{`

`maximum_a = integer;`

`}`

`}`

`else`

`{`

`char * message;`

`if (error == 1)`

`{`

`message = "Invalid integer";`

`}`

`else if (error == 8)`

`{`

`message = "Invalid octal integer";`

`}`

`else if (error == 16)`

`{`

`message = "Invalid hexadecimal integer";`

`}`

`else if (error == 255)`

`{`

`message = "Integer is too large";`

`}`

`else`

`{`

`message = "Unknown error";`

`}`

`fprintf (stderr, "%s: %s: %s\n", program_file_name, message, argument_array [argument_index]);`

`return_value = EXIT_FAILURE;`

`}`

`}`

`/* Set the print format string */`

`char format_string [] = "%lu %llu %llu\n";`

`format_string  = base_indicator;`

`format_string  = base_indicator;`

`format_string  = base_indicator;`

`if (maximum_a > 1)`

`{`

`if (minimum_a < 2)`

`{`

`minimum_a = 2;`

`}`

`/* Create an array to store some prime numbers */`

`unsigned short int upper_bound_for_primes = sqrt (maximum_a);`

`if ((upper_bound_for_primes > 2) && ((upper_bound_for_primes % 2) == 0))`

`{`

`upper_bound_for_primes --;`

`}`

`if ((upper_bound_for_primes > 3) && ((upper_bound_for_primes % 3) ==0))`

`{`

`upper_bound_for_primes -= 2;`

`}`

`unsigned short int prime_array_size = ((upper_bound_for_primes - 1) / 3) + 2;`

`unsigned short int * prime_array;`

`/* The initial size of the array of prime numbers is likely to be too large, so we can try a smaller size if we cannot allocate enough memory for the initial size */`

`while ((prime_array_size > 0) && ((prime_array = malloc (prime_array_size * sizeof (unsigned short int))) == NULL))`

`{`

`prime_array_size *= 0.9;`

`}`

`if (prime_array == NULL)`

`{`

`perror (program_file_name);`

`return_value = EXIT_FAILURE;`

`}`

`else`

`{`

`unsigned char const false = 0;`

`unsigned char const true = 1;`

`/* Find prime numbers and store them in the array */`

`unsigned short int candidate = 2;`

`unsigned short int prime_count = 0;`

`unsigned char found_enough_primes = false;`

`unsigned char step = 4;`

`unsigned char upper_bound_for_square_roots = sqrt (upper_bound_for_primes);`

`unsigned char square_root = 0;`

`/* Store the first two prime numbers */`

`if (prime_array_size > 0)`

`{`

`prime_array [prime_count] = candidate;`

`prime_count ++;`

`candidate = 3;`

`if (prime_array_size > 1)`

`{`

`prime_array [prime_count] = candidate;`

`prime_count ++;`

`candidate = 5;`

`/* Find more prime numbers */`

`while ((prime_count < prime_array_size) && (found_enough_primes == false))`

`{`

`/* We have found enough primes if we have reached the last candidate */`

`if (candidate >= upper_bound_for_primes)`

`{`

`found_enough_primes = true;`

`}`

`else if (step == 2)`

`{`

`/* Adjust the step size to skip the odd numbers greater than 5 that are equal to 3 mod 6, because they are divisible by 3 */`

`step = 4;`

`}`

`else`

`{`

`step = 2;`

`}`

`/* Calculate the square root of the candidate */`

`while ((square_root < upper_bound_for_square_roots) && (((square_root + 1) * (square_root + 1)) <= candidate))`

`{`

`square_root ++;`

`}`

`/* Test whether any of the prime numbers that are less than the square root of the canddidate divide the candidate. We need not test whether the prime numbers 2 or 3 divide the candidate, because we have only considered candidates that are not divisible by 2 or 3. */`

`unsigned char found_factor = false;`

`unsigned short int prime_counter;`

`for (prime_counter = 2; (found_factor == false) && (prime_counter < prime_count) && (prime_array [prime_counter] <= square_root); prime_counter ++)`

`{`

`if ((candidate % (prime_array [prime_counter])) == 0)`

`{`

`/* Found a factor, so the candidate is not prime */`

`found_factor = true;`

`}`

`}`

`if (found_factor == false)`

`{`

`/* Did not find a factor, so the candidate is prime and we can add it to the array */`

`prime_array [prime_count] = candidate;`

`prime_count ++;`

`}`

`candidate += step;`

`}`

`}`

`}`

`if (found_enough_primes == false)`

`{`

`/* When we have not found enough prime numbers, reduce the maximum a to test */`

`unsigned long int new_maximum_a = ((1ul * candidate) * candidate) - 1;`

`if (new_maximum_a < maximum_a)`

`{`

`fprintf (stderr, "%s: Unable to allocate enough memory to store all the necessary primes, so reducing the maximum value of a to %lu\n", program_file_name, new_maximum_a);`

`return_value = EXIT_FAILURE;`

`maximum_a = new_maximum_a;`

`}`

`}`

`else if (prime_array_size > prime_count)`

`{`

`/* When we have found enough prime numbers, try to reduce the size of the array, so that it is only as large as we need */`

`void * new_prime_array;`

`prime_array_size = prime_count;`

`new_prime_array = realloc (prime_array, prime_array_size * sizeof (unsigned short int));`

`if ((new_prime_array != NULL) && (new_prime_array != prime_array))`

`{`

`prime_array = new_prime_array;`

`}`

`}`

`/* The minimum number that is divisible by n different prime numbers is the product of the lowest n prime numbers */`

`unsigned char maximum_prime_factor_count = 0;`

`unsigned long int result = maximum_a;`

`while ((result > 1) && (maximum_prime_factor_count < prime_count))`

`{`

`result /= prime_array [maximum_prime_factor_count];`

`if (result >= 1)`

`{`

`maximum_prime_factor_count ++;`

`}`

`}`

`if (result > 1)`

`{`

`if (maximum_a < 6)`

`{`

`maximum_prime_factor_count = 1;`

`result = 0;`

`}`

`else if (maximum_a < 30)`

`{`

`maximum_prime_factor_count = 2;`

`result = 0;`

`}`

`}`

`if (result > 1)`

`{`

`fprintf (stderr, "%s: Unable to calculate the maximum number of prime factors\n", program_file_name);`

`return_value = EXIT_FAILURE;`

`}`

`else`

`{`

`/* Create an array to store the factors that are each a power of a single prime number */`

`unsigned long int * prime_factor_array = malloc (maximum_prime_factor_count * sizeof (unsigned long int));`

`if (prime_factor_array == NULL)`

`{`

`perror (program_file_name);`

`return_value = EXIT_FAILURE;`

`}`

`else`

`{`

`unsigned short int const factor_array_size = 1 << (maximum_prime_factor_count - 1);`

`unsigned long int * factor_array = malloc (factor_array_size * sizeof (unsigned long int));`

`if (factor_array == NULL)`

`{`

`perror (program_file_name);`

`return_value = EXIT_FAILURE;`

`}`

`else`

`{`

`/* Create an array to store the square of each prime number */`

`unsigned long int * squared_prime_array = malloc (prime_count * sizeof (unsigned long int));`

`if (squared_prime_array != NULL)`

`{`

`unsigned short int prime_counter;`

`for (prime_counter = 0; prime_counter < prime_count; prime_counter ++)`

`{`

`squared_prime_array [prime_counter] = (1ul * prime_array [prime_counter]) * prime_array [prime_counter];`

`}`

`}`

`/* Search for primitive triples by considering the factors of a that are prime numbers raised to a positive power. (These are prime powers.) The square of these factors are factors of a squared that have no common prime factors. We can use each of the products of these factors when they are small enough. */`

`unsigned char reached_maximum = false;`

`double boundary_factor = sqrt (2) - 1;`

`unsigned long int a;`

`for (a = minimum_a; reached_maximum == false; a ++)`

`{`

`if (a == maximum_a)`

`{`

`reached_maximum = true;`

`}`

`unsigned char const a_mod_4 = a % 4;`

`if (a_mod_4 != 2)`

`{`

`unsigned long int result;`

`unsigned short int prime_counter;`

`unsigned char prime_factor_count = 0;`

`if (a_mod_4 == 0)`

`{`

`result = a / 2;`

`prime_counter = 0;`

`}`

`else`

`{`

`result = a;`

`prime_counter = 1;`

`}`

`unsigned long int upper_bound_for_lower_factor = result * boundary_factor;`

`unsigned long long int square = (1ull * result) * result;`

`unsigned short int combination_count = 1;`

`while ((result > 1) && (prime_counter < prime_count) && (((squared_prime_array != NULL) && (squared_prime_array [prime_counter] <= result)) || ((squared_prime_array == NULL) && (((1ull * prime_array [prime_counter]) * prime_array [prime_counter]) <= result))))`

`{`

`if ((result % (prime_array [prime_counter])) == 0)`

`{`

`unsigned long int factor = prime_array [prime_counter];`

`result /= prime_array [prime_counter];`

`while ((result % (prime_array [prime_counter])) == 0)`

`{`

`factor *= prime_array [prime_counter];`

`result /= prime_array [prime_counter];`

`}`

`if (factor <= upper_bound_for_lower_factor)`

`{`

`unsigned long long int const squared_factor = (1ull * factor) * factor;`

`if (squared_factor <= upper_bound_for_lower_factor)`

`{`

`prime_factor_array [prime_factor_count] = squared_factor;`

`prime_factor_count ++;`

`combination_count *= 2;`

`}`

`}`

`}`

`prime_counter ++;`

`}`

`if ((result > 1) && (result <= upper_bound_for_lower_factor))`

`{`

`unsigned long long int const factor = (1ull * result) * result;`

`if (factor <= upper_bound_for_lower_factor)`

`{`

`prime_factor_array [prime_factor_count] = factor;`

`prime_factor_count ++;`

`combination_count *= 2;`

`}`

`}`

`unsigned char prime_factor_counter;`

`if (prime_factor_count > 1)`

`{`

`/* Sort prime factors */`

`unsigned char finished = false;`

`for (prime_factor_counter = 1; (finished == false) && (prime_factor_counter < prime_factor_count); prime_factor_counter ++)`

`{`

`unsigned char secondary_prime_factor_counter;`

`unsigned char swaps = false;`

`for (secondary_prime_factor_counter = prime_factor_counter; secondary_prime_factor_counter < prime_factor_count; secondary_prime_factor_counter ++)`

`{`

`if (prime_factor_array [secondary_prime_factor_counter - 1] > prime_factor_array [secondary_prime_factor_counter])`

`{`

`/* Swap prime factors */`

`unsigned long int const factor = prime_factor_array [secondary_prime_factor_counter - 1];`

`prime_factor_array [secondary_prime_factor_counter - 1] = prime_factor_array [secondary_prime_factor_counter];`

`prime_factor_array [secondary_prime_factor_counter] = factor;`

`swaps = true;`

`}`

`}`

`if (swaps == false)`

`{`

`finished = true;`

`}`

`}`

`}`

`/* Consider each product of the prime power factors and add those that are small enough to an array of factors */`

`unsigned short int combination_counter;`

`unsigned short int factor_count = 0;`

`for (combination_counter = 0; combination_counter < combination_count; combination_counter ++)`

`{`

`unsigned long long int factor = 1;`

`for (prime_factor_counter = 0; (factor <= upper_bound_for_lower_factor) && (prime_factor_counter < prime_factor_count); prime_factor_counter ++)`

`{`

`if ((combination_counter & (1 << prime_factor_counter)) > 0)`

`{`

`factor *= prime_factor_array [prime_factor_counter];`

`}`

`}`

`if (factor <= upper_bound_for_lower_factor)`

`{`

`factor_array [factor_count] = factor;`

`factor_count ++;`

`}`

`}`

`if (factor_count > 0)`

`{`

`unsigned short int factor_counter;`

`if (factor_count > 1)`

`{`

`/* Sort factors */`

`unsigned char finished = false;`

`for (factor_counter = 1; (finished == false) && (factor_counter < factor_count); factor_counter ++)`

`{`

`unsigned char secondary_factor_counter;`

`unsigned char swaps = false;`

`for (secondary_factor_counter = factor_counter; secondary_factor_counter < factor_count; secondary_factor_counter ++)`

`{`

`if (factor_array [secondary_factor_counter - 1] > factor_array [secondary_factor_counter])`

`{`

`/* Swap factors */`

`unsigned long int const factor = factor_array [secondary_factor_counter - 1];`

`factor_array [secondary_factor_counter - 1] = factor_array [secondary_factor_counter];`

`factor_array [secondary_factor_counter] = factor;`

`swaps = true;`

`}`

`}`

`if (swaps == false)`

`{`

`finished = true;`

`}`

`}`

`}`

`/* Use each factor to derive the values of b and c */`

`for (factor_counter = 0; factor_counter < factor_count; factor_counter ++)`

`{`

`unsigned long int lower_factor = factor_array [factor_count - factor_counter - 1];`

`unsigned long long int upper_factor = square / lower_factor;`

`unsigned long long int b, c;`

`if (a_mod_4 == 0)`

`{`

`b = upper_factor - lower_factor;`

`c = upper_factor + lower_factor;`

`}`

`else`

`{`

`b = (upper_factor - lower_factor) / 2;`

`c = (upper_factor + lower_factor) / 2;`

`}`

`printf (format_string, a, b, c);`

`}`

`}`

`}`

`}`

`if (squared_prime_array != NULL)`

`{free (squared_prime_array);`

`squared_prime_array = NULL;`

`}`

`free (factor_array);`

`factor_array = NULL;`

`}`

`free (prime_factor_array);`

`prime_factor_array = NULL;`

`}`

`}`

`}`

`free (prime_array);`

`prime_array = NULL;`

`}`

`}`

`}`

`return return_value;`

`}`