Top
高速フーリエ変換の実装は AtCoder のやつが暗唱できそうなレベルで分かりやすい。
それを元にした実装が左だが、これは誤差が怖く、整数剰余環を使った右側の実装がよく知られている。
NTT とか FMT とか言われていて、NTT という言葉の方がメジャーっぽいが個人的には有名な企業名とかぶらない FMT を推したい。
#!/usr/bin/python2
# -*- coding: utf-8 -*-
# †
from math import cos, sin, pi


def myfft(a, inv):
    n = len(a)
    if n == 1:
        return a
    m = n / 2
    a0 = myfft(a[0::2], inv)
    a1 = myfft(a[1::2], inv)
    theta = 2 * pi / n
    if inv:
        theta *= -1
    z = complex(cos(theta), sin(theta))
    pz = 1
    r = [None] * n
    for i in xrange(n):
        r[i] = a0[i%m] + pz * a1[i%m]
        pz *= z
    return r

def fft(a):
    return myfft(a, False)

def ifft(a):
    n = len(a)
    b = myfft(a, True)

    for i in xrange(n):
        b[i] /= n
    return b

def convol(a, b):
    n = 1
    while n < len(a) + len(b):
        n <<= 1
    a += [0] * (n - len(a))
    b += [0] * (n - len(b))
    a = fft(a)
    b = fft(b)
    c = [ai * bi for ai, bi in zip(a, b)]
    c = ifft(c)
    return c
          
#!/usr/bin/python2
# -*- coding: utf-8 -*-
# †
O = 3 # 原始根 \omega
M = (5 << 37) * 211 + 1

def myfmt(a, inv):
    n = len(a)
    if n == 1:
        return a
    m = n / 2
    a0 = myfmt(a[0::2], inv)
    a1 = myfmt(a[1::2], inv)
    z = pow(O, (M-1)/n, M)
    if inv:
        z = pow(z, M-2, M) # 乗法逆元を求めている
    pz = 1
    r = [None] * n
    for i in xrange(n):
        r[i] = a0[i%m] + pz * a1[i%m]
        r[i] %= M
        pz = (pz * z) % M
    return r

def fmt(a):
    return myfmt(a, False)

def ifmt(a):
    n = len(a)
    b = myfmt(a, True)
    inv = pow(n, M-2, M) # 乗法逆元を求めている
    for i in xrange(n):
        b[i] = b[i] * inv % M
    return b

def convol(a, b):
    n = 1
    while n < len(a) + len(b):
        n <<= 1
    a += [0] * (n - len(a))
    b += [0] * (n - len(b))
    a = fmt(a)
    b = fmt(b)
    c = [ai * bi % M for ai, bi in zip(a, b)]
    c = ifmt(c)
    return c
          

羃剰余の計算が重いが、必要となるパラメータのバリエーションは少ないのでメモ化が有効。
暗唱するにあたっての障壁は $\omega$ とか M とかの値をどうやって覚えとくかか…。

(ΦωΦ)<おしまい