# Linear Recurrence와 Kitamasa

Linear Recurrence는 앞의 몇 항으로 다음 항이 결정되는 수열입니다. Matrix Exponentiation으로도 풀 수 있지만, 차수 `K`가 크고 `N`이 매우 클 때는 characteristic polynomial을 이용하는 Kitamasa 방식이 더 직접적입니다.

## 문제 신호

| 문제 표현 | 접근 |
| --- | --- |
| `a_n = c1*a_{n-1} + ... + ck*a_{n-k}` | linear recurrence |
| `n`이 `10^18`처럼 크다 | fast exponentiation of polynomial |
| `K`가 수십~수천 | Kitamasa 후보 |
| 여러 n번째 항을 물어본다 | precompute/fast doubling 변형 고려 |
| 생성함수 분모가 주어진다 | Bostan-Mori 후보 |

`K`가 작으면 matrix exponentiation도 충분합니다. `K`가 커질수록 `K x K` 행렬 곱보다 polynomial reduction 관점이 유리해집니다.

## Characteristic Polynomial

점화식이 아래와 같다고 합시다.

```text
a_n = c0*a_{n-1} + c1*a_{n-2} + ... + c_{K-1}*a_{n-K}
```

그러면 companion relation은 아래입니다.

```text
x^K = c0*x^{K-1} + c1*x^{K-2} + ... + c_{K-1}
```

`x^n`을 이 관계로 계속 줄이면 `degree < K`인 polynomial이 됩니다.

```text
x^n mod P(x) = p0 + p1*x + ... + p_{K-1}*x^{K-1}
```

그럼 답은 아래처럼 계산됩니다.

```text
a_n = p0*a_0 + p1*a_1 + ... + p_{K-1}*a_{K-1}
```

## Kitamasa 구현

아래 구현은 `coeff[i]`가 `a_n`에서 `a_{n-i-1}`에 곱해지는 계수라는 convention을 사용합니다.

`K>=1`, initial.size()==coeff.size()==K, n>=0을 전제로 합니다. 공개 nth 함수는 계수를 정규화합니다. 내부 함수에는 정규화된 길이 K 벡터를 전달합니다. Kitamasa 자체는 역원이 필요 없으므로 합성수 modulus에도 적용할 수 있지만 BM은 field가 필요합니다.

