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

          Begin

            I := NMax;

            While (I > 1) And (C[I] = 0) Do I := I - 1;

            Dlina := I

          End;

     При ее разработке было использовано следующее соображение: если число не равно нулю, то количество цифр в его записи равно номеру первой цифры, отличной от нуля, если просмотр числа осуществляется от старшего разряда к младшему. Если же длинное число равно нулю, то получается, что количество цифр в его записи равно одной, что и требовалось.

     Ну  и, наконец, главная процедура, ради которой и была проделана вся  предшествующая работа. При составлении  алгоритма используется идея умножения "столбиком", хотя в нашем варианте сложение выполняется не по окончанию умножения, а по ходу его, т.е. перемножив очередные цифры, сразу же добавляем результирующую цифру в нужный разряд и формируем перенос в следующий разряд.

     (Процедура  умножения длинных чисел. A, B —  множители, C — произведение)

      Procedure Multiplication(A, B : DlChislo; Var C : DlChislo);

       Var I, J : Integer; P : Digit; VspRez : 0..99;

       Begin

         Zero(C);

         For I := 1 To Dlina(A) Do {Цикл по количеству цифр в первом числе}

          Begin

             P:= 0; {Первоначально перенос равен  нулю}

             For J:= 1 To Dlina(B) Do {Цикл по количеству цифр во втором числе}

              Begin

                VspRez := A[I] * B[J] + P + C[I + J - 1];

                C[I + J - 1]:= VspRez Mod 10; {Очередное значение цифры в разряде I + J - 1}

                P:= VspRez Div 10 {Перенос в следующий разряд}

              End;

            C[I + J]:= P {последний перенос может  быть отличен от нуля, запишем  его в пока ещё свободный  разряд}

          End

       End;

     Целиком программа представлена в Приложении 3. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

     БЫСТРОЕ_УМНОЖЕНИЕ

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

     Составим  алгоритм умножения:

     1. Найти наименьшее число Len –  степень двойки: Len . A.Size + B.Size. Для  рассматриваемых чисел Len=8.

     2. Дополнить A и B ведущими незначащими нулями до Len. А нашем примере A=(3,4,0,0,0,0,0,0), B=(5,4,3,2,1,0,0,0).

     3. Вычислить БПФ действительных  векторов на обоих массивах  цифр.

     4. Скалярно перемножить преобразованные  вектора, получив вектор размера  Len.

     5. Применить обратное преобразование Фурье, результатом которого будет свертка.

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

       Цифры больших чисел хранятся  в целочисленном формате. Поэтому  для БПФ их необходимо скопировать во временные массивы типа с плавающей точкой. Если предполагается получать результат максимальной длины N, то необходимо выделить для них память как минимум размера

     MaxLen=2k, где MaxLen – минимальная степень  двойки, большая N. Например, если  максимальный результат будет  состоять из 1000 цифр по основанию BASE, то минимальный объем памяти MaxLen=1024, так как именно такой длины БПФ будет вычисляться.

     real *LongNum1=NULL, *LongNum2=NULL;

     // Инициализацию можно делать только 1 раз за всю программу. 

     void FastMulInit(ulong Len) {

     ulong MaxLen;

     if ((Len & -Len) == Len) // Len = степень двойки

     MaxLen = Len;

     else { // иначе вычислить MaxLen – наименьшую степень 2,

     MaxLen = 2; // такую что 2MaxLen . Len

     do MaxLen*=2; while (MaxLen < Len);

     }

     LongNum1 = new real[MaxLen];

     LongNum2 = new real[MaxLen];

     }

     // Деинициализация

     void FastMulDeInit() {

     delete LongNum1;

     delete LongNum2;

     }

     Разобранная в соответствующем разделе функция RealFFT() производит преобразование “на  месте”, возвращая результирующие векторы в сжатом виде.

     Соответственно, функция скалярного произведения должна учитывать такой формат.

     // Скалярное умножение комплексных векторов в сжатом виде: LongNum1 = LongNum1*LongNum2

     void RealFFTScalar(real *LongNum1, const real *LongNum2, ulong Len) {

     Complex *CF1= (Complex*)LongNum1;

     const Complex *CF2=(Complex*)LongNum2;

     // первые два элемента - сгруппированные в одно комплексное действительные числа

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

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

     for (ulong x = 1; x < Len/2; x++) // остальные – комплексные, как им и

     CF1[x] = CF1[x]*CF2[x]; // следует быть после  ДПФ 

     }

     Сделаем более подробный разбор последнего шага.

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

     a[0] a[N/2] a[1] a[2] ….. a[N/2-1]

     b[0] b[N/2] b[1] b[2] ….. b[N/2-1]

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

     Часть вопросов, связанных с точностью, была решена в обсуждении БПФ. Однако даже при абсолютно точной тригонометрии будут накапливаться ошибки вычислений, так как операции арифметические операции не могут производиться с абсолютной точностью. Поэтому размер используемого типа данных должн быть достаточно большим, чтобы ошибка на любом знаке была меньше 0.5.

     Например, при использовании типа данных размера 1, дробь 1/3 представляется в виде 0.3. При сложении трех дробей получаем

     1/3 + 1/3 + 1/3 = (в формате чисел с плавающей  точкой) 0.3 + 0.3 + 0.3 = 0.9

     Если  же размер – две цифры, то 1/3 = 0.33,

     1/3 + 1/3 + 1/3 = 0.33 + 0.33 + 0.33 = 0.99. Точность вычислений  сильно возросла.

     Вообще  говоря, путей увеличения точности два. Один из них связан с увеличением длины используемого типа: От float к double, далее к long double, потом к double double2…

     Другой  подход более гибок. Он предполагает при фиксированном типе уменьшать  длину основания BASE. Таким образом  число станет длиннее, будет занимать больше памяти, но за счет этого увеличивается точность.

     Ограничения БПФ-умножения 

     Интересную  оценку для ошибок дал Колин Персиваль .

     Пусть требуется перемножить векторы из 2n координат с использованием БПФ для векторов с действительными координатами. Тогда из его результатов следует, что погрешность err умножения x на y оценивается сверху выражением

     err < 2n BASE2 ( .*3n + . 5 (3n+4) + .(3n+3) ) (2)

     где . - точность сложения/умножения, . - точность тригонометрических вычислений,

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

     Например, при используемом типе double(53 бита), .=2-53. Ошибки тригонометрии ограничены величиной .=./ 2 .

     Оценим  верхнюю границу ошибок (2) числом .. Приблизительно посчитав константы, получаем 2n BASE2 2-53 (11.83 n + 11.07) < .

     Для чисел длины 220 это ведет к неравенству BASE < 4100. Такова оценка худшего случая, обоснованная математически.

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

     При округлении следует смотреть на разницу  между округляемым значением  и результатом округления. Если она  меньше 0.2, то умножение, скорее всего, дает правильный результат, если больше – рекомендуется уменьшить BASE или воспользоваться другим базовым типом для хранения коэффициентов.

     После выполнения шага 5 нет готового произведения, а есть лишь свертка – результат  без переносов. Как уже говорилось при рассмотрении пирамиды умножения, значения коэффициентов свертки могут быть много больше основания, достигая 2N*BASE2. Если дополнительно вспомнить, что при обратном преобразовании Фурье происходит деление результатов выполнения функции RealFFT() на N , то максимальный размер цифры становится равен 2N2*BASE2, поэтому следует позаботиться, чтобы не произошло переполнения. В частности, не следует объявлять BASE длиннее 4х десятичных цифр.

     Последние два типа поддерживают далеко не все процессоры  

     В качестве резюме к вышесказанному, заметим, что проблемы всего три:

     1. Точность тригонометрии 

     2. Точность при вычислении БПФ 

     3. Переполнение базового типа.

     Первая  проблема решена при обсуждении БПФ. Вторая и третья решаются путем уменьшения BASE, либо увеличения базового типа. При  этом эффективность алгоритма падает, так как меньшее основание означает удлинение количества цифр, а больший базовый тип не всегда доступен.

     Следующая функция преобразует свертку Convolution длины Len в большое число C, делая  округления и выполняя переносы. В  конце выполнения переменная MaxError будет содержать максимальную ошибку округления.

     RealFFT() не производит нормализацию результата, поэтому ее необходимо сделать  здесь же.

     real MaxError;

     void CarryNormalize(real *Convolution, ulong Len, BigInt &C) {

     real inv = 1.0 / (Len/2); // коэффициент нормализации

     // ДПФ проводилось над “комплексным”  вектором 

     // в 2 раза меньшей длины 

     real RawPyramid, Pyramid, PyramidError, Carry = 0;

     short *c = C.Coef;

     ulong x;

     // В C поместим только ту часть  результата, которая туда влезает 

     // ДПФ имеет длину, равную 2k, но вектор коэффициентов

     // мог быть инициализован на  меньшее количество элементов,  не на степень 2.

     if ( Len > C.SizeMax ) Len=C.SizeMax;

     MaxError = 0;

     for (x = 0; x < Len; x++) {

     RawPyramid = Convolution[x] * inv + Carry;

     // Прибавляем 0.5, чтобы округление произошло к ближайшему целому

     Pyramid = floor(RawPyramid + 0.5);

     PyramidError = fabs(RawPyramid - Pyramid);

     if (PyramidError > MaxError)

     MaxError = PyramidError;

     Carry = floor(Pyramid / BASE); // вычисляем переносы

     c[x] = (short)(Pyramid - Carry * BASE);

     }

     // Все готово, осталось лишь установить  размер С, по первому 

     // ненулевому коэффициенту.

     do { x--; } while (c[x]==0);

     C.Size = x+1;

     }

     Теперь  можно реализовать алгоритм целиком.

     // Вычислить C = A*B, работает вызов FastMul(A, B, A)

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

     ulong x;

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

     if (!LongNum1 || !LongNum2) error("FastMul not initalized.");

     // Шаг 1

     ulong Len = 2;

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