алгоритм - Бинарный поиск по монотонной функции стандартными средствами C++ - возможен ли?


5

Есть алгоритм бинарного поиска по монотонной функции. Возможно ли его реализовать на стандартных встроенных средствах из C++ (STL)?

Например, бинарный поиск по массиву я нашёл, он реализован и реализован хорошо. Однако, про вещественный бинарный поиск ни слова.

Сам примерный алгоритм, про который идёт речь:

#define EPS 1E-9

double f(double x)
{
   ///some monotonically increasing function on [a, b], for example f(x) = x^3:
   return x*x*x;
}

double binarySearch(double C, double a, double b)
{
   double low = a, high = b;
   double mid;
   while(abs(low-high) > EPS)
   {
      mid = low + (high - low) / 2;
      if f(mid) < C
          low = mid;
      else
          high = mid;
   }
   return mid;
}
Источник
  •  116
  •  2
  • 7 янв 2019 2019-01-07 11:49:14
@VTT А так же, как я знаю, очень часто во избежание зацикливаний или слишком долгого времени работы используют в алгоритме не while, а for с захардкоженным количеством итераций. — 6 янв 20192019-01-06 19:07:02.000000
@VTT Поэтому, как я знаю, в таких случаях идёт не сравнение F(m) == need, а сравнение need - F(m) < EPS. — 6 янв 20192019-01-06 19:05:20.000000
@HolyBlackCat Да. Но даже если функция непрерывна и монотонна это не значит, что она примет все значения на интервале из-за погрешностей вычислений с плавающей точкой. Это же не математическая сферическая функция в вакууме. — 6 янв 20192019-01-06 19:00:44.000000
@VTT Похоже, подразумеваются непрерывные монотонные функции. — 6 янв 20192019-01-06 18:58:36.000000
"речь идёт о поиске аргумента функции, при котором она принимавет нужное значение" - а как вы предлагаете рассматривать ситуацию, когда значение лежит в нужном интервале, но функция его никогда не принимает? — 6 янв 20192019-01-06 18:57:33.000000

2 ответа

4

Это возможно, но не все готовые средства присутствуют в стандартной библиотеке С++.

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

  2. Нам нужен итератор произвольного доступа по "вирутальным" последовательностям, т.е. по результатам вызова функции, генерируемым "на лету", а не по физическим значениям в памяти. Такие итераторы или нечто подобное, насколько я знаю, присутствуют в Boost, но не в нынешней стандартной библиотеке С++. (Что странно, ибо концепция это весьма натуральная и востребованная.)

  3. Понятие строгого равенства не очень хорошо работает (мягко выражаясь) с плавающей арифметикой, поэтому для бинарного поиска придется использовать граничные алгоритмы: std::lower_bound, std::upper_bound, std::equal_range...

Вот пример такого итератора, написанного "на коленке" (вполне может быть, что некоторые операции "избыточны", т.е. не будут востребованными используемыми алгоритмами, но я не задавался этим вопросом, а сразу реализовал "побольше")

#include <cstdint>
#include <iterator>

template <typename F, typename ARG, std::uintmax_t PRECISION>
struct FunctionIterator
{
  using iterator_category = std::random_access_iterator_tag;
  using value_type = ARG;
  using difference_type = std::intmax_t;
  using pointer = ARG *;
  using reference = ARG;

  F f;
  std::intmax_t x;

  FunctionIterator(F f, ARG x) : f(std::move(f)), x((std::intmax_t) (x * PRECISION))
    {}

  ARG arg() const
    { return (ARG) x / PRECISION; }
  ARG operator *() const
    { return f(arg()); }

  FunctionIterator &operator++()   
    { ++x; return *this; }
  FunctionIterator &operator++(int)   
    { FunctionIterator old = *this; ++*this; return old; }
  FunctionIterator &operator +=(std::intmax_t rhs)   
    { x += rhs; return *this; }
  friend FunctionIterator operator +(const FunctionIterator &lhs, std::intmax_t rhs)
    { return FunctionIterator(lhs) += rhs; }

  FunctionIterator &operator--()   
    { --x; return *this; }
  FunctionIterator &operator--(int)   
    { FunctionIterator old = *this; --*this; return old; }
  FunctionIterator &operator -=(std::intmax_t rhs)   
    { x -= rhs; return *this; }
  friend FunctionIterator operator -(const FunctionIterator &lhs, std::intmax_t rhs)
    { return FunctionIterator(lhs) -= rhs; }

  friend std::intmax_t operator -(const FunctionIterator &lhs, const FunctionIterator &rhs)
    { return lhs.x - rhs.x; }

  friend bool operator ==(const FunctionIterator &lhs, const FunctionIterator &rhs)
    { return lhs.x == rhs.x; }
  friend bool operator !=(const FunctionIterator &lhs, const FunctionIterator &rhs)
    { return lhs.x != rhs.x; }
  friend bool operator <(const FunctionIterator &lhs, const FunctionIterator &rhs)
    { return lhs.x < rhs.x; }
};

Заполучив такой итератор, мы сможем передавать его в стандартные алгоритмы, в т.ч. для выполнения поиска.

Пример использования его для выполнения бинарного поиска на отрезке, скажем [-100, 100] c точностью 0.001

#include <algorithm>
#include <iostream>

double f(double x)
{
  return x * x * x;
}

int main()
{
  FunctionIterator<double (*)(double), double, 1000> 
    a(f, -100), 
    b(f, 100),
    c = std::lower_bound(a, b, 5.0);
  std::cout << c.arg() << " " << *c << std::endl;
}

Вывод

1.71 5.00021

Увеличив множитель до 100000 получим

1.70998 5.00004

И т.д., стараясь следить за опасностью переполнения std::intmax_t.

  • 7 янв 2019 2019-01-07 06:28:59
Да, написать такое извращение надо постараться... Наверное, такое было только во времена XT, когда в играх таблицы синусов считали заранее и в виде дробей - как два целых числа. — 7 янв 20192019-01-07 11:48:41.000000
5

Нет, в стандартной библиотеке такого нет.

Но, откровенно говоря, ваш алгоритм несколько, гм, удивляет. Тогда уж имеет смысл делать бинарный поиск для решения f(x) == 0. Ну, для обобщенности... я тут набросал на коленке :)

template<
    typename Double,
    typename Func,
    typename = std::enable_if_t<std::is_floating_point<Double>::value>
    >
    Double binEq(Double left, Double right, Double eps, Func f)
{ 
    using std::swap;
    Double x = left;

    if (left > right) swap(left,right);

    Double fl = f(left), fr = f(right);

    if (fl*fr > 0) throw std::exception("Wrong range");

    while ( right - left > eps )
    {
        x = (left + right)/Double(2.0);

        if (fl * f(x) < 0) right = x;
        else               left  = x;
    }
    return x;
}
  • 6 янв 2019 2019-01-06 19:45:16
Студия - в своем репертуаре. :) У std::exception не должно быть конструктора с параметром-строкой. (Минус не мой.) — 6 янв 20192019-01-06 21:01:45.000000
@HolyBlackCat В VC++ компилируется, я под ним ваял... Тут обработка может быть любой - просто сообщение, что решение не гарантировано, так что и пробовать не будем :( — 6 янв 20192019-01-06 19:43:18.000000
Из-за std::exception не компилируется. Нужно вместо него std::runtime_error или что-то подобное. — 6 янв 20192019-01-06 19:01:23.000000