> **코드 환경: 일반 C++17 학습용.** 헤더·STL을 허용하는 로컬 예제입니다. h-contest 제출에 옮길 때는 [공통 코드](https://h.readiz.com/learn/cpp-common-library)와 문제의 공개 API에 맞춰 필요한 부분을 바꿉니다.

```cpp compile-check
#include <vector>
using namespace std;

const long long MOD_RECURRENCE = 998244353;

long long normalizeMod(long long value) {
    value %= MOD_RECURRENCE;
    if (value < 0) {
        value += MOD_RECURRENCE;
    }
    return value;
}

vector<long long> combinePolynomial(
    const vector<long long>& left,
    const vector<long long>& right,
    const vector<long long>& coeff
) {
    int k = (int)coeff.size();
    vector<long long> temp(2 * k - 1, 0);

    for (int i = 0; i < k; ++i) {
        for (int j = 0; j < k; ++j) {
            temp[i + j] = (temp[i + j] + left[i] * right[j]) % MOD_RECURRENCE;
        }
    }

    for (int degree = 2 * k - 2; degree >= k; --degree) {
        long long value = temp[degree];
        if (value == 0) {
            continue;
        }
        for (int j = 1; j <= k; ++j) {
            temp[degree - j] = (temp[degree - j] + value * coeff[j - 1]) % MOD_RECURRENCE;
        }
    }

    temp.resize(k);
    return temp;
}

vector<long long> coefficientOfNthPower(long long n, const vector<long long>& coeff) {
    int k = (int)coeff.size();
    vector<long long> result(k, 0);
    vector<long long> base(k, 0);

    result[0] = 1;
    if (k == 1) {
        base[0] = coeff[0];
    } else {
        base[1] = 1;
    }

    while (n > 0) {
        if (n & 1LL) {
            result = combinePolynomial(result, base, coeff);
        }
        base = combinePolynomial(base, base, coeff);
        n >>= 1LL;
    }

    return result;
}

long long nthLinearRecurrence(
    const vector<long long>& initial,
    const vector<long long>& coeff,
    long long n
) {
    int k = (int)coeff.size();
    if (n < (long long)initial.size()) {
        return normalizeMod(initial[(int)n]);
    }

    vector<long long> normalizedCoeff = coeff;
    for (auto& value : normalizedCoeff) value = normalizeMod(value);
    vector<long long> weight = coefficientOfNthPower(n, normalizedCoeff);
    long long answer = 0;
    for (int i = 0; i < k; ++i) {
        answer = (answer + weight[i] * normalizeMod(initial[i])) % MOD_RECURRENCE;
    }
    return answer;
}
```

`coeff`와 `initial`의 길이는 같아야 합니다. 초기항은 `a_0..a_{K-1}` 순서입니다.

## 작은 예시

Fibonacci는 아래 점화식입니다.

```text
F_n = F_{n-1} + F_{n-2}
F_0 = 0, F_1 = 1
```

그러면 입력은 아래처럼 됩니다.

```text
initial = [0, 1]
coeff = [1, 1]
```

`x^n mod (x^2 - x - 1)`의 계수를 구한 뒤 `F_0`, `F_1`에 곱하면 `F_n`이 됩니다.

## Matrix Exponentiation과 비교

| 방식 | 시간 | 특징 |
| --- | ---: | --- |
| Matrix exponentiation | `O(K^3 log N)` | 구현 직관적, 전이 일반화 쉬움 |
| Kitamasa 기본형 | `O(K^2 log N)` | 선형 점화식 특화 |
| NTT 최적화 | `O(K log K log N)` 근처 | 구현 복잡 |
| Bostan-Mori | `O(K log K log N)` | 생성함수 분수 형태에 강함 |

문제가 단순 선형 점화식이면 Kitamasa가 깔끔합니다. 상태 전이가 sparse하거나 다른 구조가 있으면 행렬 방식이 더 읽기 쉬울 수 있습니다.

## Berlekamp-Massey와의 연결

처음 몇 항만 주어지고 점화식을 모르면 Berlekamp-Massey로 최소 선형 점화식을 추정할 수 있습니다.

```text
sequence prefix -> Berlekamp-Massey -> coeff -> Kitamasa nth term
```

다만 이 조합은 모듈러 field 위에서 동작합니다. 합성수 mod나 실수 근사 수열에는 그대로 적용하면 안 됩니다.

## 시간 복잡도

기본 구현은 polynomial 곱셈과 reduction에 `O(K^2)`가 들고, 거듭제곱에 `O(log N)`번 사용합니다.

```text
O(K^2 log N)
```

`K`가 5000 이상이면 이 구현도 부담될 수 있습니다. 그때는 NTT 기반 polynomial reduction이나 Bostan-Mori를 고려합니다.

## 상태 DP에서 점화식 찾기

길이 n에서 11을 포함하지 않는 이진 문자열 수는 a[0]=1,a[1]=2, 계수 [1,1]을 위 함수에 넣습니다. 끝 문자가 0/1인 두 상태의 선형 전이를 합치면 이 점화식이 나옵니다.

S개 상태의 고정 선형 전이 A에서 uᵀAⁿv는 Cayley-Hamilton에 의해 차수 S 이하 점화식을 가집니다. 이 상한을 알고 field 위에서 정확한 앞 2S항을 만들면 BM으로 복원할 수 있습니다. min/max는 일반 field 선형 전이가 아니지만 XOR는 GF(2) 덧셈입니다. 다항식 방정식의 생성함수가 항상 유리함수인 것은 아니며, 유한 선형 전이는 유리 생성함수로 연결됩니다.

## 대표 로컬 연습: K차 선형 점화식의 N번째 항

고정된 소수 mod `998244353`에서 아래 점화식을 따르는 수열의 `a_N`을 구합니다.

```text
a_n = c0*a_{n-1} + c1*a_{n-2} + ... + c_{K-1}*a_{n-K}
```

초기항은 `a_0, a_1, ..., a_{K-1}` 순서로 주어집니다.

### 입력

```text
K N
a_0 a_1 ... a_{K-1}
c0 c1 ... c_{K-1}
```

- `1 <= K <= 300`
- `0 <= N <= 10^18`
- 모든 항과 계수는 `0 <= value < 998244353`

### 출력

```text
a_N mod 998244353
```

### 예시

```text
2 10
0 1
1 1
```

```text
55
```

Fibonacci를 `F_0 = 0`, `F_1 = 1`, `F_n = F_{n-1} + F_{n-2}`로 둔 입력입니다.

## 손으로 따라가는 Trace

Fibonacci의 companion relation은 아래입니다.

```text
x^2 = x + 1
```

따라서 `x^N mod (x^2 - x - 1)`을 `p0 + p1*x` 꼴로 줄이면:

```text
F_N = p0*F_0 + p1*F_1
```

`N = 10`일 때 필요한 거듭제곱은 아래처럼 줄어듭니다.

| 항 | 줄이기 전 | relation 적용 후 | 계수 `[p0, p1]` |
| --- | --- | --- | --- |
| `x^1` | `x` | `x` | `[0, 1]` |
| `x^2` | `x^2` | `x + 1` | `[1, 1]` |
| `x^4` | `(x + 1)^2 = x^2 + 2x + 1` | `3x + 2` | `[2, 3]` |
| `x^8` | `(3x + 2)^2 = 9x^2 + 12x + 4` | `21x + 13` | `[13, 21]` |
| `x^10` | `x^8 * x^2 = (21x + 13)(x + 1)` | `55x + 34` | `[34, 55]` |

마지막 계수 `[34, 55]`로 `34*F_0 + 55*F_1 = 55`가 됩니다.

## 구현 기준

```cpp
#include <iostream>
#include <vector>
using namespace std;

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    int k;
    long long n;
    cin >> k >> n;

    vector<long long> initial(k), coeff(k);
    for (long long& value : initial) {
        cin >> value;
        value = normalizeMod(value);
    }
    for (long long& value : coeff) {
        cin >> value;
        value = normalizeMod(value);
    }

    cout << nthLinearRecurrence(initial, coeff, n) << '\n';
}
```

## 검증용 Case

| 입력 요약 | 기대값 | 확인 포인트 |
| --- | ---: | --- |
| `K=2`, Fibonacci, `N=0` | 0 | 초기항을 바로 반환 |
| `K=2`, Fibonacci, `N=1` | 1 | 초기항 index가 0-based |
| `K=2`, Fibonacci, `N=10` | 55 | coefficient trace와 일치 |
| `K=1`, `a_n=3a_{n-1}`, `a_0=2`, `N=4` | 162 | `k=1`에서 `x mod P(x)=c0` 처리 |
| `K=3`, tribonacci `0,0,1`, coeff `1,1,1`, `N=5` | 4 | reduction 방향 검증 |

## Stress 기준

작은 입력에서는 직접 점화식을 전개하는 naive 구현과 비교합니다.

1. `K <= 6`, `N <= 80`으로 random initial/coeff를 생성합니다.
2. naive로 `a_K`부터 `a_N`까지 순서대로 계산합니다.
3. 같은 입력을 Kitamasa 구현에 넣어 결과가 같은지 비교합니다.
4. `K=1`, `N<K`, 계수가 0인 경우, 모든 초기항이 0인 경우를 별도 deterministic case로 둡니다.

Berlekamp-Massey와 연결할 때는 BM이 찾은 coeff로 이 연습의 `nthLinearRecurrence`을 호출하고, BM에 쓰지 않은 holdout 항을 하나 더 비교해야 합니다.