Блог пользователя Maxim

Автор Maxim, 6 лет назад, По-русски

Как выглядит типичный бинарный поиск в вещественных числах? Так:

float l = 0.0f, r = 1e14f;
for (int i = 0; i < iteration_count && l + eps < r; ++i)
{
    float mid = 0.5f * (l + r);
    if (check(mid))
        r = mid;
    else
        l = mid;
}

Здесь l и r — границы, выставляемые исходя из условия, check(x) — функция, которая начиная с какого-то x возвращает true, до него — false, и нужно найти этот x. Как правило, l или 0.5 * (l + r) — ответ.

Но есть несколько минусов:

  1. Количество итераций еще нужно подобрать, балансируя между точностью ответа и временем работы.
  2. Ответ с погрешностью.

Что можно с этим сделать?

Для начала разберемся, как представляются вещественные числа. Детально можно почитать, к примеру, тут и тут. Воспользуемся тем, что если проделать операцию Y = real(int(X) + 1), где X — вещественоое число, int(X) принимает вещественное число и вовзращает целое, побитовое представление которого такое же, как у Х (то есть, возвращает тот же Х, но интерпретируемый как целое число), real(X) делает противоположное — интерпретирует целое как вещественное, то Y будет равен следующему представимому вещественному числу после Х. Грубо говоря, это аналогично Y = X + eps, где eps минимально возможный и X >= 0 и Y = X - eps, если X <= 0 (в обоих случаях включен ноль, потому что в вещественных числах есть и 0, и -0). Грубо — потому что не для всех Х одинаковый eps. Должен отметить, что есть NaN, inf, максимальное положительное или отрицательное число, но это не то, о чем я хочу говорить. Можно заметить, что для вещественных L и R, где L < R, выполняется условие int(L) < int(L) + 1 < int(L) + 2 < ... < int(R), аналогично как L < L + eps < L + 2 * eps < ... < R. Значит, можно делать бинарный поиск с границами int(L), int(R). Теперь мы умеем так:

float l_real = 0.0f, r_real = 1e14f;
uint32_t l = reinterpret_cast<uint32_t&>(l_real), r = reinterpret_cast<uint32_t&>(r_real);
while (l < r)
{
    uint32_t mid = l + (r - l) / 2;
    if (check(reinterpret_cast<float&>(mid)))
        r = mid;
    else
        l = mid + 1;
}

То же самое работает для double и uint64_t. Можно использовать int для float и long long для double. Должен отметить, что стандартом не гарантировано, что float и int по 32 бита, а double и long long — по 64 бита, но я не встречался с тем, чтоб это было не так.

Какие плюсы?

  1. Максимально возможная точность.
  2. Малое количество итераций (до 32 для float и до 64 для double).

Что важно помнить?

  • l, r и mid — какие-то не несущие смысл числа, все что важно — сохранение порядка в обоих интерпретациях.
  • l + r может переполниться, поэтому l + (r - l) / 2 вместо (l + r) / 2.
  • Есть неудобство с границами с разными знаками. Для отрицательных чисел нужно перевести значение из прямого в дополнительный код при переводе из вещественного в целое, чтоб сохранить порядок. Так же можно сделать сначала проверку в нуле, чтоб у границ был одинаковый знак, или прибавить к значениям константу. На моем опыте, задач, где границы с разными знаками, крайне мало.

Не будет лишним рассмотреть картинки с распределением вещественных чисел на коородинатное прямой. Для примера, отмечу, что на отрезке от 0 до std::numeric_limints<float>::max() медиана (она же и первый mid в бинарном поиске) будет равна 1.5. То есть, от 0 до 1.5 представимо столько же чисел, сколько и от 1.5 до максимально представимого. Для float в обычном подходе понадобится 129 итераций только для того, чтоб правая граница стала меньше 1.5, если ответ, к примеру, 0 (не стоит забывать, что обычно в задачах начальная правая граница значительно меньше). То есть, в обычном подходе представимые значения отбрасываются крайне плохо.

Пример решения задачи, используя такой подход: 45343644

Такой же подход применим и к тернарному поиску (кто хочет сделать классный пример? ^_^).

  • Проголосовать: нравится
  • +110
  • Проголосовать: не нравится

»
6 лет назад, # |
Rev. 2   Проголосовать: нравится +73 Проголосовать: не нравится

It's very interesting, but here's the way I see it.

What is better:

for (int whatever = 0; whatever < 200; whatever++)

+ simple

+ easy to understand

+ fast enough 99% of the time

- not fast enough 1% of the time

+ 200 is good enough for any problem

OR:

uint32_t l = reinterpret_cast<uint32_t&>(l_real), r = reinterpret_cast<uint32_t&>(r_real);, blahblahblah

- C++ black magic

- who knows when it will break because of some compiler issue from 2012

+ save a couple iterations

