菜猫喝水

【TiTle】

身为全电大最聪明的猫，菜猫他有 $n$ 个不同水杯，为了方便喝水，菜猫给第 $i$ 号水杯里面盛上i毫升水。

【Content】

有一天，猫口渴了。他决定喝一些水来解渴。因为活得太久，喝水在猫看来也应该是一种艺术，不同的水杯里的水一起喝可口程度不同。所以他决定在他的 $n$ 个水杯中挑出可口度最大的 $k$ 个水杯来喝水。这 $k$ 个水杯的可口度是他们的盛水量的最大公约数。

【Standard Input】

两个空格分开的正整数 $n$ 和 $k$。

【Standard Output】

一个整数，为最大的可口度。

输入样例

4 2

输出样例

2

【样例解释】

菜猫一共有 $4$ 个杯子，里面分别装了 $1$ mL，$2$ mL，$3$ mL 和 $4$ mL 水，现在要挑出两杯水，使得这些水的毫升数的最大公约数最大。可以挑选$2$ mL和 $4$ mL，这样最大公约数为 $2$，最为可口。但假如挑选的是 $1$ mL 和 $3$ mL，最大公约数仅为 $1$，不是很可口。所以可口度最大为 $2$.
对于 100% 的数据，$k$，$n$ $\leq$ $2\cdot 10^9$

【示例代码】

#include<bits/stdc++.h>
using namespace std;
const int mo = 1e9+7,maxn = 12;
int n,col[3],cha[3] = {0},id = 0,has[99][99],m,cnt = 0,tot1[999],tot2[999];
char a[maxn + 1],b[maxn + 1],t;
struct node{
    long long v[80][80];
}nil,ans,f;
node operator * (node a,node b){
    node c = nil;
    long long t ;
    for (int i = 0 ; i <= id; ++i)
        for (int j = 0 ; j <= id ; ++j){
            t = 0;
            for (int k = 0 ; k <= id ; ++k) t = (t + a.v[i][k] * b.v[k][j]) % mo;
            c.v[i][j] = t;
        }
    return c;
}
int main(){
    scanf("%s%s",a,b);  n = strlen(a);
    for (int i = 0 ; i < 3 ; ++i) scanf("%d",&col[i]);
    scanf("%d",&m);
    for (int i = 0 ; i < n ; ++i){
        cha[(b[i] - a[i] + 3) % 3] ++;
        t = a[i];
        while (t != b[i]){
            ++cnt; m -= col[t - 'a'];
            ++t;
            if (t > 'c') t -= 3;
        }
    }
    if (m < 0){puts("0"); return 0;}
    m = m / (col[0] + col[1] + col[2]) * 3 + cnt;
    
    for (int i = 0 ; i <= n ; ++i)
        for (int j = 0 ; i + j <= n ; ++j){
            has[i][j] = ++id;
            tot1[id] = i;
            tot2[id] = j;
        }
    memset(nil.v,0,sizeof(nil.v));
    ans = f = nil;
    ans.v[ has[cha[1]][cha[2]] ][0] = 1;
    int x,y;
    for (int i = 1 ; i <= id ; ++i){
        x = tot1[i]; y = tot2[i];
        if (n - x - y) f.v[i][has[x + 1][y]] = x + 1;
        if (x) f.v[i][has[x - 1][y + 1]] = y + 1;
        if (y) f.v[i][has[x][y - 1]] = n - x - y + 1;
    }
    f.v[0][0] = 1; f.v[0][1] = 1;
    for (;m;m>>=1,f=f*f) if (m &1) ans = f * ans;
    cout<<(ans.v[0][0] + ans.v[1][0]) % mo; 
    return 0;
}