Programming Praxis

2 Responses to “Two Sieving Problems”

1. Graham said

   August 27, 2013 at 11:07 PM

   My first answers are in J; it’s definitely cheating when your language has a
   “prime factors of” function.
   Problem 1:

   +/@(1=#@q:)(10x^50)+1+i.1e6%2
   

   This takes a very long time to get anywhere; the second is much faster:

   +/@(1=#@q:)1e6+1+4*1e6%4
   

   I originally had simpler answers, but decided that only working over odd
   numbers (first question) or numbers equivalent to 1 mod 4 (second question) was
   worth the uglification of my code.

   My second answers are in Ruby; I can’t imagine trying to do these in a language
   like C or C++, even with the help of multiprecision libraries like GMP or NTL.
   I went with the Miller-Rabin primality test.

   require 'mathn'
   
   def decomp n
     s, d = 0, n
     while d.even?
       s, d = s + 1, d >> 1
     end
     return s, d
   end
   
   
   def modpow n, e, m
     r = 1
     while e > 0
       if e.odd?
         r = (r * n) % m
       end
       e >>= 1
       n = (n * n) % m
     end
     return r
   end
   
   def miller_rabin n, k=42
     s, d= decomp(n - 1)
     k.times do
       a = 2 + rand(n - 4)
       x = modpow a, d, n
       next if [1, n - 1].include? x
       flag = (s - 1).times do
         x = (x * x) % n
         return false if x == 1
         break n - 1 if x == n - 1
       end
       next if flag == n - 1
       return false
     end
     return true
   end
   
   class Integer
     def prime?
       return false if self < 2
       return true if self == 2
       return false if self.even?
       return miller_rabin self
     end
   end
   
   def problem1
     (10**50 + 1 .. 10**50 + 10**6).step(2).select(&:prime?).length
   end
   
   def problem2
     (10**6 + 1 .. 2 * 10**6).step(4).select(&:prime?).length
   end
   
2. Paul said

   August 30, 2013 at 3:42 PM

   A Python version with segments.

   import time
   from math import ceil, sqrt
   
   from gmpy2 import is_prime
   
   from ma.primegmp import sieve
   
   def erase(flags, primes, q, n):
       for pk, qk in zip(primes, q):
           for i in xrange(qk, n, pk):
               flags[i] = 0
   
   def gen_sieve_seg(l, r, number_miller_rabin=25):
       """segmented sieve
          uses prime numbes <= 1000000
          checks candidates with miller rabin
       """
       primes = list(sieve(1000000))
       primes.remove(2)
       q = [(-(l + 1 + pk) // 2) % pk for pk in primes]
       n = ((r - l + 1) // 2)
       flags = [1] * n
       erase(flags, primes,q, n)
       for i, f in enumerate(flags):
           if f:
               cand = l + i * 2 + 1
               if is_prime(cand, number_miller_rabin):
                   yield cand
                   
   def gen_sieve_seg_4np1(l, r):
       """segmented sieve
          ony numbers 4*n+1 are sieved
       """
       primes = list(sieve(int(ceil(sqrt(r)))))
       primes.remove(2)
       assert not l % 4
       q = [(-(l + 1 - (pk % 4) * pk) // 4) % pk for pk in primes]
       n = (r - l) // 4
       flags = [1] * n
       erase(flags, primes,q, n)
       for i, f in enumerate(flags):
           if f:
               yield l + i * 4 + 1
               
   low = 10 ** 50
   high = low + 10 ** 6
   
   t0 = time.clock()
   print sum(1 for p in gen_sieve_seg(low, high, 5)),
   print "{:6.1f} sec".format(time.clock() - t0)
   
   t0 = time.clock()
   print sum(1 for p in gen_sieve_seg_4np1(1000000, 2000000)),
   print "{:6.3f} sec".format(time.clock() - t0)
   
   """
   8737    4.3 sec
   35241  0.044 sec
   """