Finding n'th term in a (linear) recurrence relation using matrix exponentiation

Question

I am trying to write (what I think is a beautiful piece of code) a Java program that would find the n'th term of a recursively defined sequence through efficient matrix exponentiation, in the following manner:

Given is a linear recurrence relation:

f(n) = c_1 * f(n-1) + c_2 * f(n-2)+ ... + c_k * f(n-k) = sum_{i=1}^k c_i * f(n-i)

With known integer coefficients C_k = c_{1<=i<=k} and known integer base cases B_k = b_{1<=i<=k}

For demonstration of notation, for the well known fibonacci sequence, this becomes C_2 = {1,1} and base cases B_2 = {1,1}, f(n) = f(n-1) + f(n-2).

The Java method takes as arguments the desired term index n, a mod m (because terms grow large very fast), an array of the coefficients C_k and an array of the base cases B_k. The method returns the value f(n) mod m.

The class looks like this:

public class test2 {

public static void main(String[] args) 
{
    long n = 13;
    long m = 100000000L;
    long [] baseCases = {1,1};
    long [] coefficients = {1,1};
    long v = f(n , m , coefficients , baseCases);
    
    System.out.println(v);
}

public static long f (long n , long m , long [] coefficients , long [] baseCases) //main method, taking arguments of nth term, mod m , coefficients and base cases
{
    // generate the transformation matrix {c1, c2 , c3} , {1,0,0...} , {0,1,0,0...} , {0,0,1,0....}........
    int k = coefficients.length;
    long [][] M = new long [k][k];
    
    for (int i = 0 ; i < k ; i++)
        M[0][i] = coefficients[i];
    
    for (int i = 1 ; i < k ; i++)
    {
        for (int j = 0 ; j < k ; j++)
            M[i][j] = 0;
        M[i][i - 1] = 1;
    }
    
    // generate the "bases" matrix
    long [][] B = new long [baseCases.length][baseCases.length];
    for (int i = 0 ; i < baseCases.length ; i++)
        B[0][i] = baseCases[i];
    for (int i = 1 ; i < B.length ; i++)
        for (int j = 0 ; j < B[0].length ; j ++)
            B[i][j] = 0;
                
    long v = ME(B , M , n - k + 1 , m); //the matrix raised to the n-k+1 power would yield the nth term. for fibonacci, k = 2 , so n-1 exponentioation
    
    return v;
}

public static void mult (long M1 [][] , long M2 [][] , long m) //method to multiply 2 matrices
{
    int r = Math.max(M1.length, M2.length);
    int c = Math.max(M1[0].length, M2[0].length);
    int mulTemp[][] = new int[r][c]; //temporary matrix to store values
    
    for (int i = 0 ; i < r ; i++)
    {
        for (int j = 0 ; j < c ; j++)
        {
            mulTemp[i][j] = 0;
            for (int k = 0 ; k < r ; k++)
                mulTemp[i][j] += M1[i][k] * M2[k][j] % m;
        }
    }
    
    //transferring the values back to M1
    for (int i = 0 ; i < M1.length ; i++) 
        for (int j = 0 ; j < M1[0].length ; j++)
            M1[i][j] = mulTemp[i][j]; 
}
  
public static long ME (long B[][] , long M[][] , long n , long m) //method for fast matrix exponentiation 
{ 
    if (n == 1)
    {
        long s = 0;
        for (int i = 0 ; i < B[0].length ; i++)
            s += B[0][i];
        
        return s;
    }
  
    ME(B , M , n / 2 , m);
  
    mult(M  , M , m);
  
    if ((n & 1) != 0)
        mult(B , M , m);

    long s = 0;
    for (int i = 0 ; i < B[0].length ; i++)
        s += B[0][i];
    
    return s;
}

}

The code is relatively concise and clear, and commented as necessary, but it doesn't do as intended. For example, as is currently written in the main method, I am trying to generate the fibonacci sequence with it. Fibonacci numbers are indeed generated, but not with the respective n's (for n=13 it gives 5 and it should give 233, for n=10 it gives 89 instead of 55).

I think I am mishandling the matrix methods, and also I am not using the base cases correctly, I don't know exactly how to use them in the code.. I assume there is something wrong with how I transform the base cases array into the B matrix.. If I understand theory correctly, the transformation matrix is of the form:

c1 , c2 , c3 , c4....

 1  , 0  , 0  , 0 ....

 0  , 1  , 0  , 0 ....

 0  , 0  , 1  , 0 ....

......................

............... 1  , 0 

And from there it is a matter of matrix exponentiation, using the base cases matrix inbetween.

I would like to know how to change this code to be able to deal with any linear recurrence relation, with the known given (integer) parameters as I mentioned, and not just the Fibonacci sequence (I am only using this sequence in the post because it's the easiest to understand and implement). I can easily look up the Fibonacci matrix exponentiation algorithm but it didn't give me much insight as to how to generalize this idea in an algorithm, even though I think the math is clear.

What I am really asking is how to implement this theory in a code.

Answer 1

Your construction of the matrix looks OK.

I don't know what you're trying to do with the base cases, though. The base cases are a vector. There should be no base cases matrix at all.

The matrix M transforms a vector of consecutive samples forward one step, so M*[f(n), f(n-1)...] = [f(n+1), f(n)...]. The base case vector is just B=[f(k),f(k-1)...f(1)]

To find f(n), where n>k, you just get (M^(n-k)*B)[0]