Computing even huger Fibonacci numbers in Java - follow-up

Question

See the previous and initial iteration.

I have incorporated almost all the suggestions by Peter Taylor:

• The actual method returns a BigInteger instead of String.
• The actual method can deal with Fibonacci numbers with negative indices, so there is no any need to throw IllegalArgumentException.
• The only base case that is checked is \$n = 0\$.
• fibonacciMatrix renamed to pow.
• main is now more sensible.

(The matrix power method is still recursive.)

See what I have:

import java.math.BigInteger;
public class LargeFibonacciNumbers {
 public static BigInteger fibonacci(int n) {
 if (n == 0) {
 return BigInteger.ZERO;
 }
 BigInteger[][] matrix = new BigInteger[2][2];
 matrix[0][0] = BigInteger.ONE;
 matrix[0][1] = BigInteger.ONE;
 matrix[1][0] = BigInteger.ONE;
 matrix[1][1] = BigInteger.ZERO;
 BigInteger tmp = pow(matrix, Math.abs(n))[0][1];
 if (n > 0 || ((n & 1) == 1)) {
 return tmp;
 } else {
 return tmp.negate();
 }
 }
 private static BigInteger[][] multiply(BigInteger[][] matrix1, 
 BigInteger[][] matrix2) {
 BigInteger[][] ret = new BigInteger[2][2];
 ret[0][0] = matrix1[0][0].multiply(matrix2[0][0])
 .add(matrix1[0][1].multiply(matrix2[1][0]));
 ret[0][1] = matrix1[0][0].multiply(matrix2[0][1])
 .add(matrix1[0][1].multiply(matrix2[1][1]));
 ret[1][0] = matrix1[1][0].multiply(matrix2[0][0])
 .add(matrix1[1][1].multiply(matrix2[1][0]));
 ret[1][1] = matrix1[1][0].multiply(matrix2[0][1])
 .add(matrix1[1][1].multiply(matrix2[1][1]));
 return ret;
 }
 private static BigInteger[][] pow(BigInteger[][] matrix, int n) {
 if (n == 1) {
 // End the recursion.
 return matrix;
 }
 BigInteger[][] tmp = pow(matrix, n >> 1);
 tmp = multiply(tmp, tmp);
 if ((n & 1) == 1) {
 tmp = multiply(matrix, tmp);
 }
 return tmp;
 }
 public static void main(String[] args) {
 if (args.length > 0) {
 int n;
 try {
 n = Integer.parseInt(args[0]);
 } catch (NumberFormatException ex) {
 System.err.println("Could not parse \"" + args[0] + 
 "\" as an integer.");
 return;
 }
 System.out.println(fibonacci(n));
 } else {
 System.out.println("Usage: java -jar File.jar N");
 System.out.println("Where N is the index of the Fibonacci number " +
 "to compute.");
 }
 }
}

So, what do you think?

• coderodde: @EricDuminil You are wrong. My solution performs \$\Theta(\log N)\$ multiplications of a \2ドル \times 2\$ matrices, and the one in the link \$\Theta(N)\$ additions. For instance, my version computes the millionth Fibonacci number in half a second; the one behind the link takes almost 40 seconds.
• coderodde: @EricDuminil That's a nice result!
• coderodde: @EricDuminil I have already rolled the Java port, but be my guest.
• CiaPan: By recursion (or mathematical induction) all your matrices are symmetric, so you can save approximately 1/4 of work by not calculating the item at [1][0] and using item [0][1] instead.
• CiaPan: Additionally, as you're 'manually' multiply matrices instead of using some general, three-loop routine or matrix multiplication, you can get rid of a 2×2 matrix and use a 3–items vector instead, which could save some overhead on a second dimension indexing operations.

Answer 1

Theory

As described in a related Python answer, it's possible to refactor the recursion into a for loop and to use matrix arithmetics without using matrices.

Code

