Я делаю класс BigInt как упражнение по программированию. В base-65536 он использует вектор 2-х знаковых дополнений со знаком (так, чтобы 32-битные умножения не переполнялись. Я увеличу базу, как только получу ее полностью работоспособной).
Все основные математические операции закодированы, с одной проблемой: деление болезненно медленно с основным алгоритмом, который я смог создать. (Это работает как двоичное деление для каждой цифры отношения … Я не собираюсь публиковать это, если кто-то не хочет видеть это ….)
Вместо моего медленного алгоритма я хочу использовать Ньютона-Рафсона, чтобы найти (смещенное) обратное и затем умножить (и сместить). Я думаю, у меня есть голова вокруг основы: вы даете формулу (x1 = x0 (2 — x0 * делитель)) хорошее начальное предположение, а затем после некоторого количества итераций х сходится к обратному. Эта часть кажется достаточно простой … но я сталкиваюсь с некоторыми проблемами при попытке применить эту формулу к большим целым числам:
Проблема 1:
Потому что я работаю с целыми числами … ну … я не могу использовать дроби. Это, кажется, заставляет x всегда расходиться (делитель x0 * должен быть <2 похоже?). Моя интуиция говорит мне, что должна быть некоторая модификация уравнения, которая позволила бы ему работать целыми числами (с некоторой точностью), но я действительно изо всех сил пытаюсь выяснить, что это такое. (Моя нехватка математических навыков побеждает меня здесь …) Я думаю, что мне нужно найти какое-то эквивалентное уравнение, где вместо d есть д * [база ^ somePower]? Может ли быть какое-то уравнение, как (x1 = x0 (2 — x0 * d)) что работает с целыми числами?
Проблема 2:
Когда я использую формулу Ньютона, чтобы найти обратную величину некоторых чисел, в результате получается лишь небольшая фракция ниже того, каким должен быть ответ … напр. при попытке найти обратную величину 4 (в десятичной дроби):
x0 = 0.3
x1 = 0.24
x2 = 0.2496
x3 = 0.24999936
x4 = 0.2499999999983616
x5 = 0.24999999999999999999998926258176
Если бы я представлял числа в base-10, я бы хотел получить результат 25 (и помнить, чтобы сдвинуть произведение вправо на 2). С некоторыми взаимными ответами, такими как 1/3, вы можете просто обрезать результат, если знаете, что у вас достаточно точности. Но как я могу вытащить правильный ответ из приведенного выше результата?
Извините, если это слишком расплывчато или я прошу слишком много. Я просмотрел Википедию и все исследовательские работы, которые смог найти в Google, но чувствую, что бьюсь головой о стену. Я ценю любую помощь, которую кто-нибудь может мне дать!
…
Изменить: Получил алгоритм работает, хотя он намного медленнее, чем я ожидал. Я действительно потерял много скорости по сравнению со своим старым алгоритмом, даже для чисел с тысячами цифр … Я все еще что-то упускаю. Это не проблема с умножением, которое очень быстро. (Я действительно использую алгоритм Карацубы).
Для тех, кто заинтересован, вот моя текущая итерация алгоритма Ньютона-Рафсона:
bigint operator/(const bigint& lhs, const bigint& rhs) {
if (rhs == 0) throw overflow_error("Divide by zero exception");
bigint dividend = lhs;
bigint divisor = rhs;
bool negative = 0;
if (dividend < 0) {
negative = !negative;
dividend.invert();
}
if (divisor < 0) {
negative = !negative;
divisor.invert();
}
int k = dividend.numBits() + divisor.numBits();
bigint pow2 = 1;
pow2 <<= k + 1;
bigint x = dividend - divisor;
bigint lastx = 0;
bigint lastlastx = 0;
while (1) {
x = (x * (pow2 - x * divisor)) >> k;
if (x == lastx || x == lastlastx) break;
lastlastx = lastx;
lastx = x;
}
bigint quotient = dividend * x >> k;
if (dividend - (quotient * divisor) >= divisor) quotient++;
if (negative)quotient.invert();
return quotient;
}
А вот мой (действительно уродливый) старый алгоритм, который работает быстрее:
bigint operator/(const bigint& lhs, const bigint & rhs) {
if (rhs == 0) throw overflow_error("Divide by zero exception");
bigint dividend = lhs;
bigint divisor = rhs;
bool negative = 0;
if (dividend < 0) {
negative = !negative;
dividend.invert();
}
if (divisor < 0) {
negative = !negative;
divisor.invert();
}
bigint remainder = 0;
bigint quotient = 0;
while (dividend.value.size() > 0) {
remainder.value.insert(remainder.value.begin(), dividend.value.at(dividend.value.size() - 1));
remainder.value.push_back(0);
remainder.unPad();
dividend.value.pop_back();
if (divisor > remainder) {
quotient.value.push_back(0);
} else {
int count = 0;
int i = MSB;
bigint value = 0;
while (i > 0) {
bigint increase = divisor * i;
bigint next = value + increase;
if (next <= remainder) {
value = next;
count += i;
}
i >>= 1;
}
quotient.value.push_back(count);
remainder -= value;
}
}
for (int i = 0; i < quotient.value.size() / 2; i++) {
int swap = quotient.value.at(i);
quotient.value.at(i) = quotient.value.at((quotient.value.size() - 1) - i);
quotient.value.at(quotient.value.size() - 1 - i) = swap;
}
if (negative)quotient.invert();
quotient.unPad();
return quotient;
}
Прежде всего, вы можете реализовать деление во времени O(n^2)
и с разумной константой, так что это не намного медленнее, чем наивное умножение. Однако, если вы используете Карацуба-как алгоритм, или даже FFT-основанный алгоритм умножения, то вы действительно можете ускорить свой алгоритм деления с помощью Ньютона-Рафсона.
Итерация Ньютона-Рафсона для вычисления обратной величины x
является q[n+1]=q[n]*(2-q[n]*x)
,
Предположим, мы хотим рассчитать floor(2^k/B)
где B
является положительным целым числом. без потери общности, B≤2^k
; в противном случае частное 0
, Итерация Ньютона-Рафсона для x=B/2^k
доходность q[n+1]=q[n]*(2-q[n]*B/2^k)
, мы можем изменить это как
q[n+1]=q[n]*(2^(k+1)-q[n]*B) >> k
Каждая итерация такого типа требует только целочисленных умножений и битовых сдвигов. Это сходится к floor(2^k/B)
? Не обязательно. Тем не менее, в худшем случае, это в конечном итоге чередуется между floor(2^k/B)
а также ceiling(2^k/B)
(Докажите это!). Таким образом, вы можете использовать не очень умный тест, чтобы увидеть, если вы в этом случае, и извлечь floor(2^k/B)
, (этот «не очень умный тест» должен быть намного быстрее, чем умножения в каждой итерации; однако, было бы неплохо оптимизировать эту вещь).
Действительно, расчет floor(2^k/B)
достаточно для того, чтобы рассчитать floor(A/B)
для любых натуральных чисел A,B
, принимать k
такой, что A*B≤2^k
и проверить floor(A/B)=A*ceiling(2^k/B) >> k
,
Наконец, простая, но важная оптимизация для этого подхода состоит в том, чтобы усекать умножения (то есть вычислять только старшие биты произведения) на ранних итерациях метода Ньютона-Рафсона. Причина для этого заключается в том, что результаты ранних итераций далеки от фактора, и не имеет значения выполнять их неточно. (Уточните этот аргумент и покажите, что если вы сделаете это правильно, вы можете разделить два ≤n
-битные целые во времени O(M(2n))
при условии, что вы можете умножить два ≤k
-битные целые во времени M(k)
, а также M(x)
является возрастающей выпуклой функцией).
Если я вижу это правильно, главное улучшение — выбор хорошего начального значения для x. Зная, сколько цифр в делителе, вы знаете, где должен быть старший значащий бит, как
1/x = pow(2,log2(1/x))
1/x = pow(2,-log2(x))
1/x >= pow(2,-floor(log2(x)))
floor (log2 (x)) просто является индексом самого значительного набора битов.
Ньютон-Рафсон — это алгоритм аппроксимации, не подходящий для использования в целочисленной математике. Вы получите ошибки округления, которые приведут к возникновению проблем, которые вы видите. Вы можете решить проблему с числами с плавающей запятой, а затем посмотреть, получите ли вы целое число с точностью до указанного числа цифр (см. Следующий абзац)
Что касается второй проблемы, выберите точность (количество десятичных разрядов), которую вы хотите для точности, и округлите до этой точности. Если вы выбрали двадцать цифр точности в задаче, вы округлите до 0,25. Вам просто нужно выполнить итерацию, пока требуемые цифры точности не станут стабильными. Вообще, представление иррациональных чисел на компьютере часто приводит к неточности.