math - 在有限域上实现 FFT

Question

我想使用 NTT 实现多项式乘法。我关注了Number-theoretic transform (integer DFT)而且它似乎有效。

现在我想在有限域Z_p[x]上实现多项式的乘法，其中p是任意素数。

与以前的无界情况相比，现在系数以 p 为界，这是否会改变任何内容？

特别是，原始NTT需要找到素数N作为工作模，该模大于(输入向量最大元素的大小)^2 *(输入向量的长度) + 1 这样结果就不会溢出。如果结果无论如何都会受到 p 素数的限制，那么模数可以小到什么程度？请注意，p - 1 不必采用(某个正整数)*(输入向量的长度) 的形式。

编辑:我从上面的链接复制粘贴了源代码来说明问题:

# 
# Number-theoretic transform library (Python 2, 3)
# 
# Copyright (c) 2017 Project Nayuki
# All rights reserved. Contact Nayuki for licensing.
# https://www.nayuki.io/page/number-theoretic-transform-integer-dft
#

import itertools, numbers

def find_params_and_transform(invec, minmod):
    check_int(minmod)
    mod = find_modulus(len(invec), minmod)
    root = find_primitive_root(len(invec), mod - 1, mod)
    return (transform(invec, root, mod), root, mod)

def check_int(n):
    if not isinstance(n, numbers.Integral):
        raise TypeError()

def find_modulus(veclen, minimum):
    check_int(veclen)
    check_int(minimum)
    if veclen < 1 or minimum < 1:
        raise ValueError()
    start = (minimum - 1 + veclen - 1) // veclen
    for i in itertools.count(max(start, 1)):
        n = i * veclen + 1
        assert n >= minimum
        if is_prime(n):
            return n

def is_prime(n):
    check_int(n)
    if n <= 1:
        raise ValueError()
    return all((n % i != 0) for i in range(2, sqrt(n) + 1))

def sqrt(n):
    check_int(n)
    if n < 0:
        raise ValueError()
    i = 1
    while i * i <= n:
        i *= 2
    result = 0
    while i > 0:
        if (result + i)**2 <= n:
            result += i
        i //= 2
    return result

