# Euler Solution 216

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/(45)... 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()
```