Invariant Properties

Getting an Infinite List of Primes in Java

Bear Giles | June 29, 2014

A common problem is to determine the prime factorization of a number. The brute force approach is trial division (Wikipedia, Khan Academy) but that requires a lot of wasted effort if multiple numbers must be factored.

One widely used solution is the Sieve of Eratosthenes (Wikipedia, Math World). It is easy to modify the Sieve of Eratosthenes to contain the largest prime factor of each composite number. This makes it extremely cheap to subsequently compute the prime factorization of numbers.

If we only care about primality we can either use a bitmap with the Sieve of Eratosthenes, or use the Sieve of Atkin).

(Sidenote: for clarity I’m leaving out the common optimizations that follow from the facts that a prime number is always “1 mod 2, n > 2” and “1 or 5 mod 6, n > 5”. This can substantially reduce the amount of memory required for a sieve.)

public enum SieveOfEratosthenes {
SIEVE;
private int[] sieve;
private SieveOfEratosthenes() {
// initialize with first million primes - 15485865
// initialize with first 10k primes - 104729
sieve = initialize(104729);
}
/**
* Initialize the sieve.
*/
private int[] initialize(int sieveSize) {
long sqrt = Math.round(Math.ceil(Math.sqrt(sieveSize)));
long actualSieveSize = (int) (sqrt * sqrt);
// data is initialized to zero
int[] sieve = new int[actualSieveSize];
for (int x = 2; x < sqrt; x++) {
if (sieve[x] == 0) {
for (int y = 2 * x; y < actualSieveSize; y += x) {
sieve[y] = x;
}
}
}
return sieve;
}
/**
* Is this a prime number?
*
* @FIXME handle n >= sieve.length!
*
* @param n
* @return true if prime
* @throws IllegalArgumentException
* if negative number
*/
public boolean isPrime(int n) {
if (n < 0) {
throw new IllegalArgumentException("value must be non-zero");
}
boolean isPrime = sieve[n] == 0;
return isPrime;
}
/**
* Factorize a number
*
* @FIXME handle n >= sieve.length!
*
* @param n
* @return map of prime divisors (key) and exponent(value)
* @throws IllegalArgumentException
* if negative number
*/
private Map<Integer, Integer> factorize(int n) {
if (n < 0) {
throw new IllegalArgumentException("value must be non-zero");
}
final Map<Integer, Integer> factors = new TreeMap<Integer, Integer>();
for (int factor = sieve[n]; factor > 0; factor = sieve[n]) {
if (factors.containsKey(factor)) {
factors.put(factor, 1 + factors.get(factor));
} else {
factors.put(factor, 1);
}
n /= factor;
}
// must add final term
if (factors.containsKey(n)) {
factors.put(n, 1 + factors.get(n));
} else {
factors.put(n, 1);
}
return factors;
}
/**
* Convert a factorization to a human-friendly string. The format is a
* comma-delimited list where each element is either a prime number p (as
* "p"), or the nth power of a prime number as "p^n".
*
* @param factors
* factorization
* @return string representation of factorization.
* @throws IllegalArgumentException
* if negative number
*/
public String toString(Map factors) {
StringBuilder sb = new StringBuilder(20);
for (Map.Entry entry : factors.entrySet()) {
sb.append(", ");
if (entry.getValue() == 1) {
sb.append(String.valueOf(entry.getKey()));
} else {
sb.append(String.valueOf(entry.getKey()));
sb.append("^");
sb.append(String.valueOf(entry.getValue()));
}
}
return sb.substring(2);
}
}

public enum SieveOfEratosthenes {
    SIEVE;
    
    private int[] sieve;

    private SieveOfEratosthenes() {
        // initialize with first million primes - 15485865
        // initialize with first 10k primes - 104729
        sieve = initialize(104729);
    }

    /**
     * Initialize the sieve.
     */
    private int[] initialize(int sieveSize) {
        long sqrt = Math.round(Math.ceil(Math.sqrt(sieveSize)));
        long actualSieveSize = (int) (sqrt * sqrt);

        // data is initialized to zero
        int[] sieve = new int[actualSieveSize];

        for (int x = 2; x < sqrt; x++) {
            if (sieve[x] == 0) {
                for (int y = 2 * x; y < actualSieveSize; y += x) {
                    sieve[y] = x;
                }
            }
        }

        return sieve;
    }

    /**
     * Is this a prime number?
     *
     * @FIXME handle n >= sieve.length!
     * 
     * @param n
     * @return true if prime
     * @throws IllegalArgumentException
     *             if negative number
     */
    public boolean isPrime(int n) {
        if (n < 0) {
            throw new IllegalArgumentException("value must be non-zero");
        }

        boolean isPrime = sieve[n] == 0;

        return isPrime;
    }

    /**
     * Factorize a number
     *
     * @FIXME handle n >= sieve.length!
     * 
     * @param n
     * @return map of prime divisors (key) and exponent(value)
     * @throws IllegalArgumentException
     *             if negative number
     */
    private Map<Integer, Integer> factorize(int n) {
        if (n < 0) {
            throw new IllegalArgumentException("value must be non-zero");
        }

        final Map<Integer, Integer> factors = new TreeMap<Integer, Integer>();

        for (int factor = sieve[n]; factor > 0; factor = sieve[n]) {
            if (factors.containsKey(factor)) {
                factors.put(factor, 1 + factors.get(factor));
            } else {
                factors.put(factor, 1);
            }

            n /= factor;
        }

        // must add final term
        if (factors.containsKey(n)) {
            factors.put(n, 1 + factors.get(n));
        } else {
            factors.put(n, 1);
        }

        return factors;
    }

