Арифметика многоразрядных чисел

     if ( Len < 8 ) Len = 8; // FFT работает  

     // Шаг 2. Переписываем число во временный массив и дополняем ведущими нулями

     // Порядок цифр обратный, поэтому  ведущие нули будут в конце 

     for (x = 0; x < A.Size; x++) LongNum1[x] = a[x];

     for (; x < Len; x++) LongNum1[x] = 0.;

     // Шаги 3,4.

     RealFFT(LongNum1, Len, 1);

     if (&A == &B) {

     // Если имеем дело с возведением  в квадрат – можно сэкономить  одно БПФ 

     RealFFTScalar ( LongNum1, LongNum1, Len );

     } else {

     // Иначе обработать и второе  число 

     for (x = 0; x < B.Size; x++) LongNum2[x] = B.Coef[x];

     for (; x < Len; x++) LongNum2[x] = 0.;

     RealFFT(LongNum2, Len, 1);

     RealFFTScalar ( LongNum1, LongNum2, Len );

     }

     RealFFT(LongNum1, Len, -1); // Шаг 5

     CarryNormalize( LongNum1, Len, C); // Шаг 6

     }

     По  окончании умножения в MaxError хранится наибольшая ошибка округления. Для правильности результата ей лучше быть меньше 2.

     Применение БПХ для вычисления свертки

     Умножение оперирует с действительными  числами. Как мы говорили, в этом случае лучшие результаты может дать преобразование Хартли. Сделаем более совершенный вариант быстрого умножения на его основе.

     Индексы считаются по mod N, то есть, вместо элементов  с индексом N нужно брать элементы с нулевым индексом. Все шаги алгоритма  БЫСТРОЕ_УМНОЖЕНИЕ остаются такими же, кроме шага 4, где координаты вектора следует вычислять по формуле (3).

     Существование двух стилей БПХ и БПФ(разбиения  по частоте и по времени) позволяет  сделать одно существенное улучшение  алгоритма. А именно, полностью избежать вызова FFTReOrder(), который занимает около 10% времени.

     Действительно, БПХ_В требует на вход реорганизованный вектор и дает - обычный, в то время  как БПХ_Ч работает наоборот: на входе  – обычный, а в конце работы – реорганизованный.

     Поэтому можно изменить шаги следующим образом:

     Векторы хранятся с обычным порядком индексов.

     3. Вычислить БПХ_Ч на обоих массивах  цифр.

     Вектор  на выходе в реорганизованном виде.

     4. Вычислить коэффициенты ДПХ(a . b) по формуле (3), учитывая изменения  в порядке индексов.

     Вектор  все еще в реорганизованном виде.

     5. Вычислить БПХ_В, результатом  которого будет ненормализованная  свертка 

     Окончательный вектор имеет обычный порядок  индексов.

     Все необходимые для реализации кирпичики  уже есть. Единственно, нужно научиться  применять 

     формулу (3) к векторам, “перетасованным” функцией FFTReOrdeк().

     В качестве примера рассмотрим вычисление ДПХ(a . b) для a=БПХ_Ч(a1, a2, …) , b=БПХ_Ч(b1, b2, ...) на месте вектора a.

     Элементы  с индексами 0 и N/2 можно обработать отдельно. Для них формула приобретает  особенно простой вид: a0 = a0b0 , aN/2 = aN/2bN/2.

     Для остальных элементов стоит учесть, что на вход формулы поступают  числа с двух позиций, а выход записывается на одну. Чтобы не произошло лишнее затирание, можно вычислять по два значения одновременно, используя симметричность выражения.

     Теперь  элементы вектора a с суммой индексов N будут преобразовываться попарно  и “на месте”.

     Образуется  “бабочка свертки”, выполнение которой  для k=1…N/2-1 дает необходимый результат.

     Если  порядок индексов обычный, то вычисления производятся непосредственно по формулам.

     Для 8-элементных векторов:

     Вычисление  ДПХ свертки через ДПХ векторов. Индексы идут в обычном порядке, N=8

     Интереснее  с реорганизованными векторами. Чтобы учесть обратный битовый порядок, можно выполнять сначала бабочки для 2 элементов, затем для 4, для 8.. И так далее с шагом *2. Первые два элемента имеют индексы 0, N/2 и их следует обработать отдельно.

     Step = 2

     Step = 4

     Step = 8

     Вычисление  ДПХ свертки через ДПХ векторов. Индексы в обратном битовом порядке, N=16

     // Вычисление ДПХ свертки на месте LongNum1

     // Не производится умножение   в формуле (3), поэтому элементы 

     // получаются в 2 раза больше, чем  должны быть.

     a0 a8 a4 a12 a2 a10 a6 a14 a1 a9 a5 a13 a3 a11 a7 a15

     b0 b8 b4 b12 b2 b10 b6 b14 b1 b9 b5 b13 b3 b11 b7 b15

     a0 a8 a4 a12 a2 a10 a6 a14 a1 a9 a5 a13 a3 a11 a7 a15

     b0 b8 b4 b12 b2 b10 b6 b14 b1 b9 b5 b13 b3 b11 b7 b15

     a0 a8 a4 a12 a2 a10 a6 a14 a1 a9 a5 a13 a3 a11 a7 a15

     b0 b8 b4 b12 b2 b10 b6 b14 b1 b9 b5 b13 b3 b11 b7 b15

     a0 a1 a2 a3 a4 a5 a6 a7

     b0 b1 b2 b3 b4 b5 b6 b7

     void FHTConvolution(real *LongNum1, const real *LongNum2, ulong Len) { ulong Step=2, Step2=Step*2;

     ulong x,y;

     // Индексы 0, N/2. Как и все остальное, по формуле (3) без .

     LongNum1[0] = LongNum1[0] * 2.0 * LongNum2[0];

     LongNum1[1] = LongNum1[1] * 2.0 * LongNum2[1];

     while (Step < Len) { // Step – текущий шаг вычисления бабочек

     for (x=Step,y=Step2-1;x<Step2;x+=2,y-=2) {

     // Преобразовать элементы с индексами  x, y по формуле (3) без .

     real h1p,h1m,h2p,h2m;

     real s1,d1;

     h1p=LongNum1[x];

     h1m=LongNum1[y];

     s1=h1p+h1m;

     d1=h1p-h1m;

     h2p=LongNum2[x];

     h2m=LongNum2[y];

     LongNum1[x]=(h2p*s1+h2m*d1);

     LongNum1[y]=(h2m*s1-h2p*d1);

     }

     Step*=2;

     Step2*=2;

     }

     }

     Вычисление  пары коэффициентов ДПХ свертки  таким способом требует 4 сложений и 4 умножений, в то время как для преобразования Фурье эти цифры составляют 2 комлексных умножения , то есть 8 действительных умножений и 4 сложения. Однако, комплексный вектор в 2

     раза  короче, поэтому ДПХ требует на 2 сложения больше для пары элементов  – всего на N сложений больше для всего вектора.

     Теперь  уже можно организовать быстрое  умножение с использованием функций FHT_F(), FHT_T(). Для этого нужно добавить в конец FastMulInit() вызов CreateSineTable(MaxLen) и изменить коэффициент нормализации в CarryNormalize() на inv = 0.5 / Len: длина преобразования Len, а 0.5 появляется из-за того, что элементы свертки в два раза больше, чем надо.

     // Быстрое умножение с совместным  использованием разбиения по  времени и по частоте 

     // Работает только для чисел,  результат умножения которых имеет 8 или более разрядов

     // Ограничение – следствие базиса  рекурсии в 8 элементов для  FFT_T.

     // В любом случае, для таких маленьких  чисел обычное умножение намного  лучше. 

     void FastMul(const BigInt &A,const BigInt &B, BigInt &C) {

     ulong Len = 2;

     while ( Len < A.Size + B.Size ) Len *=2;

     // (**)

     ulong x;

     const short *a=A.Coef, *b=B.Coef;

     for (x = 0; x < A.Size; x++) LongNum1[x] = a[x];

     for (; x < Len; x++) LongNum1[x] = 0.0;

     FHT_F(LongNum1,Len);

     if (a == b) {

     FHTConvolution(LongNum1,LongNum1,Len);

     } else {

     for (x = 0; x < B.Size; x++) LongNum2[x] = b[x];

     for (; x < Len; x++) LongNum2[x] = 0.0;

     FHT_F(LongNum2,Len);

     FHTConvolution(LongNum1,LongNum2,Len);

     }

     FHT_T(LongNum1,Len);

     CarryNormalize(LongNum1, Len, C);

     Сравнение “быстрого” и “школьного” умножения

     Итак, алгоритм разобран и запрограммирован. Теперь самое время сравнить два  метода.

     Для сравнения возьмем БПХ-умножение  как метод, оптимизированный нами гораздо  лучше БПФ, и функцию Mul(), реализованную  при создании класса BigInt. Умножаются два числа одинакового размера.

     Сравнительное время умножения 

     0

     500

     1000

     1500

     2000

     2500

     3000

     3500

     4000

     4500

     0 50 100 150 200 250

     длина числа (десятичные цифры)

     время

     Как видно, поведение методов принципиально  различается. Скачок в БПХ-методе возникает при переходе через степень двойки. Например, при перемножении 128-значных чисел вычисляется ДПХ 128-координатных векторов, а уже для умножения 129-значных векторов приходится работать с наименьшей степенью двух, большей 129 – с 256-значными векторами.

     Исходя  из тестовых данных, начиная со 150-300 десятичных цифр(в зависимости от компилятора эффективность кода разная) БПХ-алгоритм становится быстрее.

     Естественным  следствием из такого поведения алгоритмов является использование обычного умножения для чисел длины менее 40(возможно, 80) разрядов и применение “быстрого” метода при больших размерах.

     // вставить вместо (**) в функцию  FastMul

     if ( Len < 40 ) {

     BigInt Atmp(A), Btmp(B);

     Mul(Atmp,Btmp,C);

     return;

     }

     Временные числа Atmp, Btmp создаются, чтобы сохранить возможность запуска FastMul(A,B,A), так как Mul(A, B, A) всегда обнуляет место для результата(здесь – число A) перед процедурой.

     Конечно, рассмотренные реализации алгоритмов не являются наиболее оптимальными. Однако, судя по многочисленным тестам, хорошо запрограммированное “быстрое” умножение опережает школьный алгоритм того же качества в районе 250 цифр.

     Точность  вычислений и ее улучшения.

     Как уже говорилось, в отличие от школьного  алгоритма, “быстрое” умножение может ошибаться.

     Очень плохими векторами в этом смысле являются (BASE-1, BASE-1, …, BASE-1). При BASE=10000 БПФ-умножения  позволяло умножать числа до 4 миллионов  цифр. Оценка БПХ несколько лучше, здесь верхняя граница достигала 16 миллионов.

     Очень сильно помочь может “балансировка” исходных векторов. Координаты входных векторов лежат в от 0 до BASE-1. Перед умножением преобразуем их к интервалу [-BASE/2…BASE/2-1], а после - вернем цифры в обычный диапазон.

     Не  важно, используется ли БПХ или БПФ, опыт показывает, что при таком подходе ошибка на случайном векторе уменьшается во много раз. Если умножение работает на пределе точности, то рекомендуется использовать этот способ.

     Математически доказано, что максимальная длина  чисел в этом случае увеличивается в 4 раза. То есть, если из формулы (2) следует, что для N=220 и типа double необходимое основание BASE < 4100, то при балансировке будет BASE < 16400. Практически, с основанием 10000 можно перемножать числа даже до 230, однако это выходит за рамки математических гарантий

правильности  результата.

      Заключение

     Альтернативные  методы умножения.

     Существуют  рекурсивные алгоритмы, основанные на подходе “разделяй-и-властвуй”. В качестве примера, можно упомянуть  методы Карацубы, Тома-Кука. Они позволяют  свести умножение большого числа к умножениям двух-трех в 2 раза меньших и нескольким сложениям, уменьшая