def find_primitive_root(degree, totient, mod):
    check_int(degree)
    check_int(totient)
    check_int(mod)
    if not (1 <= degree <= totient < mod):
        raise ValueError()
    if totient % degree != 0:
        raise ValueError()
    gen = find_generator(totient, mod)
    root = pow(gen, totient // degree, mod)
    assert 0 <= root < mod
    return root

def find_generator(totient, mod):
    check_int(totient)
    check_int(mod)
    if not (1 <= totient < mod):
        raise ValueError()
    for i in range(1, mod):
        if is_generator(i, totient, mod):
            return i
    raise ValueError("No generator exists")

def is_generator(val, totient, mod):
    check_int(val)
    check_int(totient)
    check_int(mod)
    if not (0 <= val < mod):
        raise ValueError()
    if not (1 <= totient < mod):
        raise ValueError()
    pf = unique_prime_factors(totient)
    return pow(val, totient, mod) == 1 and all((pow(val, totient // p, mod) != 1) for p in pf)

def unique_prime_factors(n):
    check_int(n)
    if n < 1:
        raise ValueError()
    result = []
    i = 2
    end = sqrt(n)
    while i <= end:
        if n % i == 0:
            n //= i
            result.append(i)
            while n % i == 0:
                n //= i
            end = sqrt(n)
        i += 1
    if n > 1:
        result.append(n)
    return result

def transform(invec, root, mod):
    check_int(root)
    check_int(mod)
    if len(invec) >= mod:
        raise ValueError()
    if not all((0 <= val < mod) for val in invec):
        raise ValueError()
    if not (1 <= root < mod):
        raise ValueError()

    outvec = []
    for i in range(len(invec)):
        temp = 0
        for (j, val) in enumerate(invec):
            temp += val * pow(root, i * j, mod)
            temp %= mod
        outvec.append(temp)
    return outvec

def inverse_transform(invec, root, mod):
    outvec = transform(invec, reciprocal(root, mod), mod)
    scaler = reciprocal(len(invec), mod)
    return [(val * scaler % mod) for val in outvec]

def reciprocal(n, mod):
    check_int(n)
    check_int(mod)
    if not (0 <= n < mod):
        raise ValueError()
    x, y = mod, n
    a, b = 0, 1
    while y != 0:
        a, b = b, a - x // y * b
        x, y = y, x % y
    if x == 1:
        return a % mod
    else:
        raise ValueError("Reciprocal does not exist")

def circular_convolve(vec0, vec1):
    if not (0 < len(vec0) == len(vec1)):
        raise ValueError()
    if any((val < 0) for val in itertools.chain(vec0, vec1)):
        raise ValueError()
    maxval = max(val for val in itertools.chain(vec0, vec1))
    minmod = maxval**2 * len(vec0) + 1
    temp0, root, mod = find_params_and_transform(vec0, minmod)
    temp1 = transform(vec1, root, mod)
    temp2 = [(x * y % mod) for (x, y) in zip(temp0, temp1)]
    return inverse_transform(temp2, root, mod)

vec0 = [24, 12, 28, 8, 0, 0, 0, 0]
vec1 = [4, 26, 29, 23, 0, 0, 0, 0]

print(circular_convolve(vec0, vec1))

def modulo(vec, prime):
    return [x % prime for x in vec]

print(modulo(circular_convolve(vec0, vec1), 31))

打印:

[96, 672, 1120, 1660, 1296, 876, 184, 0]
[3, 21, 4, 17, 25, 8, 29, 0]

但是，当我将 minmod = maxval**2 * len(vec0) + 1 更改为 minmod = maxval + 1 时，它会停止工作:

[14, 16, 13, 20, 25, 15, 20, 0]
[14, 16, 13, 20, 25, 15, 20, 0]

为了按预期工作，最小的 minmod(上面链接中的 N)是多少？

Answer 1

如果您输入n整数与某个素数 q 绑定(bind)(任何 mod q 而不仅仅是质数都是相同的)您可以将其用作 max value +1但要注意，你不能将它用作素数 p对于 NTT 因为 NTT 质数 p具有特殊的性质。所有这些都在这里:

• Translation from Complex-FFT to Finite-Field-FFT

所以我们每个输入的最大值是 q-1但在任务计算期间(对 2 个 NTT 结果进行卷积)，第一层结果的大小可能会上升到 n.(q-1)但当我们对它们进行卷积时，最终 iNTT 的输入幅度将上升至:

m = n.((q-1)^2)

如果您在 NTT 上执行与 m 不同的操作等式可能会改变。

现在让我们回到 p简而言之，您可以使用任何质数 p坚持这些:

p mod n == 1
p > m

并且存在 1 <= r,L < p这样:

p mod (L-1) = 0
r^(L*i) mod p == 1 // i = { 0,n }
r^(L*i) mod p != 1 // i = { 1,2,3, ... n-1 }

如果这一切都满足的话 p是n次单位根，可用于NTT。找到这样的质数以及 r,L查看上面的链接(有找到这样的 C++ 代码)。

例如，在字符串乘法期间，我们采用 2 个字符串对它们进行 NTT 处理，然后对结果进行卷积，并 iNTT 返回结果(即两个输入大小的总和)。例如:

                                99999999999999999999999999999999
                               *99999999999999999999999999999999
----------------------------------------------------------------
9999999999999999999999999999999800000000000000000000000000000001

q = 10两个操作数都是 9^32 所以 n=32因此m = 9*9*32 = 2592找到的质数是 p = 2689 。如您所见，结果匹配，因此不会发生溢出。但是，如果我使用仍然适合所有其他条件的任何较小的素数，结果将不匹配。我专门用它来尽可能地拉伸(stretch) NTT 值(所有值都是 q-1 并且大小等于 2 的相同幂)

如果您的NTT速度很快并且n不是 2 的幂，那么您需要为每个 NTT 将零填充到最接近的更高或等于 2 的幂。但这不应该影响m作为零填充的值不应增加值的大小。我的测试证明了这一点，因此对于卷积，您可以使用:

m = (n1+n2).((q-1)^2)/2

哪里n1,n2是零填充之前的原始输入大小。

有关实现 NTT 的更多信息，您可以查看我的 C++ 版本(经过广泛优化):

• Modular arithmetics and NTT (finite field DFT) optimizations

回答您的问题:

1. 是的，您可以利用输入 mod q 的事实。但你不能使用q如p !!!

2. 您可以使用minmod = n * (maxval + 1)仅适用于单个 NTT(或 NTT 的第一层)，但当您在 ​​NTT 使用期间将它们与卷积链接时，您不能将其用于最终的 INTT 阶段!!!

但是，正如我在评论中提到的，最简单的是使用最大可能 p适合您正在使用的数据类型，并且可用于支持的所有 2 次方大小的输入。

这基本上使你的问题变得无关紧要。我能想到的唯一不可能/不需要的情况是任意精度数字，其中没有“没有”最大限制。变量 p 绑定(bind)有许多性能问题作为搜索 p确实很慢(甚至可能比 NTT 本身还慢)并且可变 p禁用了所需的模块化算术的许多性能优化，使得 NTT 非常慢。