선형 점화식

• 점화식의 종류는 무궁무진하지만 제일 기본으로 다루는 것은 선형점화식으로, a[n] = c[1]*a[n-1] + c[2]*a[n-2] + … + c[k]*a[n-k] + f(n) 형태로 표현되는 점화식이다. 이때 c[i]는 모두 상수이고, f는 산술적 함수이다.
• 선형점화식에서 f(n)=0 일 경우, 동차선형점화식(homogeneous linear recurrence relation) 이라고 부른다. 그렇지 않는 경우는 비동차선형점화식이라고 부른다
• 조합론에서 점화식을 푼다고 표현하는 것은, 점화식을 일반항으로 표현하는 것을 말한다. 간단하게 요약하면
  • 동차선형점화식에 대해서는 특성방적식을 이용하는 일반화된 풀이법이 존재한다. 이것이 가장 기본이다.
  • 비동차선형점화식의 경우는 특수해을 구하여 선형점화식으로 변형해서 푼다. 특수해 b[n]을 구하면 x[n]=a[n]-b[n]으로 새로운 수열을 정의해서 x[n]을 동차점화식으로 표현할수 있다. f(n)=d, f(n)=dn, f(n)=dn^2, f(n)=c*n^d 등등의 형태에 대해서 특수해를 구하는 법들이 알려져있지만 일반적인 방법은 없는것 같다.
  • 선형 점화식이 아닌 경우에도 이리저리 잘 치환해서 선형점화식으로 변형해서 푸는 연습을 한다
• 하지만, PS에서의 목표는 점화식을 푸는것이 아니다. 보통 원하는 것은 n번째 항, a[n]의 값을 계산하는 것인데, 이때 점화식을 풀어서 일반항을 구한뒤 대입한다는 방법은 별로 효과적이지 못하다. 예를들어 피보나치수의 일반항은 $ a_n = \frac{1}{\sqrt5} \Bigg(\bigg( \frac{1+\sqrt5}{2} \bigg)^n -\bigg( \frac{1-\sqrt5}{2} \bigg)^n \Bigg) $ 인데, 실제로 이 식에 n을 대입해서 값을 구하려고 시도하는 것은 실수 오차로 원하는 값을 얻을수 없다.
• 따라서 PS에서는 일반항을 구하는것이 아니라, n번째 항을 효율적으로 구하는 방법을 다루게 된다. 조합론에서와 마찬가지로 동차선형점화식을 해결하는 것이 기본 목표이다. 비동차선형점화식이나 비선형점화식의 경우 조합론에서 사용하는 방법과 같은 방식으로 동차선형점화식으로 변환해서 풀수 있을 것 같긴 하지만, 실제 문제에서 여기까지 필요한 경우는 거의 본적이 없다..

• DP에서 얻어진 점화식이 dp[n] = c[1]*dp[n-1] + c[2]*dp[n-2] + … + c[k]*dp[n-k] 형태일 경우. 이런 형태의 점화식은 상수계수의 선형 점화식 이라고 불린다
  • 대표적인 예로 피보나치 수열이 있다
• 일반적인 DP를 이용하면 O(n) (정확히는 O(nk)) 시간에 풀 수 있지만, 더 빠르게 푸는 방법들이 존재한다.
• 가장 이해하기 간단한 방법은 행렬의 거듭제곱을 이용해서 O(logn), (정확히는 O(k^3logn))에 푸는 방법이다. k가 작은 경우에는 이것으로 충분하다.
• k가 큰 경우를 해결하기 위해 원래 알려져 있던 방법은 키타마사법이었다. 나이브하게 구현하면 O(k^2*logn)에 동작하게 할수 있다. k가 많이 큰 경우에는 내부 로직에 포함된 다항식 곱셈과 다항식 나눗셈을 FFT를 이용해서 계산하도록 구현하면 O(k*logk*logn)에 동작하게 된다.
• (!) 하지만! 이제 키타마사법은 잊어버려도 된다. 키타마사법보다 갖은 시간복잡도를 가지면서도, 상수가 작아서 더 빠르게 돌아가는 '보스탄-모리 알고리즘' 이 2020년에 발표되었기 때문이다!
  • 빠를 뿐더러 구현하기에도 더 간단한데, 다항식 곱셈과 다항식 나눗셈이 모두 필요한 키타마사법에 비해, 보스탄-모리는 다항식 곱셈만 필요로 하기 때문이다. 다항식 곱셈만이라면 FFT를 이용해서 구하는 것도 간단하다. 다항식 나눗셈을 FFT로 구현하려면 좀더 복잡하다.
