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

#cpp #алгоритм #математика #stl

Есть алгоритм бинарного поиска по монотонной функции. Возможно ли его реализовать
на стандартных встроенных средствах из 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;
}

    

Ответы

Ответ 1

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


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


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

#include 
#include 

template 
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 
#include 

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

int main()
{
  FunctionIterator 
    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.


Ответ 2

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

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

template<
    typename Double,
    typename Func,
    typename = std::enable_if_t::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;
}