Matrix_exponential Template use for solving recurrence relation

/* Template of matrix exponential written by Rahul Kumar Singh (selfcompiler) date :19-oct-2012 time --8:41 P.M. */
/* ~~~~~~~~**@#*#@**~~~~~~~~~*/
#include<cstdio>
#include<iostream>
#include<math.h>
#include<vector>
#include<utility>
#include<map>
using namespace std;
#define K 3 
//size of matrix using 1 based index not 0'!!!!! ........
#define LIMIT 10000000000000
//maximum value of n
#define MODULO 1000000007
// modulo value
#define ll long long int
/*matrix structure */
struct MATRIX{
        ll mat[K+1][K+1];
    MATRIX()  //constructor
    {
        for(int i=0;i<=K;i++)
            for(int j=0;j<=K;j++)
                mat[i][j]=0;
    }//end constructor
    void print_matrix() //print matrix
    {
        for(int i=1;i<=K;i++)
        {
            for(int j=1;j<=K;j++)
                cout<<mat[i][j]<<" ";
            cout<<std::endl;
        }
    } //end print matrix

};//end of matrix structre
/* main program start here */
int main()
{
    struct MATRIX A;  
    ll column[K+1]={0,1,1,1};  //column vector to find the answer
    A.mat[1][1]=1,A.mat[1][2]=1,A.mat[1][3]=1;  //transformation matrix of 3 X 3 (1 based index)
    A.mat[2][1]=1,A.mat[2][2]=0;A.mat[2][3]=0;
    A.mat[3][1]=0,A.mat[3][2]=1,A.mat[3][3]=0;
    map<int,struct MATRIX> I;
    I[1]=A;  //transformation matrix
        ll next=2,p=1,p1,p2;
while(next<=LIMIT)    //precomputation of matrix of power 1,2,4,8,16,.......
{
             MATRIX object;
         p1=p2=p;
         for(int i=1;i<=K;i++)    //matrix multiplication start
         for(int j=1;j<=K;j++)
             for(int k=1;k<=K;k++)
             {
                  object.mat[i][j]+=(I[p1].mat[i][k]*I[p2].mat[k][j])%MODULO;
                  if(object.mat[i][j]>=MODULO)
                           object.mat[i][j]%=MODULO;

             }//matrix multipliaction end
         I[next]=object;
         p=next;
         next=next*2;
               
}//end of precomputation
    ll tc,n,count=0,pw2;
    cin>>tc;
    ll term[K+1]={0,1,1,1};     //r=q+w+e;   //term[1]>term[2]>term[3];
while(tc--)
{
    cin>>n;
    MATRIX obj;
        count=1;
                if(n>=4)
        {
        n=n-3;
        while(n)  //compute matrix of power n
        {


                     ll x=n&-n;  //farthest set bit of n
                          if(count)  
              {      //first encounter for initialization
                                 count=0;  
                 obj=I[x];

              }
              else
              {
                    MATRIX temp;  //make a temp object
                 for(int i=1;i<=K;i++)// matrix multipliaction
                 for(int j=1;j<=K;j++)
                     for(int k=1;k<=K;k++)
                         {
                           temp.mat[i][j]+=(I[x].mat[i][k]*obj.mat[k][j])%MODULO;
                           if(temp.mat[i][j]>=MODULO)
                           temp.mat[i][j]%=MODULO;
                         }//matrix multipliaction end
                obj=temp;
              }
              n=n-x;
        ll ans=0;
        for(int i=1;i<=K;i++)  //compute answer
        {
            ans+=(column[i]*obj.mat[1][i])%MODULO;
            if(ans>=MODULO)
                ans%=MODULO;
        }// finally computed
             cout<<ans<<std::endl;
        }
        else
        {
            cout<<term[n]<<std::endl;
        }

}

         return 0;
}

   
 // complexity O(logN)