Subtask 1:

Đặt hàm ~DP[i]~ là tổng của tích của các dãy có ~F(A) = i~. Ta có thể viết công thức quy hoạch động là: ~DP[i] = \sum_{j=1}^{i-1} DP[j] \cdot (i - j)~ ~(1)~ , khi mà thêm phần tử ~j~ vào cuối các dãy có ~F(A) = i~ sẽ nhân giá trị tích lên ~j~.

Nếu đặt ~D[i]~ là tổng tiền tố của hàm ~DP[i]~ sao cho: ~D[i] = DP[1] + DP[2] + \ldots + DP[i] = \sum_{j=1}^{i} DP[j]~. Thì ~(1)~ có thể suy diễn thành: ~DP[i] = D[1] + D[2] + \ldots + D[i-1]~

Đặt thêm hàm ~C[i]~ là tổng tiền tố của hàm ~D[i]~, ta có thể viết lại công thức như sau:

• ~DP[i] = C[i-1]~

• ~D[i] = D[i-1] + DP[i]~

• ~C[i] = C[i-1] + D[i]~

Từ đó ta chỉ cần một vòng lặp để duyệt qua các chỉ số ~i~ và tính lần lượt các giá trị ~DP[i], D[i], C[i]~.

Độ phức tạp: ~O(R)~.

Subtask 2:

Ở đây ta cần áp dụng thêm điều kiện các phần tử trong tập hợp phải chia hết cho ít nhất một phần tử thuộc ~K~. Đặt ~S = LCM(K)~.

Ở đây ta viết lại công thức ~(1)~ như sau

• ~DP[i] = \sum_{j=1}^i DP[i-j] \times j, \quad \forall j(\exists x, j | K_x)~

Ý tưởng ở đây là chia tổng tiền tố ~C[i], D[i]~ thành các thành phần đồng dư riêng biệt ~C[i][k], D[i][k], 0 \leq k < S~ để tính tổng các ~DP[j]~ sao cho ~(i-j) \equiv k \pmod{S}~.

Ta định nghĩa lại 2 hàm sau:

• ~D[i][k]~ là tổng các tiền tố ~DP[j]~ thỏa ~(i-j) \equiv k \pmod{S}~

• ~C[i][k]~ là tổng các tiền tố ~D[j][k]~, là tổng các ~DP[j]~ nhân hệ số của nó là ~(i-j+1)~ thỏa ~(i-j) \equiv k \pmod{S}~

Viết lại công thức quy hoạch động như sau:

• ~DP[i] = \sum_{u=0, (u+1)|K}^{S-1} C[i-1][u]~, thỏa điều kiện ~u+1~ chia hết cho một phần tử thuộc tập ~K~.

• ~D[i][0] = D[i-1][S-1] + DP[i]~

• ~C[i][0] = C[i-1][S-1] + D[i][0]~

• ~D[i][k] = D[i-1][k-1], \forall k=1..S-1~

• ~C[i][k] = C[i-1][k-1] + D[i][k], \forall k=1..S-1~

Chuyển từ ~i-1~ đến ~i~ thì giá trị đồng dư ~k~ cũng phải tăng lên 1, nên ~DP[i][k]~, ~C[i][k]~ chuyển trạng thái từ ~D[i-1][k-1]~, ~C[i-1][k-1]~.

Ta duyệt hết index ~i~ để tính mọi giá trị ~DP[i]~.

Độ phức tạp: ~O(R)~.

Subtask 3:

Nhắc lại công thức ở subtask 1:

• ~DP[i] = C[i-1]~

• ~D[i] = D[i-1] + DP[i]~

• ~C[i] = C[i-1] + D[i]~

Để viết hàm nhân ma trận, ta biến đổi một chút công thức truy hồi như sau:

• ~D[i] = D[i-1] + DP[i]~

• ~C[i] = C[i-1] + D[i-1] + DP[i]~

• ~DP[i+1] = C[i] = C[i-1] + D[i-1] + DP[i]~

Đổi index ở ~DP[i+1]~ để viết ma trận cho dễ

$$\left(\begin{array}{ccc} 1 & 1 & 1 \\ 1 & 1 & 0 \\ 1 & 1 & 1 \end{array} \right) \times \left(\begin{array}{c} DP[i] \\ D[i-1] \\ C[i-1] \end{array} \right) = \left(\begin{array}{c} DP[i+1] \\ D[i] \\ C[i] \end{array} \right)$$

Độ phức tạp: ~O(\log_2{R})~.

Subtask 4:

Hàm ma trận ở đây đã tăng kích thước lên thành ~420 \times 420~ theo công thức tổng quát ở subtask 2.

