Convenient and near-optimal binary search on floating point numbers

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

Автор nor, 2 месяца назад, По-английски

https://nor-blog.codeberg.page/posts/2024-03-05-floating-point-binsearch/

nor
2 месяца назад
5

Комментарии (5)

Написать комментарий?

ToxicPie9

2 месяца назад, # |

+52

the blog: *about binary search*

the code:

template <std::size_t N_BITS>
using int_least_t = std::conditional_t<
    N_BITS <= 8, std::uint8_t,
    std::conditional_t<
        N_BITS <= 16, std::uint16_t,
        std::conditional_t<
            N_BITS <= 32, std::uint32_t,
            std::conditional_t<
                N_BITS <= 64, std::uint64_t,
                std::conditional_t<N_BITS <= 128, __uint128_t, void>>>>>;

→ Ответить

nor

2 месяца назад, # ^ |

+37

Please tell Bjarne Stroustrup to add templated int_least_t to the STL instead of the concrete type definitions :(

→ Ответить

chromate00

2 месяца назад, # |

+11

"I am pretty sure that this method has been explored before" Yes.

→ Ответить

tofael1104

2 месяца назад, # |

← Rev. 2 →

Great Blog!

→ Ответить

ToxicPie9

2 месяца назад, # |

Since most problems allow some relative error in the answer, it is usually not necessary to do the full 64 iterations. For example, if the answer is positive and we can tolerate a $$$10^{-6}$$$ relative error, then we can do binary search on doubles like this:

while (get_int(a) + (1ULL << 32) < get_int(b)) {
    auto m = std::midpoint(get_int(a), get_int(b));
    if (predicate(get_float(m))) {
        a = get_float(m);
    } else {
        b = get_float(m);
    }
}

Outside of annoying edge cases (0, NaN, inf, subnormal, etc.), the condition get_int(a) + (1ULL << 32) >= get_int(b) implies the relative error of $$$a$$$ and $$$b$$$ is at most $$$2^{-20}$$$ < $$$10^{-6}$$$.

Proof

The 64 bits in an IEEE 754 double has the following format:

 sign bit
 | exponent (11 bits)
 | |           mantissa (52 bits)
 v v           v
+-+-----------+----------------------------------------------------+
|0|10000000000|1000000000000000000000000000000000000000000000000000|
+-+-----------+----------------------------------------------------+

Except the special cases mentioned, it stores the value $$$2^{\text{exponent - 1023}} \times (1 + \text{mantissa} \times 2^{-52})$$$. The above example shows $$$2^{\text{1024 - 1023}} \times (1 + 2^{51} \times 2^{-52}) = 3.0$$$.

Let $$$e_a, m_a, e_b, m_b$$$ be the exponents and mantissas of $$$a$$$ and $$$b$$$.

Case 1: $$$e_a = e_b$$$. Then we know that $$$m_b - m_a \le 2^{32}$$$, so

$$$ \begin{split} b - a &= 2^{e_b - 1023} \times (1 + m_b \times 2^{-52}) - 2^{e_a - 1023} \times (1 + m_a \times 2^{-52}) \\ &= 2^{e_a - 1023} \times (m_b - m_a) \times 2^{-52} \\ &\le 2^{e_a - 1023} \times 2^{-20} \\ &\le a \times 2^{-20}. \end{split} $$$

Case 2: $$$e_a + 1 = e_b$$$. Then we know that $$$(2^{52} + m_b) - m_a \le 2^{32}$$$, and since $$$m_a < 2^{52}$$$, $$$m_b < 2^{32}$$$, so $$$m_a \ge m_b + 2^{52} - 2^{32} > m_b + m_b \times (2^{20} - 1) = m_b \times 2^{20}$$$. Therefore

$$$ \begin{split} b - a &= 2^{(e_a + 1) - 1023} \times (1 + m_b \times 2^{-52}) - 2^{e_a - 1023} \times (1 + m_a \times 2^{-52}) \\ &= 2^{e_a - 1023} \times (2 + 2 m_b \times 2^{-52}) - 2^{e_a - 1023} \times (1 + m_a \times 2^{-52}) \\ &= 2^{e_a - 1023} \times (1 + (2m_b - m_a) \times 2^{-52}) \\ &= 2^{e_a - 1023} \times (2^{52} + 2m_b - m_a) \times 2^{-52} \\ &\le 2^{e_a - 1023} \times (2^{32} + m_b) \times 2^{-52} \\ &< 2^{e_a - 1023} \times (2^{-20} \times (2^{52} + m_a)) \times 2^{-52} \\ &= a \times 2^{-20}. \end{split} $$$

QED.

The same could work for any $$$b < 52$$$ for $$$2^{-(52 - b)}$$$ precision. How to modify nor's template to support this is an exercise for the reader.

→ Ответить