Last active
December 22, 2015 23:29
-
-
Save pjt33/6547705 to your computer and use it in GitHub Desktop.
Enumerator for practical numbers
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
package com.akshor.pjt33.math; | |
import java.util.*; | |
public class PracticalNumbers | |
{ | |
public static void main(String[] args) { | |
Generator gen = new Generator(); | |
// Print all practical numbers up to 30000 | |
for (int pr = gen.next(); pr <= 30000; pr = gen.next()) { | |
System.out.println(pr); | |
} | |
} | |
public static class Generator implements Iterator<Integer> | |
{ | |
private final PriorityQueue<Tuple> queue; | |
private final Map<Integer, Integer> nextPrime; | |
public Generator() { | |
queue = new PriorityQueue<Tuple>(); | |
queue.offer(new Tuple(1, 2, 1, 2)); | |
nextPrime = new HashMap<Integer, Integer>(); | |
nextPrime.put(2, 3); | |
} | |
@Override | |
public boolean hasNext() { | |
return !queue.isEmpty(); | |
} | |
@Override | |
public Integer next() { | |
Tuple t = queue.poll(); | |
queue.offer(new Tuple(t.n * t.p, t.p, t.sigma * (t.p_pow * t.p -1) / (t.p_pow - 1), t.p * t.p_pow)); | |
int q = primeSuccessor(t.p); | |
if (q <= 1 + t.sigma) { | |
queue.offer(new Tuple(t.n * q, q, t.sigma * (q + 1), q * q)); | |
} | |
// If it was a new prime, also consider replacing it. | |
// This is an optimisation over pushing all of the tuples for p < q <= 1 + sigma at once. | |
if (t.p_pow == t.p * t.p && q <= 1 + t.sigma / (t.p + 1)) { | |
queue.offer(new Tuple(t.n / t.p * q, q, t.sigma / (t.p + 1) * (q + 1), q*q)); | |
} | |
return t.n; | |
} | |
@Override | |
public void remove() { | |
throw new UnsupportedOperationException(); | |
} | |
private int primeSuccessor(int p) { | |
Integer q = nextPrime.get(p); | |
if (q == null) { | |
q = p + 2; | |
while (!isPrime(q)) q += 2; | |
nextPrime.put(p, q); | |
} | |
return q; | |
} | |
private boolean isPrime(int p) { | |
for (Integer q : nextPrime.keySet()) { | |
if (p % q == 0) return false; | |
} | |
return true; | |
} | |
} | |
private static class Tuple implements Comparable<Tuple> | |
{ | |
public final int n; | |
public final int p; // Largest prime factor of n | |
public final int sigma; // Sum of divisors of n | |
public final int p_pow; // p^(a+1) where p^a divides n and p^(a+1) doesn't | |
public Tuple(int n, int p, int sigma, int p_pow) { | |
this.n = n; | |
this.p = p; | |
this.sigma = sigma; | |
this.p_pow = p_pow; | |
} | |
@Override | |
public int compareTo(Tuple that) { | |
return n < that.n ? -1 : n == that.n ? 0 : 1; | |
} | |
} | |
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# Tested with Python 2.7 | |
import Queue | |
def practical_numbers(): | |
nextPrime = {2: 3} | |
# Each element of Q is a tuple (n, p, sigma(n), p^a) | |
Q = Queue.PriorityQueue() | |
Q.put((1, 2, 1, 2)) | |
while True: | |
n, p, sigma, pPow = Q.get() | |
Q.put((n*p, p, sigma*(pPow*p-1)/(pPow-1), p*pPow)) | |
if p not in nextPrime: | |
# Bertrand's postulate guarantees that there's a prime p' such that p < p' < 2p. | |
nextPrime[p] = min(c for c in xrange(p+1, 2*p) if not any(c % k == 0 for k in nextPrime.keys())) | |
q = nextPrime[p] | |
if (q < sigma+1): | |
Q.put((n*q, q, sigma*(q+1), q*q)) | |
# If it was a new prime, also consider replacing it. | |
# This is an optimisation over pushing all of the tuples for p < q <= 1 + sigma at once. | |
if (pPow == p*p and q < 2+sigma/(p+1)): | |
Q.put((n/p*q, q, sigma/(p+1)*(q+1), q*q)) | |
yield n | |
# Sample usage: | |
if __name__ == '__main__': | |
pr = practical_numbers() | |
for i in range(8): | |
print pr.next() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment