# Euler Solution 216

### From ProgSoc Wiki

Consider numbers t(n) of the form t(n) = 2n^(2)-1 with n > 1. The first such numbers are 7, 17, 31, 49, 71, 97, 127 and 161. It turns out that only 49 = 7*7 and 161 = 7*23 are not prime. For n ≤ 10000 there are 2202 numbers t(n) that are prime.

How many numbers t(n) are prime for n ≤ 50,000,000 ?

## Python by Althalus and SanguineV

Not a solution, just the (so far) most efficient isPrime test I've managed to write.

It all falls apart on this line for very large numbers:

x = a**d % n

For big d's, this step takes a very long time (ie when testing isPrime(4999999), s=1,d=2499999 - python doesn't handle x^2499999 well).... Works beautifully for numbers where d is smaller though.

**Update:** When you are doing powers modulo some number you can use the modulo at intermediate stages and reach the same result. So I often use a helper function to do the work:

# Raise "n" to the power "p", modulo some "m" # Note, works for "p"s that are positive integers def powMod(n,p,m): res = n while p > 1: res = (n * res) % m p = p - 1 return res

This saves problems with:

- limits of internal functions to deal with massive powers
- potential space issues with representing HUGE numbers

I suspect (but have not proven) that this is faster in practice than using "n**p % m" in general, for large "p"s.

I adapated the "isPrime" code to use the new function and also from some Haskell code I have been playing with. Mostly it is a rewrite/copy of the original (see further below). I also used more bitwise operations and saved on repeated arithemtic operations. I suspect it might speed things up a bit, but it isn't clear this will be fast enough to solve everything by itself.

# Uses the Miller-Rabin algorithm to test for primality # Should not be called on n < 2 or even n def isPrime(n,its=5): n1 = n ^ 1 s,d = 0,n1 # Changed this, on the basis that "not" is usuall optimised better than "==" # while d & 1 == 0: while not d & 1: s,d = s + 1, d>>1 for i in xrange(its): a = random.randrange(4,n) - 2 x = powMod(a,d,n) if x == 1 or x == n1: continue for j in xrange(s-1): x = x * x % n if x == 1: return False if x == n1: break if x ^ n1: return False return True

This seems to produce the same results as the algorithm below by Althalus (not surprising, it's mostly copied and optimised):

def isPrime(n, levels=5): An implementation of the Miller-Rabin method. n is the number to check for pseudo-primality levels specifies how many different 'a' values to test with. higher vlaue = more accuracy if n < 2: return False if n % 2 == 0: return False s,d = 0, n-1 while d%2 == 0: s,d = s+1,d/2 for z in xrange(levels): a = randrange(2,n-2) x = a**d % n if x == 1 or x == n-1: continue for y in xrange(s-1): x = x**2 % n if x == 1: return False if x == n - 1: break if x == n - 1: continue return False return True

Possibly related/helpful Euler's pseudoprimes: http://en.wikipedia.org/wiki/Euler_pseudoprime

**Re prime testing:** There are a variety of tests out there for primality, it appears that Miller-Rabin is accepted as the fastest algorithm (although AKS is more mathematically sound). Given we are testing (as per the code above) through 5 iterations, the chances of accidentally letting a prime slip through are supposedly 1/(4^{5})... about 1 in 1024. By those odds we should expect to see a lot of false positives in our results!

An interesting avenue to investigate may be the factors that do/do not appear in the potential prime numbers generated. For example, I did a test run over the first 2000 potential numbers and found that the following were prime factors of one or more of them:

[7,17,23,31,41,47,71,73,79,89,97,103,113,127,137,151,167,191,193,199 ,223,233,239,241,257,263,271,281,311,313,337,353,359,367,383,401,409,431 ,433,439,449,457,463,479,487,503,521,569,577,593,599,601,607,617,631,641 ,647,673,719,727,743,751,761,769,809,823,839,857,863,881,887,911,919,929 ,937,953,967,977,983,991,1009,1031,1033,1039,1049,1063,1087,1097,1103,1129 ,1151,1153,1193,1201,1217,1223,1231,1249,1279,1289,1297,1303,1319,1321,1327 ,1361,1399,1423,1433,1439,1447,1471,1481,1489,1511,1543,1567,1583,1601,1607 ,1663,1721,1759,1783,1801,1823,1831,1871,1873,1879,1999,2017,2081,2087,2089 ,2111,2129,2137,2153,2383,2393,2503,2521,2593,2689]

... and now I have to run, back to this later.

I'm convinced this is solvable without doing any primality testing. I can't see the pattern though.

Following snippet is a simple generator that returns the next number for t(n) = 2n^(2)-1 for n > 1. Don't think it actually HELPS, but hey, I was bored.

def getN(): n,c=1,6 while True: n,c=n+c,c+4 yield n n=getN() while True: print n.next()