Matrix_exponentiation.cpp

// general implementation of matrix exponentiation .. ..
// implementation code by adist98
#include<bits/stdc++.h>

using namespace std;

int main()
{
    unsigned long long int w;
    cin >> w;
    unsigned long long int a[w][w] = {{0}};
    unsigned long long int temp[w][w] = {{0}};
    for(int i=0; i<w; i++){
        for(int j=0; j<w; j++){
            cin >> a[i][j];
            temp[i][j] = a[i][j];
        }
    }
    int q = 1000;
    while(q--){
        unsigned long long int x,y,n;
        cin >> x >> y >> n;
        unsigned long long int result[w][w];
        for(int i=0; i<w; i++){
            for(int j=0; j<w; j++){
                if(i == j){
                    result[i][j] = 1;
                }else{
                    result[i][j] = 0;
                }
            }
        }
        while(n > 0){
            if(n%2 == 1){
                unsigned long long int temp3[w][w];
                for(int i=0; i<w; i++){
                    for(int j=0; j<w; j++){
                        temp3[i][j] = 0;
                    }

                }
                for (int i = 0; i < w; i++)
                {

                    for (int j = 0; j < w; j++)
                    {

                        for (int k = 0; k < w; k++)
                            temp3[i][j] += (result[i][k]%1000000007) *
                                         (a[k][j]%1000000007);
                            temp3[i][j] = temp3[i][j]%1000000007;
                    }

                }
                for(int i=0; i<w; i++){
                    for(int j=0; j<w; j++){
                        result[i][j] = temp3[i][j];
                    }
                }
            }
            // complete squaring operation

            unsigned long long int temp2[w][w] = {{0}};
            for(int i=0; i<w; i++){
                for(int j=0; j<w; j++){
                    temp2[i][j] = 0;
                }
            }
            for (int i = 0; i < w; i++)
            {

                for (int j = 0; j < w; j++)
                {

                    for (int k = 0; k < w; k++)
                        temp2[i][j] += (a[i][k]%1000000007) *
                                     (a[k][j]%1000000007);
                        temp2[i][j] = temp2[i][j]%1000000007;
                }

            }
            for(int i=0; i<w; i++){
                for(int j=0; j<w; j++){
                    a[i][j] = temp2[i][j];
                }
            }
            n = n/2;
        }
        cout << result[x-1][y-1]%1000000007 << endl;
        for(int i=0; i<w; i++){
            for(int j=0; j<w; j++){
                a[i][j] = temp[i][j];
            }
        }
    }

}