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:

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:


... 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():
 while True:
  yield n

while True: 
Personal tools