pi_license.gp

/*
The PI License

Finds a four digit prime number n such that n is
	- Balanced: A prime equidistant from the previous and subsequent prime numbers.
	- Odious: A number is odious if the binary representation of the number as an odd number of 1's, and is evil if the number of 1's is even.
	- Happy: Add the squares of the digits to get a new number. Continue the process and if you ever get to $1$, the original number is happy.
	- Apocalyptic: If 2^n contains the sequence 666, then n is apocalyptic.
	
Author: John Peach
Wild Peaches, 22 July 2024

Load with:
\r pi_license	

*/
\\---------------------------------------------------------------------------------------------------------------------------------------------
/*

Primes between two numbers

Paramters:
	n1,n2: Lower and upper bounds for prime list. Assumes n1 < n2.
	
Returns
	prime_list: Primes between n1 and n2


Example call:
prime_list = primes_between(1e3, 1e4);
Number of primes found: 1061

*/
primes_between(n1,n2) = {

	\\ Primes less than n2. The primepi function counts the number k of primes up to n. 
	\\ The function primes returns the first k primes.
	p2 = primes(primepi(n2));
	
	\\ Primes less than n1
	p1 = primes(primepi(n1));
	
	\\ Vector of 4-digit primes 
	prime_list = setminus(p2,p1);
	
	\\ Print the number of primes between n1 and n2
	print("Number of primes found: " length(prime_list));
	
	\\ Return the list of primes
	return (prime_list)	
}
\\---------------------------------------------------------------------------------------------------------------------------------------------
/*
Balanced primes found in ordered list of primes.
A prime is balanced if it is midway between adjacent primes.

Paramters:
	search_list: Ordered vector of primes 
	
Returns
	balanced_list: Primes midway between nearest primes

Example call:
balanced_primes = balanced(search_list);
Number of balanced primes found: 50
*/
balanced(search_list) = {

	\\ Number of primes in list
	n_primes = length(search_list);

	balanced_primes = [];
	for(k = 2, n_primes-1,
		if(search_list[k] - search_list[k-1] == search_list[k+1] - search_list[k],
			balanced_primes = concat(balanced_primes, search_list[k])
		)
	);
	
	\\ Print number of balanced primes found
	print("Number of balanced primes found: " length(balanced_primes));
	
	\\ Return the list of balanced primes
	return(balanced_primes)	
}
\\---------------------------------------------------------------------------------------------------------------------------------------------
/*
Odious numbers found in a list of integers.
A number is odious if the sum of the 1's in the binary representation is odd, and evil otherwise.

Paramters:
	search_list: List of integers
	
Returns
	odious_list: Integers with an odd number of 1's in binary


Example call:
odious_list = odious(search_list);
Number of odious integers found: 27

*/
odious(search_list) = {

	\\ Number of integers in list
	n_search = length(search_list);

	odious_list = [];
	for (k = 1, n_search,
		if (vecsum(binary(search_list[k])) % 2 == 1, odious_list = concat(odious_list, search_list[k]))
	);

	\\ Print number of odious integers found
	print("Number of odious integers found: " length(odious_list));
	
	\\ Return the list of odious primes
	return(odious_list)	
}
\\---------------------------------------------------------------------------------------------------------------------------------------------
/*
Happy primes found in list of integers. 

A number is happy if the sequence formed by summing the square of the digits ends in 1. The sequence will always end in one of ten numbers.
Ref: https://mathworld.wolfram.com/HappyNumber.html#:~:text=For%20example%2C%20starting%20with%207,(OEIS%20A007770).
Ref: https://www.numbersaplenty.com/set/happy_number/


Paramters:
	search_list: Vector of integers
	
Returns
	happy_list: Happy numbers found in search_list


Example call:
happy_list = happy(search_list);
Number of happy integers found: 7

*/
happy(search_list) = {

	\\ Number of input integers
	n_search = length(search_list);

	\\ Happy sequences always contain one of the following numbers. Zero has been omitted because the only sequence with zero is the trivial case.
	happy_endings = [1, 4, 16, 20, 37, 42, 58, 89, 145];
	happy_list = [];

	for (k = 1, n_search,
		\\ Select next element from input list
		S = search_list[k];
		
		\\ Generate the sequence, stopping when the iteration reaches one of ten numbers (0 omitted).
		\\ The vecsearch operator requires a sorted vector.
		while(!vecsearch(happy_endings, S),
			S = vecsum([x^2 | x <- digits(S)]);
		);
		
		\\ Number is happy if the sequence ends in 1.
		if (S == 1, happy_list = concat(happy_list, search_list[k]));
	);
	
	\\ Print number of happy integers found
	print("Number of happy integers found: " length(happy_list));
	
	\\ Return the list of happy integers
	return(happy_list)
}
\\---------------------------------------------------------------------------------------------------------------------------------------------
/*
Finds apocalyptic numbers from list of integers. 

A number n is apocalyptic if 2^n contains 666. The search for apocalyptic numbers requires conversion to a string.
Ref: https://mathworld.wolfram.com/HappyNumber.html#:~:text=For%20example%2C%20starting%20with%207,(OEIS%20A007770).
Ref: https://www.numbersaplenty.com/set/happy_number/


Paramters:
	search_list: Vector of integers
	
Returns
	apocalypse_list: 2 raised to a number contains the sequence "666"


Example call:
apocalypse_list = apocalyptic(search_list);
Number of apocalyptic integers found: 7
*/
apocalyptic(search_list) = {

	\\ Number of integers in input vector
	n_search = length(search_list);
	
	apocalyptic_list = [];
	for (k = 1, n_search,
		
		\\ Convert 2^p to a string
		S = Str(2^search_list[k]);
		
		\\ Search S for length 3 substrings containing "666"
		j = 1;
		\\ Flag indicating "666" found
		apocalypse_now = 0;
		while(j < length(S) - 2 && !apocalypse_now,
			if(substr(S,j,3) == "666", apocalypse_now = 1, j += 1)
		);
		if ( apocalypse_now, apocalyptic_list = concat(apocalyptic_list, search_list[k]) );
	);
	
	\\ Print number of apocalyptic integers found
	print("Number of apocalyptic integers found: " length(apocalyptic_list));
	
	\\ Return the list of apocalyptic integers
	return(apocalyptic_list)
}
\\---------------------------------------------------------------------------------------------------------------------------------------------
\\ Define a function to extract a substring beginning at n with length m
\\ Ref: https://stackoverflow.com/questions/66966205/is-there-a-better-way-to-get-substring-in-pari-gp
substr(s,n,m) = strchr(Vecsmall(s)[n..n+m-1]);
\\---------------------------------------------------------------------------------------------------------------------------------------------