- weird negative edge cases

  • »
    »
    6 лет назад, # ^ |
      Проголосовать: нравится 0 Проголосовать: не нравится

    + optimize the precision

  • »
    »
    6 лет назад, # ^ |
      Проголосовать: нравится +3 Проголосовать: не нравится

    Thanks for sharing your point of view. Let me share the way I see it then.

    + simple

    + easy to understand

    Cycle from l to r is also simple and easy to understand. Approach that in the article is also simple and easy to understand imho.

    + 200 is good enough for any problem

    32 is good enough for not less amount of problems. Let me give an example: what if we have binary search in binary search and a heavy check function? Then this approach is ~40 times (200 * 200 / (32 * 32)) faster.

    - C++ black magic

    Why did you mark it as disadvantage? Come to the dark side ;-)

    - who knows when it will break because of some compiler issue from 2012

    I do believe that the probability of that is not more than the break of simple cycle.

    - weird negative edge cases

    Not weird imho.

    • »
      »
      »
      5 лет назад, # ^ |
      Rev. 2   Проголосовать: нравится +3 Проголосовать: не нравится

      oh my god, it really works

      thanks a lot

      10 WA/TL attempts with iterative binary search bring me here

  • »
    »
    6 лет назад, # ^ |
    Rev. 3   Проголосовать: нравится +9 Проголосовать: не нравится

    C++ black magic

    who knows when it will break because of some compiler issue from 2012

    It's not even black magic, it's undefined behavior. It's already broken.

    From https://en.cppreference.com/w/cpp/language/reinterpret_cast:

    double d = 0.1;
    std::int64_t n;
    static_assert(sizeof n == sizeof d);
    // n = *reinterpret_cast<std::int64_t*>(&d); // Undefined behavior
    std::memcpy(&n, &d, sizeof d); // OK
    n = std::bit_cast<std::int64_t>(d); // also OK
    • »
      »
      »
      6 лет назад, # ^ |
        Проголосовать: нравится +10 Проголосовать: не нравится

      it's undefined behavior

      It's true because of type aliasing rule. But from the same link:

      Note that many C++ compilers relax this rule

      Also I've mentioned that there is no guarantee that float and int are 32-bit and double and long long are 64-bit. So it's already UB. But it doesn't really matter because it's almost impossible to find the compiler that misbehaves. It's individual choice between hypothetical misbehaving and more precisely/faster solution. Also it's still more probable to fail writing simple cycle as I've mentioned earlier. Actually, there is no compiler that follows the standard at 100%. If there is some discomfort — it is possible to use memset from your example.

      It's already broken

      Strange hyperbolization.

»
6 лет назад, # |
Rev. 2   Проголосовать: нравится 0 Проголосовать: не нравится
Ternary Search - Min Function

For finding the max of a unimodal function just reverse the sign.

»
5 лет назад, # |
Rev. 2   Проголосовать: нравится -8 Проголосовать: не нравится

For the regular floating-point search, instead of an iteration count, can you alternatively specify r-l > 1e-7? (Or whatever tolerance you need.)

  • »
    »
    5 лет назад, # ^ |
      Проголосовать: нравится +18 Проголосовать: не нравится

    No, it doesn't work well.

    • »
      »
      »
      5 лет назад, # ^ |
        Проголосовать: нравится +8 Проголосовать: не нравится

      Can you let me know why? It seems identical in functionality to a fixed iteration count, but easier to understand.

      • »
        »
        »
        »
        5 лет назад, # ^ |
          Проголосовать: нравится +13 Проголосовать: не нравится

        I'm not really sure why, I think it has something to do with floating point precision, but it does not work well. I remember trying to solve a couple problems like that and I was always getting WA no matter what tolerance I had, then I changed to a fixed iteration count and boom ez AC.

  • »
    »
    5 лет назад, # ^ |
      Проголосовать: нравится +18 Проголосовать: не нравится

    r — l > eps works badly on functions like y = 1000000000000x and you need to find f(x)

    f(r) — f(l) > eps works badly on functions like y = 0.00000000000001x and you need to find x

  • »
    »
    5 лет назад, # ^ |
      Проголосовать: нравится 0 Проголосовать: не нравится

    I tried to solve this problem during the 2015 ICPC Latin American regionals using r-l>EPS without any success.

    After the contest I tried using iteration count and it was accepted :(

  • »
    »
    5 лет назад, # ^ |
      Проголосовать: нравится +3 Проголосовать: не нравится

    If you are binary searching the answer and will just print $$$(L + R) / 2$$$ after the binary search is done, then you can use something like this:

    while(R - L > eps * max(1.0, L)) { ... }
    

    This checks for a relative or absolute error. It assumes that $$$0 \leq L \leq R$$$ though. A general version is quite tricky and I'm not even sure that I will get it right here:

    #define SMALLEST (L < 0 && 0 < R ? 0.0 : min(abs(L), abs(R)))
    while(R - L > eps * max(1.0, SMALLEST)) { ... }
    

    The idea is to get the minimum possible value of the answer — the minimum value of something in interval $$$[L, R]$$$.

    My advice is: just iterate a constant number of times.

    Bonus
  • »
    »
    5 лет назад, # ^ |
      Проголосовать: нравится 0 Проголосовать: не нравится

    I have also seen people do something like this:

    If every test file has a single test case and say, the time limit is around 2 seconds, then you may also write the binary search as:

    while((double)clock() / CLOCKS_PER_SEC < 1.8) {
         //Obtain mid, compute stuff and update lo and hi accordingly
    }
    

    It seems more of a magic to me, but if you can choose the constant in the while loop in a way that your execution time is within the time limit, you will get a very good precision.

»
4 года назад, # |
  Проголосовать: нравится 0 Проголосовать: не нравится

Максим спасибо!