Перевод с комплексного БПФ на конечное поле-БПФ
Добрый день!
Я пытаюсь разработать алгоритм NTT, основанный на наивной рекурсивной реализации FFT, которую я уже имею.
Рассмотрим следующий код (coefficientsдлина, пусть будет m, это точная степень двух):
/// <summary>
/// Calculates the result of the recursive Number Theoretic Transform.
/// </summary>
/// <param name="coefficients"></param>
/// <returns></returns>
private static BigInteger[] Recursive_NTT_Skeleton(
IList<BigInteger> coefficients,
IList<BigInteger> rootsOfUnity,
int step,
int offset)
{
// Calculate the length of vectors at the current step of recursion.
// -
int n = coefficients.Count / step - offset / step;
if (n == 1)
{
return new BigInteger[] { coefficients[offset] };
}
BigInteger[] results = new BigInteger[n];
IList<BigInteger> resultEvens =
Recursive_NTT_Skeleton(coefficients, rootsOfUnity, step * 2, offset);
IList<BigInteger> resultOdds =
Recursive_NTT_Skeleton(coefficients, rootsOfUnity, step * 2, offset + step);
for (int k = 0; k < n / 2; k++)
{
BigInteger bfly = (rootsOfUnity[k * step] * resultOdds[k]) % NTT_MODULUS;
results[k] = (resultEvens[k] + bfly) % NTT_MODULUS;
results[k + n / 2] = (resultEvens[k] - bfly) % NTT_MODULUS;
}
return results;
}
Работало для сложных БПФ (заменить BigInteger со сложным числовым типом (у меня был свой)). Это не работает здесь, хотя я изменил процедуру нахождения примитивных корней единства надлежащим образом.
Предположительно, проблема заключается в следующем: rootsOfUnity Переданный параметр изначально содержал только первую половину mСложные корни единства в следующем порядке:
omega^0 = 1, omega^1, omega^2, ..., omega^(n/2)
Этого было достаточно, потому что на этих трех строках кода:
BigInteger bfly = (rootsOfUnity[k * step] * resultOdds[k]) % NTT_MODULUS;
results[k] = (resultEvens[k] + bfly) % NTT_MODULUS;
results[k + n / 2] = (resultEvens[k] - bfly) % NTT_MODULUS;
Первоначально я использовал тот факт, что на любом уровне рекурсии (для любого n а также i), сложный корень единства -omega^(i) = omega^(i + n/2),
Однако это свойство явно не выполняется в конечных полях. Но есть ли аналог, который позволил бы мне вычислить только первую половину корней?
Или я должен продлить цикл с n/2 в n и предварительно вычислить все mкорни единства?
Может быть, есть другие проблемы с этим кодом?..
Заранее большое спасибо!
3 ответа
Недавно я тоже хотел реализовать NTT для быстрого умножения вместо DFFT. Прочитайте много запутанных вещей, разные буквы везде и нет простого решения, а также мои знания в области конечных полей устарели, но сегодня я наконец-то понял (после 2 дней попыток и аналогии с коэффициентами DFT), так что вот мои идеи для NTT:
вычисление
X(i) = sum(j=0..n-1) of ( Wn^(i*j)*x(i) );где
X[]NTT трансформируетсяx[]размераnгдеWnэто основа NTT. Все вычисления выполняются по целочисленной арифметике по модулюmod pникаких комплексных чисел нигде нет.Важные ценности
Wn = r ^ L mod pявляется основой для NTTWn = r ^ (p-1-L) mod pявляется основой для INTTRn = n ^ (p-2) mod pмасштабирует мультипликативную константу для INTT~(1/n)pпремьер, чтоp mod n == 1а такжеp>max'maxмаксимальное значение x[i] для NTT или X[i] для INTTr = <1,p)L = <1,p)а также делитp-1r,Lдолжны быть объединены такr^(L*i) mod p == 1еслиi=0или жеi=nr,Lдолжны быть объединены такr^(L*i) mod p != 1если0 < i < nmax'является максимальным значением подрезультата и зависит отnи тип вычисления. Для одного (я)NTT этоmax' = n*maxно для свертки двухnразмерные векторы этоmax' = n*max*maxи т.д. См. Реализация БПФ над конечными полями для получения дополнительной информации об этом.функциональная комбинация
r,L,pотличается для разныхnэто важно, вы должны пересчитать или выбрать параметры из таблицы перед каждым уровнем NTT (
nвсегда половина предыдущей рекурсии).
Вот мой код C++, который находит r,L,p параметры (нужна модульная арифметика, которая не включена, вы можете заменить ее на (a+b)%c,(ab)%c,(a*b)%c,... но в этом случае остерегайтесь переполнений, особенно для modpow а также modmul) Код не оптимизирован, но есть способы значительно ускорить его. Кроме того, таблица простых чисел довольно ограничена, поэтому либо используйте SoE, либо любой другой алгоритм для получения простых чисел до max' для того, чтобы работать безопасно.
DWORD _arithmetics_primes[]=
{
2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,
179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,
419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,
661,673,677,683,691,701,709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,
947,953,967,971,977,983,991,997,1009,1013,1019,1021,1031,1033,1039,1049,1051,1061,1063,1069,1087,1091,1093,1097,1103,1109,1117,1123,1129,1151,
0}; // end of table is 0, the more primes are there the bigger numbers and n can be used
// compute NTT consts W=r^L%p for n
int i,j,k,n=16;
long w,W,iW,p,r,L,l,e;
long max=81*n; // edit1: max num for NTT for my multiplication purposses
for (e=1,j=0;e;j++) // find prime p that p%n=1 AND p>max ... 9*9=81
{
p=_arithmetics_primes[j];
if (!p) break;
if ((p>max)&&(p%n==1))
for (r=2;r<p;r++) // check all r
{
for (l=1;l<p;l++)// all l that divide p-1
{
L=(p-1);
if (L%l!=0) continue;
L/=l;
W=modpow(r,L,p);
e=0;
for (w=1,i=0;i<=n;i++,w=modmul(w,W,p))
{
if ((i==0) &&(w!=1)) { e=1; break; }
if ((i==n) &&(w!=1)) { e=1; break; }
if ((i>0)&&(i<n)&&(w==1)) { e=1; break; }
}
if (!e) break;
}
if (!e) break;
}
}
if (e) { error; } // error no combination r,l,p for n found
W=modpow(r, L,p); // Wn for NTT
iW=modpow(r,p-1-L,p); // Wn for INTT
и вот мои медленные реализации NTT и INTT (я еще не дошел до быстрого NTT,INTT), они оба успешно протестированы с умножением Шёнхаге-Штрассена.
//---------------------------------------------------------------------------
void NTT(long *dst,long *src,long n,long m,long w)
{
long i,j,wj,wi,a,n2=n>>1;
for (wj=1,j=0;j<n;j++)
{
a=0;
for (wi=1,i=0;i<n;i++)
{
a=modadd(a,modmul(wi,src[i],m),m);
wi=modmul(wi,wj,m);
}
dst[j]=a;
wj=modmul(wj,w,m);
}
}
//---------------------------------------------------------------------------
void INTT(long *dst,long *src,long n,long m,long w)
{
long i,j,wi=1,wj=1,rN,a,n2=n>>1;
rN=modpow(n,m-2,m);
for (wj=1,j=0;j<n;j++)
{
a=0;
for (wi=1,i=0;i<n;i++)
{
a=modadd(a,modmul(wi,src[i],m),m);
wi=modmul(wi,wj,m);
}
dst[j]=modmul(a,rN,m);
wj=modmul(wj,w,m);
}
}
//---------------------------------------------------------------------------
dst это целевой массив src это исходный массив n размер массива m это модуль (p) w является основой (Wn)
надеюсь, это кому-то поможет. Если я что-то забыл, пожалуйста, напишите...
[edit1: быстрый NTT/INTT]
Наконец-то мне удается заставить быстрый NTT / INTT работать. Был немного сложнее, чем обычный FFT:
//---------------------------------------------------------------------------
void _NFTT(long *dst,long *src,long n,long m,long w)
{
if (n<=1) { if (n==1) dst[0]=src[0]; return; }
long i,j,a0,a1,n2=n>>1,w2=modmul(w,w,m);
// reorder even,odd
for (i=0,j=0;i<n2;i++,j+=2) dst[i]=src[j];
for ( j=1;i<n ;i++,j+=2) dst[i]=src[j];
// recursion
_NFTT(src ,dst ,n2,m,w2); // even
_NFTT(src+n2,dst+n2,n2,m,w2); // odd
// restore results
for (w2=1,i=0,j=n2;i<n2;i++,j++,w2=modmul(w2,w,m))
{
a0=src[i];
a1=modmul(src[j],w2,m);
dst[i]=modadd(a0,a1,m);
dst[j]=modsub(a0,a1,m);
}
}
//---------------------------------------------------------------------------
void _INFTT(long *dst,long *src,long n,long m,long w)
{
long i,rN;
rN=modpow(n,m-2,m);
_NFTT(dst,src,n,m,w);
for (i=0;i<n;i++) dst[i]=modmul(dst[i],rN,m);
}
//---------------------------------------------------------------------------
[Edit3]
Я оптимизировал свой код (в 3 раза быстрее, чем код выше), но все же я не удовлетворен этим, поэтому я начал новый вопрос с ним. Там я оптимизировал свой код еще дальше (примерно в 40 раз быстрее, чем код выше), поэтому его скорость почти такая же, как у FFT с плавающей запятой того же размера. Ссылка на это здесь:
Чтобы преобразовать БПФ Кули-Тьюки (комплексное) в модульный арифметический подход, т. Е. NTT, необходимо заменить комплексное определение на омега. Чтобы подход был чисто рекурсивным, вам также необходимо пересчитать омега для каждого уровня на основе текущего размера сигнала. Это возможно, потому что мин. подходящий модуль уменьшается при перемещении вниз по дереву вызовов, поэтому модуль, используемый для корня, подходит для нижних уровней. Кроме того, поскольку мы используем один и тот же модуль, при перемещении по дереву вызовов может использоваться один и тот же генератор. Кроме того, для обратного преобразования вы должны предпринять дополнительный шаг, чтобы получить пересчитанный омега a и вместо этого использовать как омега: b = a ^ -1 (с помощью обратной операции по модулю). В частности, b = invMod(a, N) st b * a == 1 (mod N), где N - выбранный простой модуль.
Переписывание выражения с использованием омеги с использованием периодичности все еще работает в модульной арифметической области. Вам также нужно найти способ определить модуль (простое число) для задачи и действительный генератор.
Мы отмечаем, что ваш код работает, хотя это не MWE. Мы расширили его, используя здравый смысл, и получили правильный результат для приложения полиномиального умножения. Вы просто должны предоставить правильные значения омега, возведенные в определенные силы.
Хотя ваш код работает, как и во многих других источниках, вы удваиваете интервалы для каждого уровня. Это не приводит к рекурсии, которая столь же чиста, хотя; это оказывается идентичным пересчету омега на основе текущего размера сигнала, потому что мощность для определения омега обратно пропорциональна размеру сигнала. Повторим еще раз: вдвое меньший размер сигнала подобен квадратному омеге, что похоже на удвоение мощности омеги (что можно сделать для удвоения расстояния). Хорошая вещь о подходе, который имеет дело с пересчетом омеги, состоит в том, что каждая подзадача является более четко завершенной сама по себе.
Есть статья, в которой показана некоторая математика для модульного подхода; это статья Бактира и Сунара от 2006 года. См. статью в конце этого поста.
Вам не нужно расширять цикл с n / 2 до n.
Так что, да, некоторые источники, которые говорят, что просто опускают другое определение омеги для модульного арифметического подхода, охватывают множество деталей.
Другая проблема заключается в том, что важно признать, что размер сигнала должен быть достаточно большим, чтобы не допустить переполнения результирующего сигнала во временной области, если мы выполняем свертку. Кроме того, может быть полезно найти определенные реализации для возведения в степень, зависящие от модуля, которые бывают быстрыми, даже если мощность довольно велика.
Рекомендации
- Бактир и Сунар - Достижение эффективного полиномиального умножения в полях Ферма с использованием быстрого преобразования Фурье (2006)
Вы должны убедиться, что корни единства действительно существуют. В R есть только 2 корня единицы: 1 и -1, поскольку только для них x^n=1 может быть истинным.
В C у вас бесконечно много корней из единицы: w=exp(2*pi*i/N) является примитивным N-м корнем из единицы и все w^k при 0<=k
Теперь к вашей проблеме: вы должны убедиться, что кольцо, в котором вы работаете, предлагает то же свойство: достаточно корней единства.
Шенхаге и Штрассен ( http://en.wikipedia.org/wiki/Sch%C3%B6nhage%E2%80%93Strassen_algorithm) используют целые числа по модулю 2^N+1. У этого кольца достаточно корней единства. 2^N == -1 - 2-й корень единицы, 2^(N/2) - 4-й корень единицы и так далее. Кроме того, эти корни единства имеют то преимущество, что они представляют собой степени двойки и могут быть реализованы в виде двоичных сдвигов (с последующей операцией по модулю, которая сводится к сложению / вычитанию).
Я думаю, что QuickMul ( http://www.cs.nyu.edu/exact/doc/qmul.ps) работает по модулю 2^N-1.