矩乘快速幂+高精度

\(f_i = f_{i - 1} + f_{i - 3}\)

#include <iostream>
#include <cstdio>
#include <cstring>

using namespace std;

const int N = 500;
struct INT
{
    int a[N];
    
    INT()
    {
        memset(a, 0, sizeof(a));
    }

    void rd()
    {
        char s[N];
        scanf("%s", s);
        a[0] = strlen(s);
        for(int i = 1; i <= a[0]; i++)
            a[i] = s[a[0] - i] - '0';
        return;
    }

    void print(char c)
    {
        if(!a[0]) cout << 0;
        for(int i = a[0]; i >= 1; i--)
            cout << a[i];
        cout << c;
        return;
    }

    INT operator + (const INT &b) const
    {
        INT c;
        c.a[0] = max(a[0], b.a[0]);
        for(int i = 1; i <= c.a[0]; i++)
        {
            c.a[i] += a[i] + b.a[i];
            c.a[i + 1] += c.a[i] / 10;
            c.a[i] %= 10;
        }
        if(c.a[c.a[0] + 1]) c.a[0]++;
        if(c.a[0] > 100) c.a[0] = 100;
        return c;
    }

    bool operator < (const INT &b) const
    {
        if(a[0] != b.a[0]) return a[0] < b.a[0];
        for(int i = a[0]; i >= 1; i--)
            if(a[i] != b.a[i]) return a[i] < b.a[i];
        return false;
    }
    
    INT operator - (const int &b) const
    {
        INT t;
        memcpy(t.a, a, sizeof(a));
        t.a[1] -= b;
        for(int i = 1; i < t.a[0]; i++)
            if(t.a[i] < 0) t.a[i] += 10, t.a[i + 1]--;
            else break;
        return t;
    }

    INT operator * (const INT &b) const
    {
        INT c;
        for(int i = 1; i <= a[0]; i++)
        {
            int w = 0;
            for(int j = 1; j <= b.a[0]; j++)
            {
                c.a[i + j - 1] += a[i] * b.a[j] + w;
                w = c.a[i + j - 1] / 10;
                c.a[i + j - 1] %= 10;
            }
            c.a[i + b.a[0]] += w;
        }
        c.a[0] = a[0] + b.a[0];
        if(!c.a[c.a[0]]) c.a[0]--;
        if(c.a[0] > 100) c.a[0] = 100;
        return c;
    }

    INT operator / (const int &b) const
    {
        INT c;
        memcpy(c.a, a, sizeof(a));
        for(int i = a[0]; i > 1; i--)
        {
            c.a[i - 1] += (c.a[i] % b) * 10;
            c.a[i] /= b;
        }
        c.a[1] /= b;
        while(c.a[0] && !c.a[c.a[0]]) c.a[0]--;
        return c;
    }

    void operator = (int b)
    {
        while(b) a[++a[0]] = b % 10, b /= 10;
    }
} n;

struct Matrix
{
    INT f[5][5];

    void print()
    {
        for(int i = 1; i <= 3; i++)
        {
            for(int j = 1; j <= 3; j++)
                f[i][j].print(' ');
            cout << endl;
        }
    }
    friend Matrix operator * (Matrix a, Matrix b)
    {
        Matrix res;
        for(int i = 1; i <= 3; i++)
            for(int j = 1; j <= 3; j++)
                for(int k = 1; k <= 3; k++)
                    res.f[i][j] = res.f[i][j] + (a.f[i][k] * b.f[k][j]);
        return res;
    }

    friend Matrix operator ^ (Matrix a, INT p)
    {
        Matrix b;
        b.f[1][1] = 1, b.f[2][2] = 1, b.f[3][3] = 1;
        while(p.a[0])
        {
            if(p.a[1] & 1) b = b * a;
            a = a * a;
            p = p / 2;
        }
        return b;
    }
} A, B;

int main()
{
    ios :: sync_with_stdio(false);
    n.rd();
    A.f[1][4] = 1, A.f[1][1] = 2, A.f[1][2] = 1, A.f[1][3] = 1;
    B.f[1][1] = 1, B.f[1][2] = 1, B.f[2][3] = 1, B.f[3][1] = 1;

    if(n.a[0] == 1 && n.a[1] <= 3) A.f[1][4 - n.a[1]].print('\n');
    else A = A * (B ^ (n - 3)), A.f[1][1].print('\n');
    return 0;
}