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.
My first answers are in J; it’s definitely cheating when your language has a
“prime factors of” function.
Problem 1:
This takes a very long time to get anywhere; the second is much faster:
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 endA 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 """