• 보스탄-모리로 구현할때, 대충 k가 60 이하면 FFT를 안쓰는게 더 빠르고, 그 이상이면 FFT를 쓰는것이 더 빠른것 같다
• 아래는 보스탄-모리 법을 알기 전에 작성했던 키타마사법 관련 내용인데.. 지우기는 뭐해서 보관해둔다

.

• 행렬의 거듭제곱을 이용해서 O(logn)에 푸는 방법이 있다.
• 인접 3항간의 점화식, dp[n] = a*dp[n-1] + b*dp[n-2] 에 대한 경우를 예로 들면
  • $ \begin{bmatrix}dp[n] \\ dp[n-1] \end{bmatrix}=\begin{bmatrix}a & b \\1 & 0 \end{bmatrix}\begin{bmatrix}dp[n-1] \\ dp[n-2] \end{bmatrix} $ 으로 쓸 수 있으니까
  • $ \begin{bmatrix}dp[n] \\ dp[n-1] \end{bmatrix}=\begin{bmatrix}a & b \\1 & 0 \end{bmatrix}^{n-1}\begin{bmatrix}dp[1] \\ dp[0] \end{bmatrix}$ 이 된다
• 중간에 있는 행렬의 (n-1)승을 구하는 것은, 거듭제곱의 빠른 계산 (Exponentiation by squaring)을 이용해서 O(logn)에 구할 수 있다.
• 더 많은 인접항에 대한 점화식도 같은 방식으로 풀 수 있다.
• 그리고 상수항이 포함된 선형비동차점화식, dp[n] = c[1]*dp[n-1] + c[2]*dp[n-2] + … + c[k]*a[n-k] + c0 도, 행렬로 모델링이 가능하기 때문에, 이 방식을 적용할 수 있다.
  • 아래에서 설명할 키타마사법이나 보스탄-모리알고리즘으로는 이 형태를 바로 풀지 못하므로 선형동차점화식으로 변형해주는 과정이 필요한데 (약간 확신이 없다. 확인 필요;), 행렬곱 방식은 그럴 필요가 없다.
• 코드
  • teflib.combinatorics.linear_homogeneous_recurrence
• 같은 점화식에 대해서 여러개의 n값을 처리해야 할 경우에는, M^2, M^4, M^8, … 행렬을 매번 만드는 대신, 한번만 만들어놓고 쿼리마다 재사용하는 편이 더 효율적이다. (teflib에 구현되어 있지는 않다)

• 위의 식에서 k가 커지면, 행렬 계산에 드는 시간도 고려해야 한다. kxk 행렬을 단순히 곱하는 데에는 O(k^3)이므로, 시간복잡도가 O(k^3 * logn)이다.
  • 그냥 행렬을 안쓰고 계산하면 O(k*n)
• 키타마사법은 이것을 O(k^2 * logn)으로, 더 나아가서는 O(klogk * logn)으로 줄일 수 있는 방법을 제공한다.
• 키타마사법의 기본은 A[n]이 A[n-1], …, A[n-k]의 선형 결합으로 표현되었던 것을, A[k], A[k-1], …, A[1]의 선형 결합으로 바꿔서 표현하는 것이다. 즉 아래와 같은 형태로 변환되는 새로운 계수들 d의 값을 구하고, 처음 k항의 값을 대입하면 바로 n번째 항의 값을 계산할 수 있다.
  • $ A[n] = c[1]A[n-1] + c[2]A[n-2] + ... + c[k]A[n-k] \\ = ... \\ = d[1]A[k] + d[2]A[k-1] + ... + d[k]A[1] $