import java.math.BigInteger;
public class LargeFibonacciNumbers
{
 public static BigInteger fibonacci(int n) {
 BigInteger f_n = BigInteger.ZERO;
 BigInteger f_n_plus_1 = BigInteger.ONE;
 BigInteger two = BigInteger.valueOf(2);
 BigInteger f_n_plus_1_squared, f_2n, f_n_squared, f_2n_plus_1;
 for (int i = Integer.SIZE - Integer.numberOfLeadingZeros(n); i >= 0; i--) {
 f_n_squared = f_n.multiply(f_n);
 f_n_plus_1_squared = f_n_plus_1.multiply(f_n_plus_1);
 f_2n = f_n.multiply(f_n_plus_1).multiply(two).subtract(f_n_squared);
 f_2n_plus_1 = f_n_squared.add(f_n_plus_1_squared);
 if ((n >> i & 1) == 1) {
 f_n = f_2n_plus_1;
 f_n_plus_1 = f_2n.add(f_2n_plus_1);
 } else {
 f_n = f_2n;
 f_n_plus_1 = f_2n_plus_1;
 }
 }
 return f_n;
 }
 public static void main(String[] args) {
 int n = 10000000;
 long start = System.nanoTime();
 System.out.println(fibonacci(n).bitLength());
 long finish = System.nanoTime();
 System.out.println((finish - start) / 1000000 + " ms");
 }
}

Comments

bitLength() isn't defined for ints, so this method has been used instead.

Performance

This code calculates the 10-millionth Fibonacci number in less than 1.5s on my computer. In comparison, the matrix version needs almost 4s.

Note that the test code only displays the bit length of the result, not the result itself. fibonacci(1E7) has more than 2 million digits, and it takes much longer to display it in the console than to calculate it.

Answer 2

(I don't intend to go into computing Fibonacci numbers fast (Takahashi, Daisuke: "A fast algorithm for computing large Fibonacci numbers"), but keep the focus on the code presented.)
Which is not doc-commented.
With pow() cleaned up (rename tmp to squared would be misleading, squareRoot not much better, toBeSquared long/abstruse - fubar), multiply(m1, m2) comes to the fore:

• uses fixed dimensions
• hard coded "iteration" of indices (highly repetitive)
• iteration order suboptimal regarding caching
• representing 2(n)-dimensional matrices as what n-dimensional array Java has to offer, it isn't trivial to traverse the matrix in "Gray code order"

Short of implementing the Strassen algorithm, you can try to approach sweet oblivion in (even a ×ばつ2) matrix multiplication by reversing the iteration order of columns/rows (so that elements used last in the previous iteration get used first in the current one) and picking an evaluation order where each dot product reuses one column or row from the previous dot product.
The ease of getting all those indexes wrong is another reason to factor out dot():

static final BigInteger
 ONE = BigInteger.ONE,
 ZERO = BigInteger.ZERO,
 FIBONACCI[][] = {
 { ONE, ONE },
 { ONE, ZERO }
 };
/** Computes the dot product for <code>result[first][second]</code>.
 * @param result
 * @param left source of elements with first index <code>first</code>
 * @param right source of elements with second index <code>second</code>
 * @return dot product
 */
private static BigInteger dot(BigInteger[][] result,
 BigInteger[][] left,
 BigInteger[][] right,
 int first, int second) {
// trying to decide what to present to CR resulted in TLE
 boolean up = 0 == ((first ^ second) & 1);
 BigInteger sum = ZERO;
 for (int i = up ? 0 : result.length - 1,
 beyond = up ? result.length : -1,
 step = up ? 1 : -1 ;
 i != beyond ; i += step)
 sum = sum.add(left[first][i].multiply(right[i][second]));
 return result[first][second] = sum; }
/** Multiplies matrices <code>left</code> and <code>right</code>*/
private static BigInteger[][] multiply(BigInteger[][] left, 
 BigInteger[][] right) {
 BigInteger[][] ret = new BigInteger[2][2];
 dot(ret, left, right, 0, 0);
 dot(ret, left, right, 0, 1);
 dot(ret, left, right, 1, 1);
 dot(ret, left, right, 1, 0);
 return ret;
}
/** Computes the <code>n</code>th Fibonacci number. */
public static BigInteger fibonacci(long n) {
 if (n == 0)
 return ZERO;
 BigInteger tmp = pow(FIBONACCI, Math.abs(n))[0][1];
 return 0 < n || (n & 1) == 1 ? tmp : tmp.negate();
}

The for-setup in dot() looks a reason to try

 if (up)
 <for-up>
 accumulator()
 else
 <for-down>
 accumulator()