Convenient and near-optimal binary search on floating point numbers

nor's blog

By nor, 2 months ago, In English

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

nor
2 months ago
5

Comments (5)

Write comment?

ToxicPie9

2 months ago, # |

+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>>>>>;

→ Reply

nor

2 months ago, # ^ |

+37

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

→ Reply

chromate00

2 months ago, # |

+11

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

→ Reply

tofael1104

2 months ago, # |

← Rev. 2 →

Great Blog!

→ Reply

ToxicPie9

2 months ago, # |

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.

→ Reply