• 여기서 d[i]의 값들을 구하기 위한 중간과정이 필요한데.. 처음에 A[n]을 A[n-1]~A[n-k]의 선형결합으로 표현한 식에서, A[n-1]=c[1]A[n-2]+…+c[k]A[n-k-1]을 대입해주면 A[n]을 A[n-2]~A[n-k-1]의 선형결합으로 바꿀수 있다. 이것을 n번 반복하면 이론적으로 A[n]을 A[k], A[k-1], …, A[1]의 선형결합으로 바꿀수 있다. 하지만 이렇게 하면 결국 n번의 연산이 필요하니 이것을 효율적으로 하는 방법을 생각해보자
• 효율적으로 계산할수 있는 형태로 바꾸는 방법을 설명하는 방식이 크게 두가지가 있는데
  • 첫번째는 그냥 진짜로 식을 잘 정리하는 것이다. 과정을 따라가는게 아주 어려운것은 아니지만, 기억하기에는 불편하기는 하다.. 이 방식의 설명은 https://casterian.net/archives/254을 참고.
  • 두번째는, 위에서 처럼 대입해서 차수를 한단계식 낮추는 것이, 다항식의 나눗셈을 하는 과정과 완전히 동일하다는 점을 이용하는 것이다. https://justicehui.github.io/hard-algorithm/2021/03/13/kitamasa/ 에 잘 설명이 되어있는데, 저 과정을 실제로 따라가보면, 위의 방법은 $x^n$ 을 $f(x) = x^k - c_1\cdot x^{k-1} - c_2\cdot x^{k-2} - ... - c_k\cdot x^0 $ 으로 나눈 나머지를 구하는 것과 동일하다. 이렇게 접근하면 이해도 편하고 기억도 쉬우니 이쪽으로 생각하자.
• 그럼 결국 $x^n mod f(x)$ 을 계산하면 새로운 계수들 d[i]를 구할 수 있다. x^n mod f(x) 는 거듭제곱의 빠른 계산 (Exponentiation by squaring) 방법으로 logn번의 곱셈과 나머지 연산으로 계산 가능하다. (헷갈리면 그냥 정수에 대해서 a^n % b 를 거듭제곱의 빠른 계산 (Exponentiation by squaring)으로 구한다고 생각하자).
• f(x)의 차수가 k이고 이걸로 계속 나머지를 취해주니까, 이 과정에서 곱셈은 k차 다항식끼리, 나눗셈은 2k차 다항식을 k차 다항식으로 나누는 형식이 된다. 그냥 나이브하게 다항식 곱셈과 나눗셈을 구현하면 각각 O(k^2)에 처리된다. 총 시간복잡도는 O(k^2 * logn)
• 그리고, 고속 푸리에 변환 (Fast Fourier Transform, FFT) 를 이용해서 다항식 곱셈과 나눗셈을 구현하면 각각 O(klogk)에 처리가능하고, 총 시간복잡도는 O(klogklogn). 이것을 따로 '고속 키타마사법' 이라고 부르는 경우도 있다.
• '키타마사법'이라는 이름은 PS 커뮤니티에서만 불리는 명칭인것 같다. 일본의 PS 커뮤니티에서 유래된 것 같은데, 그래서 그 'きたまさ'라는 사람이 대체 누구인지는 이리저리 찾아봤지만 알아내는 데에 실패했다. 대신 알게된 것은 이 방법이 처음으로 제안된 것은 Fiduccia의 1985년 논문인 An Efficient Formula for Linear Recurrences 이라는 것. 정확히 부르려면 Fiduccia's algorithm 이라고 부르는 것이 올바를 것 같지만 그러기에는 이미 늦은듯.
• 코드는 대충 이런식이다,

  • def linear_recurrence_kitamasa(coef, seed, n, mod):
    
        def polymul(a, b):
            ret = [0] * (len(a) + len(b) - 1)
            for i, a_i in enumerate(a):
                for j, b_j in enumerate(b):
                    ret[i + j] += a_i * b_j
            return [x % MOD for x in ret]
    
        def polymod(a, b):
            ret = a[:]
            for i in reversed(range(len(b) - 1, len(a))):
                k = i - len(b) + 1
                for j, b_j in enumerate(b):
                    ret[k + j] -= ret[i] * b_j
                ret = [x % MOD for x in ret]
            return ret[:len(b)]
    
        d = [1]
        xn = [0, 1]
        f = [-c for c in reversed(coef)] + [1]
        while n:
            if n & 1:
                d = polymod(polymul(d, xn), f)
            n >>= 1
            xn = polymod(polymul(xn, xn), f)
    
        return sum(s_i * d_i for s_i, d_i in zip(seed, d)) % mod

  • 다항식 곱셈과 나눗셈 함수를 각각 나이브하게 구현해서 사용한 버전으로, 저 함수들은 FFT를 사용해서 구현할수 있다.