    /**
     * Convert a factorization to a human-friendly string. The format is a
     * comma-delimited list where each element is either a prime number p (as
     * "p"), or the nth power of a prime number as "p^n".
     * 
     * @param factors
     *            factorization
     * @return string representation of factorization.
     * @throws IllegalArgumentException
     *             if negative number
     */
    public String toString(Map factors) {
        StringBuilder sb = new StringBuilder(20);

        for (Map.Entry entry : factors.entrySet()) {
            sb.append(", ");

            if (entry.getValue() == 1) {
                sb.append(String.valueOf(entry.getKey()));
            } else {
                sb.append(String.valueOf(entry.getKey()));
                sb.append("^");
                sb.append(String.valueOf(entry.getValue()));
            }
        }

        return sb.substring(2);
    }
}

This code has a major weakness – it will fail if the requested number is out of range. There is an easy fix – we can dynamically resize the sieve as required. We use a Lock to ensure multithreaded calls don’t get the sieve in an intermediate state. We need to be careful to avoid getting into a deadlock between the read and write locks.

private final ReadWriteLock lock = new ReentrantReadWriteLock();
/**
* Initialize the sieve. This method is called when it is necessary to grow
* the sieve.
*/
private void reinitialize(int n) {
try {
lock.writeLock().lock();
// allocate 50% more than required to minimize thrashing.
initialize((3 * n) / 2);
} finally {
lock.writeLock().unlock();
}
}
/**
* Is this a prime number?
*
* @param n
* @return true if prime
* @throws IllegalArgumentException
* if negative number
*/
public boolean isPrime(int n) {
if (n < 0) {
throw new IllegalArgumentException("value must be non-zero");
}
if (n > sieve.length) {
reinitialize(n);
}
boolean isPrime = false;
try {
lock.readLock().lock();
isPrime = sieve[n] == 0;
} finally {
lock.readLock().unlock();
}
return isPrime;
}
/**
* Factorize a number
*
* @param n
* @return map of prime divisors (key) and exponent(value)
* @throws IllegalArgumentException
* if negative number
*/
private Map<Integer, Integer> factorize(int n) {
if (n < 0) {
throw new IllegalArgumentException("value must be non-zero");
}
final Map<Integer, Integer> factors = new TreeMap<Integer, Integer>();
try {
if (n > sieve.length) {
reinitialize(n);
}
lock.readLock().lock();
for (int factor = sieve[n]; factor > 0; factor = sieve[n]) {
if (factors.containsKey(factor)) {
factors.put(factor, 1 + factors.get(factor));
} else {
factors.put(factor, 1);
}
n /= factor;
}
} finally {
lock.readLock().unlock();
}
// must add final term
if (factors.containsKey(n)) {
factors.put(n, 1 + factors.get(n));
} else {
factors.put(n, 1);
}
return factors;
}

    private final ReadWriteLock lock = new ReentrantReadWriteLock();

    /**
     * Initialize the sieve. This method is called when it is necessary to grow
     * the sieve.
     */
    private void reinitialize(int n) {
        try {
            lock.writeLock().lock();
            // allocate 50% more than required to minimize thrashing.
            initialize((3 * n) / 2);
        } finally {
            lock.writeLock().unlock();
        }
    }

    /**
     * Is this a prime number?
     * 
     * @param n
     * @return true if prime
     * @throws IllegalArgumentException
     *             if negative number
     */
    public boolean isPrime(int n) {
        if (n < 0) {
            throw new IllegalArgumentException("value must be non-zero");
        }

        if (n > sieve.length) {
            reinitialize(n);
        }

        boolean isPrime = false;
        try {
            lock.readLock().lock();
            isPrime = sieve[n] == 0;
        } finally {
            lock.readLock().unlock();
        }

        return isPrime;
    }

    /**
     * Factorize a number
     * 
     * @param n
     * @return map of prime divisors (key) and exponent(value)
     * @throws IllegalArgumentException
     *             if negative number
     */
    private Map<Integer, Integer> factorize(int n) {
        if (n < 0) {
            throw new IllegalArgumentException("value must be non-zero");
        }

        final Map<Integer, Integer> factors = new TreeMap<Integer, Integer>();

        try {
            if (n > sieve.length) {
                reinitialize(n);
            }

            lock.readLock().lock();
            for (int factor = sieve[n]; factor > 0; factor = sieve[n]) {
                if (factors.containsKey(factor)) {
                    factors.put(factor, 1 + factors.get(factor));
                } else {
                    factors.put(factor, 1);
                }

                n /= factor;
            }
        } finally {
            lock.readLock().unlock();
        }

        // must add final term
        if (factors.containsKey(n)) {
            factors.put(n, 1 + factors.get(n));
        } else {
            factors.put(n, 1);
        }

        return factors;
    }

Iterable<Integer> and foreach loops

In the real world it’s often easier to use a foreach loop (or explicit Iterator) than to probe a table item by item. Fortunately it’s easy to create an iterator that’s built on top of our self-growing sieve.

