The below function only works if $$$n$$$ is even. One simple way to handle the odd case is to break the expression as $$$a^n = a \cdot a^{n-1}$$$. We can get the answer for $$$a^{n-1}$$$ because it's even ans whether $$$n$$$ is even or odd, in $$$\leq 2$$$ steps $$$n$$$ definitely gets halved which ensures that the time complexity remains $$$O(log(n))$$$.

int binexpo(int a,int n)
{
    if (n == 1) return a;
    if (n % 2) return a * binexpo(a, n - 1); // this line happen only #SetBits times as on every step n is halving
    int L = binexpo(a, n / 2);
    int R = L;
    return R * L;    
}

Naive Approach ($$$O(log(n)^{2})$$$)

int binexpo(int x,int n){/* above code goes here */}

int GP(int r, int n)
{
    if (n == 1)
        return 1;

    if (n % 2)
        return 1 + r * GP(r, n - 1);

    int L = GP(r, n / 2);
    int R = binexpo(r, n / 2) * L;
    return L + R;
}

Now, applying our optimization, we compute { $$$ { r^n ,1 + r \dots + r^{n-1} } $$$ } together for a final time complexity of $$$O(logn)$$$

// {r^n, 1 + r + r^2 ... r^(n-1)}
pair<int, int> GP(int r, int n)
{
    if (n == 1)
        return {r, 1};

    if (n & 1) {
        auto [x, y] = GP(r, n - 1);
        return {r * x, 1 + r * y};
    }

    auto [x, L] = GP(r, n / 2);
    int R = x * L;
    return {x * x, R + L};
}

$$$R = (\frac{n}{2} + 0) \cdot r^{\frac{n}{2}} + (\frac{n}{2} + 1) \cdot r^{\frac{n}{2}+1} + \dots + (\frac{n}{2} + i) \cdot r^{\frac{n}{2}+i} + \dots + (\frac{n}{2} + \frac{n}{2}-1) \cdot r^{\frac{n}{2} + \frac{n}{2}-1} $$$

$$$\implies R = \frac{n}{2} \cdot ( r^{\frac{n}{2}} + r^{\frac{n}{2}+1} + \dots + r^{\frac{n}{2}+i} + \dots + r^{n-1} ) + r^{\frac{n}{2}} \cdot \underbrace{ \left( 0 \cdot r^{0} + 1 \cdot r^{1} + \dots + i \cdot r^{i} + \dots + (\frac{n}{2}-1) \cdot r^{\frac{n}{2}-1}\right)}_\text{L} $$$

$$$ \implies R = r^{\frac{n}{2}} \cdot \left( \frac{n}{2} \cdot \underbrace{( 1 + r^{1} + \dots + r^{i} + \dots + r^{\frac{n}{2}-1} )}_\text{GP series} + L \right) $$$

We shall call this GP series as $$$L_{0}$$$, whose sum we can already calculate efficiently. Hence,

$$$ R = r^{\frac{n}{2}} \cdot \left( \frac{n}{2} \cdot L_{0} + L \right) $$$

$$$ \implies R = f(L_{0},L) = r^{\frac{n}{2}} \cdot \left( \frac{n}{2} \cdot L_{0} + L \right) $$$

// {r^n , 1+r+r^2 ... + r^(n-1) , 0 + r + 2r^2 + ... + (n-1)r^(n-1) }
array<int, 3> APGP(int r, int n)
{
    if (n == 1)
        return {r, 1, 0};

    if (n & 1) {
        auto [x, L0, L] = APGP(r, n - 1);
        return {r * x, 1 + r * L0, r * L0 + r * L};
    }

    auto [x, L0, L] = APGP(r, n / 2);
    int R0 = x * L0;
    int R = x * (n / 2 * L0 + L);
    return {x * x, L0 + R0, L + R};
}

Case 1: for $$$m=0$$$ , which is simple GP series

$$$R_{0,\frac{n}{2}} = r^{\frac{n}{2}} \cdot \left( \sum\limits_{j=0}^{0} \binom{0}{j} \left( \frac{n}{2}\right) ^{0-j} \cdot L_{j,\frac{n}{2}} \right) = r^{\frac{n}{2}} \cdot L_{0,\frac{n}{2}}$$$

Case 2: for $$$m=1$$$ , which is AP-GP series

