Реализация БПФ над конечными полями

Question

Я хотел бы реализовать умножение полиномов с использованием NTT. Я следовал за Теоретико-числовым преобразованием (целое число DFT), и это, кажется, работает.

Теперь я хотел бы реализовать умножение многочленов над конечными полями Z_p[x] где p произвольное простое число.

Изменит ли это что-нибудь, что коэффициенты теперь ограничены pпо сравнению с прежним неограниченным случаем?

В частности, в оригинальном NTT требуется найти простое число N в качестве рабочего модуля, который больше, чем (magnitude of largest element of input vector)^2 * (length of input vector) + 1 так что результат никогда не переполняется. Если результат будет ограничен этим p в любом случае, насколько мал модуль может быть? Обратите внимание, что p - 1 не должно быть формы (some positive integer) * (length of input vector),

Редактировать: я скопировал исходный код по ссылке выше, чтобы проиллюстрировать проблему:

# 
# 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 (любой mod q не просто премьер будет такой же) Вы можете использовать его как max value +1 но будьте осторожны, вы не можете использовать его как премьер p для NTT, потому что премьер NTT p имеет особые свойства. Все они здесь:

• Перевод с комплексного БПФ на конечное поле-БПФ

поэтому наше максимальное значение каждого входа 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, Как вы видите результат совпадения, переполнение не происходит. Однако, если я использую любое меньшее простое число, которое все еще соответствует всем другим условиям, результат не будет совпадать. Я использовал это специально, чтобы максимально растянуть значения NTT (все значения q-1 и размеры равны одинаковой степени 2)

Если ваш NTT быстрый и n не является степенью 2, тогда вам нужно обнулить площадку до ближайшей более высокой или равной степени 2 для каждого NTT. Но это не должно влиять на m значение в качестве нуля не должно увеличивать величину значений. Мое тестирование доказывает это, поэтому для свертки вы можете использовать:

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

где n1,n2 являются исходными размерами сырья перед нулевой

Для получения дополнительной информации о реализации NTT вы можете проверить мой в C++ (значительно оптимизирован):

• Модульная арифметика и NTT (конечное поле DFT) оптимизации

Итак, чтобы ответить на ваши вопросы:

1. да, вы можете воспользоваться тем, что вход mod q но вы не можете использовать q как p !!!

2. Ты можешь использовать minmod = n * (maxval + 1) только для одного NTT (или первого уровня NTT), но, поскольку вы объединяете их в цепочку во время использования NTT, вы не можете использовать это для заключительного этапа INTT!!!

Однако, как я уже упоминал в комментариях, проще всего использовать максимально возможное p он соответствует типу данных, который вы используете, и может использоваться для всех типов поддерживаемых входных данных.

Что в основном делает ваш вопрос неактуальным. Единственный случай, который я могу придумать, где это невозможно / желательно, - это числа с произвольной точностью, где нет максимального ограничения. Существует много проблем с производительностью, связанных с переменной p как поиск p очень медленно (может быть даже медленнее, чем само NTT), а также переменная p отключает многие оптимизации производительности модульной арифметики, необходимые для замедления NTT.