/**
* @see java.util.List#get(int)
*
* We can use a cache of the first few (1000? 10,000?) primes
* for improved performance.
*
* @param n
* @return nth prime (starting with 2)
* @throws IllegalArgumentException
* if negative number
*/
public Integer get(int n) {
if (n < 0) {
throw new IllegalArgumentException("value must be non-zero");
}
Iterator<Integer> iter = iterator();
for (int i = 0; i < n; i++) {
iter.next();
}
return iter.next();
}
/**
* @see java.util.List#indexOf(java.lang.Object)
*/
public int indexOf(Integer n) {
if (!isPrime(n)) {
return -1;
}
int index = 0;
for (int i : sieve) {
if (i == n) {
return index;
}
index++;
}
return -1;
}
/**
* @see java.lang.Iterable#iterator()
*/
public Iterator<Integer> iterator() {
return new EratosthenesListIterator();
}
public ListIterator<Integer> listIterator() {
return new EratosthenesListIterator();
}
/**
* List iterator.
*
* @author Bear Giles <bgiles@coyotesong.com>
*/
static class EratosthenesListIterator extends AbstractListIterator<Integer> {
int offset = 2;
/**
* @see com.invariantproperties.projecteuler.AbstractListIterator#getNext()
*/
@Override
protected Integer getNext() {
while (true) {
offset++;
if (SIEVE.isPrime(offset)) {
return offset;
}
}
// we'll always find a value since we dynamically resize the sieve.
}
/**
* @see com.invariantproperties.projecteuler.AbstractListIterator#getPrevious()
*/
@Override
protected Integer getPrevious() {
while (offset > 0) {
offset--;
if (SIEVE.isPrime(offset)) {
return offset;
}
}
// we only get here if something went horribly wrong
throw new NoSuchElementException();
}
}
}

   /**
     * @see java.util.List#get(int)
     *
     * We can use a cache of the first few (1000? 10,000?) primes
     * for improved performance.
     *
     * @param n
     * @return nth prime (starting with 2)
     * @throws IllegalArgumentException
     *             if negative number
     */
    public Integer get(int n) {
        if (n < 0) {
            throw new IllegalArgumentException("value must be non-zero");
        }

        Iterator<Integer> iter = iterator();
        for (int i = 0; i < n; i++) {
            iter.next();
        }

        return iter.next();
    }

    /**
     * @see java.util.List#indexOf(java.lang.Object)
     */
    public int indexOf(Integer n) {
        if (!isPrime(n)) {
            return -1;
        }

        int index = 0;
        for (int i : sieve) {
            if (i == n) {
                return index;
            }
            index++;
        }
        return -1;
    }
   /**
     * @see java.lang.Iterable#iterator()
     */
    public Iterator<Integer> iterator() {
        return new EratosthenesListIterator();
    }

    public ListIterator<Integer> listIterator() {
        return new EratosthenesListIterator();
    }

    /**
     * List iterator.
     *
     * @author Bear Giles <bgiles@coyotesong.com>
     */
    static class EratosthenesListIterator extends AbstractListIterator<Integer> {
        int offset = 2;

        /**
         * @see com.invariantproperties.projecteuler.AbstractListIterator#getNext()
         */
        @Override
        protected Integer getNext() {
            while (true) {
                offset++;
                if (SIEVE.isPrime(offset)) {
                    return offset;
                }
            }
 
            // we'll always find a value since we dynamically resize the sieve.
        }

        /**
         * @see com.invariantproperties.projecteuler.AbstractListIterator#getPrevious()
         */
        @Override
        protected Integer getPrevious() {
            while (offset > 0) {
                offset--;
                if (SIEVE.isPrime(offset)) {
                    return offset;
                }
            }

            // we only get here if something went horribly wrong
            throw new NoSuchElementException();
        }
    }
}

IMPORTANT: The code

for (int prime : SieveOfEratosthenes.SIEVE) { ... }

for (int prime : SieveOfEratosthenes.SIEVE) { ... }

is essentially an infinite loop. It will only stop once the JVM exhausts the heap space when allocating a new sieve.

In practice this means that the maximum prime we can maintain in our sieve is around 1 GB. That requires 4 GB with 4-byte ints. If we only care about primality and use a common optimization that 4 GB can hold information on 64 GB values. For simplicity we can call this 9-to-10 digit numbers (base 10).

What if we put our sieve on a disk?

There is no reason why the sieve has to remain in memory. Our iterator can quietly load values from disk instead of an in-memory cache. A 4 TB disk, probably accessed in raw mode, would seem to bump the size of our sieve to 14-to-15 digit numbers (base 10). In fact it will be a bit less because we’ll have to double the size of our primitive types from int to long, and then probably to an even larger format.

More! More! More!

We can dramatically increase the effective size of our sieve by noting that we only have to compute sqrt(n) to initialize a sieve of n values. We can flip that and say that a fully populated sieve of n values can be used to populate another sieve of n² values. In this case we’ll want to only populate a band, not the full n² sieve. Our in-memory sieve can now cover values up to roughly 40 digit numbers (base 10), and the disk-based sieve jumps to as much as 60 digit numbers (base 10), minus the space require for the larger values.

There is no reason why this approach can’t be taken even further – use a small sieve to bootstrap a larger transient sieve and use it, in turn, to populate an even larger sieve.

But how long will this take?

Aye, there’s the rub. The cost to initialize a sieve of n values is O(n²). You can use various tweaks to reduce the constants but at the end of the day you’re visiting every node once (O(n)), and then visiting some rolling value proportional to n beyond each of those points. For what it’s worth this is a problem where keeping the CPU’s cache architecture could make a big difference.

In practical terms any recent system should be able to create a sieve containing the first million primes within a few seconds. Bump the sieve to the first billion primes and the time has probably leapt to a week, maybe a month if limited JVM heap space forces us to use the disk heavily. My gut instinct is that it will take a server farm months to years to populate a TB disk

Why bother?

For most of us the main takeaway is a demonstration of how to start a collection with a small seed, say a sieve with n = 1000, and transparently grow it as required. This is easy with prime numbers but it isn’t a huge stretch to imagine the same approach being used with, oh, RSS feeds. We’re used to thinking of Iterators as some boring aspect of Collections but in fact they give us a lot of flexibility when used as part of an Iterable.

There is also a practical reason for a large prime sieve – factoring large numbers. There are several good algorithms for factoring large numbers but they’re expensive – even “small” numbers may take months or years on a server farm. That’s why the first step is always doing trial division with “small” primes – something that may take a day by itself.

Source Code

