при реализации этой задачи столкнулся с проблемой медленных вычислений
Сперва решил ускорить log10 - не помогло, хоть и вычисляет быстрее чем log10 из стандартной библиотеки - (см. сравнение в коцне вопроса), но использование таблицы степеней 10 оказалось не сильно хорошей идей, хотя бин поиск должен был повлиять...
Вобщем, не могу понять что ещё и главное как можно оптимизировать кроме вычисления десятичного логарифма
#include
#include
#include
#include
const long double _10_POWERS[40] =
{
1e+0, 1e+1, 1e+2, 1e+3, 1e+4, 1e+5, 1e+6, 1e+7, 1e+9, 1e+10,
1e+11, 1e+12, 1e+13, 1e+13, 1e+14, 1e+15, 1e+16, 1e+17, 1e+18, 1e+19,
1e+20, 1e+21, 1e+22, 1e+23, 1e+24, 1e+25, 1e+26, 1e+27, 1e+28, 1e+29,
1e+30, 1e+31, 1e+32, 1e+33, 1e+34, 1e+35, 1e+36, 1e+37, 1e+38, 1e+39
};
static inline uint32_t log10_fast(long double x)
{
//uint32_t res = 0;
int l = 0, r = 40 - 1;
while (l <= r)
{
int mid = l + ((r -l) >> 1);
if ( x >= _10_POWERS[mid] && x < _10_POWERS[mid + 1] )
{ return mid; }
if (x >= _10_POWERS[mid])
l = mid;
else
r = mid;
}
return 0;
};
uint32_t compute(int n, std::vector>& a)
{
long double x = 0.0f;
uint32_t s = 0;
const long double _2_96 = pow(2, 96);
const long double _2_64 = pow(2, 64);
const long double _2_32 = pow(2, 32);
for (int i = 0; i < n; ++i)
{
for (int j = i + 1; j < n; ++j)
{
x = _2_96 * (a[i][0] ^ a[j][0]);
x += _2_64 * (a[i][1] ^ a[j][1]);
x += _2_32 * (a[i][2] ^ a[j][2]);
x += (a[i][3] ^ a[j][3]);
s += log10_fast(x);
}
}
return 2 * s;
}
int main(int argc, char const *argv[])
{
int n;
std::cin >> n;
std::vector> a(n, std::vector(4));
for (int i = 0; i < n; ++i)
{
for (int j = 0; j < 4; ++j)
{
std::cin >> a[i][j];
}
}
std::cout << compute(n, a) << '
';
return 0;
}
Какие есть идеи по поводу оптимизаций тут ?
Может есть другой более быстрый способ вычислить log10 ? Или может дело не в логарифме ?
p.s.
сравнение log10 и log10_fast
uint32_t s = 0;
high_resolution_clock::time_point t1 = high_resolution_clock::now();
for (unsigned i = 0; i < 1e+8; ++i)
{
s += log10( static_cast(rand()) );
}
high_resolution_clock::time_point t2 = high_resolution_clock::now();
duration dur = duration_cast( t2 - t1 );
std::cout << dur.count() << '
'; // 6.374 sec
s = 0;
t1 = high_resolution_clock::now();
for (unsigned i = 0; i < 1e+8; ++i)
{
s += log10_fast( static_cast(rand()) );
}
t2 = high_resolution_clock::now();
dur = duration_cast( t2 - t1 );
std::cout << dur.count() << '
'; // 5.907 sec
Ответ
Мне кажется, что главная проблема заключается в том, что вы много раз вызываете функцию log10. Давайте распишем сумму логарифмов как логарифм произведения:
Если чисел немного, так что их произведение не вызовет переполнение, то можно так и посчитать.
Другая проблема заключается в использовании дробных чисел, они не так быстро перемножаются, как целые. Хочу предложить приближённое решение в целых числах. Заметим, что если, например, a[i][1] xor a[j][1] не равно нулю, то a[i][3] xor a[j][3] и a[i][4] xor a[j][4] можно не считать, так как их добавка к xor'у будет очень маленькой. Рассмотрим следующий алгоритм:
Для каждой пары (i, j), приближённо считаем A_i xor A_j в виде x * 2^k, где x, k --- некоторые целые числа, причём 0 <= x < 2^32
Перемножаем полученные значения следующим образом: (x1 * 2^k1) * (x2 * 2^k2) = (x1 * x2) * 2^(k1+k2) = x * 2^(k1+k2), причём 0 <= x <= 2^64. Представляем x в виде x=y * 2^m', причём 0 <= y < 2^32. Итак, (x1 * 2^k1) * (x2 * 2^k2) = y * 2^(k1+k2+m)
По факту мы посчитали не десятичный логарифм, а двоичный, чтобы получить десятичный логарифм, нужно домножить на log_10(2).
Собственно, код:
#include
#include
#include
#include
#include "/home/dima/C++/debug.h"
using namespace std;
double compute(const vector &a) {
static const uint64_t two_power_32 = 1ull << 32;
int n = a.size();
// текущий накопленный результат равен value * 2^power_index
uint64_t value = 1;
int power_index = 0;
for (int i = 0; i < n; ++i) {
for (int j = 0; j < i; ++j) {
uint64_t xor1 = a[i][0] ^a[j][0];
uint64_t xor2 = a[i][1] ^a[j][1];
uint64_t value_current;
if (xor1 == 0) {
value_current = xor2;
} else if (xor1 >= two_power_32) {
value_current = xor1;
power_index += 64;
} else {
assert(0 <= xor1 && xor1 < two_power_32);
value_current = (xor1 << 32) + (xor2 >> 32);
power_index += 32;
}
while (value_current >= two_power_32) {
value_current /= 2;
++power_index;
}
assert(0 <= value_current && value_current < two_power_32);
assert(0 <= value && value < two_power_32);
value *= value_current;
while (value >= two_power_32) {
value /= 2;
++power_index;
}
}
}
// result = log10(value * 2^power_index)
// result = log10(value) + log10(2^power_index)
// result = log10(value) + power_index * log10(2)
double result = log10(value) + power_index * log10(2);
return result * 2;
}
int main() {
freopen("input.txt", "r", stdin);
int n;
cin >> n;
vector a(n);
for (int i = 0; i < n; ++i) {
for (int j = 0; j < 2; ++j) {
uint32_t ai1, ai2;
cin >> ai1 >> ai2;
a[i][j] = (uint64_t(ai1) << 32) + ai2;
}
}
cout << compute(a) << endl;
return 0;
}
К сожалению, я не сравнивал производительность, но я искренне верю, что это работает быстрее, чем n^2 раз вычислять log10
Обновление: я тут потестировал, при n=5000 моя реализация чуть медленнее вашей оригинальной. Всё дело в этих циклах:
while (value >= two_power_32) {
value /= 2;
++power_index;
}
Их можно переписать разными способами, вот вариант для GCC:
static const uint64_t two_power_32 = 1ull << 32;
inline void divide_until_less_then_two_power_32(uint64_t &value, int &power_index) {
// Эта функция эквивалентна этим строчкам:
// while (value >= two_power_32) {
// value /= 2;
// ++power_index;
// }
if (value < two_power_32) {
return;
}
int power_index_delta = 32 - __builtin_clzll(value);
power_index += power_index_delta;
value >>= power_index_delta;
assert(0 <= value && value < two_power_32);
}
Полный код:
#include
#include
#include
#include
using namespace std;
#include "/home/dima/C++/debug.h"
static const uint64_t two_power_32 = 1ull << 32;
inline void divide_until_less_then_two_power_32(uint64_t &value, int &power_index) {
// Эта функция эквивалентна этим строчкам:
// while (value >= two_power_32) {
// value /= 2;
// ++power_index;
// }
if (value < two_power_32) {
return;
}
int power_index_delta = 32 - __builtin_clzll(value);
power_index += power_index_delta;
value >>= power_index_delta;
assert(0 <= value && value < two_power_32);
}
double compute(const vector &a) {
int n = a.size();
// текущий накопленный результат равен value * 2^power_index
uint64_t value = 1;
int power_index = 0;
for (int i = 0; i < n; ++i) {
for (int j = 0; j < i; ++j) {
uint64_t xor1 = a[i][0] ^a[j][0];
uint64_t xor2 = a[i][1] ^a[j][1];
uint64_t value_current;
if (xor1 == 0) {
value_current = xor2;
} else if (xor1 >= two_power_32) {
value_current = xor1;
power_index += 64;
} else {
assert(0 <= xor1 && xor1 < two_power_32);
value_current = (xor1 << 32) + (xor2 >> 32);
power_index += 32;
}
divide_until_less_then_two_power_32(value_current, power_index);
assert(0 <= value_current && value_current < two_power_32);
assert(0 <= value && value < two_power_32);
value *= value_current;
divide_until_less_then_two_power_32(value, power_index);
}
}
// result = log10(value * 2^power_index)
// result = log10(value) + log10(2^power_index)
// result = log10(value) + power_index * log10(2)
double result = log10(value) + power_index * log10(2);
return result * 2;
}
int main() {
freopen("input.txt", "r", stdin);
int n;
cin >> n;
vector a(n);
for (int i = 0; i < n; ++i) {
for (int j = 0; j < 2; ++j) {
uint32_t ai1, ai2;
cin >> ai1 >> ai2;
a[i][j] = (uint64_t(ai1) << 32) + ai2;
}
}
cout << compute(a) << endl;
return 0;
}
Если я всё правильно посчитал, то эта версия работает в два раза быстрее.
Обновление 2: исправил ошибку (добавил строчку power_index += 32;)