Ở subtask 2:

• ~DP[i] = \sum_{u=0, (u+1)|K}^{S-1} C[i-1][u]~ thỏa điều kiện ~u+1~ chia hết cho một phần tử thuộc tập ~K~

• ~D[i][0] = D[i-1][S-1] + DP[i]~

• ~C[i][0] = C[i-1][S-1] + D[i][0]~

• ~D[i][k] = D[i-1][k-1], \forall k = 1..S-1~

• ~C[i][k] = C[i-1][k-1] + D[i][k], \forall k = 1..S-1~

• ~S[i] = S[i-1] + \sum_{u=0, (u+1)|K}^{S-1} C[i-1][u]~

Viết lại công thức như sau:

• ~DP[i] = \sum_{u=0, (u+1) | K}^{S-1} C[i-1][u],~ thỏa điều kiện ~u+1~ chia hết cho một phần tử thuộc tập ~K~

• ~D[i][0] = D[i-1][S-1] + \sum_{u=0, (u+1) | K}^{S-1} C[i-1][u]~

• ~C[i][0] = C[i-1][S-1] + D[i-1][S-1] + \sum_{u=0, (u+1) | K}^{S-1} C[i-1][u]~

• ~D[i][k] = D[i-1][k-1], \quad \forall k = 1, \ldots, S-1~

• ~C[i][k] = C[i-1][k-1] + D[i-1][k-1], \quad \forall k = 1, \ldots, S-1~

• ~S[i] = S[i-1] + \sum_{u=0, (u+1) | K}^{S-1} C[i-1][u]~

$$X \times \left(\begin{array}{c} S[i-1] \\ DP[i-1] \\ \\ D[i-1][0] \\ C[i-1][0] \\ D[i-1][1] \\ C[i-1][1] \\ \vdots \\ D[i-1][79] \\ C[i-1][79] \\ \end{array} \right) = \left( \begin{array}{c} S[i] \\ DP[i] \\ \\ D[i][0] \\ C[i][0] \\ D[i][1] \\ C[i][1] \\ \vdots \\ D[i][79] \\ C[i][79] \\ \end{array} \right)$$

Lưu ý:

• Cần phải nháp và code tay khá kĩ cho ma trận

• Có thể đổi sang ~S[i-2] = S[i-1] + DP[i-1]~ để bớt đi một nhân tử ~DP[i]~

Độ phức tạp: ~O([LCM(K_1, \ldots ,K_m)]^3 \times \log_2{R})~

#include<bits/stdc++.h>
#define endl '\n'
#define fi first
#define se second
using namespace std;
typedef long double ld;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> pii;
typedef pair<ld,ld> ldld;
typedef pair<ll,ll> llll;

const int MOD=1e9+7;

ll l,r;
int n,s=1;
int k[121];

namespace sub1
{
    bool ok()
    {
        return s==1&&r<=1000000;
    }
    const int N=1e6+1;
    ll dp[N],d[N],c[N];
    void solve()
    {
        c[0]=1;
        for(int i=1;i<=r;i++)
        {
            dp[i]=c[i-1];
            d[i]=(d[i-1]+dp[i])%MOD;
            c[i]=(c[i-1]+d[i])%MOD;
        }
        cout << (c[r]-c[l-1]+MOD)%MOD << endl;
    }
}

namespace sub2
{
    bool ok()
    {
        return r<=1000000;
    }
    const int N=1e6+1;
    const int S=121;
    ll dp[2],d[2][S],c[2][S];
    void solve()
    {
        memset(dp,0,sizeof(dp));
        memset(d,0,sizeof(d));
        memset(c,0,sizeof(c));
        vector<int> vu;
        for(int i=1;i<=n;i++) for(int j=k[i];j<=s;j+=k[i]) vu.push_back(j);
        sort(vu.begin(),vu.end());
        vu.erase(unique(vu.begin(),vu.end()),vu.end());
        c[0][0]=d[0][0]=1;
        ll ans=0;
        for(int i=1;i<=r;i++)
        {
            int cc=(i&1);
            dp[cc]=0;
            for(int u:vu) (dp[cc]+=c[cc^1][u-1])%=MOD;
            d[cc][0]=(d[cc^1][s-1]+dp[cc])%MOD;
            c[cc][0]=(c[cc^1][s-1]+d[cc][0])%MOD;
            for(int k=1;k<s;k++)
            {
                d[cc][k]=d[cc^1][k-1];
                c[cc][k]=(c[cc^1][k-1]+d[cc][k])%MOD;
            }
            if(i>=l) (ans+=dp[cc])%=MOD;
        }
        cout << ans << endl;
    }
}