The good news is that I have published the source code for this… and the bad news is it’s part of ongoing doodling when I’m doing Project Euler problems. (There are no solutions here – it’s entirely explorations of ideas inspired by the problems. So the code is a little rough and should not be used to decide whether or not to bring me in for an interview (unless you’re impressed): http://github.com/beargiles/projecteuler.

Comments: No Comments »
Categories: java, math
Comments rss
Trackback

Fibonacci and Lucas Sequences

Bear Giles | June 24, 2014

This posts touches on three of my favorite topics – math, transferring knowledge through experience (tutorial unit tests) and the importance of research.

Most developers are aware of the Fibonacci sequence, mostly through job interviews.

To briefly recap the series is defined as

F(n) = F(n-1) + F(n-2), n > 2
F(1) = F(2) = 1

There’s a variant definition

F(n) = F(n-1) + F(n-2), n > 1
F(1) = 1
F(0) = 0

There are four well-known solutions to the white-board question “write code to calculate F(n)”.

Recursion – you need to mention this to show that you’re comfortable with recursion but you must also mention that it’s a Really Bad Idea since it requires O(2ⁿ) time and space stack since you double the work for each n.

Recursion with memoization – this can be a good approach if you point out it’s a good generalization. Basically it’s recursion but you maintain a cache (the memoization) so you only need to make the recursive call once – subsequent recursive calls just look up the cached value. This is a flexible technique since it can be used for any pure recursive function. (That is, a recursive function that depends solely on its inputs and has no side effects.) The first calls require O(n) time, stack and heap space. I don’t recall if it matters if you do the recursive call on the smaller or larger value first.

If you have a persistent cache subsequent calls require O(1) time and stack space and O(n) heap space.

Iteration – if you can’t cache the values (or just want to efficiently initialize a cache) you can use an iterative approach. It requires O(n) time but only O(1) stack and heap space.

Direct approximation – finally there is a well-known approximation using φ, or a variant using sqrt(5). It is O(1) for time, stack space, and heap space. It’s a good approach if you 1) use a lookup table for the smallest values and 2) make sure n is not too big.

The last point is often overlooked. The approximation only works as long as you don’t exceed the precision of your floating point number. F(100,000) should be good. F(1,000,000,000,000) may not be. The iterative approach isn’t practical with numbers this large.

Research

Did you know there’s two other solutions with performance O(lg(n)) (per Wikipedia) in time and space? (I’m not convinced it’s O(lg(n)) since it’s not a divide-and-conquer algorithm – the two recursive calls do not split the initial work between them – but with memoization it’s definitely less than O(n). I suspect but can’t quickly prove it’s O(lg²(n)).)

Per Wikipedia we know

F(2n-1) = F²(n) + F²(n-1)
F(2n) = F(n)(F(n) + 2F(n-1))

It is straightforward to rewrite this as a recursive method for F(n).

There is another property that considers three cases – F(3n-2), F(3n-1) and F(3n). See the code for details.

These sites provide many additional properties of the Fibonacci and related Lucas sequences. Few developers will ever need to know these properties but in those rare cases an hour of research can save days of work.

Implementation

We can now use our research to implement suitable methods for the Fibonacci and Lucas sequences.

Fibonacci calculation

(This code does not show an optimization using direct approximation for uncached values for sufficiently small n.)

/**
* Get specified Fibonacci number.
* @param n
* @return
*/
@Override
public BigInteger get(int n) {
if (n < 0) {
throw new IllegalArgumentException("index must be non-negative");
}
BigInteger value = null;
synchronized (cache) {
value = cache.get(n);
if (value == null) {
int m = n / 3;
switch (n % 3) {
case 0:
value = TWO.multiply(get(m).pow(3))
.add(THREE.multiply(get(m + 1)).multiply(get(m))
.multiply(get(m - 1)));
break;
case 1:
value = get(m + 1).pow(3)
.add(THREE.multiply(get(m + 1)
.multiply(get(m).pow(2))))
.subtract(get(m).pow(3));
break;
case 2:
value = get(m + 1).pow(3)
.add(THREE.multiply(get(m + 1).pow(2)
.multiply(get(m))))
.add(get(m).pow(3));
break;
}
cache.put(n, value);
}
}
return value;
}

    /**
     * Get specified Fibonacci number.
     * @param n
     * @return
     */
    @Override
    public BigInteger get(int n) {
        if (n < 0) {
            throw new IllegalArgumentException("index must be non-negative");
        }

        BigInteger value = null;

        synchronized (cache) {
            value = cache.get(n);

            if (value == null) {
                int m = n / 3;

                switch (n % 3) {
                case 0:
                    value = TWO.multiply(get(m).pow(3))
                               .add(THREE.multiply(get(m + 1)).multiply(get(m))
                                         .multiply(get(m - 1)));

                    break;

                case 1:
                    value = get(m + 1).pow(3)
                                .add(THREE.multiply(get(m + 1)
                                                        .multiply(get(m).pow(2))))
                                .subtract(get(m).pow(3));

                    break;

                case 2:
                    value = get(m + 1).pow(3)
                                .add(THREE.multiply(get(m + 1).pow(2)
                                                        .multiply(get(m))))
                                .add(get(m).pow(3));

                    break;
                }

                cache.put(n, value);
            }
        }

        return value;
    }

Fibonacci Iterator

/**
* ListIterator class.
* @author bgiles
*/
private static final class FibonacciIterator extends ListIterator {
private BigInteger x = BigInteger.ZERO;
private BigInteger y = BigInteger.ONE;
public FibonacciIterator() {
}
public FibonacciIterator(int startIndex, FibonacciNumber fibonacci) {
this.idx = startIndex;
this.x = fibonacci.get(idx);
this.y = fibonacci.get(idx + 1);
}
protected BigInteger getNext() {
BigInteger t = x;
x = y;
y = t.add(x);
return t;
}
protected BigInteger getPrevious() {
BigInteger t = y;
y = x;
x = t.subtract(x);
return x;
}
}

    /**
     * ListIterator class.
     * @author bgiles
     */
    private static final class FibonacciIterator extends ListIterator {
        private BigInteger x = BigInteger.ZERO;
        private BigInteger y = BigInteger.ONE;

        public FibonacciIterator() {
        }

        public FibonacciIterator(int startIndex, FibonacciNumber fibonacci) {
            this.idx = startIndex;
            this.x = fibonacci.get(idx);
            this.y = fibonacci.get(idx + 1);
        }

        protected BigInteger getNext() {
            BigInteger t = x;
            x = y;
            y = t.add(x);

            return t;
        }

        protected BigInteger getPrevious() {
            BigInteger t = y;
            y = x;
            x = t.subtract(x);

            return x;
        }
    }