• Bostan-Mori algorithm
• 키타마사법과 마찬가지로 k차 선형 점화식의 n번째 항을 O(k^2 * logn) 또는 O(klogk * logn)에 구할 수 있는 방법이다.
• 비교적 최근 (2020년)에 나온 방법이라서, 아직도 키타마사법에 비해 덜 알려진 편이지만, 실행속도나 구현난이도 면에서 모두 키타마사법보다 우위에 있다.
• 참고자료
  • https://blog.myungwoo.kr/149,
  • https://infossm.github.io/blog/2021/12/21/Bostan-Mori_Algorithm/
  • https://arxiv.org/abs/2008.08822 (원본논문)
• 설명은 위의 링크들을 참조. 전명우 블로그쪽의 설명이 간명한 편이라 처음 이해하기에는 그쪽이 좋다. 삼소멤 블로그의 글은 좀더 방대한데, FFT를 사용해서 처리할 때 효율적으로 하는 방법에 대한 설명들이 있고, 키타마사법과의 속도 비교 벤치마크까지 있다.
• 그냥 기본적으로 구현하면 상당히 간단하다. 대충 이런식이 된다

  • def linear_recurrence(coef, seed, n, mod):
        d = len(coef)
        convolution_func = _naive_convolution if d <= FFT_THRE else fft.multiply
        q = [1] + [(-x) % mod for x in coef]
        p = [x % mod for x in convolution_func(seed, q)[:d]]
        while n:
            r = q[:]
            r[1::2] = [mod - x for x in q[1::2]]
            p = [x % mod for x in convolution_func(p, r)[n & 1 :: 2]]
            q = [x % mod for x in convolution_func(q, r)[::2]]
            n >>= 1
        return p[0]

• d가 작으면 그냥 convolution을 나이브하게 구현해서 O(k^2logn)에 푸는것이 낫고, d가 크면 FFT를 이용한 convolution으로 O(klogklogn)에 계산하는 것이 효율적이다. 기준은 대충 60정도로 잡았다 (실험적으로 얻은 것은 아니고 다른 코드에 적혀있던 값을 참고했다)
• 나이브하게 컨볼루션을 계산하는 코드는 이런식이다

  • def _naive_convolution(a, b):
        ret = [0] * (len(a) + len(b) - 1)
        for i, a_i in enumerate(a):
            for j, b_j in enumerate(b):
                ret[i + j] += a_i * b_j
        return ret

• FFT를 쓰게 될 경우 직접 구현하게 된다면 거의 항상 NTT를 써서 구현해야 한다. 구하는 값이 커지므로 문제 자체에서 모듈러를 취한 값을 구하도록 출제되기 때문이다. (문제에서 요구하는 모듈러가 NTT 친화적인 소수가 아니라면 좀 이야기가 복잡해진다..). 그리고 NTT를 구현할때에도 DFT doubling 등의 테크닉을 써서 더욱 최적화 하는 것도 가능한 것 같다. 하지만, 우리는 어차피 Teflib FFT를 사용할 것이기 때문에 이런 것들은 신경쓸 필요가 없다. 오히려 이 경우에 귀찮은 것은 음수 계수를 포함한 다항식의 곱셈을 Teflib FFT로 계산하기가 번거롭다는 것인데, 어차피 모듈러를 취한 값을 구해야 한다면 이 문제는 깔끔하게 처리된다.
• 실제로 구현해서 실행한 결과는 k가 클때는 물론이고, k=2나 3처럼 작은 경우에도 행렬 거듭제곱을 이용한 방법보다 빠르게 동작했다. 그래서 결국 이전까지 썼왔던 행렬 거듭제곱법의 구현체인 teflib의 linear_homogeneous_recurrence 함수는 이 더이상 사용할 필요가 없어졌다. 내부 구현을 바꾸려다가, 아예 이 기회에 함수 이름도 좀더 짧게 줄이기로 했다. linear_homogeneous_recurrence 대신에 linear_recurrence 로 새로 함수를 만들었다.

• teflib.combinatorics.linear_recurrence
• teflib.combinatorics.linear_recurrence_large

