Как избежать переполнения в expr. A * B — C * D
Мне нужно вычислить выражение, которое выглядит так:
A*B - C*Dгде их типы: signed long long int A, B, C, D;
Каждое число может быть действительно большим (не выходя за пределы своего типа). В то время как A*B может вызвать переполнение, в то же время выражение A*B - C*D может быть действительно маленьким. Как я могу вычислить это правильно?
Например: MAX * MAX - (MAX - 1) * (MAX + 1) == 1, где MAX = LLONG_MAX - n а n — некоторое натуральное число.
Решение
Это кажется слишком тривиальным, я думаю.
Но A*B это тот, который может переполниться.
Вы могли бы сделать следующее, не теряя точности
A*B - C*D = A(D+E) - (A+F)D
= AD + AE - AD - DF
= AE - DF
^smaller quantities E & F
E = B - D (hence, far smaller than B)
F = C - A (hence, far smaller than C)
Это разложение может быть сделано дальше.
Как отметил @Gian, во время операции вычитания может потребоваться осторожность, если тип длинный без знака.
Например, в случае, если у вас есть вопрос, это займет всего одну итерацию,
MAX * MAX - (MAX - 1) * (MAX + 1)
A B C D
E = B - D = -1
F = C - A = -1
AE - DF = {MAX * -1} - {(MAX + 1) * -1} = -MAX + MAX + 1 = 1
Другие решения
Простейшим и наиболее общим решением является использование представления, которое не может переполниться, либо с помощью библиотеки длинных целых чисел (например, http://gmplib.org/) или представление с использованием структуры или массива и реализация некоторого типа длинного умножения (то есть разделение каждого числа на две 32-битные половины и выполнение умножения, как показано ниже:
(R1 + R2 * 2^32 + R3 * 2^64 + R4 * 2^96) = R = A*B = (A1 + A2 * 2^32) * (B1 + B2 * 2^32)
R1 = (A1*B1) % 2^32
R2 = ((A1*B1) / 2^32 + (A1*B2) % 2^32 + (A2*B1) % 2^32) % 2^32
R3 = (((A1*B1) / 2^32 + (A1*B2) % 2^32 + (A2*B1) % 2^32) / 2^32 + (A1*B2) / 2^32 + (A2*B1) / 2^32 + (A2*B2) % 2^32) %2^32
R4 = ((((A1*B1) / 2^32 + (A1*B2) % 2^32 + (A2*B1) % 2^32) / 2^32 + (A1*B2) / 2^32 + (A2*B1) / 2^32 + (A2*B2) % 2^32) / 2^32) + (A2*B2) / 2^32
Предполагая, что конечный результат умещается в 64 бита, вам на самом деле не нужно большинство битов R3 и ни одного из R4
Обратите внимание, что это не является стандартным, поскольку оно основано на переполнении со знаком. (GCC имеет флаги компилятора, которые позволяют это.)
Но если вы просто сделаете все расчеты в long longРезультат применения формулы непосредственно:
(A * B - C * D) будет точным, пока правильный результат вписывается в long long,
Вот обходной путь, который опирается только на определяемое реализацией поведение приведения целого числа без знака к целому числу со знаком. Но сегодня можно ожидать, что это сработает практически во всех системах.
(long long)((unsigned long long)A * B - (unsigned long long)C * D)
Это приводит к вводу unsigned long long где поведение переполнения гарантированно будет изменено стандартом. Возвращение к целому числу со знаком в конце — это часть, определяемая реализацией, но сегодня она будет работать практически во всех средах.
Если вам нужно более педантичное решение, я думаю, что вы должны использовать «длинную арифметику»
Это должно работать (я думаю):
signed long long int a = 0x7ffffffffffffffd;
signed long long int b = 0x7ffffffffffffffd;
signed long long int c = 0x7ffffffffffffffc;
signed long long int d = 0x7ffffffffffffffe;
signed long long int bd = b / d;
signed long long int bdmod = b % d;
signed long long int ca = c / a;
signed long long int camod = c % a;
signed long long int x = (bd - ca) * a * d - (camod * d - bdmod * a);
Вот мой вывод:
x = a * b - c * d
x / (a * d) = (a * b - c * d) / (a * d)
x / (a * d) = b / d - c / a
now, the integer/mod stuff:
x / (a * d) = (b / d + ( b % d ) / d) - (c / a + ( c % a ) / a )
x / (a * d) = (b / d - c / a) - ( ( c % a ) / a - ( b % d ) / d)
x = (b / d - c / a) * a * d - ( ( c % a ) * d - ( b % d ) * a)
Вы могли бы рассмотреть вычисление как наибольший общий фактор для всех ваших значений, а затем разделить их на этот фактор, прежде чем выполнять арифметические операции, а затем снова умножить. Это предполагает, что такой фактор существует, однако (например, если A, B, C а также D оказываются относительно простыми, у них не будет общего фактора).
Точно так же вы могли бы рассмотреть возможность работы с логарифмическими масштабами, но это будет немного страшно, учитывая числовую точность.
Если результат помещается в long long int, тогда выражение A * B-C * D нормально, поскольку оно выполняет арифметический мод 2 ^ 64 и даст правильный результат. Проблема в том, чтобы узнать, подходит ли результат в long long int. Чтобы обнаружить это, вы можете использовать следующий трюк с использованием double:
if( abs( (double)A*B - (double)C*D ) > MAX_LLONG )
Overflow
else
return A*B-C*D;
Проблема этого подхода заключается в том, что вы ограничены точностью мантиссы двойных чисел (54 бита?), Поэтому вам нужно ограничить произведения A * B и C * D до 63 + 54 бита (или, возможно, чуть меньше).
E = max(A,B,C,D)
A1 = A -E;
B1 = B -E;
C1 = C -E;
D1 = D -E;
затем
A*B - C*D = (A1+E)*(B1+E)-(C1+E)(D1+E) = (A1+B1-C1-D1)*E + A1*B1 -C1*D1
Вы можете написать каждое число в массиве, каждый элемент является цифрой, и выполнить вычисления как многочлены. Возьмите получившийся полином, который является массивом, и вычислите результат, умножив каждый элемент массива на 10 до степени позиции в массиве (первая позиция является наибольшей, а последняя — нулем).
Число 123 может быть выражено как:
123 = 100 * 1 + 10 * 2 + 3
для которого вы просто создаете массив [1 2 3],
Вы делаете это для всех чисел A, B, C и D, а затем умножаете их на полиномы. Получив полученный многочлен, вы просто восстанавливаете число по нему.