Lucas calculation

/**
* Get specified Lucas number.
* @param n
* @return
*/
public BigInteger get(int n) {
if (n < 0) {
throw new IllegalArgumentException("index must be non-negative");
}
BigInteger value = null;
synchronized (cache) {
value = cache.get(n);
if (value == null) {
value = Sequences.FIBONACCI.get(n + 1)
.add(Sequences.FIBONACCI.get(n - 1));
cache.put(n, value);
}
}
return value;
}

    /**
     * Get specified Lucas number.
     * @param n
     * @return
     */
    public BigInteger get(int n) {
        if (n < 0) {
            throw new IllegalArgumentException("index must be non-negative");
        }

        BigInteger value = null;

        synchronized (cache) {
            value = cache.get(n);

            if (value == null) {
                value = Sequences.FIBONACCI.get(n + 1)
                                           .add(Sequences.FIBONACCI.get(n - 1));
                cache.put(n, value);
            }
        }

        return value;
    }

Lucas iterator

/**
* ListIterator class.
* @author bgiles
*/
private static final class LucasIterator extends ListIterator {
private BigInteger x = TWO;
private BigInteger y = BigInteger.ONE;
public LucasIterator() {
}
public LucasIterator(int startIndex, LucasNumber lucas) {
idx = startIndex;
this.x = lucas.get(idx);
this.y = lucas.get(idx + 1);
}
protected BigInteger getNext() {
BigInteger t = x;
x = y;
y = t.add(x);
return t;
}
protected BigInteger getPrevious() {
BigInteger t = y;
y = x;
x = t.subtract(x);
return x;
}
}

   /**
     * ListIterator class.
     * @author bgiles
     */
    private static final class LucasIterator extends ListIterator {
        private BigInteger x = TWO;
        private BigInteger y = BigInteger.ONE;

        public LucasIterator() {
        }

        public LucasIterator(int startIndex, LucasNumber lucas) {
            idx = startIndex;
            this.x = lucas.get(idx);
            this.y = lucas.get(idx + 1);
        }

        protected BigInteger getNext() {
            BigInteger t = x;
            x = y;
            y = t.add(x);

            return t;
        }

        protected BigInteger getPrevious() {
            BigInteger t = y;
            y = x;
            x = t.subtract(x);

            return x;
        }
    }

Education

What is the best way to educate other developers about the existence of these unexpected relationships? Code, of course!

What is the best way to educate other developers about the existence of code that demonstrates these relationships? Unit tests, of course!

It is straightforward to write unit tests that simultaneous verify our implementation and inform other developers about tricks they can use to improve their code. The key is to provide a link to additional information.

Fibonacci Sequence

