【NOI2011】兔农（循环节）

我居然没看题解瞎搞出来了？

题解：

---

不难想到找到每次减1的位置，然后减去它对最终答案的贡献。

假设有一个地方是\(x,1(mod~k)\)

那么减了1后就变成了\(x,0\)。

然后可以推到\(x,0,x,x\)。

可以看做以\(x,x\)为开头，做新的序列。

设y为x在mod k意义下的逆元，那么下次1的地方就是原斐波拉契序列中第一次出现y的位置。

如果没有逆元那就结束了。

原斐波拉契序列的非循环长度是\(O(k)\)级别的，这个可以通过随机序列第一次出现相同元素来理解，为在mo意义下应该是可以看做随机的，值域大小是\(O(k^2)\)，那么期望长度就是\(O(\sqrt {k^2})\)。

所以斐波拉契序列做个3k左右然后预处理每一个元素第一次出现的位置。

然后x也是有循环的，因为x<k，所以肯定是\(O(k)\)的，那么就可以一起做了。

最后就是个等比矩阵求和，由于不一定有逆，所以要用分治法去求等比矩阵和。

由于循环开始之前和最后结尾的地方都要判，所以有点复杂。

Code：

---

#include<bits/stdc++.h>
#define fo(i, x, y) for(int i = x, B = y; i <= B; i ++)
#define ff(i, x, y) for(int i = x, B = y; i <  B; i ++)
#define fd(i, x, y) for(int i = x, B = y; i >= B; i --)
#define ll long long
#define pp printf
#define hh pp("\n")
using namespace std;

ll n; int k, mo;

const int N = 3e6 + 5;

int f[N], fi[N];
int p[N], p0, us[N];

int gcd(int x, int y) { return (!y ? x : gcd(y, x % y));}
void exgcd(int a, int b, int &x, int &y) {
    if(!b) {x = a, y = 0; return;}
    exgcd(b, a % b, y, x); y -= (a / b) * x;
}
int qni(int p, int q) {
    int x, y;
    exgcd(p, q, x, y);
    x = (x % q + q) % q;
    return x;
}

int l = -1, r;

void work() {
    int x = 1;
    while(1) {
        if(gcd(x, k) > 1) break;
        int y = qni(x, k);
        if(fi[y]) {
            p0 ++;
            p[p0] = p[p0 - 1] + fi[y];
            if(us[y]) {
                l = us[y]; r = p0;
                break;
            }
            us[y] = p0;
            x = (ll) f[fi[y] - 1] * x % k;
        } else break;
    }
}

struct jz {
    ll a[2][2];
    jz() {
        a[0][0] = a[1][1] = 1;
        a[0][1] = a[1][0] = 0;
    }
};

jz operator * (jz a, jz b) {
    jz c;
    fo(i, 0, 1) fo(j, 0, 1) 
        c.a[i][j] = (a.a[i][0] * b.a[0][j] + a.a[i][1] * b.a[1][j]) % mo;
    return c;
}

jz operator + (jz a, jz b) {
    fo(i, 0, 1) fo(j, 0, 1) a.a[i][j] = (a.a[i][j] + b.a[i][j]) % mo;
    return a;
}

jz operator - (jz a, jz b) {
    fo(i, 0, 1) fo(j, 0, 1) a.a[i][j] = (a.a[i][j] - b.a[i][j] + mo) % mo;
    return a;
}

jz ksm(jz x, ll y) {
    jz s;
    for(; y; y /= 2, x = x * x)
        if(y & 1) s = s * x;
    return s;
}

jz z;

jz ksb(jz x, ll y) {
    if(y == 1) return x;
    if(y & 1) return ksb(x, y - 1) * x + x;
    jz a = ksb(x, y / 2);
    return a + a * ksm(x, y / 2);
}
jz calc(jz x, ll y) {
    jz s = jz();
    if(y > 0) s = s + ksb(x, y);
    return s;
}

int main() {
    scanf("%lld %d %d", &n, &k, &mo);
    f[1] = f[2] = 1;
    fo(i, 3, 3000000) {
        f[i] = (f[i - 2] + f[i - 1]) % k;
        if(!fi[f[i]]) fi[f[i]] = i;
    }
    work();
    z.a[0][1] = z.a[1][0] = z.a[1][1] = 1; z.a[0][0] = 0;
    ll ans = ksm(z, n).a[1][0];
    if(l == -1) {
        fo(i, 1, p0) if(p[i] <= n)
            ans = (ans - ksm(z, n - p[i]).a[1][1] + mo) % mo;
    } else {
        fo(i, 1, l) if(p[i] <= n)
            ans = (ans - ksm(z, n - p[i]).a[1][1] + mo) % mo;
        if(p[l] <= n) {
            ll v = (n - p[l]) / (p[r] - p[l]);
            fo(i, l + 1, r) {
                ll n2 = (ll) p[i] + v * (p[r] - p[l]);
                if(n2 <= n) ans = (ans - ksm(z, n - n2).a[1][1] + mo) % mo;
            }
            if(v > 0) {
                ll st = v * (p[r] - p[l]) + p[l];
                jz sa; memset(sa.a, 0, sizeof sa.a);
                fo(i, l + 1, r) sa = sa + ksm(z, p[r] - p[i]);
                jz c = ksm(z, n - st) * sa * calc(ksm(z, p[r] - p[l]), v - 1);
                ans = (ans - c.a[1][1] + mo) % mo;
            }
        }
    }
    pp("%lld\n", ans);
}

转载于:https://www.cnblogs.com/coldchair/p/11385378.html