Курсова робота - арифметика багаторозрядних чисел
Вступ
Для безлічі додатків надаються процесором базових типів цілком вистачає. Однак, зустрічається багато задач, вихідні дані яких занадто великі. Число з 1000 цифр не поміститься ні в один регістр. Тому комп'ютерне подання таких чисел і операції над ними доводиться реалізовувати самостійно.
При цьому час виконання зовнішнього алгоритму, що використовує такі числа, дуже сильно залежить від ефективності їх реалізації. Наприклад, якщо оцінка часу визначається O (n2) множеннями, то використання для цієї операції в два рази більше швидкого алгоритму дає
прискорення в 4 рази. Тому ми будемо, мабуть, найбільш серйозно зацікавлені не просто в правильних, але можливо більш ефективних алгоритмах, які при необхідності можна реалізувати на пристойному рівні.
БИСТРОЕ_УМНОЖЕНІЕ
Ми поставили собі за мету написати не тільки, як відбувається множення над довгими числами, а й встановити, як його зробити більш швидким.
Складемо алгоритм множення:
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; // таку що 2MaxLen. Len
do MaxLen * = 2; while (MaxLen
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
>
Зробимо більш детальний розбір останнього кроку.
Всі обчислення відбуваються в форматі чисел з плаваючою точкою, використовують ірраціональні числа, тому результат буде не набором цілих чисел, а, скоріше, наближенням до нього. Для отримання відповіді кожну координату згортки необхідно округлити до найближчого цілого числа.
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
// Додаємо 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
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
// Крок 2. Переписуємо число під тимчасовий масив і доповнюємо провідними нулями
// Порядок цифр зворотний, тому провідні нулі будуть в кінці
for (x = 0; x