NYG的序列拆分 详解（递推矩阵快速幂）

9.19
思路：
题目描述有点恼火，其实就是列出所有合法的数（由n个[l,r]的数组成），然后把数分为连续的若干段让每一段都是一个等比数列，计算方案总数。
std的题解：
算法一
直接按定义计算，O(rnn!)，期望得分10。
算法二
r <= 2，可以推递推式，期望得分10。
算法三
注意到题目实际上是这样一个问题：将若干个等比数列拼成一个长度为n的序列，有多少种方案。由于数列中的数都是整数，除去公比为1的情况，每个等比数列的长度都不会很长，最长为log r。
这样我们可以直接求出所有可能的等比数列，不妨枚举其首项和公比，即可在O(r^2 log r)时间内找出所有等比数列，接下来O(n)递推即可，期望得分40。
算法四
矩阵乘法优化递推，复杂度降为O(log n log3 r + r2 log r)，期望得分50。
算法五
对于r <= 10^5, 可能存在一些复杂度不够优的找等比数列的方法，矩阵乘法优化的情况下可获得70分。
算法六
复杂度瓶颈在于找等比数列，考虑优化。
由于公比一定是有理数，不妨设其为x/y，由于公比大于1的等比数列的个数跟小于1的数量是完全相同的，不妨设x > y，设等比数列长为n，则an = x^(n-1)/y^(n-1)*a1，不妨令t = an^(x1-n) =
a1^y(1-n)，注意到一旦x,y确定，由于an = t x^(n-1) <= r; a1 = t y^(n-1) >= l，合法的t的个数为
对于n = 1; n = 2个数特殊处理，接下来再枚举x,y即可，注意可能会算重，还需要保证gcd(x,y)=1。注意到此时由于n>=3，我们只需要在sqrt(r)内枚举x,y即可。
总的复杂度为O(log n log3 r + r log r)，事实上还达不到这个复杂度。

接下来我来说说矩阵到底是这么构造的（std嫌问题太蠢啦根本不讲）【摊手
f(i)表示n==i时的ans。
显然有f(i) = f(i-1) * cnt1 + f(i-2) * cnt2 + f(i-3) * cnt3……
考虑构造矩阵：
首先肯定有：
cnt1
cnt2
cnt3
cnt4
…
…
0
（姑且叫做矩阵1）
最后要把我们的ans也就是f(n)放在0,0的位置（我们把ans矩阵称作矩阵2）。
而0,0是由矩阵2的第一行乘上矩阵1的第一列（反过来也可以，转一转矩阵也无所谓啦）。
所以矩阵2的第一行也就出来啦
ai,a(i-1) ,a(i-2) ,a(i-3) ,a(i-4),…1, 0, 0…0
要让矩阵2一直保证这种状态，就要对矩阵1做出一些改变了。
cnt1 1 0 0…0
cnt2 0 1 0…0
cnt3 0 0 1…0
cnt4 0 0 0…0
.… 0 0 0…0
.… 0 0 0…1
.0 0 0 0…0
好像就这么愉快的结束了一样。。。可是好像还忽略的一个问题，公比为一的我们还没处理呢。
所以刚才的f(i)是不包括公比为一的情况的，那么真正的f(i) = f(i-1) * cnt1 + f(i-2) * cnt2 + f(i-3) * cnt3…… + (r-l+1) （加上公比为一且长度为i的个数）。
所以矩阵1就变为了这样：
cnt1 1 0 0…0
cnt2 0 1 0…0
cnt3 0 0 1…0
cnt4 0 0 0…0
.… 0 0 0…0
.… 0 0 0…1
r-l+1 0 0 0…0
所以矩阵2就变为了这样：
ai,a(i-1) ,a(i-2) ,a(i-3) ,a(i-4),…1, 0, 0…1
所以矩阵1就变为了这样：
cnt1 1 0 0…0 0
cnt2 0 1 0…0 0
cnt3 0 0 1…0 0
cnt4 0 0 0…0 0
.… 0 0 0…0 0
.… 0 0 0…1 0
.0 0 0 0…0 0
r-l+1 0 0 0…0 1
我们可以发现最开始的矩阵2就是：
a1, a0, 0, 0, 0…1
而a1就是cnt1，a0就是1，这不就跟矩阵1的第一行一样一样的吗？
那么直接把矩阵1^n不就是ans矩阵啦！（反正矩阵2第一行下面是什么又不碍事）。
于是就愉快的做完啦:-）

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#include <cmath>
#include <ctime>
#define LL long long
#define Set(i, v) memset(i, v, sizeof i)
using namespace std;

const int mod = 1e9 + 7;
const int Mx = 24;//(24可过，不动脑子就开35)

struct Matrix{
    const static int N = 30;
    int M[N][N];
    Matrix(int v = 0){/*单位矩阵*/
        Set(M, 0);
        if(v == 1) for(int i=0; i<N; ++i) M[i][i] = 1;
    }
    Matrix operator * (const Matrix& B) const{
        Matrix C;
        for(int k=0; k<N; ++k)
            for(int i=0; i<N; ++i)
                if(M[i][k])
                    for(int j=0; j<N; ++j)
                        C.M[i][j] = (C.M[i][j] + 1LL * M[i][k] * B.M[k][j]) % mod;
        return C;
    }
    Matrix operator ^ (LL idx){
        Matrix ret(1)/*单位矩阵*/, A = *this;
        while(idx){
            if(idx & 1) ret = ret * A;
            A = A * A;
            idx >>= 1;
        }
        return ret;
    }
};

int l, r;
LL n;
LL cnt[Mx + 10];

int gcd(int aa, int bb){ return !bb ? aa : gcd(bb, aa % bb);}

LL work(){
    int R = sqrt(r) + 1;//枚举的范围 
    for(int i=0; i<=Mx; ++i) cnt[i] = 0;
    cnt[1] = r - l + 1;
    cnt[2] = 1LL * (r - l + 1) * (r - l) % mod;//单独处理一个数或两个数的情况 
    for(int i=1; i<=R; ++i) 
        for(int j=1; j<i; ++j) //因为之后要*2，姑且先不算公比为1的情况 
            if(gcd(i, j) == 1){//防止算重 
                int x = i, y = j;
                for(int k=2; 1LL*x*i<=r; ++k)//k = n - 1;
                    x *= i, y *= j, cnt[k+1] += max(r / x - (l - 1) / y, 0);
            }//an/x^(n-1) = a1/y^(n-1) => an/x = a1/y
            //cnt[k+1] = (r/(x/y) - l + 1) / y = r/x - (l-1)/y （需要乘出来是整数） 
    for(int i=3; i<=Mx; ++i) cnt[i] = (cnt[i] << 1) % mod;
    //公比大于1的等比数列的个数跟小于1的数量是完全相同的，所以*2 
    Matrix A;
    for(int i=1; i<=Mx; ++i) A.M[i - 1][0] = cnt[i];
    for(int i=1; i<Mx; ++i) A.M[i - 1][i] = 1;
    A.M[Mx][Mx] = A.M[0][Mx] = 1;
    A.M[Mx][0] = r - l + 1;//考虑公比为一的情况 
    Matrix ans;
    ans.M[0][0] = 1;
    return (ans * (A ^ n)).M[0][0];
}

int main(){
    freopen ("excellent.in", "r", stdin);
    freopen ("excellent.out", "w", stdout);
    int T;
    scanf("%d", &T);
    while(T--){
        scanf("%I64d%d%d", &n, &l, &r);
        printf("%I64d\n", work());
    }
    //printf("%.4lf\n", 1.0 * clock() / CLOCKS_PER_SEC);
    return 0;
}