Euler Solution 216

From ProgSoc Wiki

Jump to: navigation, search

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()
Personal tools