$$$R_{1,\frac{n}{2}} = r^{\frac{n}{2}} \cdot \left( \sum\limits_{j=0}^{1} \binom{1}{j} \left( \frac{n}{2}\right) ^{1-j} \cdot L_{j,\frac{n}{2}} \right) = r^{\frac{n}{2}} \cdot \left( \left( \frac{n}{2}\right) ^{1} \cdot L_{0,\frac{n}{2}} + L_{1,\frac{n}{2}} \right)$$$

Dimensions are : $$$ (m+1,1) = (m+1,m+2) \cdot (m+2,1) $$$.

We will return the array $$$(S_{n,0},S_{n,1}, \dots , S_{n,m}, r^n)$$$.

Dimension are: $$$ (m+1,1) = (m+1,m+1) \cdot (m+1,1) $$$.

We will return the array $$$(L_{0, n/2} + R_{0, n/2},L_{1, n/2} + R_{1, n/2}, \dots , L_{m, n/2} + R_{m, n/2}, r^n)$$$.

1. Precalculation of $$$nCr$$$ takes $$$O(m^2)$$$.

2. During our computations, whenever we have even $$$n$$$ we need to compute $$$(\frac{n}{2})^j$$$ for each $$$j \leq m$$$, which takes $$$O(m \cdot log(n))$$$.

Total complexity $$$O(m^2 + m \cdot log(n) + m^2 \cdot log(n))$$$ asymptotically remains same $$$O(m^2 \cdot log(n))$$$ .

using int64 = long long;
int MOD = 1e9+7;
const int M=1002;
void add(int&x,int y){
    x+=y;
    if(x>=MOD)x-=MOD;
}
void mul(int&x,int y){
    x = 1ll * x * y % MOD;
}

int ncr[M][M]{0};
int nCr(int nn,int rr){
 return ncr[nn][rr];   
}
// Precalculate calculate n choose r in O(M^2)
void nCrinit(){
    ncr[0][0]=1;
    for(int n=0;n<M;n++)ncr[n][n]=1,ncr[n][0]=1;
    for(int n=1;n<M;n++){
        for(int r=1;r<n;r++){
            add(ncr[n][r],ncr[n-1][r-1]);
            add(ncr[n][r],ncr[n-1][r]);
        }
    }
}

// Precalculation of powers of n in O(M * logN)
using vt = vector<int>;
map<int64,vt> ni;
int n_i(int64 n,int i){
    return ni[n][i];
}
void n_iinit(int64 n){
    while(n>0){
        if(n&1){
            n--;
            continue;
        }
        vt x(M); x[0]=1;
        for(int i=1;i<M;i++)x[i] = 1ll * x[i-1] * (n/2 % MOD) % MOD;
        ni[n/2] = x;
        n>>=1;
    }
    ni[1] = vector<int>(M,1);
}

// Function to generate the series {S(n,0),S(n,1),...,S(n,m),r^n}
vector<int> Series(int r,int64 n,int m)
{
    if(n==1){
        vector<int> a(m+2,0);
        a[0]=1;
        a[m+1]=r;
        return a;
    }
    if(n&1){
        // Recursive case for odd n
        vector<int>a = Series(r,n-1,m);
        vector<int> S(m+2,0);
        add(S[m+1],1ll * r * a[m+1]%MOD);
        add(S[0],1);
        add(S[0],1ll*r*a[0]%MOD);
        for(int i=1;i<=m;i++){
            for(int j=0;j<=i;j++){
                add(S[i], 1ll * nCr(i,j) * a[j] % MOD);
            }
            mul(S[i],r);
        }
        return S;
    }

    // Recursive case for even n
    vector<int> L = Series(r,n/2,m);
    vector<int> R(m+1,0);
    for(int i=0;i<=m;i++){
        for(int j=0;j<=i;j++){
            int h = 1;
            mul(h,nCr(i,j));
            mul(h,n_i(n/2,i-j));
            mul(h,L[m+1]);
            mul(h,L[j]);
            add(R[i],h);
        }
    }
    vector<int> S(m+2,0);
    S[m+1] = 1ll * L[m+1] * L[m+1] % MOD;
    for(int i=0;i<=m;i++){
        // S[i]=(MOD + R[i] + L[i])%MOD;
        add(S[i],R[i]);
        add(S[i],L[i]);
    }
    return S;
}