4

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

Сперва решил ускорить 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
13
  • Меня терзают смутные сомнения насчет точности вычислений... Ведь может же быть такой вариант, что число окажется больше/меньше границы, а принято будет как граничное...
    – Harry
    22 мая 2017 в 4:02
  • @Harry больше меньше какой границы ? если речь о числе x то во всяком случае точность, я думаю, что здесь не принципиальна так как логарифм должен давать одно и тоже число для целого набора чисел например лежащих в диапазоне 10^5 до 99^5 - число 5, а что касается быстроты операций с целыми числами над числами с плавающей точкой то тут да, понятно быстрее, но это уже другой вопрос
    – ampawd
    22 мая 2017 в 10:06
  • Ну сами подумайте - представим, что у вас в битах - ровно в какой-то там степени. Вы уверены, что логарифм не вернет значение на 1 меньше, например?...
    – Harry
    22 мая 2017 в 10:08
  • @Harry это зависит в какую сторону округлять я думаю, мне непонятна какая тут связь между точностями вычислений и производительностью ?
    – ampawd
    22 мая 2017 в 10:11
  • @Harry во всяком случае с моими "не точными вычислениями" неправильных результатов я не увидел запуская это решение на сайте
    – ampawd
    22 мая 2017 в 10:15

2 ответа 2

5

Мне кажется, что главная проблема заключается в том, что вы много раз вызываете функцию 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;)

6
  • Произведение может запросто нарушить точность. Впрочем, поскольку ее нет и в исходном решении, то...
    – Harry
    22 мая 2017 в 4:25
  • идеи, конечно хорошие, касательно оптимизаций, но вот запуская последний ваш вариант почему то выдаётся 26.0014 , когда правильный ответ 22 для теста 2 0 0 0 2324 0 2332 0 0
    – ampawd
    22 мая 2017 в 13:13
  • @ampawd, да, в ответе действительно была ошибка. Я исправил, теперь на этом тесте результат 45.2673, что очень близко к правильному ответу, который, кстати не равен 44.
    – diralik
    22 мая 2017 в 14:06
  • 44 будет если брать целочисленное значение логарифма округлённое вниз и умноженное на 2 поэтому это тоже верный ответ, а там вещественный логарифм считается
    – ampawd
    22 мая 2017 в 14:10
  • спасибо за идеи, старание
    – ampawd
    22 мая 2017 в 14:13
4

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

Я бы делал так - как минимум для повышения точности (потому что ваш 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

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

3
  • на самом деле я уже решил проблему, дело было в функции log10_fast там цикл зацикливался когда x == 0 отсюда и тормоза бесконечные - так то эта функция значительно повысила производительность, на серверах этих всего 0.358 сек, а с стардартным log10 по времени вылетало, но что-то мне подсказывает что если использоваться вашу таблицу можно ещё выиграть во времени :)
    – ampawd
    22 мая 2017 в 17:43
  • 1
    @ampawd Еще раз - главное, что оно обеспечит ТОЧНОСТЬ - см. мой дополненный ответ: в нем конкретные примеры неприятностей. То, что вы на них в своих расчетах не попадали - не моя вина :)
    – Harry
    22 мая 2017 в 18:43
  • ну такое возможно да, не спорю
    – ampawd
    22 мая 2017 в 19:03

Ваш ответ

By clicking “Отправить ответ”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Всё ещё ищете ответ? Посмотрите другие вопросы с метками или задайте свой вопрос.