оптимизировать вычисления

305
24 мая 2017, 08:09

при реализации этой задачи столкнулся с проблемой медленных вычислений

Сперва решил ускорить log10 - не помогло, хоть и вычисляет быстрее чем log10 из стандартной библиотеки - (см. сравнение в коцне вопроса), но использование таблицы степеней 10 оказалось не сильно хорошей идей, хотя бин поиск должен был повлиять...

Вобщем, не могу понять что ещё и главное как можно оптимизировать кроме вычисления десятичного логарифма

#include <iostream>
#include <vector>
#include <math.h>
#include <stdint.h>
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<std::vector<uint32_t>>& 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<std::vector<uint32_t>> a(n, std::vector<uint32_t>(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) << '\n';
    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<long double>(rand()) );
}     
high_resolution_clock::time_point t2 = high_resolution_clock::now();
duration<double> dur = duration_cast<milliseconds>( t2 - t1 );
std::cout << dur.count() << '\n';  // 6.374 sec

s = 0;
t1 = high_resolution_clock::now();
for (unsigned i = 0; i < 1e+8; ++i)
{
    s += log10_fast( static_cast<long double>(rand()) );
}    
t2 = high_resolution_clock::now(); 
dur = duration_cast<milliseconds>( t2 - t1 );
std::cout << dur.count() << '\n'; // 5.907 sec
Answer 1

Мне кажется, что главная проблема заключается в том, что вы много раз вызываете функцию log10. Давайте распишем сумму логарифмов как логарифм произведения:

Если чисел немного, так что их произведение не вызовет переполнение, то можно так и посчитать.

Другая проблема заключается в использовании дробных чисел, они не так быстро перемножаются, как целые. Хочу предложить приближённое решение в целых числах. Заметим, что если, например, a[i][1] xor a[j][1] не равно нулю, то a[i][3] xor a[j][3] и a[i][4] xor a[j][4] можно не считать, так как их добавка к xor'у будет очень маленькой. Рассмотрим следующий алгоритм:

  1. Для каждой пары (i, j), приближённо считаем A_i xor A_j в виде x * 2^k, где x, k --- некоторые целые числа, причём 0 <= x < 2^32.
  2. Перемножаем полученные значения следующим образом: (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)
  3. По факту мы посчитали не десятичный логарифм, а двоичный, чтобы получить десятичный логарифм, нужно домножить на log_10(2).

Собственно, код:

#include <iostream>
#include <vector>
#include <cassert>
#include <cmath>
#include "/home/dima/C++/debug.h"
using namespace std;
double compute(const vector<uint64_t[2]> &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<uint64_t[2]> 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 <iostream>
#include <vector>
#include <cassert>
#include <cmath>
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<uint64_t[2]> &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<uint64_t[2]> 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;)

Answer 2

У меня не готовый код, но...

Я бы делал так - как минимум для повышения точности (потому что ваш double никак не отловит точное значение, например, того же 1039). Сделал бы 128-разрядное число - как 

unsigned long a[4];

Далее, вот такую табличку степеней 10 в таком точном представлении -

unsigned long p10[39][4] =
{
    { 0x00000001, 0x00000000, 0x00000000, 0x00000000 },   // 1
    { 0x0000000a, 0x00000000, 0x00000000, 0x00000000 },   // 10
    { 0x00000064, 0x00000000, 0x00000000, 0x00000000 },   // 100
    { 0x000003e8, 0x00000000, 0x00000000, 0x00000000 },   // 1000
    { 0x00002710, 0x00000000, 0x00000000, 0x00000000 },   // 10000
    { 0x000186a0, 0x00000000, 0x00000000, 0x00000000 },   // 100000
    ...
    { 0x00000000, 0x098a2240, 0x5a86c47a, 0x4b3b4ca8 },   // 10000000000....0
    { 0x00000000, 0x5f655680, 0x8943acc4, 0xf050fe93 },   // 10000000000....00

(полностью - здесь: http://vpaste.net/RljZF). Ну, или если удобнее - то наоборот, от старшего к младшему. Далее я бы написал простую функцию сравнения вот таких 128-разрядных чисел - очень просто, начиная со старшего - и вперед...

И логарифмировал бы без всяких переводов в doubleы - аналогично вашему поиску при взятии логарифма. Причем надо еще померить, дает ли что-то при таком небольшом количестве бинарный поиск или нет. Можно поиграться, начиная с поиска по старшему элементу. При несовпадении - сразу определяется логарифм, при совпадении - переходим к следующему и так далее... Написать у вас, я думаю, проблем не составит.

Плюсы - точность, не используется арифметика с плавающей точкой.

Update о точности...

Рассмотрим значение 0x00000000 0x00000000 0x00038d7e 0xa4c67FFE, т.е. число 999999999999998.

Очевидно, что значение его логарифма, floor до целого - 14.

Теперь вычисляем ваше значение - 

double x = pow(2,32)*0x00038d7e + 0xa4c67FFE;

VC++ 2015 дает для

double x = pow(2,32)*0x00038d7e + 0xa4c67FFE;
printf("%.10lf\n",x);
printf("%.10lf\n",log10(x));
int l = log10(x);
printf("%d\n",l);

следующие результаты:

999999999999998.0000000000
15.0000000000
15

Вы можете возразить, что логарифм вы считаете не так... но проверьте сами, что число 0x00000000 0x00000000 0x0DE0B6B3 0xA763FFF8 - т.е. 999999999999999992 - даст при вычислении вашим способом - с умножением на pow(2,...) - число 

double x == 1000000000000000000.00000

Так что значения логарифмов у вас все равно для некоторых чисел окажутся неверными.

READ ALSO
Как сделать блок схему по функции кода С++? [требует правки]
c++

Как сделать блок схему по функции кода С++? [требует правки]

 310
Какой алгоритм использовать?
c++

Какой алгоритм использовать?

Есть n городов, которые необходимо соединить дорогами, так, чтобы можно было добраться из любого города в любой другой (напрямую или через...

 258
Максимальный размер QString в Qt
c++

Максимальный размер QString в Qt

Сколько вообще символов в QString вставить? Я смотрю там есть какой то лимит или я незнаю, короче часть символов просто режется при выводе или...

 681
Передача файлов через TCP
c++

Передача файлов через TCP

Использую TCP-клиент и TCP-серверПередаю строку на сервер

 227