• 타일 채우기 2: 차수가 2차인 선형 점화식. 행렬 거듭제곱을 이용해도 충분히 풀수 있다.
• 15572: 차수가 최대 1,000인 선형 점화식. O(k^2logn) 에 풀어야 한다
• 13725: 차수가 최대 30,000인 선형 점화식. O(klogklogn) 에 풀어야 한다

• 전이를 행렬형태로 표현하면, 위에서 설명한 행렬 거듭제곱을 이용한 풀이를 적용할 수 있다.

• 벌리캠프-매시 알고리즘을 통해서 선형 동차 점하식을 찾아낼수 있다.
• Berlekamp의 통일된 한국어 표기가 없어서 다양한 표기로 알려졌고, PS쪽에는 그중에서도 벌레캠프라는 발음이 널리 퍼져버려서 어떻게 사람 이름이 벌레캠프? 이런 뻘글이 유행하기도 했다ㅠㅠ 한국어 위키에서는 벌러캠프, 백준 태그 기준으로는 벌리캠프로 표기되어있고, 바둑의 끝내기에 관련된 연구때문에 바둑관련 기사에도 가끔 등장하는데 그쪽에서는 벌캄프가 많이 쓰인다. 여기에서는 벌리캠프로 쓰기로 하겠다.
• 알고리즘 설명에 관해서는 구사과님의 글 이 거의 완벽한 레퍼런스이다. 배경지식/기본 사용법/활용법 등에 대해서 자세히 정리가 되어있고, '내부 구조는 나도 모른다', '코드 그냥 복붙해서 쓰면 된다' 라는 문장으로 인해서, 모든 사람들이 '구사과도 원리 이해 못하고 쓴대'라는 완벽한 근거로 템플릿만 복붙해서 쓰는 것을 마음 편히 할수 있게 되었다.
• 이 알고리즘이 해주는 것은, 주어진 점화식의 처음 항들로부터 해당되는 최소 길이의 선형 동차 점화식을 찾는 것이다. 시간복잡도는 n개의 초기항을 넘겨주었을 경우 O(n^2)이다 (정확히는 O(nk+nlogmod) 이라고 한다).
• 알고리즘을 사용할때 몇가지 조건이 있는데
  • 찾으려는 점화식의 차수가 k라면, 최소 3k개의 초기항을 주어야 한다
  • 모듈러가 소수여야 한다
• 이 알고리즘은 활용처가 꽤 많은데, 생각보다 많은 식들이 선형 동차 점화식으로 표현될수 있기 때문이다.
  • 상수항을 추가로 갖는 선형비동차점화식도 선형동차점화식으로 변환 가능
  • 연립선형점화식도 선형동차점화식으로 변환가능
  • 일반적으로 상태 전이가 행렬로 표현가능 하면 선형동차점화식으로 표현가능
  • k차 다항식 f(x)를 k차 선형동차점화식으로 표현 가능. f(n)이 점화식의 n번째 항이 된다.
    • 물론 라그랑주 인터폴레이션으로 f(x)를 바로 찾아서 f(n)을 계산할 수도 있지만, 따로 코드 짜기가 귀찮으면 점화식으로 접근해도 된다
• 구현
  • 구사과님의 가르침에 따라, 나도 구현은 기존 코드를 복붙해서 쓰도록 했다. 구현된 기존 코드들은 찾아보면 많이 나온다.
  • 약간의 최적화와, 편의성 개선을 추가했다. 찾아진 점화식의 크기가 넣어준 초기항 개수의 1/3보다 큰경우에는 점화식을 신뢰할수 없으므로 예외를 발생시킨다든가..
  • 코드 텦립에 구현.
    • teflib.combinatorics.find_linear_recurrence
  • 코드를 사용하는 방식은 문제에 따라 다르다. 점화식이 입력값에 무관하게 일정한 경우에는, 로컬에서 알고리즘을 돌려서 구한 점화식을 그냥 코드에 박아 넣는것이 가능하다. 입력에 따라 점화식이 달라져서, 런타임에서 직접 점화식을 찾아야 하는 경우에는 이런식으로 쓴다.

    • seeds = compute_seeds()
      coefs = combinatorics.find_linear_recurrence(seeds, MOD)
      answer = combinatorics.linear_recurrence(coefs, seeds, n, MOD)