public class FibonacciNumberTest extends AbstractRecurrenceSequenceTest {
private static final BigInteger MINUS_ONE = BigInteger.valueOf(-1);
/**
* Constructor
*/
public FibonacciNumberTest() throws NoSuchMethodException {
super(FibonacciNumber.class);
}
/**
* Get number of tests to run.
*/
@Override
public int getMaxTests() {
return 300;
}
/**
* Verify the definition is properly implemented.
*
* @return
*/
@Test
@Override
public void verifyDefinition() {
for (int n = 2; n < getMaxTests(); n++) {
BigInteger u = seq.get(n);
BigInteger v = seq.get(n - 1);
BigInteger w = seq.get(n - 2);
Assert.assertEquals(u, v.add(w));
}
}
/**
* Verify initial terms.
*/
@Test
@Override
public void verifyInitialTerms() {
verifyInitialTerms(Arrays.asList(ZERO, ONE, ONE, TWO, THREE, FIVE, EIGHT));
}
/**
* Verify that every third term is even and the other two terms are odd.
* This is a subset of the general divisibility property.
*
* @return
*/
@Test
public void verifyEvenDivisibility() {
for (int n = 0; n < getMaxTests(); n += 3) {
Assert.assertEquals(ZERO, seq.get(n).mod(TWO));
Assert.assertEquals(ONE, seq.get(n + 1).mod(TWO));
Assert.assertEquals(ONE, seq.get(n + 2).mod(TWO));
}
}
/**
* Verify general divisibility property.
*
* @return
*/
@Test
public void verifyDivisibility() {
for (int d = 3; d < getMaxTests(); d++) {
BigInteger divisor = seq.get(d);
for (int n = 0; n < getMaxTests(); n += d) {
Assert.assertEquals(ZERO, seq.get(n).mod(divisor));
for (int i = 1; (i < d) && ((n + i) < getMaxTests()); i++) {
Assert.assertFalse(ZERO.equals(seq.get(n + i).mod(divisor)));
}
}
}
}
/**
* Verify the property that gcd(F(m), F(n)) = F(gcd(m,n)). This is a
* stronger statement than the divisibility property.
*/
@Test
public void verifyGcd() {
for (int m = 3; m < getMaxTests(); m++) {
for (int n = m + 1; n < getMaxTests(); n++) {
BigInteger gcd1 = seq.get(m).gcd(seq.get(n));
int gcd2 = BigInteger.valueOf(m).gcd(BigInteger.valueOf(n))
.intValue();
Assert.assertEquals(gcd1, seq.get(gcd2));
}
}
}
/**
* Verify second identity (per Wikipedia): sum(F(i)) = F(n+2)-1
*/
@Test
public void verifySecondIdentity() {
BigInteger sum = ZERO;
for (int n = 0; n < getMaxTests(); n++) {
sum = sum.add(seq.get(n));
Assert.assertEquals(sum, seq.get(n + 2).subtract(ONE));
}
}
/**
* Verify third identity (per Wikipedia): sum(F(2i)) = F(2n+1)-1 and
* sum(F(2i+1)) = F(2n)
*/
@Test
public void verifyThirdIdentity() {
BigInteger sum = ZERO;
for (int n = 0; n < getMaxTests(); n += 2) {
sum = sum.add(seq.get(n));
Assert.assertEquals(sum, seq.get(n + 1).subtract(ONE));
}
sum = ZERO;
for (int n = 1; n < getMaxTests(); n += 2) {
sum = sum.add(seq.get(n));
Assert.assertEquals(sum, seq.get(n + 1));
}
}
/**
* Verify fourth identity (per Wikipedia): sum(iF(i)) = nF(n+2) - F(n+3) + 2
*/
@Test
public void verifyFourthIdentity() {
BigInteger sum = ZERO;
for (int n = 0; n < getMaxTests(); n++) {
sum = sum.add(BigInteger.valueOf(n).multiply(seq.get(n)));
BigInteger x = BigInteger.valueOf(n).multiply(seq.get(n + 2))
.subtract(seq.get(n + 3)).add(TWO);
Assert.assertEquals(sum, x);
}
}
/**
* Verify fifth identity (per Wikipedia): sum(F(i)^2) = F(n)F(n+1)
*/
public void verifyFifthIdentity() {
BigInteger sum = ZERO;
for (int n = 0; n < getMaxTests(); n += 2) {
BigInteger u = seq.get(n);
BigInteger v = seq.get(n + 1);
sum = sum.add(u.pow(2));
Assert.assertEquals(sum, u.multiply(v));
}
}
/**
* Verify Cassini's Identity - F(n-1)F(n+1) - F(n)^2 = -1^n
*/
@Test
public void verifyCassiniIdentity() {
for (int n = 2; n < getMaxTests(); n += 2) {
BigInteger u = seq.get(n - 1);
BigInteger v = seq.get(n);
BigInteger w = seq.get(n + 1);
BigInteger x = w.multiply(u).subtract(v.pow(2));
Assert.assertEquals(ONE, x);
}
for (int n = 1; n < getMaxTests(); n += 2) {
BigInteger u = seq.get(n - 1);
BigInteger v = seq.get(n);
BigInteger w = seq.get(n + 1);
BigInteger x = w.multiply(u).subtract(v.pow(2));
Assert.assertEquals(MINUS_ONE, x);
}
}
/**
* Verify doubling: F(2n-1) = F(n)^2 + F(n-1)^2 and F(2n) =
* F(n)(F(n-1)+F(n+1)) = F(n)(2*F(n-1)+F(n).
*/
@Test
public void verifyDoubling() {
for (int n = 1; n < getMaxTests(); n++) {
BigInteger u = seq.get(n - 1);
BigInteger v = seq.get(n);
BigInteger w = seq.get(n + 1);
BigInteger x = v.multiply(v).add(u.pow(2));
Assert.assertEquals(seq.get((2 * n) - 1), x);
x = v.multiply(u.add(w));
Assert.assertEquals(seq.get(2 * n), x);
x = v.multiply(v.add(TWO.multiply(u)));
Assert.assertEquals(seq.get(2 * n), x);
}
}
/**
* Verify tripling.
*/
@Test
public void verifyTripling() {
for (int n = 1; n < getMaxTests(); n++) {
BigInteger u = seq.get(n - 1);
BigInteger v = seq.get(n);
BigInteger w = seq.get(n + 1);
BigInteger x = TWO.multiply(v.pow(3))
.add(THREE.multiply(v).multiply(u).multiply(w));
Assert.assertEquals(seq.get(3 * n), x);
x = w.pow(3).add(THREE.multiply(w).multiply(v.pow(2)))
.subtract(v.pow(3));
Assert.assertEquals(seq.get((3 * n) + 1), x);
x = w.pow(3).add(THREE.multiply(w.pow(2)).multiply(v)).add(v.pow(3));
Assert.assertEquals(seq.get((3 * n) + 2), x);
}
}
}

public class FibonacciNumberTest extends AbstractRecurrenceSequenceTest {
    private static final BigInteger MINUS_ONE = BigInteger.valueOf(-1);

    /**
     * Constructor
     */
    public FibonacciNumberTest() throws NoSuchMethodException {
        super(FibonacciNumber.class);
    }

    /**
     * Get number of tests to run.
     */
    @Override
    public int getMaxTests() {
        return 300;
    }

    /**
     * Verify the definition is properly implemented.
     *
     * @return
     */
    @Test
    @Override
    public void verifyDefinition() {
        for (int n = 2; n < getMaxTests(); n++) {
            BigInteger u = seq.get(n);
            BigInteger v = seq.get(n - 1);
            BigInteger w = seq.get(n - 2);
            Assert.assertEquals(u, v.add(w));
        }
    }

    /**
     * Verify initial terms.
     */
    @Test
    @Override
    public void verifyInitialTerms() {
        verifyInitialTerms(Arrays.asList(ZERO, ONE, ONE, TWO, THREE, FIVE, EIGHT));
    }

    /**
     * Verify that every third term is even and the other two terms are odd.
     * This is a subset of the general divisibility property.
     *
     * @return
     */
    @Test
    public void verifyEvenDivisibility() {
        for (int n = 0; n < getMaxTests(); n += 3) {
            Assert.assertEquals(ZERO, seq.get(n).mod(TWO));
            Assert.assertEquals(ONE, seq.get(n + 1).mod(TWO));
            Assert.assertEquals(ONE, seq.get(n + 2).mod(TWO));
        }
    }

