Challenge 2 results

DonMcCormack · May 21, 2020 05:31 PM

Note: A challenge entry from @David_Burnham was submitted prior to the deadline but was quarantined by my email system. I only discovered it after the original blog post was published. I am updating the results to reflect changes based on the new entry.

Thanks for everyone who submitted code for Challenge 2. Scoring was challenge, since some code worked better for smaller values and some for larger values. The table below gives the results for each part of the challenge. There were three runs for isPrime, 1000 5-digit random odd numbers not divisible by 5, 8- and 9- digit odd numbers not divisible by 5, and 12- to 15- digit odd numbers not divisible by 5. Results were based on the final two conditions. Points were divided between the challenges, 1.5 point for 1^st, 1 for 2^nd, and 0.5 for 3^rd. Times for the fastest submissions are also given (in seconds).

Place	getPrimes(1)	getPrimes(2)	isPrime(Small)	isPrime(Big)
1	@David_Burnham (0.1070)	@David_Burnham (0.0930)	@leelaim (0.8756)	@landon (16.89)
2	@leelaim	@P_Anderson	@David_Burnham	@David_Burnham
3	@landon	@MathStatChem	@P_Anderson	@P_Anderson

@leelaim isPrime(Small) time was as fast as any by a JMP staff member (+3 pts.). @landon isPrime(Big) script uses the Parallel Assign function to improve speed. Two of the entries cached their results for getPrimes using Namespace. This approach is helpful when the results need to be accessed more than once (run times are based on uncached results, however). All four of these entries have been included below.

Here is the leaderboard after two challenges:

Name	Points	Entries
@P_Anderson	27.5	2
@leelaim	26.5	2
@MathStatChem	24	2
@David_Burnham	23	2
@landon	12	1

Challenge 3 will be posted Wednesday, June 1.

Challenge 2 asked for finding primes two possible ways. Given a positive integer value n, you had to find: 1) n primes or 2) all primes under n. Let's start with 2) since it is somewhat easier to implement. The simplest methods for finding primes is trial division. Starting with 2 as a candidate and an empty set of primes, check if anything from the current set of primes divides the candidate number evenly. If not, add it to the set of primes and continue to the next candidate. Since there is nothing in the prime set to begin, 2 is added. Checking 3, it does not divide anything in {2} and is added to the prime set. The number 4, however does divide something in {2,3} and is not added. An example implementation is shown below.

	getPrimes2 = Function( {n}, {primeSet,i,j},
		primeSet = [];
		For(i=2,i<=n,i++,
			For(j=2,j<=N Row(primeSet),j++,
				If(Mod(i,j)==0, Break());
			);
			If(j>N Row(primeSet),primeSet = primeSet |/ i);
		);
		primeSet;
	);

There are a few advantages to this approach: It’s easy to understand, the code is concise and, it can be easily adapted for both prime finding challenges. Unfortunately, it’s extremely slow. Finding the primes below 100,000 took me 19.43 seconds.

A slightly more complex approach is to use a sieve. Rather than testing a number into primality, sieves test numbers out of primality. The oldest and best known is the Sieve of Eratosthenes. Start by assuming everything is prime. Beginning with 2 as the first prime, find all numbers that divide it evenly and remove them from the prime set. Continue with the next number still in the prime set. Below is a non-optimized version.

	getPrimes2 = Function( {n}, {primeSet,i,j},
		primeSet    = J( n, 1, 1 );
		primeSet[1] = 0;
		For( i = 2, i <= n, i++,
			For( j = i, j <= n / i, j++,
				primeSet[j * i] = 0
			)
		);
		Loc( primeSet );
	);

The most obvious advantage is speed. My time for finding all primes below 100,000 was 0.2811 sec, about 70X faster, and 3.4025 sec for 1,000,000. The code is still concise and easy to understand. A minor disadvantage is that you need to create a vector to track non-prime numbers (that is, you can’t grow the list by concatenation -- which we will see later is not necessarily an advantage). This makes the problem of finding n primes more difficult because you don’t know exactly how large a vector you will need.

Below are two optimized versions. The first starts by removing all multiples of two, then only loops over the odd numbers from 3 onward. Sqrt(n) is used for the upper bound of the outer loop because all non-primes less than that will have already been sieved out. For example, to find all primes under 100, you only need to sieve to 10. The sieved-out numbers for any prime greater than 10 will be the product of that prime times a number less than 10, which will have been removed with some other prime number (for example, 11*5 is removed with the sieve for 5, 11*9 with the sieve for 3, etc.). These improvements cut the runtime by more than half to 1.3545. The second shows the advantage of using JSL’s matrix operations for the inner loop, with times of 0.01625 for 1,000,000 and 2.6694 for 100,000,000, a near 100-fold improvement!

	getPrimes2 = Function( {n}, {primeSet,i,j},
		primeSet = J( n, 1, 1 );
		primeSet[4::n::2] = 0;
		primeSet[1] = 0;
		For( i = 3, i <= Sqrt( n ), i+=2,
			For( j = i, j <= n / i, j++,
				primeSet[j * i] = 0
			)
		);
		Loc( primeSet );
	);

	getPrimes2 = Function( {n}, {primeSet,i,j,upper},
		primeSet = J( n, 1, 1 );
		primeSet[4::n::2] = 0;
		primeSet[1] = 0;
		i = 3;
		upper = Sqrt(n);
		For(i=3,i<=upper,i+=2,
			If(primeSet[i] == 1,
				primeSet[(i*2)::n::i] = 0;
			);
		);
		Loc( primeSet );
	);

Can you do better? The Sieve of Atkin is theoretically more efficient and was implemented in Primes Atkins example.JSL and Primes Atkins (optimized).JSL. The first gives a non-optimized and easier to understand implementation (2.236 sec for 1,000,000). Better than the non-optimized version of Eratosthenes. The optimized version, however, couldn’t beat optimized Eratosthenes (4.914 sec for 100,000,000). In addition, it is much longer and more complex, programmatically and theoretically.

What about reading from a table of numbers? While this is not in the spirit of writing an algorithm, I didn’t say you couldn’t. @Craige_Hales provided a solution using this approach. Given the problem specification of n not to exceed 1e8, the size of the script to contain the compressed values took over 4 million characters of code, 4.6MB of disk space and about 5.5 seconds to load, regardless of the size of n. Slower than the other solutions. While neither the fastest or the shortest, I have included Craige's script because of his unique approach to the problem, use of compression, and Parallel Assign.

Finding n primes is a more challenging problem. Assuming you use a sieve, you’ll need to create a vector of values up to the nth prime. To do this, an upper bound for that value needs to be estimated. Underestimation should be avoided as it may require rerunning part or all of the algorithm to find unsieved values. Estimate value of nth prime.jmp shows 6 different estimates. Est1 gives a case where values are underestimated for small n. Many estimates only apply to n greater than a specific value and care needs to be taken to provide alternative estimates in these cases. Est2 is good for n ≥ 13 and Est3 improves on that for larger values. Unfortunately, good bounds are rarely simple. Est4 was used in Primes fast submission.jsl. While it doesn’t do as good a job and underestimates for very large values of n, it is simple and easy to remember. Est5 was derived using Nonlinear with the actual values. Est6 improves on Est5 by slightly perturbing the values so that none of the values in the table are underestimated. It is not guaranteed underestimates don’t exist between 100 and 1e8. The other table columns give the under/overestimates for each technique.

Determining primality could follow a similar path to generating primes. Given a number x, you could let i = 2. If Mod(n,i) == 0 stop, otherwise increment i. Continue until Mod(n,i) == 0 or i == n. This code below benefits from only checking odd i after the initial test at i = 2 and stopping when i > = Ceiling(Sqrt(n)).

	isPrime = Function({inVals},{i,j,outVec,nextVal,uBound},
		outVec = J(N Row(inVals),1,-1);
		For(i=1,i<=N Row(inVals),i++,
			nextVal = inVals[i];
			If(Mod(nextVal,2) == 0, 
				outVec[i] == 0;
				
			,//ELSE
				j = 3;
				uBound = Ceiling(Sqrt(nextVal));
				While(j < uBound & Mod(nextVal,j) != 0, j+=2);
				If(j >= uBound, outVec[i] = 1,outVec[i] = 0);
			)
		);
		outVec;
	);

The script above takes about 1.4 seconds to check all 1,000 eight- and nine-digit numbers in the table Table of big numbers.jmp (Small).

To improve this further, you can generate a vector of prime numbers (up to Ceiling(Sqrt(n))) and iterate through them.

	getPrimes = Function( {n}, {primeSet,i,j,upper},
		primeSet = J( n, 1, 1 );
		primeSet[4::n::2] = 0;
		primeSet[1] = 0;
		i = 3;
		upper = Sqrt(n);
		For(i=3,i<=upper,i+=2,
			If(primeSet[i] == 1,
				primeSet[(i*2)::n::i] = 0;
			);
		);
		Loc( primeSet );
	);
	
	isPrime = Function({inVals},{i,j,outVec,primeVec,nextVal},
		outVec = J(N Row(inVals),1,-1);
		primeVec = getPrimes(Ceiling(Sqrt(Max(inVals))));
		For(i=1,i<=N Row(inVals),i++,
			nextVal = inVals[i];
			j=1;
			While(j < N Row(primeVec) & Mod(nextVal,primeVec[j]) != 0, j++);
			If(j >= N Row(primeVec), outVec[i]=1,outVec[i]=0);
		);
		outVec;
	);

This cuts the time by more than half to under 0.6 seconds.

Two general notes to improve speed: In situations where a matrix value needs to be used more than once, it is faster to assign the value to a variable, using the variable instead of the indexed matrix. While this may potentially make the code longer, it will save time. Try this by replacing nextVal in the code above with its indexed matrix value (inVals[i]). Second, Loc and matrix concatenation slow as the size of the matrix grows. If you can work with a sorted matrix, Loc Sorted performs a binary search and is much faster. See Speed testing matrix operations.jsl for additional examples on matrix referencing and concatenation.

Trying to find primes in the column of large numbers (Big) is more problematic. To find the first 10 takes about 1.9 seconds. At this point, you have to rely on a specialized algorithm for improvements. Miller-Rabin (MR) is a commonly used and relatively easily implemented technique for testing primality. A deterministic version is given in millerRabin1.jsl. MR relies on modular exponentiation and care must be taken to keep from incurring overflow errors. By itself, MR takes 5.28 seconds to work through all 1000 Big numbers. It may surprise you that MR is slower than the first isPrime algorithm above (2.34 seconds) at working through the Small numbers. This is likely because it becomes more computationally efficient as numbers increase in size, suggesting an improvement: test to a certain size prime using 2) above. If the number is prime, stop. Otherwise, switch to MR. This approach can be found in millerRabin2.jsl. Some tuning was performed to find the ideal number of primes to check. This cut the time in half, to 2.33 seconds.