Two Sieving Problems

August 27, 2013

The count of distinct factors of n is a function known as ω(n) in number theory. Again we use a sieve, but this time the sieve counts each factor instead of clearing a flag:

(define (omega n k)
  (let ((sieve (make-vector (+ n 1) 0)) (zs (list)))
    (do ((p 2 (+ p 1))) ((< n p))
      (when (zero? (vector-ref sieve p))
        (do ((i p (+ i p))) ((< n i))
          (vector-set! sieve i (+ (vector-ref sieve i) 1)))))
    (do ((p 2 (+ p 1))) ((< n p) (reverse zs))
      (when (= (vector-ref sieve p) k) (set! zs (cons p zs))))))

The sieve adds 1 for each new factor, and starts at p instead of p2 because it has to find composites as well as primes.

You can see the three programs, with their supporting code, at http://programmingpraxis.codepad.org/gDCKy2xH.

About these ads

Pages: 1 2 3

2 Responses to “Two Sieving Problems”

  1. Graham said

    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

    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
    """
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 637 other followers

%d bloggers like this: