Умножение и деление с округлением вверх
Как посчитать значение математического выражения?
ceil[a*b:c]
, где *, : - операции умножения и деления в множестве вещественных чисел, ceil - округление вверх.
Пример: ceil[2*3:5] = ceil[6:5] = ceil[1.2] = 2.
Ограничения в коде: a, b, c типа unsigned long long.
Примеры в коде:
assert(mul_div_ceil(2ull, 3ull, 5ull) == 2ull);
const auto max = std::numeric_limits<unsigned long long>::max();
assert(mul_div_ceil(2ull, max, 2ull) == max);
assert(mul_div_ceil(max, 2ull, 2ull) == max);
Проблемы: переполнение при наивном умножении (в смысле получения неправильного результата, а не с точки зрения языка), потеря точности при наивном "перевод в double -> деление -> умножение -> перевод в unsigned long long".
Для простого деления с округлением вверх (ceil[a:b]) это можно сделать так:
(a % b == 0) ? a / b : a / b + 1
Ответы (7 шт):
Для простого деления все делается проще -
(a+b-1)/b
Вероятность переполнения при умножении можно попробовать уменьшить, находя НОД и сокращая на него - для a и c и для b и c.
Как избежать этих неприятностей вообще без использования чего-то подлиннее - не знаю, мне кажется, что это невозможно.
Надо сначала поделить a или b на c, а после умножить на оставшийся множитель. Тогда переполнения не будет.
Может быть, не самое красивое решение, но рабочее. Решил просто символьную арифметику использовать
#include <vector>
typedef unsigned long long ull;
ull mul_div_ceil(ull a, ull b, ull c) {
//Преобразование в массив
std::vector<int> num_a,num_b;
while (a) {
num_a.push_back(a % 10);
a /= 10;
}
while (b) {
num_b.push_back(b % 10);
b /= 10;
}
//Умножение
std::vector<int> c_num(num_a.size() + num_b.size());
for (size_t i = 0; i < num_a.size(); ++i)
for (int j = 0, carry = 0; j < (int)num_b.size() || carry; ++j) {
long long cur = c_num[i + j] + num_a[i] * 1ll * (j < (int)num_b.size() ? num_b[j] : 0) + carry;
c_num[i + j] = int(cur % 10);
carry = int(cur / 10);
}
while (c_num.size() > 1 && c_num.back() == 0)
c_num.pop_back();
//Деление
ull carry = 0;
for (int i = (int)c_num.size() - 1; i >= 0; --i) {
ull cur = c_num[i] + carry * 1ll * 10;
c_num[i] = ull(cur / c);
carry = ull(cur % c);
}
while (c_num.size() > 1 && c_num.back() == 0)
c_num.pop_back();
//Обратное преобразование
ull answer = 0;
for (int i = (int)c_num.size() - 1; i >= 0; --i)
answer = answer * 10 + c_num[i];
return carry ? answer + 1 : answer;
}
Вот частичное решение для случаев, когда a,b - большие, c - маленькое (не превышает корня из UNSIGNED_LONG_LONG_MAX).
Для начала - немного математики:
a = a1 * c + a2, где a1 = a / c (целочисленное деление), a2 = a % c b = b1 * c + b2, где b1 = b / c (целочисленное деление), b2 = a % c
Распишем и упростим выражение ab / c:
ab / c = (a1c + a2) * (b1c + b2) / c = (a1b1cc + a1b2c + a2b1c + a2b2) / c = a1b1c + a1b2 + a2b1 + a2*b2/c.
Первые три слагаемых - целые числа. Если их реально посчитать без переполнения, то задача решаема. В противном случае задача не решаема, и ничего вы с ней не сделаете.
Последнее слагаемое гарантированно считается без переполнения, потому что и a2, и b2 меньше, чем c, и поэтому a2*b2 < UNSIGNED_LONG_LONG_MAX.
Теперь код:
function getResult(unsigned long long a, unsigned long long b, unsigned long long c)
{
unsigned long long a1 = a / c;
unsigned long long a2 = a % c;
unsigned long long b1 = b / c;
unsigned long long b2 = b % c;
return a1*b1*c + a1*b2 + a2*b1 + ceil(a2*b2/c);
}
Один множитель можно уменьшать пополам и результат умножения удвоить. Всё реализуется рекурсией. Тип divtype представляет целое число плюс числитель дробной части. Делитель подразумевается с.
// g++-9 -Wall -Wextra -Wpedantic -std=c++11 -Os mul_div_ceil.cpp -o mul_div_ceil
# include <cassert>
# include <limits>
# include <utility>
# include <iostream>
// div + rem / c
// rem < c
struct divtype {
unsigned long long int div ;
unsigned long long int rem ;
} ;
// (x+y/c) * 2 => x*2 + y*2 / c + перенос
divtype mul2 ( divtype a , unsigned long long int c ) {
if ( c - a . rem > a . rem )
return { a . div * 2 , a . rem * 2 } ;
return { a . div * 2 + 1 , a . rem - ( c - a . rem ) } ; }
// (x+y/c) + b/c => x + (y+b) / c + перенос
divtype add ( divtype a , unsigned long long int b ,
unsigned long long int c ) {
if ( c - a . rem > b )
return { a . div , a . rem + b } ;
return { a . div + 1 , b - ( c - a . rem ) } ; }
// a*b/c == ([a/2]*b/c)*2 + (a%2 ? b/c : 0)
// a <= b <= c/2
divtype mul_div_ceil_recursion ( unsigned long long int a ,
unsigned long long int b ,
unsigned long long int c ) {
if ( a == 1 )
return { 0 , b } ;
divtype a2 = mul_div_ceil_recursion ( a / 2 , b , c ) ;
a2 = mul2 ( a2 , c ) ;
if ( a % 2 )
a2 = add ( a2 , b , c ) ;
return a2 ; }
unsigned long long int mul_div_ceil (
unsigned long long int a ,
unsigned long long int b ,
unsigned long long int c ) {
// a == ac * c + ar
unsigned long long int ac = a / c ;
unsigned long long int ar = a % c ;
// b == bc * c + br
unsigned long long int bc = b / c ;
unsigned long long int br = b % c ;
// a * b / c ==
// (ac * c + ar) * (bc * c + br) / c ==
// ac * bc * c + ac * br + ar * bc + ar * br / c
unsigned long long int result = ac * bc * c + ac * br + ar * bc ;
if ( ar && br ) {
// ceil ( ar * br / c ) , ar < c , br < c
// Оптимизация если числитель большой
// ar * br / c == br - ( c - ar ) * br / c
bool flagceil = true ;
if ( c - ar < ar ) {
ar = c - ar ;
result += br ;
flagceil = false ; }
// ar * br / c == ar - ar * ( c - br ) / c
if ( c - br < br ) {
br = c - br ;
if ( flagceil ) {
result += ar ;
flagceil = false ; }
else {
result -= ar ;
flagceil = true ; } }
if ( br < ar )
std :: swap ( ar , br ) ;
divtype d = mul_div_ceil_recursion ( ar , br , c ) ;
if ( flagceil ) {
if ( d . rem )
++ d . div ;
result += d . div ; }
else
result -= d . div ; }
std :: cout << a << " * " << b << " / " << c << " = " << result << std :: endl ;
return result ; }
int main() {
assert(mul_div_ceil(2ull, 3ull, 5ull) == 2ull);
const auto max = std::numeric_limits<unsigned long long>::max();
assert(mul_div_ceil(2ull, max, 2ull) == max);
assert(mul_div_ceil(max, 2ull, 2ull) == max);
assert(mul_div_ceil(3ull, max-3, max) == 3);
assert(mul_div_ceil(max-3, 3ull, max) == 3);
assert(mul_div_ceil(3ull, max, max-3) == 4); }
Результат:
> ./mul_div_ceil
2 * 3 / 5 = 2
2 * 18446744073709551615 / 2 = 18446744073709551615
18446744073709551615 * 2 / 2 = 18446744073709551615
3 * 18446744073709551612 / 18446744073709551615 = 3
18446744073709551612 * 3 / 18446744073709551615 = 3
3 * 18446744073709551615 / 18446744073709551612 = 4
При умножении двух 64-битных чисел в общем случае результат должен уметь хранить 128 бит.
Если тип данных uint128 недоступен, то его можно эмулировать, написав дополнительный код.
#include <assert.h>
#include <limits>
// хранит 128-битное число вида: high_part * 2^64 + low_part
struct double_ull
{
unsigned long long high_part;
unsigned long long low_part;
// конструирование из 64-битного unsigned long long
double_ull(unsigned long long num)
{
high_part = 0;
low_part = num;
}
// конструирование 128-битного числа из старшей и младшей 64-битных частей
double_ull(unsigned long long high, unsigned long long low)
{
high_part = high;
low_part = low;
}
// оператор сдвига числа влево
double_ull operator << (unsigned short int count_bits) const
{
return double_ull
{
(high_part << count_bits) | (low_part >> (64 - count_bits)),
low_part << count_bits
};
}
// оператор сложения с unsigned long long
double_ull operator + (unsigned long long right) const
{
return (*this + double_ull{right});
}
// оператор сложения с другим double_ull
double_ull operator + (double_ull right) const
{
unsigned long long h_part = high_part + right.high_part;
unsigned long long l_part = low_part + right.low_part;
// перенос единички в high_part при переполнении в low_part
if ((l_part < low_part) || (l_part < right.low_part))
{
h_part++;
}
return double_ull { h_part, l_part };
}
// оператор умножения
double_ull operator * (unsigned long long right) const
{
// разбиваем старшую и младшую части числа на куски в виде a1 * 2^96 + a2 * 2^64 + a3 * 2^32 + a4
unsigned long long a1 = high_part >> 32;
unsigned long long a2 = high_part & 0xFFFFFFFFull;
unsigned long long a3 = low_part >> 32;
unsigned long long a4 = low_part & 0xFFFFFFFFull;
// разбиваем right на куски в виде r1 * 2^32 + r2
unsigned long long r1 = right >> 32;
unsigned long long r2 = right & 0xFFFFFFFFull;
// умножение a * r принимает вид:
// (a1 * 2^96 + a2 * 2^64 + a3 * 2^32 + a4) * (r1 * 2^32 + r2) =
// = a1*r1 * 2^128 + (a1*r2 + a2*r1) * 2^96 + (a2*r2 + a3*r1) * 2^64 + (a3*r2 + a4*r1) * 2^32 + a4*r2
// самый старший кусок (a1*r1 * 2^128) отбросим - это переполнение
// перемножаем куски чисел, конвертируем произведения в double_ull
auto a1r2 = double_ull { a1 * r2 } << 96;
auto a2r1 = double_ull { a2 * r1 } << 96;
auto a2r2 = double_ull { a2 * r2 } << 64;
auto a3r1 = double_ull { a3 * r1 } << 64;
auto a3r2 = double_ull { a3 * r2 } << 32;
auto a4r1 = double_ull { a4 * r1 } << 32;
auto a4r2 = double_ull { a4 * r2 };
return a1r2 + a2r1 + a2r2 + a3r1 + a3r2 + a4r1 + a4r2;
}
// оператор целочисленного деления
double_ull operator / (unsigned long long right) const
{
// деление принимает вид: (m1 * 2^64 + m2) / right = m1 * 2^64 / c + m2 / right;
// выполняем деление по кускам, накапливая остаток
double_ull result { low_part / right };
double_ull rem { low_part % right };
// каждая единица, суммирующаяся в "high_part": 2^64 = (2^64 - 1) + 1 = max + 1
// деление high_part выглядит как-то так:
// high_part * 2^64 / right = high_part * ( max + 1 ) / right = high_part * max / right + high_part / right
if (high_part > 0)
{
constexpr auto max = std::numeric_limits<unsigned long long>::max();
unsigned long long d1 = max / right;
unsigned long long r1 = max % right;
// первое слагаемое правой части
result = result + (double_ull{high_part} *d1);
rem = rem + (double_ull{high_part} * r1);
// второе слагаемое правой части
result = result + (high_part / right);
rem = rem + (high_part % right);
// ещё раз поделим остаток (в нём могло накопиться число большее, чем "right")
result = result + (rem / right);
}
return result;
}
// умножение двух unsigned long long
static double_ull mul(unsigned long long a, unsigned long long b)
{
return double_ull{a} * b;
}
// деление с округлением вверх
unsigned long long div_ceil(unsigned long long c)
{
return ((*this + (c - 1)) / c).low_part;
}
};
unsigned long long mul_div_ceil(unsigned long long a, unsigned long long b, unsigned long long c)
{
return double_ull::mul(a, b).div_ceil(c);
}
int main()
{
assert(mul_div_ceil(2ull, 3ull, 5ull) == 2ull);
constexpr auto max = std::numeric_limits<unsigned long long>::max();
assert(mul_div_ceil(2ull, max, 2ull) == max);
assert(mul_div_ceil(max, 2ull, 2ull) == max);
assert(mul_div_ceil(max / 2, max / 2, max) == (max / 4 + 1));
assert(mul_div_ceil(max, max, max) == max);
assert(mul_div_ceil(max / 3 * 2, max / 3 * 2, max / 3 * 2) == (max / 3 * 2));
return 0;
}
Спасибо AlexGlebe за указание на ошибку.
Код исправлен.