Арифметика многоразрядных чисел
Арифметика
многоразрядных чисел
Оглавление
Введение…………………………………………………………
Представление
числа в любой системе
Операции
над длинными числами..………………………………………….....
Сложение многоразрядных чисел……………………………………………….6
Вычитание многоразрядных чисел ……………………………………………...7
Умножение многоразрядных чисел ……………………………………………8
Быстрое
умножение………………………………………………………
Сравнение “быстрого” и “школьного” умножения…...………………………22
Точность
вычислений и её улучшение………………………………………...23
Заключение……………………………………………………
Литература……………………………………………………
Приложение…..……………………………………………
Введение
Для множества приложений предоставляемых процессором базовых типов вполне хватает. Однако, встречается много задач, исходные данные которых слишком велики. Число из 1000 цифр не поместится ни в один регистр. Поэтому компьютерное представление таких чисел и операции над ними приходится реализовывать самостоятельно.
При этом время выполнения внешнего алгоритма, использующего такие числа, очень сильно зависит от эффективности их реализации. Например, если оценка времени определяется O(n2) умножениями, то использование для этой операции в два раза более быстрого алгоритма дает
ускорение
в 4 раза. Поэтому мы будем, пожалуй, наиболее
серьезно заинтересованы не просто в правильных,
но возможно более эффективных алгоритмах,
которые при необходимости можно реализовать
на приличном уровне.
Представление числа в любой системе счисления
Обычно, неотрицательное целое число N длины n представляется в виде N = a0 + a1*BASE + a2*BASE2 + ... + an-1 BASEn-1,
где BASE – основание системы счисления, все коэффициенты 0 . ai < BASE.
Например, число в этой интерпретации будет иметь вид
1234510 = 5 + 4*10 + 3*102 + 2*103 + 1*104 (BASE=10).
Длинное
число хранится в массиве, где i-й
элемент соответствует
В качестве примера, рассмотрим массив для : : (5, 4, 3, 2, 1).
Как видно, цифры хранятся в обратном порядке. Это – некая “заготовка на будущее”: дело в том, что реализации алгоритмов при этом имеют более естественный вид.
Такое представление N является частным случаем многочлена n-й степени P(x) = a0 + a1*x + a2*x2 + ... + an-1 xn-1, который также удобно хранить в виде массива коэффициентов. Поэтому большинство основных операций над числами при соответствующем упрощении алгоритмов работают для произвольных многочленов (сложение, умножение и т.п).
Знак числа, как и место десятичной точки, можно запомнить в отдельной переменной и учитывать при выполнении операций. Далее мы будем оперировать лишь с целыми неотрицательными числами.
Основание системы счисления BASE обычно зависит от максимального размера базового типа данных на компьютере, и выбирается, исходя из следующих соображений:
1.
Основание должно подходить
2. BASE должно быть как можно больше, чтобы уменьшить размер представления длинного числа и увеличить скорость операций с ними, но достаточно малого размера, чтобы все операции с коэффициентами использовали базовый тип данных.
3. Для удобства можно выбрать BASE как степень 10 (вывод информации, отладка).
Операции над «длинными числами»
Рассмотрим достаточно популярную в программировании задачу на работу с "длинными" числами. Реально с "астрономическими" или "микроскопическими" числами приходится сталкиваться не так уж и часто. Тем не менее, упражнения, рассматриваемые в этой публикации, могут послужить хорошей тренировкой в области программирования и занять достойное место в классах с углубленным изучением информатики или на кружках по программированию. Алгоритмы, представленные ниже, записаны на Turbo Pascal, версия 7.0 и на Си++. При желании или необходимости они могут легко быть адаптированы к любой другой программной среде. Диапазон представления целых чисел (Integer, Word, LongInt) ограничен, о чем не раз уже говорилось (впрочем, для действительных величин это замечание тоже актуально). Поэтому при решении задач всегда приходится действовать с оглядкой, — как бы не допустить возникновения ошибки выхода за диапазон или переполнения. Например, вычисляя факториал (n! = 1 * 2 * 3 * … * n), в диапазоне представления величин типа Integer удастся правильно получить только 7! = 5040, а в диапазоне представления типа LongInt — 12! = 479001600. Для больших значений, конечно, можно использовать действительные типы данных, но это уже не гарантирует точного результата. Поэтому полезно для получения точных значений при действиях с многозначными числами разработать другие способы представления таких чисел, алгоритмы выполнения арифметических и других операций, процедуры ввода и вывода результатов и т.д.
Числа, для представления которых в стандартных компьютерных типах данных не хватает количества двоичных разрядов, называются "длинными". Реализация арифметических операций над такими "длинными" числами получила название "длинной арифметики". Наиболее естественным способом представления многозначного числа является запись каждого его разряда в виде отдельного элемента линейного массива (или списка, где память под цифру будет отводиться по мере надобности, в то время как в массиве приходится заранее задавать максимальное количество элементов в нем). Пусть (для удобства дальнейших действий) разряды "длинного" числа при записи в массив нумеруются с единицы, начиная с разряда единиц, т.е. цифра из разряда единиц — элемент массива с номером один, цифра из разряда десятков — элемент массива с номером два и т.д. Определим для работы с "длинными" неотрицательными числами тип данных:
Но
если мы будем реализовывать
Сложение многоразрядных(длинных) чисел
Практически каждый умеет складывать длинные числа, используя для этого листок бумаги и карандаш. Чем мы теперь займемся – это перенесем имеющееся у нас понимание на язык, понимаемый компьютером.
Схема для сложения очень проста: складываем цифры слева направо (цифры хранятся в обратном порядке). Если зафиксировано переполнение (т.е при сложении получена цифра, большая максимально возможной в данной системе счисления), то происходит “перенос” в следующий разряд.
Найдём сумму чисел А и В. Будем складывать элементы массивов с одинаковыми номерами и хранить полученные результаты, как элементы массива С (алгоритм работы с такими массивами напоминает сложение в столбик). Если на каком-то шаге получим С[i] , то на месте C[i] оставляем остаток от деления C[i] на 10 (количество единиц данного разряда), а к элементу C[i+1] прибавим 1 (перенос 1 в следующий разряд).
Если число А состоит из N цифр, а число В – из M цифр, то число С имеет либо max(N,M) цифр, либо max(N,M)+1 цифру. Обозначим через К величину max(N,M)+1.
Перед выполнением сложения следует дополнить незначащими нулями число с меньшим количеством цифр так, чтобы количество цифр уравнялось. Этого можно достигнуть и обнулив массивы А и В до ввода цифр.
Алгоритм сложения двух многоразрядных чисел представлен в Приложении 1.
После выполнения данного алгоритма в первых К элементах массива С будет храниться сумма чисел А и В, С[1]- цифра младшего разряда. Последняя проверка позволяет убрать незначащий нуль из старшего разряда числа С.
Вычитание многоразрядных чисел
Вычитание осуществляется аналогично, с той лишь разницей, что осуществляется перенос “заимствования”. Мы работаем с положительными числами, поэтому если вычитаемое большое по размеру – происходит выход.
Если длины одинаковы, но одно больше другого - это приведет к тому, что в конце процесса останется неиспользованное заимствование, а результат будет дополнением до BASEB.Size.
Для определённости будем считать, что А>B. Если нет, то необходимо поменять числа местами, а результат сделать отрицательным.
Если число А состоит из N цифр, тогда число В не может иметь более чем N цифр. Разность будет содержать не более чем N цифр.
Найдём разность чисел А и В. Будем вычитать элементы массива В из элементов массива А с соответствующими номерами. Полученные результаты будем хранить в массиве С. Если на каком-то шаге получим C[i]<0, то к элементу C[i] прибавляем 10, а от элемента C[i+1] отнимаем 1(забираем 1 из следующего разряда).
Алгоритм вычитания двух многоразрядных чисел представлен в Приложении 2.
Умножение многоразрядных чисел
Наиболее естественным способом представления многозначного числа является запись каждого его разряда в виде отдельного элемента линейного массива (или списка, где память под цифру будет отводиться по мере надобности, в то время как в массиве приходится заранее задавать максимальное количество элементов в нем). Пусть (для удобства дальнейших действий) разряды "длинного" числа при записи в массив нумеруются с единицы, начиная с разряда единиц, т.е. цифра из разряда единиц — элемент массива с номером один, цифра из разряда десятков — элемент массива с номером два и т.д. Определим для работы с "длинными" неотрицательными числами тип данных:
Const MNax = 2000;
Type Digit = 0..9;
DlChislo = Array[1..Nmax] Of Digit;
Для
решения поставленной задачи необходимо
уметь выполнять следующие
1) ввод "длинного" числа;
2) собственно умножение двух "длинных" чисел;
3) вывод "длинного" числа;
4) определение количества цифр в записи числа.
Каждую из подзадач реализуем в виде отдельной подпрограммы. Начнем с ввода. Ввести большое число целесообразно в виде строки, а в дальнейшем преобразовать в массив цифр. В процедуре учтен указанный выше способ размещения "длинного" числа в массиве, т.е. с точки зрения пользователя число записывается как бы в обратном порядке.
(Процедура
преобразования длинного числа,
Procedure Translate(S : String; Var A : DlChislo; Var OK : Boolean);
Var I : Word;
Begin
Zero(A); I := Length(S); OK := True;
While (I >= 1) And OK Do
Begin
If S[I] In ['0'..'9']
Then A[Length(S)- I + 1]:= Ord(S[I]) - 48
Else OK := False; I := I - 1
End
End;
В процедуре вызывается подпрограмма Zero(A), назначение которой — запись нуля в каждый разряд длинного числа. Вот текст этой процедуры:
(Процедура обнуления длинного числа)
Procedure Zero(Var A : DlChislo);
Var I : Integer;
Begin
For I := 1 To NMax Do A[I] := 0;
End;
Таким образом, длинное число записано в массив, где впереди (в качестве элементов с большими номерами) стоят незначащие нули. При выполнении действий и выводе ответа они не учитываются.
Сейчас
разработаем функцию
(Функция определения количества цифр в записи длинного числа)
Function Dlina(C : DlChislo) : Integer;
Var I : Integer;
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.
Скалярно перемножить
5. Применить обратное преобразование Фурье, результатом которого будет свертка.
6.
Преобразовать свертку в
Цифры больших чисел хранятся
в целочисленном формате.
MaxLen=2k, где MaxLen – минимальная степень двойки, большая N. Например, если максимальный результат будет состоять из 1000 цифр по основанию BASE, то минимальный объем памяти MaxLen=1024, так как именно такой длины БПФ будет вычисляться.
real *LongNum1=NULL, *LongNum2=NULL;
//
Инициализацию можно делать
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;
}
Разобранная
в соответствующем разделе
Соответственно, функция скалярного произведения должна учитывать такой формат.
//
Скалярное умножение
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;
- Арифметика многочленов
- Арифметическая прогрессия
- Арифметические действия над многозначными числами
- Арифметические задачи в 5-6 классах
- Арифметические задачи как средство развития у детей логического мышления
- Арифметические и логические операции в компьютере
- Арифметические и логические основы вычислительной техники
- Аристократический отдых в Греции, на островах Ионического и Эгейского морей
- Аристотель
- Аристотель, его вклад в развитие науки
- Аристотель и Фома Аквинский
- Аристотель и Фома Аквинский
- Аристотель о государстве
- Арифметикалық өрнектеу