    /**
     * Verify general divisibility property.
     *
     * @return
     */
    @Test
    public void verifyDivisibility() {
        for (int d = 3; d < getMaxTests(); d++) {
            BigInteger divisor = seq.get(d);

            for (int n = 0; n < getMaxTests(); n += d) {
                Assert.assertEquals(ZERO, seq.get(n).mod(divisor));

                for (int i = 1; (i < d) && ((n + i) < getMaxTests()); i++) {
                    Assert.assertFalse(ZERO.equals(seq.get(n + i).mod(divisor)));
                }
            }
        }
    }

    /**
     * Verify the property that gcd(F(m), F(n)) = F(gcd(m,n)). This is a
     * stronger statement than the divisibility property.
     */
    @Test
    public void verifyGcd() {
        for (int m = 3; m < getMaxTests(); m++) {
            for (int n = m + 1; n < getMaxTests(); n++) {
                BigInteger gcd1 = seq.get(m).gcd(seq.get(n));
                int gcd2 = BigInteger.valueOf(m).gcd(BigInteger.valueOf(n))
                                     .intValue();
                Assert.assertEquals(gcd1, seq.get(gcd2));
            }
        }
    }

    /**
     * Verify second identity (per Wikipedia): sum(F(i)) = F(n+2)-1
     */
    @Test
    public void verifySecondIdentity() {
        BigInteger sum = ZERO;

        for (int n = 0; n < getMaxTests(); n++) {
            sum = sum.add(seq.get(n));
            Assert.assertEquals(sum, seq.get(n + 2).subtract(ONE));
        }
    }

    /**
     * Verify third identity (per Wikipedia): sum(F(2i)) = F(2n+1)-1 and
     * sum(F(2i+1)) = F(2n)
     */
    @Test
    public void verifyThirdIdentity() {
        BigInteger sum = ZERO;

        for (int n = 0; n < getMaxTests(); n += 2) {
            sum = sum.add(seq.get(n));
            Assert.assertEquals(sum, seq.get(n + 1).subtract(ONE));
        }

        sum = ZERO;

        for (int n = 1; n < getMaxTests(); n += 2) {
            sum = sum.add(seq.get(n));
            Assert.assertEquals(sum, seq.get(n + 1));
        }
    }

    /**
     * Verify fourth identity (per Wikipedia): sum(iF(i)) = nF(n+2) - F(n+3) + 2
     */
    @Test
    public void verifyFourthIdentity() {
        BigInteger sum = ZERO;

        for (int n = 0; n < getMaxTests(); n++) {
            sum = sum.add(BigInteger.valueOf(n).multiply(seq.get(n)));

            BigInteger x = BigInteger.valueOf(n).multiply(seq.get(n + 2))
                                     .subtract(seq.get(n + 3)).add(TWO);
            Assert.assertEquals(sum, x);
        }
    }

    /**
     * Verify fifth identity (per Wikipedia): sum(F(i)^2) = F(n)F(n+1)
     */
    public void verifyFifthIdentity() {
        BigInteger sum = ZERO;

        for (int n = 0; n < getMaxTests(); n += 2) {
            BigInteger u = seq.get(n);
            BigInteger v = seq.get(n + 1);
            sum = sum.add(u.pow(2));
            Assert.assertEquals(sum, u.multiply(v));
        }
    }

    /**
     * Verify Cassini&#039;s Identity - F(n-1)F(n+1) - F(n)^2 = -1^n
     */
    @Test
    public void verifyCassiniIdentity() {
        for (int n = 2; n < getMaxTests(); n += 2) {
            BigInteger u = seq.get(n - 1);
            BigInteger v = seq.get(n);
            BigInteger w = seq.get(n + 1);

            BigInteger x = w.multiply(u).subtract(v.pow(2));
            Assert.assertEquals(ONE, x);
        }

        for (int n = 1; n < getMaxTests(); n += 2) {
            BigInteger u = seq.get(n - 1);
            BigInteger v = seq.get(n);
            BigInteger w = seq.get(n + 1);

            BigInteger x = w.multiply(u).subtract(v.pow(2));
            Assert.assertEquals(MINUS_ONE, x);
        }
    }

    /**
     * Verify doubling: F(2n-1) = F(n)^2 + F(n-1)^2 and F(2n) =
     * F(n)(F(n-1)+F(n+1)) = F(n)(2*F(n-1)+F(n).
     */
    @Test
    public void verifyDoubling() {
        for (int n = 1; n < getMaxTests(); n++) {
            BigInteger u = seq.get(n - 1);
            BigInteger v = seq.get(n);
            BigInteger w = seq.get(n + 1);

            BigInteger x = v.multiply(v).add(u.pow(2));
            Assert.assertEquals(seq.get((2 * n) - 1), x);

            x = v.multiply(u.add(w));
            Assert.assertEquals(seq.get(2 * n), x);

            x = v.multiply(v.add(TWO.multiply(u)));
            Assert.assertEquals(seq.get(2 * n), x);
        }
    }

    /**
     * Verify tripling.
     */
    @Test
    public void verifyTripling() {
        for (int n = 1; n < getMaxTests(); n++) {
            BigInteger u = seq.get(n - 1);
            BigInteger v = seq.get(n);
            BigInteger w = seq.get(n + 1);

            BigInteger x = TWO.multiply(v.pow(3))
                              .add(THREE.multiply(v).multiply(u).multiply(w));
            Assert.assertEquals(seq.get(3 * n), x);

            x = w.pow(3).add(THREE.multiply(w).multiply(v.pow(2)))
                 .subtract(v.pow(3));
            Assert.assertEquals(seq.get((3 * n) + 1), x);

            x = w.pow(3).add(THREE.multiply(w.pow(2)).multiply(v)).add(v.pow(3));
            Assert.assertEquals(seq.get((3 * n) + 2), x);
        }
    }
}

