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;}