namespace sub3
{
    bool ok()
    {
        return s==1;
    }
    vector<vector<ll>> mul(vector<vector<ll>> &a,vector<vector<ll>> &b)
    {
        vector<vector<ll>> c(a.size(),vector<ll>(b[0].size(),0));
        for(int i=0;i<(int)a.size();i++) for(int j=0;j<(int)b[0].size();j++) for(int k=0;k<(int)b.size();k++) (c[i][j]+=a[i][k]*b[k][j]%MOD)%=MOD;
        return c;
    }
    vector<vector<ll>> base,dp;
    ll calc(ll x)
    {
        base.clear();dp.clear();
        base.resize(3);dp.resize(3);
        base[0]={1,1,1};
        base[1]={1,1,0};
        base[2]={1,1,1};
        dp[0]={1};
        dp[1]={0};
        dp[2]={1};
        while(x)
        {
            if(x&1) dp=mul(base,dp);
            base=mul(base,base);
            x>>=1;
        }
        return dp[0][0];
    }
    void solve()
    {
        cout << (calc(r)-calc(l-1)+MOD)%MOD << endl;
    }
}

namespace sub4
{
    vector<vector<ll>> mul(vector<vector<ll>> &a,vector<vector<ll>> &b)
    {
        vector<vector<ll>> c(a.size(),vector<ll>(b[0].size(),0));
        for(int i=0;i<(int)a.size();i++) for(int j=0;j<(int)b[0].size();j++) for(int k=0;k<(int)b.size();k++) (c[i][j]+=a[i][k]*b[k][j])%=MOD;
        return c;
    }
    vector<int> vu;
    vector<vector<ll>> base,dpr,dpl;
    ll calc(ll l,ll r)
    {
        /**
        s[i-1]
        dp[i-1]
        d[i-1][0]
        c[i-1][0]
        d[i-1][1]
        c[i-1][1]
        ...
        d[i-1][119]
        c[i-1][119]
        **/
        /// base
        base.clear();dpl.clear();dpr.clear();
        base.resize(2*s+2);dpr.resize(2*s+2);
        for(int i=0;i<2*s+2;i++) base[i].resize(2*s+2,0);

        base[0][0]=1;
        base[2][2*s]++;
        base[3][2*s]++;base[3][2*s+1]++;
        for(int u:vu)
        {
            base[0][2*u+1]++;
            base[1][2*u+1]++;
            base[2][2*u+1]++;
            base[3][2*u+1]++;
        }
        for(int i=1;i<=s-1;i++)
        {
            base[2+2*i][2*i]++;
            base[2+2*i+1][2*i]++;
            base[2+2*i+1][2*i+1]++;
        }
        /// dp
        dpr[0]={0};
        dpr[1]={0};
        dpr[2]={1};
        dpr[3]={1};
        for(int i=4;i/2<=s;i++) dpr[i]={0};
        dpl=dpr;
//        for(int i=0;i<2*s+2;i++)
//        {
//            for(int j=0;j<2*s+2;j++) cout << base[i][j] << ' ';
//            cout << " | " << dp[i][0] << endl;
//        }
        l--;
        while(r)
        {
            if(l&1) dpl=mul(base,dpl);
            if(r&1) dpr=mul(base,dpr);
            base=mul(base,base);
            r>>=1;l>>=1;
        }
//        cout << dp[0][0] << ' ' << dp[1][0] << endl;
        return (dpr[0][0]-dpl[0][0]+MOD)%MOD;
    }
    void solve()
    {
        vu.clear();
        for(int i=1;i<=n;i++) for(int j=k[i];j<=s;j+=k[i]) vu.push_back(j);
        sort(vu.begin(),vu.end());
        vu.erase(unique(vu.begin(),vu.end()),vu.end());
        cout << calc(l,r) << endl;
//        cout << (calc(r)-calc(l-1)+MOD)%MOD << endl;
    }
}

void solve()
{
    /// Input
    cin >> l >> r >> n;
    s=1;
    for(int i=1;i<=n;i++)
    {
        cin >> k[i];
        s=s/__gcd(s,k[i])*k[i];
    }
    /// Solve
    if(sub1::ok()) sub1::solve();
    else if(sub2::ok()) sub2::solve();
    else if(sub3::ok()) sub3::solve();
    else sub4::solve();
}

int main()
{
    freopen("a.inp","r",stdin);
    freopen("a.out","w",stdout);
    ios_base::sync_with_stdio(0);cin.tie(0);cout.tie(0);
    int t=1;
    //cin >> t;
    while(t--)
    solve();
    return 0;
}