Lucas Sequence

public class LucasNumberTest extends AbstractRecurrenceSequenceTest {
private static final FibonacciNumber fibonacci = new FibonacciNumber();
/**
* Constructor
*/
public LucasNumberTest() throws NoSuchMethodException {
super(LucasNumber.class);
}
/**
* Get number of tests to run.
*/
@Override
public int getMaxTests() {
return 300;
}
/**
* Verify the definition is properly implemented.
*
* @return
*/
@Test
@Override
public void verifyDefinition() {
for (int n = 2; n < getMaxTests(); n++) {
BigInteger u = seq.get(n);
BigInteger v = seq.get(n - 1);
BigInteger w = seq.get(n - 2);
Assert.assertEquals(u, v.add(w));
}
}
/**
* Verify initial terms.
*/
@Test
@Override
public void verifyInitialTerms() {
verifyInitialTerms(Arrays.asList(TWO, ONE, THREE, FOUR, SEVEN, ELEVEN,
BigInteger.valueOf(18), BigInteger.valueOf(29)));
}
/**
* Verify Lucas properties.
*/
@Test
public void verifyLucas() {
// L(n) = F(n-1) + F(n+1)
for (int n = 2; n < getMaxTests(); n++) {
Assert.assertEquals(seq.get(n),
fibonacci.get(n - 1).add(fibonacci.get(n + 1)));
}
}
/**
* F(2n) = L(n)F(n)
*/
@Test
public void verifyLucas2() {
for (int n = 2; n < getMaxTests(); n++) {
Assert.assertEquals(fibonacci.get(2 * n),
seq.get(n).multiply(fibonacci.get(n)));
}
}
/**
* F(n) = (L(n-1)+ L(n+1))/5
*/
@Test
public void verifyLucas3() {
for (int n = 2; n < getMaxTests(); n++) {
Assert.assertEquals(FIVE.multiply(fibonacci.get(n)),
seq.get(n - 1).add(seq.get(n + 1)));
}
}
/**
* L(n)^2 = 5 F(n)^2 + 4(-1)^n
*/
@Test
public void verifyLucas4() {
for (int n = 2; n < getMaxTests(); n += 2) {
Assert.assertEquals(seq.get(n).pow(2),
FIVE.multiply(fibonacci.get(n).pow(2)).add(FOUR));
}
for (int n = 1; n < getMaxTests(); n += 2) {
Assert.assertEquals(seq.get(n).pow(2),
FIVE.multiply(fibonacci.get(n).pow(2)).subtract(FOUR));
}
}
}

public class LucasNumberTest extends AbstractRecurrenceSequenceTest {
    private static final FibonacciNumber fibonacci = new FibonacciNumber();

    /**
     * Constructor
     */
    public LucasNumberTest() throws NoSuchMethodException {
        super(LucasNumber.class);
    }

    /**
     * Get number of tests to run.
     */
    @Override
    public int getMaxTests() {
        return 300;
    }

    /**
     * Verify the definition is properly implemented.
     *
     * @return
     */
    @Test
    @Override
    public void verifyDefinition() {
        for (int n = 2; n < getMaxTests(); n++) {
            BigInteger u = seq.get(n);
            BigInteger v = seq.get(n - 1);
            BigInteger w = seq.get(n - 2);
            Assert.assertEquals(u, v.add(w));
        }
    }

    /**
     * Verify initial terms.
     */
    @Test
    @Override
    public void verifyInitialTerms() {
        verifyInitialTerms(Arrays.asList(TWO, ONE, THREE, FOUR, SEVEN, ELEVEN,
                BigInteger.valueOf(18), BigInteger.valueOf(29)));
    }

    /**
     * Verify Lucas properties.
     */
    @Test
    public void verifyLucas() {
        // L(n) = F(n-1) + F(n+1)
        for (int n = 2; n < getMaxTests(); n++) {
            Assert.assertEquals(seq.get(n),
                fibonacci.get(n - 1).add(fibonacci.get(n + 1)));
        }
    }

    /**
     *  F(2n) = L(n)F(n)
     */
    @Test
    public void verifyLucas2() {
        for (int n = 2; n < getMaxTests(); n++) {
            Assert.assertEquals(fibonacci.get(2 * n),
                seq.get(n).multiply(fibonacci.get(n)));
        }
    }

    /**
     * F(n) = (L(n-1)+ L(n+1))/5
     */
    @Test
    public void verifyLucas3() {
        for (int n = 2; n < getMaxTests(); n++) {
            Assert.assertEquals(FIVE.multiply(fibonacci.get(n)),
                seq.get(n - 1).add(seq.get(n + 1)));
        }
    }

    /**
     * L(n)^2 = 5 F(n)^2 + 4(-1)^n
     */
    @Test
    public void verifyLucas4() {
        for (int n = 2; n < getMaxTests(); n += 2) {
            Assert.assertEquals(seq.get(n).pow(2),
                FIVE.multiply(fibonacci.get(n).pow(2)).add(FOUR));
        }

        for (int n = 1; n < getMaxTests(); n += 2) {
            Assert.assertEquals(seq.get(n).pow(2),
                FIVE.multiply(fibonacci.get(n).pow(2)).subtract(FOUR));
        }
    }
}

Conclusion

Obviously developers rarely need to compute Fibonacci numbers unless they’re working on Project Euler problems or at a job interview. This code isn’t going to have direct utility.

At the same time it’s a powerful demonstration of the value of investing an hour or two in research even if you’re sure you already know everything you need to know. You probably don’t need BigInteger implementation but some people might consider the O(lg(n)) approach preferable to the estimate using powers of φ, or could make good use of the relationships discussed on the MathWorld and Wikipedia pages.