// read_ints.cpp
// Чтение последовательности целых чисел неизвестного размера в динамический массив,
// организуемый средствами стандартной библиотеки C (malloc, realloc, free).
#include <cstdlib>
#include <iostream>
using namespace std;

// Возвращает число прочитанных чисел.
size_t read_ints(int **ints, size_t max_size = 1 << 20)
{
  size_t cap = max_size < 16 ? max_size : 16; // capacity
  auto a = (int*)malloc(cap * sizeof(int));
  *ints = a;

  if (!a)
    return 0;

  size_t i = 0; // size
  for (int x;; ++i)
  {
    if (i == cap)
    {
      if (cap == max_size)
        break;

      cap *= 2;
      if (max_size < cap)
        cap = max_size;

      if (auto r = (int*)realloc(a, cap * sizeof(int)))
        a = r;
      else
        break;
    }

    if (cin >> x)
      a[i] = x;
    else
      break;
  }

  return i;
}

/////////////////////////////////////////////////////////////////////
int main()
{
  while (true)
  {
    size_t mxs = 0;
    cin >> mxs;

    int *a = nullptr;
    size_t sz = read_ints(&a, mxs);
    cin.clear();
    cin.ignore(cin.rdbuf()->in_avail());
    cin.sync();

    for (size_t i = 0; i < sz; ++i)
      cout << a[i] << '\n';
    free(a);
  }
  
  return 0;
}

// pint.cpp
// Примитивный способ округления до целых.
// В каком диапазоне значений x эта функция работает корректно?
float pint(float x)
{
  x += 12582912.f; // 1.5f * (1 << 23)
  x –= 12582912.f;
  return x;
}

// ppow.cpp
// Эффективный алгоритм возведения в целочисленную степень.
// Логарифмическое число умножений.
// std::pow(Float, Integer) реализуется аналогично.

// Неотрицательный показатель.
float ppow(float x, unsigned n)
{
  float s = 1.f;
  for (; n != 0; n /= 2, x *= x)
    if (n % 2 != 0)
      s *= x;
  return s;
}

// Отрицательный или неотрицательный показатель.
float ppow(float x, int n)
{
  return n < 0 ? 1.f / ppow(x, -n) 
    : ppow(x, n);
}

// hypot.cpp
// Testing different hypot implementations (in float).
// V.0: standard implementation;
// V.1: just sqrt(v.v);
// V.2: usual non-overflowing code with division;
// V.3: (float)sqrt(double(v).double(v)).
#include <utility>
#include <cmath>

inline float hypot_v0(float x, float y)
{
  return std::hypotf(x, y);
}

inline float hypot_v0(float x, float y, float z)
{
  return std::hypotf(std::hypotf(x, y), z);
}

inline float hypot_v1(float x, float y)
{
  return std::sqrt(x * x + y * y);
}

inline float hypot_v1(float x, float y, float z)
{
  return std::sqrt(x * x + y * y + z * z);
}

inline float hypot_v2(float x, float y)
{
  x = std::abs(x);
  y = std::abs(y);
  
  const float
    min_ = x < y? x: y,
    max_ = x < y? y: x;

  if (max_ == 0.f || std::isinf(max_))
    return max_;
  
  const float q = min_ / max_;
  return max_ * std::sqrt(1.f + q * q);
}

inline float hypot_v2(float x, float y, float z)
{
  x = std::abs(x);
  y = std::abs(y);
  z = std::abs(z);

  if (y < x)
    std::swap(x, y);
  if (z < y)
    std::swap(z, y);

  // Now z is max.
  if (z == 0.f || std::isinf(z))
    return 0.f;

  const float qx = x / z, qy = y / z;
  return z * std::sqrt(1.f + (qx * qx + qy * qy));
}

inline float hypot_v3(float x, float y)
{
  const double x_ = x, y_ = y;
  return (float)std::sqrt(x_ * x_ + y_ * y_);
}

inline float hypot_v3(float x, float y, float z)
{
  const double x_ = x, y_ = y, z_ = z;
  return (float)std::sqrt(x_ * x_ + y_ * y_ + z_ * z_);
}


////////////////////////////////////////////////////////////////////////////////////////////////////
// Testing
#include <cfloat>
#include <cstdint>
#include <iostream>
#include <chrono>
#include <vector>

/// Perform action and measure the time difference using Clock clock.
template <class Clock, class Action>
double timing(Action && action)
{
  const auto t0 = Clock::now();
  action();
  const std::chrono::duration<double> p = Clock::now() - t0;
  return p.count();
}

/// Perform timings and output the results.
template <class Clock = std::chrono::high_resolution_clock>
class Timings
{
  std::vector<std::pair<std::string, double>> results;

public:
  Timings() = default;

  template <class Action>
  Timings(std::string msg, Action && action)
  {
    (*this)(std::move(msg), std::forward<Action>(action));
  }

  template <class Action>
  Timings& operator()(std::string msg, Action && action)
  {
    results.emplace_back(std::move(msg),
      timing<Clock>(std::forward<Action>(action)));
    return *this;
  }

  /// Output the accumulated results.
  std::ostream& operator()(std::ostream & os) const
  {
    using namespace std;
    if (results.empty())
      return os << "No results" << endl;

    const auto npt = results.front().second;
    os << results.front().first << npt << "s\n";
    
    for (auto p = begin(results), p_end = end(results); ++p != p_end;)
      os << p->first << p->second << "s, " << npt / p->second << " times\n";
    return os << endl;
  }
};


std::uint32_t ulp_dist(float a, float b)
{
  const auto
    a_ = *reinterpret_cast<const uint32_t*>(&a),
    b_ = *reinterpret_cast<const uint32_t*>(&b);
  return a_ < b_? b_ - a_: a_ - b_;
}


int main()
{
  using namespace std;
  using Clock = chrono::high_resolution_clock;
  using Clock_period = Clock::period;
  cout << (double)Clock_period::den / Clock_period::num << " tps\n";
  cout.precision(4);
  
  // 2D versions.
  const int N = 12'000'000;
  vector<float> src(N);
  for (int i = 0; i < N; i += 2)
  {
    const double phi = 1e-5 * i;
    src[i] = (float)sin(phi);
    src[i + 1] = (float)cos(phi);
  }
  
  float s[4]{};
  #define HT(n)\
    float sum = 0.f;\
    for (int i = 0; i < N; i += 2)\
      sum += hypot_v##n(src[i], src[i+1]);\
    s[n] = sum;\
  // end define
  
  cout << "2D versions\n";
  Timings<Clock>
    ("v0: ", [&]{ HT(0) })
    ("v1: ", [&]{ HT(1) })
    ("v2: ", [&]{ HT(2) })
    ("v3: ", [&]{ HT(3) })
    ("v0: ", [&]{ HT(0) })
    (cout);
    
  cout << '\n';
  
  Timings<Clock>
    ("v0: ", [&]{ HT(0) })
    ("v3: ", [&]{ HT(3) })
    ("v0: ", [&]{ HT(0) })
    ("v2: ", [&]{ HT(2) })
    ("v1: ", [&]{ HT(1) })
    (cout);
    
  cout << '\n';
  for (int i = 0; i < 4; ++i)
    cout << '\t' << s[i];
  cout << '\n';
  
  // 3D versions.
  for (int i = 0; i < N; i += 3)
  {
    const double phi = 1e-5 * i;
    const double psi = 1.5e-5 * i;
    src[i] = (float)sin(phi);
    
    const auto cos_phi = cos(phi);
    src[i + 1] = (float)(cos_phi * sin(psi));
    src[i + 2] = (float)(cos_phi * cos(psi));
  }
  
  #undef HT
  #define HT(n)\
    float sum = 0.f;\
    for (int i = 0; i < N; i += 3)\
      sum += hypot_v##n(src[i], src[i+1], src[i+2]);\
    s[n] = sum;\
  // end define
  
  cout << endl << "3D versions\n";
  Timings<Clock>
    ("v0: ", [&]{ HT(0) })
    ("v1: ", [&]{ HT(1) })
    ("v2: ", [&]{ HT(2) })
    ("v3: ", [&]{ HT(3) })
    ("v0: ", [&]{ HT(0) })
    (cout);
    
  cout << '\n';
  
  Timings<Clock>
    ("v0: ", [&]{ HT(0) })
    ("v3: ", [&]{ HT(3) })
    ("v0: ", [&]{ HT(0) })
    ("v2: ", [&]{ HT(2) })
    ("v1: ", [&]{ HT(1) })
    (cout);
    
  cout << '\n';
  for (int i = 0; i < 4; ++i)
    cout << '\t' << s[i];
  cout << '\n';
  
  // Compare quality of hypot_v1, hypot_v2 against hypot_v3.
  float fixed[] { 1e-15f, 1e-6f, 1e-2f, 1.f, 1e+2f, 1e+6f, 1e+15f };
  for (auto y: fixed)
  {
    cout << "\nFor y = " << y << endl;
    
    float max_d1x = 0.f, max_d2x = 0.f;
    int i = 0;
    uint32_t max_ulp1 = 0, max_ulp2 = 0;
    
    const auto x_max = std::sqrt(FLT_MAX);
    for (float x = FLT_MIN; x != x_max; x = std::nextafter(x, FLT_MAX))
    {
      const float
        v1 = hypot_v1(x, y),
        v2 = hypot_v2(x, y),
        v3 = hypot_v3(x, y);
        
      const auto
        d1 = ulp_dist(v1, v3),
        d2 = ulp_dist(v2, v3);
        
      if (max_ulp1 < d1)
      {
        max_ulp1 = d1;
        max_d1x = x;
      }
      
      if (max_ulp2 < d2)
      {
        max_ulp2 = d2;
        max_d2x = x;
      }
      
      if (++i == 1 << 24)
      {
        cout.put('.').flush();
        i = 0;
      }
    }
    
    cout.precision(10);
    cout << "\nmax_d1 = " << max_ulp1 << " at x = " << max_d1x;
    cout << "\nmax_d2 = " << max_ulp2 << " at x = " << max_d2x;
  }
  
  cout << endl;
  cin.ignore();
  return 0;
}

// bisection_recursive.cpp
// Решаем уравнение f(x) = 0, x in [a, b].
// Простейший вариант реализации метода бисекции:
// рекурсия, заданная где-то (зафиксированная) функция f.

float f(float);
float solve(float a, float b, float fa, float fb)
{
  if (fa == 0.f) return a;
  if (fb == 0.f) return b;
  const auto m = 0.5f * (a + b), fm = f(m);
  if (m == a || m == b) return m;
  if (signbit(fa) != signbit(fm))
    return solve(a, m, fa, fm);
  if (signbit(fb) != signbit(fm))
    return solve(m, b, fm, fb);
  return m;
}

// Удобный "вход" в рекурсию.
float solve(float a, float b)
{ return solve(a, b, f(a), f(b)); }

// bisection_recursive_functor.cpp
// Решаем уравнение f(x) = 0, x in [a, b].
// Функцию передаём в стиле C++ (как функциональный объект).
// Простейший вариант реализации метода бисекции:
// рекурсия, заданная где-то (зафиксированная) функция f.

template <class F>
float solve(F &f, float a, float b, float fa, float fb)
{
  if (fa == 0.f) return a;
  if (fb == 0.f) return b;
  const auto m = 0.5f * (a + b), fm = f(m);
  if (m == a || m == b) return m;
  if (signbit(fa) != signbit(fm))
    return solve(f, a, m, fa, fm);
  if (signbit(fb) != signbit(fm))
    return solve(f, m, b, fm, fb);
  return m;
}

// Удобный "вход" в рекурсию.
template <class F>
float solve(F f, float a, float b)
{ return solve(f, a, b, f(a), f(b)); }

// Сложный вопрос: передавать f по значению или по ссылке?
// Проще избавиться от рекурсии и необходимости передавать f повторно в рекурсивный вызов.
// (См. примеры далее.)

// nsolve.cpp
#include <cassert>
#include <cstddef>
#include <cmath>
#include <climits>

typedef float(*Function)(float);
struct Local_solve_task
{
  Function function;             // solving function(x) = 0
  float function_tolerance;      // maximal tolerable error on function value (divergence from zero)
  float interval_tolerance;      // interval span lower bound
  float left_bound, right_bound; // interval containing a root
  float left_value, right_value; // values of function in points left_bound and right_bound
};

typedef bool(*Interval_check)(Local_solve_task&); // returns true to stop the search
typedef bool(*Local_solve)(Local_solve_task&, Interval_check);


/// Uniform grid search.
/// Precomputes task.left_value and task.right_value.
/// @return total number of report invocations
std::size_t global_solve
  (
    Local_solve_task &task,
    Interval_check interval_check,
    Local_solve lsolve,
    float step
  )
{
  assert(step != 0);
  const auto a = task.left_bound, b = task.right_bound;
  const std::size_t max_step = (long)((b - a) / step);
  assert(max_step != 0);
  assert(max_step < LONG_MAX);

  std::size_t invocations = 0;

  task.left_value = task.function(task.left_bound);
  for (std::size_t i = 1; i <= max_step; ++i)
  {
    const auto next_bound = i == max_step ? b : a + i * step;
    const auto next_value = task.function(next_bound);

    task.right_bound = next_bound;
    task.right_value = next_value;

    if (std::signbit(task.left_value) != std::signbit(task.right_value))
    {
      ++invocations;
      if (lsolve(task, interval_check))
        break;
    }

    task.left_bound = next_bound;
    task.left_value = next_value;
  }

  return invocations;
}


// Метод деления пополам (бисекция).
bool bisection(Local_solve_task &task, Interval_check interval_check)
{
  const auto f = task.function;
  const auto
    itol = task.interval_tolerance,
    ftol = task.function_tolerance;

  auto
    a = task.left_bound,
    b = task.right_bound,
    fa = task.left_value,
    fb = task.right_value;

  while (true)
  {
    if (std::abs(fa) <= ftol
     || std::abs(fb) <= ftol
     || std::abs(b - a) <= itol)
      break;

    const auto m = 0.5f * (a + b);
    if (m == a || m == b)
      break;

    const auto fm = f(m);
    if (std::signbit(fa) != std::signbit(fm))
    {
      b = m;
      fb = fm;
    }
    else
    {
      assert(std::signbit(fb) != std::signbit(fm));
      a = m;
      fa = fm;
    }
  }

  task.left_bound = a;
  task.right_bound = b;
  task.left_value = fa;
  task.right_value = fb;
  return interval_check(task);
}


// Вилочный метод секущих (метод ложного положения).
bool bracket_secant(Local_solve_task &task, Interval_check interval_check)
{
  const auto f = task.function;
  const auto
    itol = task.interval_tolerance,
    ftol = task.function_tolerance;

  auto
    a = task.left_bound,
    b = task.right_bound,
    fa = task.left_value,
    fb = task.right_value;

  while (true)
  {
    if (std::abs(fa) <= ftol
     || std::abs(fb) <= ftol
     || std::abs(b - a) <= itol)
      break;

    const auto m = (a*fb - b*fa) / (fb - fa); // the only difference with bisection
    if (m == a || m == b)
      break;

    const auto fm = f(m);
    if (std::signbit(fa) != std::signbit(fm))
    {
      b = m;
      fb = fm;
    }
    else
    {
      assert(std::signbit(fb) != std::signbit(fm));
      a = m;
      fa = fm;
    }
  }

  task.left_bound = a;
  task.right_bound = b;
  task.left_value = fa;
  task.right_value = fb;
  return interval_check(task);
}


// Модифицированный вилочный метод секущих.
bool bracket_secant_illinois(Local_solve_task &task, Interval_check interval_check)
{
  const auto f = task.function;
  const auto
    itol = task.interval_tolerance,
    ftol = task.function_tolerance;

  auto
    a = task.left_bound,
    b = task.right_bound,
    fa = task.left_value,
    fb = task.right_value;

  for (int to_left = 0, to_right = 0;;)
  {
    if (std::abs(fa) <= ftol
     || std::abs(fb) <= ftol
     || std::abs(b - a) <= itol)
      break;

    const auto m =
      to_left == 2?  (to_left = 0,   (a*fb - 0.5f*b*fa) / (fb - 0.5f*fa))
    : to_right == 2? (to_right = 0, ( 0.5f*a*fb - b*fa) / (0.5f*fb - fa))
    : (a*fb - b*fa) / (fb - fa);

    if (m == a || m == b)
      break;

    const auto fm = f(m);
    if (std::signbit(fa) != std::signbit(fm))
    {
      b = m;
      fb = fm;
      to_right = 0;
      ++to_left;
    }
    else
    {
      assert(std::signbit(fb) != std::signbit(fm));
      a = m;
      fa = fm;
      to_left = 0;
      ++to_right;
    }
  }

  task.left_bound = a;
  task.right_bound = b;
  task.left_value = fa;
  task.right_value = fb;
  return interval_check(task);
}


// Метод Риддера.
bool ridders(Local_solve_task &task, Interval_check interval_check)
{
  const auto f = task.function;
  const auto
    itol = task.interval_tolerance,
    ftol = task.function_tolerance;

  auto
    a = task.left_bound,
    b = task.right_bound,
    fa = task.left_value,
    fb = task.right_value;

  while (true)
  {
    if (std::abs(fa) <= ftol
     || std::abs(fb) <= ftol
     || std::abs(b - a) <= itol)
      break;

    const auto m = 0.5f * (a + b);
    if (m == a || m == b)
      break;

    const auto
      fm = f(m),
      add = (m - a) * fm / std::sqrt(fm * fm - fa * fb),
      p = m - (fa < fb? add: -add),
      fp = f(p);

    if (std::signbit(fp) != std::signbit(fm))
    {
      if (p < m)
      {
        a = p;
        b = m;
        fa = fp;
        fb = fm;
      }
      else
      {
        a = m;
        b = p;
        fa = fm;
        fb = fp;
      }
    }
    else if (std::signbit(fp) != std::signbit(fa))
    {
      b = p;
      fb = fp;
    }
    else
    {
      assert(std::signbit(fp) != std::signbit(fb));
      a = p;
      fa = fp;
    }
  }

  task.left_bound = a;
  task.right_bound = b;
  task.left_value = fa;
  task.right_value = fb;
  return interval_check(task);
}


// Метод Риддера + метод секущих.
bool ridders_secant(Local_solve_task &task, Interval_check interval_check)
{
  const auto f = task.function;
  const auto
    itol = task.interval_tolerance,
    ftol = task.function_tolerance;

  auto
    a = task.left_bound,
    b = task.right_bound,
    fa = task.left_value,
    fb = task.right_value;

  while (true)
  {
    if (std::abs(fa) <= ftol
     || std::abs(fb) <= ftol
     || std::abs(b - a) <= itol)
      break;

    const auto m = (a*fb - b*fa) / (fb - fa);
    if (m == a || m == b)
      break;

    const auto
      fm = f(m),
      add = (m - a) * fm / std::sqrt(fm * fm - fa * fb),
      p = m - (fa < fb? add: -add),
      fp = f(p);

    if (std::signbit(fp) != std::signbit(fm))
    {
      if (p < m)
      {
        a = p;
        b = m;
        fa = fp;
        fb = fm;
      }
      else
      {
        a = m;
        b = p;
        fa = fm;
        fb = fp;
      }
    }
    else if (std::signbit(fp) != std::signbit(fa))
    {
      b = p;
      fb = fp;
    }
    else
    {
      assert(std::signbit(fp) != std::signbit(fb));
      a = p;
      fa = fp;
    }
  }

  task.left_bound = a;
  task.right_bound = b;
  task.left_value = fa;
  task.right_value = fb;
  return interval_check(task);
}


// Метод Брента-Деккера.
// https://en.wikipedia.org/wiki/Brent%27s_method
bool brent_dekker(Local_solve_task &task, Interval_check interval_check)
{
  const auto f = task.function;
  const auto
    itol = task.interval_tolerance,
    ftol = task.function_tolerance;

  auto
    a = task.left_bound,
    b = task.right_bound,
    fa = task.left_value,
    fb = task.right_value;

  if (std::abs(fa) < std::abs(fb))
  {
    auto t = a; a = b; b = t;
    t = fa; fa = fb; fb = t;
  }

  auto c = a, fc = fa, d = c;
  for (bool flag = true;;)
  {
    if (std::abs(fa) <= ftol
     || std::abs(fb) <= ftol
     || std::abs(b - a) <= itol)
      break;

    if (std::abs(fc) <= ftol)
    {
      a = c;
      fa = fc;
      break;
    }

    auto m =
      (fa != fc && fb != fc)?
        (a*fb*fc / ((fa - fb)*(fa - fc)) + b*fa*fc / ((fb - fa)*(fb - fc)) + c*fa*fb / ((fc - fa)*(fc - fb)))
      : (a*fb - b*fa) / (fb - fa); // secant fallback

    if // possible bisection fallback
      (
        (a < b && (b <= m || 4 * m < (3 * a + b)))
     || (b < a && (m <= b || (3 * a + b) < 4 * m))
     ||  flag && 2 * std::abs(m - b) >= std::abs(b - c)
     || !flag && 2 * std::abs(m - b) >= std::abs(d - c)
     ||  flag && std::abs(b - c) < itol // delta?
     || !flag && std::abs(c - d) < itol // delta?
      )
    {
      m = 0.5f * (a + b);
      flag = true;
    }
    else
      flag = false;

    if (m == a || m == b)
      break;

    const auto fm = f(m);
    d = c;
    c = b;
    fc = fb;

    if (std::signbit(fa) != std::signbit(fm))
    {
      b = m;
      fb = fm;
    }
    else
    {
      assert(std::signbit(fb) != std::signbit(fm));
      a = m;
      fa = fm;
    }
  }

  if (a < b)
  {
    task.left_bound = a;
    task.right_bound = b;
    task.left_value = fa;
    task.right_value = fb;
  }
  else
  {
    task.left_bound = b;
    task.right_bound = a;
    task.left_value = fb;
    task.right_value = fa;
  }
  return interval_check(task);
}


////////////////////////////////////////////////////////////////////////////////
// Testing.
////////////////////////////////////////////////////////////////////////////////
#include <iostream>
using namespace std;

// The counter, which is used to directly measure efficiency of different methods.
size_t f_invocations;

float f(float x)
{
  ++f_invocations;
  return 0.5f * sin(x) + sin(1.5f * x + 2.f)
    - 1.5f * sin(5.f / 17.f * x + 1.f) - atan(0.5f * x + 5.f);
}

bool print_interval(Local_solve_task &task)
{
  cout << '[' << task.left_bound << ", " << task.right_bound << "], bound values: ";
  cout << task.left_value << ", " << task.right_value << endl;
  return false; // proceed
}

int main()
{
  Local_solve_task task{ f };
  Local_solve methods[]
  {
    bisection,
    bracket_secant,
    bracket_secant_illinois,
    ridders,
    ridders_secant,
    brent_dekker
  };

  for (auto method : methods)
  {
    cout << '\n';
    task.left_bound = -20.f;
    task.right_bound = 20.f;
    f_invocations = 0;
    const auto n = global_solve(task, print_interval, method, 0.1f);
    cout << "Roots found: " << n << ", f invocations: " << f_invocations << endl;
  }

  cin.ignore();
  return 0;
}

// nsolve_functor.cpp
// nsolve + шаблоны и функторы.
#include <cassert>
#include <cstddef>
#include <cmath>
#include <limits>
#include <utility>

template <class Argument>
inline Argument half(Argument x)
{
  return Argument(0.5) * x;
}

template <class Argument, class Value>
struct Local_solve_state
{
  Argument interval_tolerance;      // interval span lower bound
  Argument left_bound, right_bound; // interval containing a root
  Value function_tolerance;         // maximal tolerable error on function value (divergence from zero)
  Value left_value, right_value;    // values of function in points left_bound and right_bound
};


/// Uniform grid search.
/// Precomputes task.left_value and task.right_value.
/// @return total number of report invocations
template <
  class Function,
  class Argument,
  class Value,
  class Interval_check,
  class Local_solve>
std::size_t global_solve
  (
    Function function,
    Local_solve_state<Argument, Value> &task,
    Interval_check interval_check,
    Local_solve lsolve,
    Argument step
  )
{
  using std::signbit;

  assert(step != 0);
  const auto a = task.left_bound, b = task.right_bound;
  const std::size_t max_step = (long)((b - a) / step);
  assert(max_step != 0);
  assert(max_step < std::numeric_limits<long>::max());

  std::size_t invocations = 0;

  task.left_value = function(task.left_bound);
  for (std::size_t i = 1; i <= max_step; ++i)
  {
    const auto next_bound = i == max_step ? b : a + i * step;
    const auto next_value = function(next_bound);

    task.right_bound = next_bound;
    task.right_value = next_value;

    if (signbit(task.left_value) != signbit(task.right_value))
    {
      ++invocations;
      if (lsolve(function, task, interval_check))
        break;
    }

    task.left_bound = next_bound;
    task.left_value = next_value;
  }

  return invocations;
}


// Метод деления пополам (бисекция).
struct Bisection
{
  template <
    class Function,
    class Argument,
    class Value,
    class Interval_check>
  bool operator()
    (
      Function f,
      Local_solve_state<Argument, Value> &task,
      Interval_check interval_check
    ) const
  {
    using std::abs;
    using std::signbit;
    const auto itol = task.interval_tolerance;
    const auto ftol = task.function_tolerance;

    auto a = task.left_bound, b = task.right_bound;
    auto fa = task.left_value, fb = task.right_value;

    while (true)
    {
      if (abs(fa) <= ftol
       || abs(fb) <= ftol
       || abs(b - a) <= itol)
        break;

      const auto m = half(a + b);
      if (m == a || m == b)
        break;

      const auto fm = f(m);
      if (signbit(fa) != signbit(fm))
      {
        b = m;
        fb = fm;
      }
      else
      {
        assert(signbit(fb) != signbit(fm));
        a = m;
        fa = fm;
      }
    }

    task.left_bound = a;
    task.right_bound = b;
    task.left_value = fa;
    task.right_value = fb;
    return interval_check(task);
  }
} const bisection;


// Вилочный метод секущих (метод ложного положения).
struct Bracket_secant
{
  template <
    class Function,
    class Argument,
    class Value,
    class Interval_check>
  bool operator()
    (
      Function f,
      Local_solve_state<Argument, Value> &task,
      Interval_check interval_check
    ) const
  {
    using std::abs;
    using std::signbit;
    const auto itol = task.interval_tolerance;
    const auto ftol = task.function_tolerance;

    auto a = task.left_bound, b = task.right_bound;
    auto fa = task.left_value, fb = task.right_value;

    while (true)
    {
      if (abs(fa) <= ftol
       || abs(fb) <= ftol
       || abs(b - a) <= itol)
        break;

      const Argument m = (a*fb - b*fa) / (fb - fa); // the only difference with bisection
      if (m == a || m == b)
        break;

      const auto fm = f(m);
      if (signbit(fa) != signbit(fm))
      {
        b = m;
        fb = fm;
      }
      else
      {
        assert(signbit(fb) != signbit(fm));
        a = m;
        fa = fm;
      }
    }

    task.left_bound = a;
    task.right_bound = b;
    task.left_value = fa;
    task.right_value = fb;
    return interval_check(task);
  }
} const bracket_secant;


// Модифицированный вилочный метод секущих.
struct Bracket_secant_illinois
{
  template <
    class Function,
    class Argument,
    class Value,
    class Interval_check>
  bool operator()
    (
      Function f,
      Local_solve_state<Argument, Value> &task,
      Interval_check interval_check
    ) const
  {
    using std::abs;
    using std::signbit;
    const auto itol = task.interval_tolerance;
    const auto ftol = task.function_tolerance;

    auto a = task.left_bound, b = task.right_bound;
    auto fa = task.left_value, fb = task.right_value;

    for (int to_left = 0, to_right = 0;;)
    {
      if (abs(fa) <= ftol
       || abs(fb) <= ftol
       || abs(b - a) <= itol)
        break;

      const Argument m =
          to_left  == 2? (to_left = 0,  (a*fb - 0.5f*b*fa) / (fb - 0.5f*fa))
        : to_right == 2? (to_right = 0, (0.5f*a*fb - b*fa) / (0.5f*fb - fa))
        : (a*fb - b*fa) / (fb - fa);

      if (m == a || m == b)
        break;

      const auto fm = f(m);
      if (signbit(fa) != signbit(fm))
      {
        b = m;
        fb = fm;
        to_right = 0;
        ++to_left;
      }
      else
      {
        assert(signbit(fb) != signbit(fm));
        a = m;
        fa = fm;
        to_left = 0;
        ++to_right;
      }
    }

    task.left_bound = a;
    task.right_bound = b;
    task.left_value = fa;
    task.right_value = fb;
    return interval_check(task);
  }
} const bracket_secant_illinois;


// Метод Риддера.
struct Ridders
{
  template <
    class Function,
    class Argument,
    class Value,
    class Interval_check>
  bool operator()
    (
      Function f,
      Local_solve_state<Argument, Value> &task,
      Interval_check interval_check
    ) const
  {
    using std::abs;
    using std::signbit;
    using std::sqrt;
    const auto itol = task.interval_tolerance;
    const auto ftol = task.function_tolerance;

    auto a = task.left_bound, b = task.right_bound;
    auto fa = task.left_value, fb = task.right_value;

    while (true)
    {
      if (abs(fa) <= ftol
       || abs(fb) <= ftol
       || abs(b - a) <= itol)
        break;

      const auto m = half(a + b);
      if (m == a || m == b)
        break;

      const auto fm = f(m);
      const Argument
        add = (m - a) * fm / sqrt(fm * fm - fa * fb),
        p = m - (fa < fb? add: -add);

      const auto fp = f(p);

      if (signbit(fp) != signbit(fm))
      {
        if (p < m)
        {
          a = p;
          b = m;
          fa = fp;
          fb = fm;
        }
        else
        {
          a = m;
          b = p;
          fa = fm;
          fb = fp;
        }
      }
      else if (signbit(fp) != signbit(fa))
      {
        b = p;
        fb = fp;
      }
      else
      {
        assert(signbit(fp) != signbit(fb));
        a = p;
        fa = fp;
      }
    }

    task.left_bound = a;
    task.right_bound = b;
    task.left_value = fa;
    task.right_value = fb;
    return interval_check(task);
  }
} const ridders;



// Метод Брента-Деккера.
// https://en.wikipedia.org/wiki/Brent%27s_method
struct Brent_Dekker
{
  template <
    class Function,
    class Argument,
    class Value,
    class Interval_check>
  bool operator()
    (
      Function f,
      Local_solve_state<Argument, Value> &task,
      Interval_check interval_check
    ) const
  {
    using std::abs;
    using std::signbit;
    const auto itol = task.interval_tolerance;
    const auto ftol = task.function_tolerance;

    auto a = task.left_bound, b = task.right_bound;
    auto fa = task.left_value, fb = task.right_value;

    if (abs(fa) < abs(fb))
    {
      std::swap(a, b);
      std::swap(fa, fb);
    }

    auto c = a, d = a;
    auto fc = fa;

    for (bool flag = true;;)
    {
      if (abs(fa) <= ftol
       || abs(fb) <= ftol
       || abs(b - a) <= itol)
        break;

      if (abs(fc) <= ftol)
      {
        a = c;
        fa = fc;
        break;
      }

      Argument m =
        (fa != fc && fb != fc)?
          (a*fb*fc / ((fa - fb)*(fa - fc)) + b*fa*fc / ((fb - fa)*(fb - fc)) + c*fa*fb / ((fc - fa)*(fc - fb)))
        : (a*fb - b*fa) / (fb - fa); // secant fallback

      if // possible bisection fallback
        (
          (a < b && (b <= m || 4 * m < (3 * a + b)))
       || (b < a && (m <= b || (3 * a + b) < 4 * m))
       || ( flag && 2 * abs(m - b) >= abs(b - c))
       || (!flag && 2 * abs(m - b) >= abs(d - c))
       || ( flag && abs(b - c) < itol) // delta?
       || (!flag && abs(c - d) < itol) // delta?
        )
      {
        m = half(a + b);
        flag = true;
      }
      else
      {
        flag = false;
      }

      if (m == a || m == b)
        break;

      const auto fm = f(m);
      d = c;
      c = b;
      fc = fb;

      if (signbit(fa) != signbit(fm))
      {
        b = m;
        fb = fm;
      }
      else
      {
        assert(signbit(fb) != signbit(fm));
        a = m;
        fa = fm;
      }
    }

    if (a < b)
    {
      task.left_bound = a;
      task.right_bound = b;
      task.left_value = fa;
      task.right_value = fb;
    }
    else
    {
      task.left_bound = b;
      task.right_bound = a;
      task.left_value = fb;
      task.right_value = fa;
    }
    return interval_check(task);
  }
} const brent_dekker;


////////////////////////////////////////////////////////////////////////////////
// Testing.
////////////////////////////////////////////////////////////////////////////////
#include <iostream>
using namespace std;

// The counter, which is used to directly measure efficiency of different methods.
size_t f_invocations;

float f(float x)
{
  ++f_invocations;
  return 0.5f * sin(x) + sin(1.5f * x + 2.f)
    - 1.5f * sin(5.f / 17.f * x + 1.f) - atan(0.5f * x + 5.f);
}

template <class Argument, class Value>
bool print_interval(Local_solve_state<Argument, Value> &task)
{
  cout << '[' << task.left_bound << ", " << task.right_bound << "], bound values: ";
  cout << task.left_value << ", " << task.right_value << endl;
  return false; // proceed
}

template <class Method>
void test(Method method)
{
  cout << '\n';

  Local_solve_state<float, float> task;
  task.left_bound = -20.f;
  task.right_bound = 20.f;
  f_invocations = 0;

  const auto n = global_solve(f, task, print_interval<float, float>, method, 0.1f);
  cout << "Roots found: " << n << ", f invocations: " << f_invocations << endl;
}


int main()
{
#define TEST(m) \
  cout << '\n' << #m; \
  test(m)

  TEST(bisection);
  TEST(bracket_secant);
  TEST(bracket_secant_illinois);
  TEST(ridders);
  TEST(brent_dekker);

  cin.ignore();
  return 0;
}

// nmin.cpp
#include <cassert>
#include <cstddef>
#include <cmath>
#include <climits>

typedef float(*Function)(float);
struct Local_min_task
{
  Function function;             // solving function(x) = 0
  float interval_tolerance;      // interval span lower bound
  float left_bound, right_bound, min_candidate; // interval containing a root + some "mid" point, which is the last minimum candidate
  float left_value, right_value, min_value;     // values of function in points left_bound, right_bound and min_candidate
};

typedef bool(*Interval_check)(Local_min_task&); // returns true to stop the search
typedef bool(*Local_minimize)(Local_min_task&, Interval_check);


/// Uniform grid search.
/// Precomputes task.left_value and task.right_value.
/// @return total number of report invocations
std::size_t global_minimize
  (
    Local_min_task &task,
    Interval_check interval_check,
    Local_minimize lmin,
    float step
  )
{
  assert(step != 0);
  const auto a = task.left_bound, b = task.right_bound;
  const std::size_t max_step = (long)((b - a) / step);
  assert(max_step != 0);
  assert(max_step < LONG_MAX);

  std::size_t invocations = 0;
  task.min_value = task.left_value = task.function(task.left_bound);
  task.min_candidate = task.left_bound;

  Local_min_task best_result = task;
  for (std::size_t i = 1; i <= max_step; ++i)
  {
    const auto next_bound = i == max_step ? b : a + i * step;
    const auto next_value = task.function(next_bound);

    task.right_bound = next_bound;
    task.right_value = next_value;

    ++invocations;
    const bool finish = lmin(task, interval_check);
    if (task.min_value < best_result.min_value)
      best_result = task;
    if (finish)
      break;

    task.left_bound = next_bound;
    task.left_value = next_value;
  }

  task = best_result;
  return invocations;
}


// Троичный поиск.
bool ternary_minimize(Local_min_task &task, Interval_check interval_check)
{
  const auto f = task.function;
  const auto itol = task.interval_tolerance;
  auto
    a = task.left_bound,
    b = task.right_bound,
    fa = task.left_value,
    fb = task.right_value;

  while (true)
  {
    if (b - a <= itol)
    {
      const auto m = 0.5f * (a + b);
      const auto fm = f(m);
      if (fa < fm)
      {
        task.min_candidate = a;
        task.min_value = fa;
      }
      else
      {
        task.min_candidate = m;
        task.min_value = fm;
      }
      break;
    }

    const auto
      third = (b - a) * (1.f/3.f),
      m1 = a + third,
      m2 = b - third;

    if (m1 == a || m2 == b || m2 <= m1)
    {
      task.min_candidate = a;
      task.min_value = fa;
      break;
    }

    const auto
      fm1 = f(m1),
      fm2 = f(m2);

    if (fm1 < fm2)
    {
      b = m2;
      fb = fm2;
    }
    else if (fm2 < fm1)
    {
      a = m1;
      fa = fm1;
    }
    else
    {
      assert(fm1 == fm2);
      a = m1;
      b = m2;
      fa = fm1;
      fb = fm2;
    }
  }

  if (fb < task.min_value)
  {
    task.min_candidate = b;
    task.min_value = fb;
  }

  task.left_bound = a;
  task.right_bound = b;
  task.left_value = fa;
  task.right_value = fb;
  return interval_check(task);
}


// Метод золотого сечения.
bool golden_ratio_minimize(Local_min_task &task, Interval_check interval_check)
{
  // Golden ratio factor.
  static const auto phi = static_cast<float>(0.5 * (std::sqrt(5.0) - 1.));
  const auto f = task.function;
  const auto itol = task.interval_tolerance;
  auto
    a = task.left_bound,
    b = task.right_bound,
    c = b - (b - a) * phi,
    d = a + (b - a) * phi,
    fc = f(c),
    fd = f(d);

  while (true)
  {
    if (d - c <= itol || c <= a || b <= d)
    {
      const auto m = 0.5f * (a + b);
      const auto fm = f(m);
      if (fc < fm)
      {
        task.min_candidate = c;
        task.min_value = fc;
      }
      else
      {
        task.min_candidate = m;
        task.min_value = fm;
      }
      break;
    }

    if (fc < fd)
    {
      b = d;
      d = c;
      fd = fc;
      c = b - (b - a) * phi;
      fc = f(c);
    }
    else if (fd < fc)
    {
      a = c;
      c = d;
      fc = fd;
      d = a + (b - a) * phi;
      fd = f(d);
    }
    else
    {
      assert(fc == fd);
      a = c;
      b = d;
      c = b - (b - a) * phi;
      d = a + (b - a) * phi;
      fc = f(c);
      fd = f(d);
    }
  }

  if (fd < task.min_value)
  {
    task.min_candidate = d;
    task.min_value = fd;
  }

  task.left_bound = c;
  task.right_bound = d;
  task.left_value = fc;
  task.right_value = fd;
  return interval_check(task);
}


// Brent's method applied to minimization.
// Borrowed from Boost (see boost.org, brent_find_minima function): boost/math/tools/minima.hpp.
bool brent_minimize(Local_min_task &task, Interval_check interval_check)
{
  const auto f = task.function;
  const auto tolerance = std::fmax(task.interval_tolerance, std::ldexp(1.f, -23));
  auto
    min = task.left_bound,
    max = task.right_bound;

  // Original comments below.
  static const auto golden = 0.3819660f;  // golden ratio, don't need too much precision here!
  float
    x,  // minima so far
    w,  // second best point
    v,  // previous value of w
    u,  // most recent evaluation point
    delta,  // The distance moved in the last step
    delta2, // The distance moved in the step before last
    fu, fv, fw, fx,  // function evaluations at u, v, w, x
    mid, // midpoint of min and max
    fract1, fract2;  // minimal relative movement in x

  x = w = v = max;
  fw = fv = fx = task.right_value;
  delta2 = delta = 0;

  while (true)
  {
    // get midpoint
    mid = (min + max) / 2;
    // work out if we're done already:
    fract1 = tolerance * std::abs(x) + tolerance / 4;
    fract2 = 2 * fract1;
    if (std::abs(x - mid) <= (fract2 - (max - min) / 2))
      break;

    if (std::abs(delta2) > fract1)
    {
      // try and construct a parabolic fit:
      auto
        r = (x - w) * (fx - fv),
        q = (x - v) * (fx - fw),
        p = (x - v) * q - (x - w) * r;

      q = 2 * (q - r);
      if (q > 0)
        p = -p;
      q = std::abs(q);

      auto td = delta2;
      delta2 = delta;

      // determine whether a parabolic step is acceptible or not:
      if ((std::abs(p) >= std::abs(q * td / 2)) || (p <= q * (min - x)) || (p >= q * (max - x)))
      {
        // nope, try golden section instead
        delta2 = (x >= mid) ? min - x : max - x;
        delta = golden * delta2;
      }
      else
      {
        // whew, parabolic fit:
        delta = p / q;
        u = x + delta;
        if (((u - min) < fract2) || ((max- u) < fract2))
          delta = (mid - x) < 0 ? -std::abs(fract1) : std::abs(fract1);
      }
    }
    else
    {
      // golden section:
      delta2 = (x >= mid) ? min - x : max - x;
      delta = golden * delta2;
    }

    // update current position:
    u = (std::abs(delta) >= fract1) ? (x + delta) : (delta > 0 ? (x + std::abs(fract1)) : (x - std::abs(fract1)));
    fu = f(u);

    if (fu <= fx)
    {
      // good new point is an improvement!
      // update brackets:
      if(u >= x)
        min = x;
      else
        max = x;
      // update control points:
      v = w;
      w = x;
      x = u;
      fv = fw;
      fw = fx;
      fx = fu;
    }
    else
    {
      // Oh dear, point u is worse than what we have already,
      // even so it *must* be better than one of our endpoints:
      if (u < x)
        min = u;
      else
        max = u;
      if ((fu <= fw) || (w == x))
      {
        // however it is at least second best:
        v = w;
        w = u;
        fv = fw;
        fw = fu;
      }
      else if ((fu <= fv) || (v == x) || (v == w))
      {
        // third best:
        v = u;
        fv = fu;
      }
    }
  }

  task.min_candidate = x;
  task.min_value = fx;
  return interval_check(task);
}


////////////////////////////////////////////////////////////////////////////////
// Testing.
////////////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <iterator>
#include <algorithm>
using namespace std;

// The counter, which is used to directly measure efficiency of different methods.
size_t f_invocations;
// "Global" minimum position and value.
float global_min, global_min_position;

float f(float x)
{
  ++f_invocations;
  return 0.5f * sin(x) + sin(1.5f * x + 2.f)
    - 1.5f * sin(5.f / 17.f * x + 1.f) - atan(0.5f * x + 5.f);
}

bool update_global_min(Local_min_task &task)
{
  if (task.min_value < global_min)
  {
    global_min = task.min_value;
    global_min_position = task.min_candidate;
  }
  return false; // proceed
}

int main()
{
  cout.precision(9);
  Local_min_task task{ f };

  static const int METHODS = 3;
  Local_minimize methods[METHODS]
  {
    ternary_minimize,
    golden_ratio_minimize,
    brent_minimize
  };

  const char *name[METHODS]
  {
    "Ternary search",
    "Golden ratio search",
    "Brent"
  };

  float gmin[METHODS], gminpos[METHODS];
  for (int i = 0; i < METHODS; ++i)
  {
    cout << name[i] << '\n';
    task.left_bound = -20.f;
    task.right_bound = 20.f;
    f_invocations = 0;
    global_min = 100.f;
    global_min_position = -100.f;
    const auto n = global_minimize(task, update_global_min, methods[i], 0.05f);
    cout << "Steps: " << n << ", f invocations: " << f_invocations << '\n';
    cout << "Global min = " << (gmin[i] = global_min) << ", at " << (gminpos[i] = global_min_position) << '\n' << endl;
  }

  auto best_i = min_element(begin(gmin), end(gmin)) - gmin;
  cout << "With min = " << gmin[best_i] << " winners are: ";
  for (int i = 0; i < METHODS; ++i)
    if (gmin[i] == gmin[best_i])
      cout << name[i] << ", ";
  cout << "\b\b." << endl;
  cin.ignore();
  return 0;
}

// nmin_functor.cpp
#include <cassert>
#include <cstddef>
#include <cmath>
#include <limits>

template <class Argument>
inline Argument half(Argument x)
{
  return Argument(0.5) * x;
}

template <class Argument>
inline Argument third(Argument x)
{
  return Argument(1./3.) * x;
}

template <class Argument, class Value>
struct Local_min_state
{
  // interval span lower bound
  Argument interval_tolerance;
  // interval containing a root + some "mid" point, which is the last minimum candidate
  Argument left_bound, right_bound, min_candidate;
  // values of function in points left_bound, right_bound and min_candidate
  Value left_value, right_value, min_value;
};


/// Uniform grid search.
/// Precomputes task.left_value and task.right_value.
/// @return total number of report invocations
template <
  class Function,
  class Argument,
  class Value,
  class Interval_check,
  class Local_minimize>
std::size_t global_minimize
  (
    Function f,
    Local_min_state<Argument, Value> &task,
    Interval_check interval_check,
    Local_minimize lmin,
    Argument step
  )
{
  step /= 2;
  assert(step != 0);

  const auto a = task.left_bound, b = task.right_bound;
  const std::size_t max_step = (long)((b - a) / step);

  assert(max_step != 0);
  assert(max_step < std::numeric_limits<long>::max());

  std::size_t invocations = 0;

  Argument x = a, y = a + step, z;
  Value fx = f(x), fy = f(y), fz;
  auto descending = fy < fx;
  task.left_bound = x;
  task.right_bound = y;
  task.left_value = fx;
  task.right_value = fy;
  task.min_candidate = descending? y: x;
  task.min_value = descending? fy: fx;

  auto best_result = task;
  for (std::size_t i = 2; i <= max_step; ++i)
  {
    z = i == max_step? b: a + i * step;
    fz = f(z);

    if (descending && fy < fz)
    {
      // A minimum is bracketed.
      task.left_bound = x;
      task.right_bound = z;
      task.left_value = fx;
      task.right_value = fz;

      ++invocations;
      const bool finish = lmin(f, task, interval_check);
      if (task.min_value < best_result.min_value)
        best_result = task;

      if (finish)
        break;

      descending = false;
    }
    else
    {
      descending = fz < fy;
    }

    x = y;
    y = z;
    fx = fy;
    fy = fz;
  }

  task = best_result;
  return invocations;
}


// Троичный поиск.
template <
  class Function,
  class Argument,
  class Value,
  class Interval_check>
bool ternary_minimize
  (
    Function f,
    Local_min_state<Argument, Value> &task,
    Interval_check interval_check
  )
{
  const auto itol = task.interval_tolerance;
  auto a = task.left_bound, b = task.right_bound;
  auto fa = task.left_value, fb = task.right_value;

  while (true)
  {
    if (b - a <= itol)
    {
      const auto m = half(a + b);
      const auto fm = f(m);
      if (fa < fm)
      {
        task.min_candidate = a;
        task.min_value = fa;
      }
      else
      {
        task.min_candidate = m;
        task.min_value = fm;
      }
      break;
    }

    const auto
      th = third(b - a),
      m1 = a + th,
      m2 = b - th;

    if (m1 == a || m2 == b || m2 <= m1)
    {
      task.min_candidate = a;
      task.min_value = fa;
      break;
    }

    const auto fm1 = f(m1), fm2 = f(m2);
    if (fm1 < fm2)
    {
      b = m2;
      fb = fm2;
    }
    else if (fm2 < fm1)
    {
      a = m1;
      fa = fm1;
    }
    else
    {
      assert(fm1 == fm2);
      a = m1;
      b = m2;
      fa = fm1;
      fb = fm2;
    }
  }

  if (fb < task.min_value)
  {
    task.min_candidate = b;
    task.min_value = fb;
  }

  task.left_bound = a;
  task.right_bound = b;
  task.left_value = fa;
  task.right_value = fb;
  return interval_check(task);
}


// Метод золотого сечения.
template <
  class Function,
  class Argument,
  class Value,
  class Interval_check>
bool golden_ratio_minimize
  (
    Function f,
    Local_min_state<Argument, Value> &task,
    Interval_check interval_check
  )
{
  using std::sqrt;
  // Golden ratio factor: (sqrt(5) - 1)/2.
  static const auto
    golden = static_cast<Argument>(half(sqrt(Argument(5)) - Argument(1)));
  const auto
    itol = task.interval_tolerance;
  auto
    a = task.left_bound,
    b = task.right_bound,
    c = b - (b - a) * golden,
    d = a + (b - a) * golden;
  auto
    fc = f(c),
    fd = f(d);

  while (true)
  {
    if (d - c <= itol || c <= a || b <= d)
    {
      const auto m = half(a + b);
      const auto fm = f(m);
      if (fc < fm)
      {
        task.min_candidate = c;
        task.min_value = fc;
      }
      else
      {
        task.min_candidate = m;
        task.min_value = fm;
      }
      break;
    }

    if (fc < fd)
    {
      b = d;
      d = c;
      fd = fc;
      c = b - (b - a) * golden;
      fc = f(c);
    }
    else if (fd < fc)
    {
      a = c;
      c = d;
      fc = fd;
      d = a + (b - a) * golden;
      fd = f(d);
    }
    else
    {
      assert(fc == fd);
      a = c;
      b = d;
      c = b - (b - a) * golden;
      d = a + (b - a) * golden;
      fc = f(c);
      fd = f(d);
    }
  }

  if (fd < task.min_value)
  {
    task.min_candidate = d;
    task.min_value = fd;
  }

  task.left_bound = c;
  task.right_bound = d;
  task.left_value = fc;
  task.right_value = fd;
  return interval_check(task);
}


// Brent's method applied to minimization.
// Borrowed from Boost (see boost.org, brent_find_minima function): boost/math/tools/minima.hpp.
template <
  class Function,
  class Argument,
  class Value,
  class Interval_check>
bool brent_minimize
  (
    Function f,
    Local_min_state<Argument, Value> &task,
    Interval_check interval_check
  )
{
  using std::fmax;
  using std::ldexp;
  using std::abs;

  const auto
    tolerance = static_cast<Argument>(
      fmax(task.interval_tolerance, static_cast<Argument>(ldexp(1., -23))));
  auto
    min = task.left_bound,
    max = task.right_bound;

  // Golden ratio factor: (sqrt(5) - 1)/2.
  static const auto
    golden = static_cast<Argument>(half(sqrt(Argument(5)) - Argument(1)));

  // Original Boost authors comments below.
  Argument
    x,  // minima so far
    w,  // second best point
    v,  // previous value of w
    u,  // most recent evaluation point
    delta,  // The distance moved in the last step
    delta2, // The distance moved in the step before last
    fract1, fract2,  // minimal relative movement in x
    mid; // midpoint of min and max

  Value
    fu, fv, fw, fx;  // function evaluations at u, v, w, x

  x = w = v = max;
  fw = fv = fx = task.right_value;
  delta2 = delta = 0;

  while (true)
  {
    // get midpoint
    mid = half(min + max);
    // work out if we're done already:
    fract1 = tolerance * abs(x) + tolerance * Argument(0.25);
    fract2 = 2 * fract1;
    if (abs(x - mid) <= (fract2 - half(max - min)))
      break;

    if (abs(delta2) > fract1)
    {
      // try and construct a parabolic fit:
      auto
        r = (x - w) * (fx - fv),
        q = (x - v) * (fx - fw),
        p = (x - v) * q - (x - w) * r;

      q = 2 * (q - r);
      if (q > 0)
        p = -p;
      q = abs(q);

      auto td = delta2;
      delta2 = delta;

      // determine whether a parabolic step is acceptible or not:
      if ((abs(p) >= abs(q * half(td))) || (p <= q * (min - x)) || (p >= q * (max - x)))
      {
        // nope, try golden section instead
        delta2 = (x >= mid) ? min - x : max - x;
        delta = golden * delta2;
      }
      else
      {
        // whew, parabolic fit:
        delta = p / q;
        u = x + delta;
        if (((u - min) < fract2) || ((max- u) < fract2))
          delta = (mid - x) < 0 ? -abs(fract1) : abs(fract1);
      }
    }
    else
    {
      // golden section:
      delta2 = (x >= mid) ? min - x : max - x;
      delta = golden * delta2;
    }

    // update current position:
    u = (abs(delta) >= fract1) ?
        (x + delta)
      : (delta > 0 ? (x + abs(fract1)) : (x - abs(fract1)));

    fu = f(u);

    if (fu <= fx)
    {
      // good new point is an improvement!
      // update brackets:
      if(u >= x)
        min = x;
      else
        max = x;
      // update control points:
      v = w;
      w = x;
      x = u;
      fv = fw;
      fw = fx;
      fx = fu;
    }
    else
    {
      // Oh dear, point u is worse than what we have already,
      // even so it *must* be better than one of our endpoints:
      if (u < x)
        min = u;
      else
        max = u;
      if ((fu <= fw) || (w == x))
      {
        // however it is at least second best:
        v = w;
        w = u;
        fv = fw;
        fw = fu;
      }
      else if ((fu <= fv) || (v == x) || (v == w))
      {
        // third best:
        v = u;
        fv = fu;
      }
    }
  }

  task.min_candidate = x;
  task.min_value = fx;
  return interval_check(task);
}


////////////////////////////////////////////////////////////////////////////////
// Testing.
////////////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <iterator>
#include <algorithm>
using namespace std;

// The counter, which is used to directly measure efficiency of different methods.
size_t f_invocations;
// "Global" minimum position and value.
float global_min, global_min_position;

float f(float x)
{
  ++f_invocations;
  return 0.5f * sin(x) + sin(1.5f * x + 2.f)
    - 1.5f * sin(5.f / 17.f * x + 1.f) - atan(0.5f * x + 5.f);
}

bool update_global_min(Local_min_state<float, float> &task)
{
  if (task.min_value < global_min)
  {
    global_min = task.min_value;
    global_min_position = task.min_candidate;
  }
  return false; // proceed
}

int main()
{
  cout.precision(9);
  Local_min_state<float, float> task
  { // Interval tolerance:
    std::sqrt(std::numeric_limits<float>::epsilon())
  };

  static const int METHODS = 3;
  using F = float(*)(float);
  using IC = bool(*)(Local_min_state<float, float>&);

  bool (*methods[METHODS])(F, Local_min_state<float, float>&, IC)
  {
    ternary_minimize<F, float, float, IC>,
    golden_ratio_minimize<F, float, float, IC>,
    brent_minimize<F, float, float, IC>
  };

  const char *name[METHODS]
  {
    "Ternary search",
    "Golden ratio search",
    "Brent"
  };

  float gmin[METHODS], gminpos[METHODS];
  for (int i = 0; i < METHODS; ++i)
  {
    cout << name[i] << '\n';
    task.left_bound = -20.f;
    task.right_bound = 20.f;
    f_invocations = 0;
    global_min = 100.f;
    global_min_position = -100.f;
    const auto n = global_minimize(f, task, update_global_min, methods[i], 0.05f);
    cout << "Minima brackets found: " << n << ", f invocations: " << f_invocations << '\n';
    cout << "Global min = " << (gmin[i] = global_min) << ", at " << (gminpos[i] = global_min_position) << '\n' << endl;
  }

  auto best_i = min_element(begin(gmin), end(gmin)) - gmin;
  cout << "With min = " << gmin[best_i] << " winners are: ";
  for (int i = 0; i < METHODS; ++i)
    if (gmin[i] == gmin[best_i])
      cout << name[i] << ", ";
  cout << "\b\b." << endl;
  cin.ignore();
  return 0;
}

// bsearch.cpp
#include <iostream>
#include <iterator>
using namespace std;

// Двоичный поиск в массиве.
const int* bsearch(int item, const int *begin, const int *end)
{
  while (1 < end - begin)
  {
    const auto mid = begin + (end - begin) / 2;
    if (item < *mid)
      end = mid;
    else
      begin = mid;
  }
  return begin;
}

int main()
{
  int a[] { 0, 1, 2, 2, 2, 3, 4, 5, 6, 8, 40 };

  for (int i: a)
  {
    const auto pos = bsearch(i, begin(a), end(a));
    cout << "a[" << (pos - a) << "] == " << *pos << '\n';
  }

  auto pos = bsearch(7, begin(a), end(a));
  cout << "a[" << (pos - a) << "] == " << *pos << '\n';

  pos = bsearch(9, begin(a), end(a));
  cout << "a[" << (pos - a) << "] == " << *pos << '\n';

  pos = bsearch(-1, begin(a), end(a));
  cout << "a[" << (pos - a) << "] == " << *pos << '\n';

  pos = bsearch(100, begin(a), end(a));
  cout << "a[" << (pos - a) << "] == " << *pos << '\n';

  cin.ignore();
  return 0;
}

// eqrange.cpp
#include <iostream>
#include <iterator>
using namespace std;

// Двоичный поиск нижней границы.
const int* lbound(const int *begin, const int *end, int item)
{
  while (1 < end - begin)
  {
    const auto mid = begin + (end - begin) / 2;
    if (*mid < item)
      begin = mid + 1;
    else
      end = mid;
  }

  return begin == end || *begin < item? end: begin;
}

// Двоичный поиск верхней границы.
const int* ubound(const int *begin, const int *end, int item)
{
  while (1 < end - begin)
  {
    const auto mid = begin + (end - begin) / 2;
    if (item < *mid)
      end = mid;
    else
      begin = mid;
  }

  return begin == end || !(item < *begin)? end: begin;
}


struct Const_range
{
  const int *b, *e;
  auto begin() const { return b; }
  auto end() const { return e; }
};

// Поиск диапазона эквивалентных (равных) элементов.
Const_range eqrange(const int *begin, const int *end, int item)
{
  const auto lb = lbound(begin, end, item);
  return { lb, ubound(lb, end, item) };
}


int main()
{
  int a[] { 0, 1, 2, 2, 2, 3, 4, 5, 6, 8, 40 };

  for (int x : eqrange(begin(a), end(a), 1))
    cout << x << ' ';

  cout << endl;
  for (int x : eqrange(begin(a), end(a), 2))
    cout << x << ' ';

  cout << endl;
  for (int x : eqrange(begin(a), end(a), 7))
    cout << x << ' ';

  cin.ignore();
  return 0;
}

// vll_stack.cpp
// void (intrusive) linked list
#include <cassert>
#include <cstdlib>
using namespace std;

// Type-agnostic (generic) code.
// Pure C.

void* list_bind(void *head, void *tail)
{
  assert(head);
  *(void**)head = tail;
  return head;
}

void* list_tail(void *head)
{
  assert(head);
  return *(void**)head;
}

void* head_value(void *head)
{
  return ((void**)head + 1);
}

void* alloc_head(size_t bytes)
{
  return calloc(1, bytes + sizeof(void*));
}

void free_head(void **head)
{
  void *tail = list_tail(*head);
  free(*head);
  *head = tail;
}

void free_list(void *head)
{
  while (head)
    free_head(&head);
}


// Use this to organize a stack.
// Uses C++ references, nullptr and IO.
// Still quite dirty at the edge (or behind?) of UB.
#include <iostream>

struct Int_stack_node
{
  void *link;
  int value;
};

Int_stack_node*& push(Int_stack_node *&stack, int value)
{
  auto new_head = alloc_head(sizeof(int));
  *(int*)head_value(new_head) = value;
  return stack = (Int_stack_node*)list_bind(new_head, stack);
}

int pop(Int_stack_node *&stack)
{
  assert(stack);
  const auto result = *(int*)head_value(stack);
  free_head((void**)&stack);
  return result;
}

void destroy(Int_stack_node *&stack)
{
  free_list(stack);
  stack = nullptr;
}

int main()
{
  Int_stack_node *stack = nullptr;
  for (int i; cin >> i;)
    push(stack, i);

  cin.clear();
  while (stack)
    cout << '\n' << pop(stack);

  cout << endl;
  return 0;
}

// vbst.cpp
// void (intrusive) binary (search) tree
#include <cassert>
#include <cstdlib>
using namespace std;

// Type-agnostic (generic) code.
// Pure C.

void* tree_bind(void *parent, void *left, void *right)
{
  assert(parent);
  void **child = (void**)parent;

  child[0] = left;
  child[1] = right;
  return parent;
}

void* tree_left(void *parent)
{
  assert(parent);
  return ((void**)parent)[0];
}

void* tree_right(void *parent)
{
  assert(parent);
  return ((void**)parent)[1];
}

void* parent_value(void *parent)
{
  return ((void**)parent + 2);
}

void* alloc_tree_node(size_t bytes)
{
  return calloc(1, bytes + 2 * sizeof(void*));
}

void free_tree_node(void *parent)
{
  free(parent);
}

// May crash due to stack overflow!
void free_tree(void *parent)
{
  if (parent)
  {
    free_tree(tree_left(parent));
    free_tree(tree_right(parent));
    free_tree_node(parent);
  }
}

// In-order tree traversal.
// May crash due to stack overflow!
void tree_inorder(void *parent, void (*visit)(void*))
{
  if (parent)
  {
    tree_inorder(tree_left(parent), visit);
    visit(parent_value(parent));
    tree_inorder(tree_right(parent), visit);
  }
}


// Use this to organize a tree.
// Uses C++ references, nullptr and IO.
// Still quite dirty at the edge (or behind?) of UB.
#include <iostream>

struct Bst_node
{
  void *left, *right;
  int value;
};

void destroy(Bst_node *&tree)
{
  free_tree(tree);
  tree = nullptr;
}

Bst_node*& insert(Bst_node *&tree, int value)
{
  auto new_node = alloc_tree_node(sizeof(int));
  *(int*)parent_value(new_node) = value;

  if (!tree)
  {
    tree = (Bst_node*)new_node;
  }
  else // Binary search.
  {
    for (Bst_node *parent = tree;;)
    {
      // A reference to a pointer (this would be a "double-star" in C).
      void *&next = value < parent->value?
        parent->left: parent->right;

      if (!next)
      {
        next = new_node;
        break;
      }

      parent = (Bst_node*)next;
    }
  }

  return tree;
}


void print_int(void *data)
{
  cout << *(int*)data << ' ';
}

int main()
{
  Bst_node *tree = nullptr;
  for (int i; cin >> i;)
    insert(tree, i);

  cin.clear();
  tree_inorder(tree, print_int);

  cout << endl;
  return 0;
}

// trie_inorder.cpp
// Simplistic trie.
#include <cstddef>
#include <cassert>
#include <iostream>
#include <iterator>
#include <string>
#include <fstream>  // вывод в файл: ofstream
using namespace std;

/// Условное значение "позиция не найдена".
const size_t NOT_FOUND = -1;

/// Найти позицию нулевого символа в массиве символов.
size_t find_null(const char *word, size_t word_len)
{
  for (size_t i = 0; i < word_len; ++i)
    if (word[i] == '\0')
      return i;
  return NOT_FOUND;
}

/// Найти позицию заданного символ в массиве символов.
size_t find_char(const char *word, size_t word_len, char ch)
{
  for (size_t i = 0; i < word_len && word[i] != '\0'; ++i)
    if (word[i] == ch)
      return i;
  return NOT_FOUND;
}

/// Максимально возможное число потомков одного узла префиксного дерева.
const size_t MAX_TRIE_CHILDREN = 8;
/// Узел префиксного дерева.
struct Trie_node
{
  Trie_node *children[MAX_TRIE_CHILDREN];
  size_t label[MAX_TRIE_CHILDREN];
  char prefix[MAX_TRIE_CHILDREN];

  /// "Конструктор по умолчанию": вызывается при создании узла.
  Trie_node() { prefix[0] = '\0'; }

  /// Поиск свободной позиции в массиве потомков.
  size_t find_empty_slot()
  {
    return find_null(prefix, MAX_TRIE_CHILDREN);
  }
  /// Поиск номера потомка, отвечающего очередному символу.
  size_t find_char(char ch) const
  {
    return ::find_char(prefix, MAX_TRIE_CHILDREN, ch);
  }

  /// Рекурсивная операция вставки распознаваемой строки с заданной меткой lbl.
  void add_word(const char *word, size_t lbl)
  {
    if (*word == '\0')
      return;

    auto index = find_char(*word);
    if (index == NOT_FOUND)
    {
      index = find_empty_slot();
      if (index == NOT_FOUND)
      {
        assert(!"Trie: impossible to insert the next child: the node is full.");
        return;
      }

      prefix[index] = *word;
      if (index + 1 < MAX_TRIE_CHILDREN)
        prefix[index + 1] = '\0'; // The next empty slot.

      children[index] = new Trie_node;
      label[index] = NOT_FOUND; // Unassigned label at first.
    }

    if (word[1] == '\0')
      label[index] = lbl;
    else
      children[index]->add_word(word + 1, lbl);
    // Аналогично find_word: рекурсия -- хвостовой вызов,
    // как здесь заменить рекурсию на цикл?
  }

  /// Извлечение метки, отвечающей заданному слову
  /// (с которой оно было ранее вставлено add_word).
  size_t find_word(const char *word) const
  {
    if (*word == '\0')
      return NOT_FOUND;

    auto index = find_char(*word);
    if (index == NOT_FOUND)
      return NOT_FOUND;

    // Это был последний ненулевой символ слова
    // (== последний символ си-строки)?
    if (word[1] == '\0')
      return label[index];

    // Нет -- продолжим с "хвостом" слова (рекурсия -- хвостовой вызов).
    // (Как заменить этот вызов на цикл?)
    return children[index]->find_word(word + 1);
  }

  /// Обход поддерева слева-направо.
  void inorder(
    void (*open)(),
    void (*visit)(char, size_t),
    void (*close)())
  {
    if (prefix[0] == '\0')
      return;

    open();
    for (size_t i = 0; i < MAX_TRIE_CHILDREN && prefix[i] != '\0'; ++i)
    {
      visit(prefix[i], label[i]);
      children[i]->inorder(open, visit, close);
    }
    close();
  }
};

/// Префиксное дерево.
struct Trie
{
  Trie_node *root;

  void add_word(const char *word, size_t label)
  {
    assert(word);
    if (!root)
      root = new Trie_node;
    root->add_word(word, label);
  }

  size_t find_word(const char *word) const
  {
    assert(word);
    if (!root)
      return NOT_FOUND;
    return root->find_word(word);
  }

  void inorder(
    void (*open)(),
    void (*visit)(char ch, size_t label),
    void (*close)())
  {
    if (root)
      root->inorder(open, visit, close);
  }
};

/// Вспомогательный код для "распечатки" дерева.
ofstream file("test.cpp"); // Файл для вывода результата.
char tab[7]; // Строка текущих отступов.
bool start;  // Флаг "ожидается первый потомок".

void print_open()
{
  file << '\n' << tab << "(\n";
  auto pos = find_null(tab, sizeof(tab) - 1);
  if (pos != NOT_FOUND)
    tab[pos] = '\t';
  start = true;
}

void print_visit(char ch, size_t label)
{
  file << (start? tab: ", ") << ch;
  if (label != NOT_FOUND)
    file << " [" << label << ']';
  start = false;
}

void print_close()
{
  auto pos = find_null(tab, sizeof(tab));
  if (pos != 0)
    tab[pos-1] = '\0';
  file << '\n' << tab << ")\n";
  start = true;
}


int main()
{
  const char *words[]
  {
    "water",
    "was",
    "while",
    "when",
    "where",
    "why",
    "what",
    "whilst",
    "then",
    "there",
    "these",
    "today",
    "to",
    "the"
  };

  Trie trie{};

  // Поместить все строки из words в префиксное дерево.
  // В качестве метки -- индекс строки в массиве.
  for (size_t i = 0, iend = end(words) - begin(words); i < iend; ++i)
    trie.add_word(words[i], i);

  file << "trie\n{";
  // Вывести префиксное дерево в файл вместе с метками.
  trie.inorder(print_open, print_visit, print_close);
  file << "}\n";
  
  // Выводить метки слов, которые вводит пользователь.
  for (string word; cin >> word;)
  {
    const auto label = trie.find_word(word.c_str()); // token(word.c_str())
    if (label == NOT_FOUND)
      cout << "<not found>\n";
    else
      cout << label << ": \"" << words[label] << "\"\n";
  }
  
  return 0;
}

// csorts_i32.cpp
// Поразрядная сортировка для int32.
#include <cstddef>
#include <cassert>
#include <cstring>
#include <cstdlib>
#include <cstdint>
#include <utility>
#include <algorithm>
using namespace std;


////////////////////////////////////////////////////////////////////////////////
// Варианты поразрядной сортировки.

// От младших к старшим.
void lsd_radix_sort(int32_t *from, int32_t *to)
{
  assert(to - from < INT_MAX);

  // 11-bit digits: 32/11 = 3 digits.
  static const unsigned
    digit_bits = 11,
    last_digit_bits = 10;

  static const int
    buckets = 1 << digit_bits,
    last_buckets = 1 << last_digit_bits;

  static const uint32_t
    mask = buckets - 1,
    last_mask = last_buckets - 1;

  using I = int;
  const I items = to - from;
  if (items < buckets)
    return sort(from, to);

  const auto bfrom = (uint32_t*)from;
  const auto buf = new uint32_t[items];

  // Three statically allocated sets of offsets:
  // all digits are counted in one scan.
  I counters[2*buckets + last_buckets]{};
  I *const offs_1 = counters;
  I *const offs_2 = offs_1 + buckets;
  I *const offs_3 = offs_2 + buckets;

  // Run #0: compute offs_1, offs_2, offs_3.
  for (I i = 0; i < items; ++i)
  {
    const auto item = bfrom[i];
    offs_1[item & mask]++;
    offs_2[(item >> digit_bits) & mask]++;
    offs_3[((item >> 2*digit_bits) - last_buckets/2) & last_mask]++;
  }

  for (I i = 0, s1 = 0, s2 = 0; i < buckets; ++i)
  {
    const auto
      t1 = offs_1[i] + s1,
      t2 = offs_2[i] + s2;

    offs_1[i] = s1;
    offs_2[i] = s2;

    s1 = t1;
    s2 = t2;
  }

  for (I i = 0, s = 0; i < last_buckets; ++i)
  {
    const auto t = offs_3[i] + s;
    offs_3[i] = s;
    s = t;
  }

  // Additional optimization for better adaptivity:
  // Runs may be disabled if only one element of
  // the corresponding offs_x array is actually non-zero.
  // This may be checked in the summing loops above.

  // Run #1: first scatter.
  for (I i = 0; i < items; ++i)
  {
    const auto item = bfrom[i];
    buf[offs_1[item & mask]++] = item;
  }

  // Run #2: second scatter.
  for (I i = 0; i < items; ++i)
  {
    const auto item = buf[i];
    bfrom[offs_2[(item >> digit_bits) & mask]++] = item;
  }

  // Run #3: final scatter.
  for (I i = 0; i < items; ++i)
  {
    const auto item = bfrom[i];
    buf[offs_3[((item >> 2*digit_bits) - last_buckets/2) & last_mask]++] = item;
  }

  memcpy(bfrom, buf, items * sizeof(uint32_t));
  delete[] buf;
}


////////////////////////////////////////////////////////////////////////////////
// От старших к младшим: младший (самый правый) разряд.
void msd_radix_sort_2(uint32_t *from, uint32_t *to, uint32_t *buf)
{
  // 11-bit digits: 32/11 = 3 digits.
  static const unsigned digit_bits = 11;
  static const int buckets = 1 << digit_bits;
  static const uint32_t mask = buckets - 1;

  using I = int;
  const I items = to - from;
  if (items < buckets)
  {
    sort(from, to);
    memcpy(buf, from, items * sizeof(uint32_t));
    return;
  }

  // Three statically allocated sets of offsets:
  // all digits are counted in one scan.
  I offs[buckets]{};

  // Compute offs.
  for (I i = 0; i < items; ++i)
  {
    const auto item = from[i];
    offs[item & mask]++;
  }

  for (I i = 0, s = 0; i < buckets; ++i)
  {
    const auto t = offs[i] + s;
    offs[i] = s;
    s = t;
  }

  // Scatter.
  for (I i = 0; i < items; ++i)
  {
    const auto item = from[i];
    buf[offs[item & mask]++] = item;
  }
}


// От старших к младшим: средний разряд.
void msd_radix_sort_1(uint32_t *from, uint32_t *to, uint32_t *buf)
{
  // 11-bit digits: 32/11 = 3 digits.
  static const unsigned shift = 11, digit_bits = 11;
  static const int buckets = 1 << digit_bits;
  static const uint32_t mask = buckets - 1;

  using I = int;
  const I items = to - from;
  if (items < buckets)
  {
    sort(from, to);
    memcpy(buf, from, items * sizeof(uint32_t));
    return;
  }

  // Three statically allocated sets of offsets:
  // all digits are counted in one scan.
  I offs[buckets]{};

  // Compute offs.
  for (I i = 0; i < items; ++i)
  {
    const auto item = from[i];
    offs[(item >> shift) & mask]++;
  }

  for (I i = 0, s = 0; i < buckets; ++i)
  {
    const auto t = offs[i] + s;
    offs[i] = s;
    s = t;
  }

  // Scatter.
  for (I i = 0; i < items; ++i)
  {
    const auto item = from[i];
    buf[offs[(item >> shift) & mask]++] = item;
  }

  // Sort blocks.
  for (I i = 0, prev = 0; i < buckets; ++i)
  {
    msd_radix_sort_2(buf + prev, buf + offs[i], from + prev);
    prev = offs[i];
  }
}


// От старших к младшим: старший разряд.
void msd_radix_sort_0(uint32_t *from, uint32_t *to, uint32_t *buf)
{
  // 11-bit digits: 32/11 = 3 digits.
  static const unsigned shift = 22, digit_bits = 10;
  static const int buckets = 1 << digit_bits;
  static const uint32_t mask = buckets - 1;

  using I = int;
  const I items = to - from;

  // Three statically allocated sets of offsets:
  // all digits are counted in one scan.
  I offs[buckets]{};

  // Compute offs.
  for (I i = 0; i < items; ++i)
  {
    const auto item = from[i];
    offs[((item >> shift) - buckets/2) & mask]++;
  }

  for (I i = 0, s = 0; i < buckets; ++i)
  {
    const auto t = offs[i] + s;
    offs[i] = s;
    s = t;
  }

  // Scatter.
  for (I i = 0; i < items; ++i)
  {
    const auto item = from[i];
    buf[offs[((item >> shift) - buckets/2) & mask]++] = item;
  }

  // Sort blocks.
  for (I i = 0, prev = 0; i < buckets; ++i)
  {
    msd_radix_sort_1(buf + prev, buf + offs[i], from + prev);
    prev = offs[i];
  }
}


// От старших к младшим: основная процедура.
void msd_radix_sort(int32_t *from, int32_t *to)
{
  const auto items = to - from;
  assert(items < INT_MAX);
  if (items < 2048)
    return sort(from, to);

  const auto buf = new uint32_t[items];
  msd_radix_sort_0((uint32_t*)from, (uint32_t*)to, buf);
  memcpy(from, buf, items * sizeof(int32_t));
  delete[] buf;
}


////////////////////////////////////////////////////////////////////////////////
// Testing.

void std_stable_sort(int32_t *a, int32_t *ae)
{
  stable_sort(a, ae);
}


#include <ctime>
#include <cmath>
#include <iostream>
#include <random>
#include <iterator>
#include <vector>

double elapsed(clock_t t)
{
  return double(clock() - t) / CLOCKS_PER_SEC;
}


void do_tests(const vector<int32_t> &master)
{
  const auto mm = minmax_element(begin(master), end(master));
  cout << "items: " << master.size()*1e-6 << "M\t";
  cout << "min: " << *mm.first << "\tmax: " << *mm.second << '\n';

  auto sorted = master;
  auto t = clock();
  sort(begin(sorted), end(sorted));
  double std_time = elapsed(t);
  cout << "\nstd::sort\t" << std_time << "s\n\n";

  struct {
    void(*sort_alg)(int32_t*, int32_t*);
    const char *name;
  } algs[]
  {
    { std_stable_sort,     "std stable" },
    { lsd_radix_sort,      "lsd radix " },
    { msd_radix_sort,      "msd radix " },
  };

  unsigned lines = 0;
  for (auto &task : algs)
  {
    auto test = master;
    t = clock();
    task.sort_alg(test.data(), test.data() + test.size());
    double this_time = elapsed(t);
    cout << task.name << '\t' << this_time << "s,\t";
    cout << std_time / this_time << '\t';
    cout << (test == sorted) << endl;

    if (++lines == 6)
    {
      lines = 0;
      cout << '\n';
    }
  }
}


int main()
{
  cout.precision(3);
  mt19937_64 rng(14456);
  vector<int32_t> master(1 << 24);

  cout << "Uniform: non-negatives\n";
  {
    uniform_int_distribution<int32_t> vg;
    for (auto &val : master)
      val = vg(rng);
  }

  do_tests(master);
  cout << endl;

  cout << "Uniform: full range\n";
  {
    uniform_int_distribution<int32_t> vg(INT32_MIN, INT32_MAX);
    for (auto &val : master)
      val = vg(rng);
  }

  do_tests(master);
  cout << endl;
  return 0;
}

// csorts_f64.cpp
// Блочная и поразрядная сортировки для массивов значений типа double.
// Два примера с разными распределениями.
#include <cstddef>
#include <cassert>
#include <cstring>
#include <cstdlib>
#include <cstdint>
#include <utility>
#include <algorithm>
using namespace std;


////////////////////////////////////////////////////////////////////////////////
// Варианты блочной сортировки.

// Нерекурсивная блочная сортировка.
void nonrec_bucket_sort(double *from, double *to)
{
  using I = ptrdiff_t;

  // Размер блока -- 32кб.
  static const I bucket_size = 32768 / sizeof(double);
  const I items = to - from;
  if (items < 2 * bucket_size)
    return sort(from, to);

  // Определить минимум и максимум.
  // g++ векторизует этот цикл, если указан параметр -ffast-math.
  auto minv = *from, maxv = *from;
  for (I i = 1; i < items; ++i)
  {
    minv = from[i] < minv ? from[i] : minv;
    maxv = maxv < from[i] ? from[i] : maxv;
  }

  if (minv == maxv)
    return;

  // Количество блоков.
  const auto buckets = 1 + items / bucket_size;
  // Коэффициент аффинной функции, вычисляющей номер блока.
  const auto factor = (buckets - 1) / (maxv - minv);

  // Выделить вспомогательную память:
  // буфер для копирования значений: items * sizeof(double),
  // счётчики блоков: buckets * sizeof(I).
  assert(items < PTRDIFF_MAX / (2 * sizeof(double)));
  const I storage_size = items * sizeof(double) + buckets * sizeof(I);
  const auto storage = malloc(storage_size);
  if (!storage)
    return sort(from, to);

  const auto buf = (double*)storage;
  const auto off = (I*)(buf + items); // Размеры блоков ("смещения" = offsets)
  memset(off, 0, buckets * sizeof(I)); // Обнулить размеры.

  // Вычислить размеры блоков.
  for (I i = 0; i < items; ++i)
    off[(I)(factor * (from[i] - minv))]++;

  // Вычислить смещения и заполнить начальные значения позиций записи.
  for (I i = 0, s = 0; i < buckets; ++i)
  {
    const auto t = off[i] + s;
    off[i] = s;
    s = t;
  }

  // Раскидать значения по блокам.
  for (I i = 0; i < items; ++i)
    buf[off[(I)(factor * (from[i] - minv))]++] = from[i];

  // Отсортировать блоки.
  for (I i = 0, prev = 0; i < buckets; ++i)
  {
    sort(buf + prev, buf + off[i]);
    prev = off[i];
  }

  // Скопировать конкатенацию блоков в исходный массив.
  memcpy(from, buf, items * sizeof(double));
  free(storage);
}




// Рекурсивная блочная сортировка.
// Использует две функции: внешнюю, создающую буфер, и
// внутреннюю -- со статическим массивом счётчиков и ограниченной глубиной рекурсии.

// Внутренняя функция рекурсивной блочной сортировки.
void rec_bucket_sort_inner(double *from, double *to, double *buf, int depth)
{
  using I = ptrdiff_t;

  // Размер блока -- 32кб.
  static const I bucket_size = 32768 / sizeof(double);
  static const I max_buckets = 2048; // Максимальное количество блоков.

  const I items = to - from;
  if (0 == depth-- || items <= 2 * bucket_size)
    return sort(from, to);

  // Определить минимум и максимум.
  // g++ векторизует этот цикл, если указан параметр -ffast-math.
  auto minv = *from, maxv = *from;
  for (I i = 1; i < items; ++i)
  {
    minv = from[i] < minv ? from[i] : minv;
    maxv = maxv < from[i] ? from[i] : maxv;
  }

  if (minv == maxv)
    return;

  // Количество блоков.
  const auto possible_buckets = 1 + items / bucket_size;
  const auto buckets = possible_buckets < max_buckets ? possible_buckets : max_buckets;

  // Коэффициент аффинной функции, вычисляющей номер блока.
  const auto factor = (buckets - 1) / (maxv - minv);

  I off[max_buckets]{};

  // Вычислить размеры блоков.
  for (I i = 0; i < items; ++i)
    off[(I)(factor * (from[i] - minv))]++;

  // Вычислить начальные значения позиций записи.
  for (I i = 0, s = 0; i < buckets; ++i)
  {
    const auto t = off[i] + s;
    off[i] = s;
    s = t;
  }

  // Раскидать значения по блокам.
  for (I i = 0; i < items; ++i)
    buf[off[(I)(factor * (from[i] - minv))]++] = from[i];

  // Скопировать конкатенацию блоков в исходный массив.
  memcpy(from, buf, items * sizeof(double));

  // Отсортировать блоки "на месте" (уже в исходном массиве).
  for (I i = 0, prev = 0; i < buckets; ++i)
  {
    rec_bucket_sort_inner(from + prev, from + off[i], buf, depth);
    prev = off[i];
  }
}


// Внешняя функция рекурсивной блочной сортировки.
void rec_bucket_sort(double *from, double *to)
{
  const auto items = to - from;
  // Выделить вспомогательную память: буфер для копирования значений.
  const auto buf = new double[items];
  // Вызвать "вложенную" рекурсивную функцию.
  rec_bucket_sort_inner(from, to, buf, 4);
  delete[] buf;
}


// Рекурсивная блочная сортировка: челночное использование буфера и основного массива.
// Использует две функции: внешнюю, создающую буфер, и
// внутреннюю -- со статическим массивом счётчиков и ограниченной глубиной рекурсии.

// Внутренняя функция рекурсивной блочной сортировки.
void rec_bucket_sort_inner2(double *from, double *to, double *buf, int depth)
{
  using I = ptrdiff_t;

  // Размер блока -- 32кб.
  static const I bucket_size = 32768 / sizeof(double);
  static const I max_buckets = 2048; // Максимальное количество блоков.

  const I items = to - from;
  if (0 == depth || items <= 2 * bucket_size)
  {
    sort(from, to);
    if (depth % 2)
      memcpy(buf, from, items * sizeof(double));
    return;
  }

  // Определить минимум и максимум.
  // g++ векторизует этот цикл, если указан параметр -ffast-math.
  auto minv = *from, maxv = *from;
  for (I i = 1; i < items; ++i)
  {
    minv = from[i] < minv ? from[i] : minv;
    maxv = maxv < from[i] ? from[i] : maxv;
  }

  if (minv == maxv)
    return;

  // Количество блоков.
  const auto possible_buckets = 1 + items / bucket_size;
  const auto buckets = possible_buckets < max_buckets ? possible_buckets : max_buckets;

  // Коэффициент аффинной функции, вычисляющей номер блока.
  const auto factor = (buckets - 1) / (maxv - minv);

  I off[max_buckets]{};

  // Вычислить размеры блоков.
  for (I i = 0; i < items; ++i)
    off[(I)(factor * (from[i] - minv))]++;

  // Вычислить начальные значения позиций записи.
  for (I i = 0, s = 0; i < buckets; ++i)
  {
    const auto t = off[i] + s;
    off[i] = s;
    s = t;
  }

  // Раскидать значения по блокам.
  for (I i = 0; i < items; ++i)
    buf[off[(I)(factor * (from[i] - minv))]++] = from[i];

  // Отсортировать блоки "на месте" (уже в исходном массиве).
  --depth;
  for (I i = 0, prev = 0; i < buckets; ++i)
  {
    rec_bucket_sort_inner2(buf + prev, buf + off[i], from + prev, depth);
    prev = off[i];
  }
}


// Внешняя функция рекурсивной блочной сортировки.
void rec_bucket_sort2(double *from, double *to)
{
  const auto items = to - from;
  // Выделить вспомогательную память: буфер для копирования значений.
  const auto buf = new double[items];
  // Вызвать "вложенную" рекурсивную функцию.
  rec_bucket_sort_inner2(from, to, buf, 4);
  delete[] buf;
}


////////////////////////////////////////////////////////////////////////////////
// Варианты поразрядной сортировки.

// От младших к старшим.
void lsd_radix_sort(double *from, double *to)
{
  static_assert(sizeof(double) == sizeof(uint64_t));
  // 11-bit digits: 64/11 = 6 digits.
  static const unsigned digit_bits = 11, buckets = 1 << digit_bits;
  static const uint64_t mask = buckets - 1;

  using I = size_t;
  const I items = to - from;
  if (items < buckets)
    return sort(from, to);

  const auto bfrom = (uint64_t*)from;
  const auto buf = new uint64_t[items];

  // Two statically allocated sets of offsets:
  // one in use, other being computed while the scan.
  I counters[2 * buckets]{};
  I *const offs_1 = counters;
  I *const offs_2 = offs_1 + buckets;

  // Run #0: negatives-positives transform, compute offs_1.
  // Negatives-positives problem: inverse negatives, flip sign bit of positives.
  for (I i = 0; i < items; ++i)
  {
    const auto raw = bfrom[i];
    const auto item = raw ^ ((0 - (raw >> 63)) | (uint64_t(1) << 63));
    bfrom[i] = item;
    offs_1[item & mask]++;
  }

  for (I i = 0, s = 0; i < buckets; ++i)
  {
    const auto t = offs_1[i] + s;
    offs_1[i] = s;
    s = t;
  }

  // Run #1: first scatter, compute offs_2.
  for (I i = 0; i < items; ++i)
  {
    const auto item = bfrom[i];
    buf[offs_1[item & mask]++] = item;
    offs_2[(item >> digit_bits) & mask]++;
  }

  for (I i = 0, s = 0; i < buckets; ++i)
  {
    offs_1[i] = 0;
    const auto t = offs_2[i] + s;
    offs_2[i] = s;
    s = t;
  }

  // Run #2: second scatter, compute offs_1.
  for (I i = 0; i < items; ++i)
  {
    const auto item = buf[i];
    bfrom[offs_2[(item >> digit_bits) & mask]++] = item;
    offs_1[(item >> 2 * digit_bits) & mask]++;
  }

  // Now more copy-paste en route.
  for (I i = 0, s = 0; i < buckets; ++i)
  {
    offs_2[i] = 0;
    const auto t = offs_1[i] + s;
    offs_1[i] = s;
    s = t;
  }

  // Run #3: third scatter, compute offs_2.
  for (I i = 0; i < items; ++i)
  {
    const auto item = bfrom[i];
    buf[offs_1[(item >> 2 * digit_bits) & mask]++] = item;
    offs_2[(item >> 3 * digit_bits) & mask]++;
  }

  for (I i = 0, s = 0; i < buckets; ++i)
  {
    offs_1[i] = 0;
    const auto t = offs_2[i] + s;
    offs_2[i] = s;
    s = t;
  }

  // Run #4: fourth scatter, compute offs_1.
  for (I i = 0; i < items; ++i)
  {
    const auto item = buf[i];
    bfrom[offs_2[(item >> 3 * digit_bits) & mask]++] = item;
    offs_1[(item >> 4 * digit_bits) & mask]++;
  }

  for (I i = 0, s = 0; i < buckets; ++i)
  {
    offs_2[i] = 0;
    const auto t = offs_1[i] + s;
    offs_1[i] = s;
    s = t;
  }

  // Run #5: fifth scatter, compute offs_2.
  for (I i = 0; i < items; ++i)
  {
    const auto item = bfrom[i];
    buf[offs_1[(item >> 4 * digit_bits) & mask]++] = item;
    offs_2[(item >> 5 * digit_bits) & mask]++;
  }

  for (I i = 0, s = 0; i < buckets; ++i)
  {
    const auto t = offs_2[i] + s;
    offs_2[i] = s;
    s = t;
  }

  // Run #6: final scatter, negatives-positives inverse transform.
  for (I i = 0; i < items; ++i)
  {
    const auto item = buf[i];
    const auto raw = item ^ (((item >> 63) - 1) | (uint64_t(1) << 63));
    bfrom[offs_2[item >> 5 * digit_bits]++] = raw;
  }

  delete[] buf;
}


// От младших к старшим: вариант 2.
void lsd_radix_sort2(double *from, double *to)
{
  static_assert(sizeof(double) == sizeof(uint64_t));
  // 11-bit digits: 64/11 = 6 digits.
  static const unsigned digit_bits = 11, buckets = 1 << digit_bits;
  static const uint64_t mask = buckets - 1;

  using I = size_t;
  const I items = to - from;
  if (items < buckets)
    return sort(from, to);

  const auto bfrom = (uint64_t*)from;
  const auto buf = new uint64_t[items];

  // Two statically allocated sets of offsets:
  // one in use, other being computed while the scan.
  I counters[2 * buckets]{};
  const auto offs_1 = &counters[0];
  const auto offs_2 = &counters[buckets];

  // Run #0: negatives-positives transform, compute offs_1 & offs_2.
  // Negatives-positives problem: inverse negatives, flip sign bit of positives.
  for (I i = 0; i < items; ++i)
  {
    const auto raw = bfrom[i];
    const auto item = raw ^ ((0 - (raw >> 63)) | (uint64_t(1) << 63));
    bfrom[i] = item;
    offs_1[item & mask]++;
    offs_2[(item >> digit_bits) & mask]++;
  }

  for (I i = 0, s1 = 0, s2 = 0; i < buckets; ++i)
  {
    const auto
      t1 = offs_1[i] + s1,
      t2 = offs_2[i] + s2;

    offs_1[i] = s1;
    offs_2[i] = s2;
    s1 = t1;
    s2 = t2;
  }

  // Run #1: first scatter.
  for (I i = 0; i < items; ++i)
  {
    const auto item = bfrom[i];
    buf[offs_1[item & mask]++] = item;
  }

  // Run #2: second scatter.
  for (I i = 0; i < items; ++i)
  {
    const auto item = buf[i];
    bfrom[offs_2[(item >> digit_bits) & mask]++] = item;
  }

  // Compute offs_1 & offs_2 for the next two runs.
  memset(counters, 0, sizeof(counters));
  for (I i = 0; i < items; ++i)
  {
    const auto item = bfrom[i];
    offs_1[(item >> 2 * digit_bits) & mask]++;
    offs_2[(item >> 3 * digit_bits) & mask]++;
  }

  for (I i = 0, s1 = 0, s2 = 0; i < buckets; ++i)
  {
    const auto
      t1 = offs_1[i] + s1,
      t2 = offs_2[i] + s2;

    offs_1[i] = s1;
    offs_2[i] = s2;
    s1 = t1;
    s2 = t2;
  }

  // Run #3: third scatter.
  for (I i = 0; i < items; ++i)
  {
    const auto item = bfrom[i];
    buf[offs_1[(item >> 2 * digit_bits) & mask]++] = item;
  }

  // Run #4: fourth scatter.
  for (I i = 0; i < items; ++i)
  {
    const auto item = buf[i];
    bfrom[offs_2[(item >> 3 * digit_bits) & mask]++] = item;
  }

  // Compute offs_1 & offs_2 for the next two runs.
  memset(counters, 0, sizeof(counters));
  for (I i = 0; i < items; ++i)
  {
    const auto item = bfrom[i];
    offs_1[(item >> 4 * digit_bits) & mask]++;
    offs_2[(item >> 5 * digit_bits) & mask]++;
  }

  for (I i = 0, s1 = 0, s2 = 0; i < buckets; ++i)
  {
    const auto
      t1 = offs_1[i] + s1,
      t2 = offs_2[i] + s2;

    offs_1[i] = s1;
    offs_2[i] = s2;
    s1 = t1;
    s2 = t2;
  }

  // Run #5: fifth scatter.
  for (I i = 0; i < items; ++i)
  {
    const auto item = bfrom[i];
    buf[offs_1[(item >> 4 * digit_bits) & mask]++] = item;
  }

  // Run #6: final scatter, negatives-positives inverse transform.
  for (I i = 0; i < items; ++i)
  {
    const auto item = buf[i];
    const auto raw = item ^ (((item >> 63) - 1) | (uint64_t(1) << 63));
    bfrom[offs_2[item >> 5 * digit_bits]++] = raw;
  }

  delete[] buf;
}


// От старших к младшим: вложенные разряды (рекурсивная функция).
void msd_radix_sort_inner(uint64_t *from, uint64_t *to, uint64_t *buf, unsigned shift)
{
  // 11-bit digits: 64/11 = 6 digits.
  static const unsigned digit_bits = 11, buckets = 1 << digit_bits;
  static const uint64_t mask = buckets - 1;

  using I = size_t;
  const I items = to - from;
  if (items < 2 * buckets)
  {
    sort(from, to);
    // Negative-positive inverse transform.
    const auto dest = shift / digit_bits % 2 ? from : buf;
    for (I i = 0; i < items; ++i)
    {
      const auto item = from[i];
      const auto raw = item ^ (((item >> 63) - 1) | (uint64_t(1) << 63));
      dest[i] = raw;
    }

    return;
  }

  // Compute offsets.
  I offs[buckets]{};
  for (I i = 0; i < items; ++i)
    offs[(from[i] >> shift) & mask]++;

  // Prepare positions.
  for (I i = 0, s = 0; i < buckets; ++i)
  {
    const auto t = offs[i] + s;
    offs[i] = s;
    s = t;
  }

  // Scatter.
  if (shift != 0)
  {
    for (I i = 0; i < items; ++i)
    {
      const auto item = from[i];
      buf[offs[(item >> shift) & mask]++] = item;
    }

    // Sort blocks.
    shift -= digit_bits;
    for (I i = 0, prev = 0; i < buckets; ++i)
    {
      msd_radix_sort_inner(buf + prev, buf + offs[i], from + prev, shift);
      prev = offs[i];
    }
  }
  else
  {
    for (I i = 0; i < items; ++i)
    {
      const auto item = from[i];
      const auto raw = item ^ (((item >> 63) - 1) | (uint64_t(1) << 63));
      buf[offs[item & mask]++] = raw;
    }
  }
}


// От старших к младшим: внешняя функция.
void msd_radix_sort(double *from, double *to)
{
  static_assert(sizeof(double) == sizeof(uint64_t));
  // 11-bit digits: 64/11 = 6 digits.
  static const unsigned
    digit_bits = 11,
    shift = 55,
    buckets = 1 << (sizeof(double) * CHAR_BIT - shift); // 1 << 9

  using I = size_t;
  const I items = to - from;
  if (items < buckets)
    return sort(from, to);

  const auto bfrom = (uint64_t*)from;
  const auto buf = new uint64_t[items];

  // Offsets.
  I offs[buckets]{};

  // Run #0: negatives-positives transform, compute offs.
  // Negatives-positives problem: inverse negatives, flip sign bit of positives.
  for (I i = 0; i < items; ++i)
  {
    const auto raw = bfrom[i];
    const auto item = raw ^ ((0 - (raw >> 63)) | (uint64_t(1) << 63));
    bfrom[i] = item;
    offs[item >> shift]++;
  }

  for (I i = 0, s = 0; i < buckets; ++i)
  {
    const auto t = offs[i] + s;
    offs[i] = s;
    s = t;
  }

  // Run #1: first scatter.
  for (I i = 0; i < items; ++i)
  {
    const auto item = bfrom[i];
    buf[offs[item >> shift]++] = item;
  }

  // Sort recursively.
  for (I i = 0, prev = 0; i < buckets; ++i)
  {
    msd_radix_sort_inner(buf + prev, buf + offs[i], bfrom + prev, shift - digit_bits);
    prev = offs[i];
  }

  delete[] buf;
}


////////////////////////////////////////////////////////////////////////////////
// Testing.

void std_stable_sort(double *a, double *ae)
{
  stable_sort(a, ae);
}


#include <ctime>
#include <cmath>
#include <iostream>
#include <random>
#include <iterator>
#include <vector>

double elapsed(clock_t t)
{
  return double(clock() - t) / CLOCKS_PER_SEC;
}


void do_tests(const vector<double> &master)
{
  const auto mm = minmax_element(begin(master), end(master));
  cout << "items: " << master.size()*1e-6 << "M\t";
  cout << "min: " << *mm.first << "\tmax: " << *mm.second << '\n';

  auto sorted = master;
  auto t = clock();
  sort(begin(sorted), end(sorted));
  double std_time = elapsed(t);
  cout << "\nstd::sort\t" << std_time << "s\n\n";

  struct {
    void(*sort_alg)(double*, double*);
    const char *name;
  } algs[]
  {
    { std_stable_sort,     "std stable" },
    { nonrec_bucket_sort,  "bucket    " },
    { rec_bucket_sort,     "rbucket   " },
    { rec_bucket_sort2,    "rbucket2  " },
    { lsd_radix_sort,      "lsd radix " },
    { lsd_radix_sort2,     "lsd radix2" },
    { msd_radix_sort,      "msd radix " },
  };

  unsigned lines = 0;
  for (auto &task : algs)
  {
    auto test = master;
    t = clock();
    task.sort_alg(test.data(), test.data() + test.size());
    double this_time = elapsed(t);
    cout << task.name << '\t' << this_time << "s,\t";
    cout << std_time / this_time << '\t';
    cout << (test == sorted) << endl;

    if (++lines == 4)
    {
      lines = 0;
      cout << '\n';
    }
  }
}


int main()
{
  cout.precision(3);
  mt19937_64 rng(14456);
  vector<double> master(1 << 24);

  cout << "Normal(0, 10) + x2\n";
  {
    normal_distribution<double> vg(0., 10.);

    double x = -10., dx = 20. / master.size();
    for (auto &val : master)
    {
      val = vg(rng) + x * x;
      x += dx;
    }
  }

  do_tests(master);
  cout << endl;

  cout << "------------------------\n"
    "Uniform in [-1e+6, 1e+6]\n";
  {
    uniform_real_distribution<double> vg(-1e+6, 1e+6);
    for (auto &val : master)
      val = vg(rng);
  }

  do_tests(master);
  cout << endl;
  return 0;
}

// csorts_s8.cpp
// Поразрядная сортировка 7-битных (ANSI) си-строк.
// Лексикографический порядок, индуцированный порядком кодов символов.
// (Впрочем, поведение strcmp, используемого для "досортировки" малых
//  блоков, не так просто и зависит от настройки локализации.)
#include <cstddef>
#include <cstring>
#include <algorithm>

// Стандартная сортировка.
void std_sort(const char **from, const char **to)
{
  std::sort(from, to,
    [](const char *a, const char *b)
    { return std::strcmp(a, b) < 0; });
}

// Стандартная устойчивая сортировка.
void std_stable_sort(const char **from, const char **to)
{
  std::stable_sort(from, to,
    [](const char *a, const char *b)
    { return std::strcmp(a, b) < 0; });
}


// Внешняя функция.
void ansi_msd_radix_sort(const char **from, const char **to)
{
  // Объявление внутренней (рекурсивной) функции сортировки.
  void ansi_msd_radix_sort_inner
    (const char **from, const char **to, const char **buf, int depth);

  // Маленькие наборы сортируем стандартной сортировкой.
  const auto items = to - from;
  if (items < 256)
    return std_sort(from, to);

  // Создать буфер и отсортировать рекурсивно.
  const auto buf = new const char*[items];
  ansi_msd_radix_sort_inner(from, to, buf, 0);
  delete[] buf;
}


// Внутренняя функция.
void ansi_msd_radix_sort_inner
  (const char **from, const char **to, const char **buf, int depth)
{
  static const int max_depth = 10;
  static const unsigned buckets = 128;
  using I = std::size_t;

  const I items = to - from;
  // Если блок маленький, или достигнута макс. глубина
  // -- "досортировать" стандартной функцией "хвосты" и закончить.
  if (items < buckets || depth == max_depth)
  {
    std_sort(from, to);

    // Скопировать элементы "обратно", поправив указатели.
    const auto dest = depth % 2? buf: from;
    for (I i = 0; i < items; ++i)
      dest[i] = from[i] - depth;

    return;
  }

  // Вычислить размеры блоков.
  I offs[buckets]{};
  for (I i = 0; i < items; ++i)
  {
    const unsigned char bucket = *from[i];
    offs[bucket]++;
  }

  // Вычислить стартовые позиции.
  for (I i = 0, s = 0; i < buckets; ++i)
  {
    const auto t = offs[i] + s;
    offs[i] = s;
    s = t;
  }

  // Разбросать по блокам в буфер и
  // сдвинуть указатели на непустые строки на один символ вперёд.
  for (I i = 0; i < items; ++i)
  {
    const auto item = from[i];
    const unsigned char bucket = *item;
    buf[offs[bucket]++] = item + (bucket != 0);
  }

  // Вернуть указатели на "пустые" строки в исходный массив и на исходные позиции.
  {
    const auto dest = depth % 2? buf: from;
    for (I i = 0, sz = offs[0]; i < sz; ++i)
      dest[i] = buf[i] - depth;
  }

  // Обработать рекурсивно непустые строки.
  ++depth;
  for (I i = 1; i < buckets; ++i)
    ansi_msd_radix_sort_inner
      (buf + offs[i-1], buf + offs[i], from + offs[i-1], depth);
}


////////////////////////////////////////////////////////////////////////////////
// Тестирование.
#include <ctime>
#include <vector>
#include <string>
#include <random>
#include <iostream>
using namespace std;


double elapsed(clock_t t)
{
  return double(clock() - t) / CLOCKS_PER_SEC;
}


vector<string> random_string_vector(size_t sz)
{
  // Инициализация генератора псевдослучайных чисел.
  mt19937_64 rng(19'08'1991);
  uniform_int_distribution<char> cg('\x20', '\x7e'); // неестественное распределение символов (равномерное).
  binomial_distribution<unsigned> lg(40, 0.5); // распределение длин строк.

  // Результирующий вектор строк.
  vector<string> strings(sz);
  for (auto &s: strings)
  {
    // Размер строки.
    s.resize(min(150u, lg(rng)));
    // Коды символов строки.
    for (auto &ch: s)
      ch = cg(rng);
  }

  return strings;
}


int main()
{
  cout.precision(3);
  cout << "An example of random strings being sorted:\n\n";
  for (auto &s: random_string_vector(12))
    cout << s << '\n';
  cout << endl;

  const auto storage = random_string_vector(1 << 22);
  vector<const char*> master(storage.size());
  for (size_t i = 0, sz = storage.size(); i < sz; ++i)
    master[i] = storage[i].c_str();

  auto sorted = master;
  auto t = clock();
  std_sort(sorted.data(), sorted.data() + sorted.size());
  double std_time = elapsed(t);
  cout << "\nstd::sort\t" << std_time << "s\n\n";

  struct {
    void (*sort_alg)(const char**, const char**);
    const char *name;
  } algs[]
  {
    { std_stable_sort,     "std stable" },
    { ansi_msd_radix_sort, "msd radix " },
  };

  unsigned lines = 0;
  for (auto &task : algs)
  {
    auto test = master;
    t = clock();
    task.sort_alg(test.data(), test.data() + test.size());
    double this_time = elapsed(t);
    cout << task.name << '\t' << this_time << "s,\t";
    cout << std_time / this_time << '\t';
    cout << (test == sorted) << endl;

    if (++lines == 4)
    {
      lines = 0;
      cout << '\n';
    }
  }

  return 0;
}

// subshashs.cpp
// substrings hashing search
#include <cstdint>
#include <cstring>
using namespace std;

/// Вычислить хэш префикса s.
uint32_t prefix_hash(const char *s)
{
  // Расширение нулями.
  char bytes[sizeof(uint32_t)]{};
  strncpy(bytes, s, sizeof(uint32_t));

  // Интерпретировать как число.
  uint32_t result;
  memcpy(&result, bytes, sizeof(uint32_t));
  return result;
}

/// Выбрать возможный префикс из набора по хэшу.
/// Возвращает успех выборки.
bool select_prefix(uint32_t hash, const char **prefix)
{
  static const char *words[]
  {
    "delta",
    "alpha",
    "beta",
    "gamma",
  };

  static const uint32_t hashes[]
  {
    0x746c'6564,
    0x6870'6c61,
    0x6174'6562,
    0x6d6d'6167,
  };

  if (hashes[hash % 4] != hash)
    return false;

  *prefix = words[hash % 4];
  return true;
}


/// Найти все вхождения искомых строк (alpha, beta, gamma, delta).
void find_all_matches(const char *text,
  void (*report_match)(const char *prefix, size_t pos))
{
  const char *prefix = nullptr;
  for (size_t pos = 0; *text != '\0'; ++text, ++pos)
  {
    if (select_prefix(prefix_hash(text), &prefix)
        && strncmp(text, prefix, strlen(prefix)) == 0)
      report_match(prefix, pos);
  }
}


////////////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <iomanip>

const char *text = "There are some alphas and betas or gammas but probably no deltas, bita is not a proper spelling.\n";
int main()
{
  cout << "Text:\n" << text;
  find_all_matches(text,
    [](const char *prefix, size_t pos)
    {
      cout << "\nfound: " << prefix << " at " << pos << '\n' << (text + pos);
    });

  return 0;
}

// rabinkarp.cpp
// Rabin-Karp algorithm
// алгоритм поиска подстроки Рабина-Карпа
// кольцевой хэш -- аналог вычисления значения числа, записанного в системе счисления с базой RK_HASH_RADIX
#include <cstring> // memcmp
using namespace std;

/// Стандарт, увы, не определяет целочисленных вариантов pow.
size_t ipow(size_t x, size_t y)
{
  size_t p = 1;
  while (y != 0)
  {
    p *= y & 1? x: 1;
    x *= x;
    y >>= 1;
  }

  return p;
}


/// Основание системы счисления для rk_hash.
const size_t RK_HASH_RADIX = 39989; // Простое число.

/// Вычислить хэш строки.
size_t rk_hash(const char *text, size_t text_size)
{
  size_t s = 0; // сумма
  while (text_size--)
  {
    s *= RK_HASH_RADIX;
    s += (unsigned char)*text++;
  }
  return s;
}

/// Выполнить поиск подстроки алгоритмом Рабина-Карпа.
size_t rk_search(const char *where, size_t where_len,
                 const char *what,  size_t what_len)
{
  if (where_len < what_len)
    return where_len;
  if (what_len == 0)
    return 0;

  // Вычислить хэш what.
  const auto
    what_hash = rk_hash(what, what_len),
    radix_pow = ipow(RK_HASH_RADIX, what_len - 1);

  // Подготовить хэш окна поиска, обработав префикс where.
  auto window_hash = rk_hash(where, what_len);

  // Идти по строке, сдвигая окно и проверяя хэш.
  // i и j -- индексы крайних символов окна.
  for (size_t i = what_len, j = 0;; ++i, ++j)
  {
    if (what_hash == window_hash
     && memcmp(where + j, what, what_len) == 0)
      return j;

    if (i == where_len)
      return where_len;

    // сдвинуть окно и поправить хэш
    const unsigned char
      old_byte = where[j],
      new_byte = where[i];

    ((window_hash -= old_byte * radix_pow)
                  *= RK_HASH_RADIX)
                  += new_byte;
  }
}


////////////////////////////////////////////////////////////////////////////////
// Тестирование.
#include <iostream>
#include <string>

int main()
{
  string line, sub;
  while (true)
  {
    getline(cin, line);
    getline(cin, sub);
    if (!cin)
      return 0;

    for (size_t p = 0; p < line.size(); p += sub.size())
    {
      p += rk_search(line.c_str() + p, line.size() - p, sub.c_str(), sub.size());
      if (p < line.size())
        cout << p << '\n';
    }
  }
}

// rabinkarp2.cpp
// Поиск всех вхождений заданного набора подстрок в тексте.
// Организует простейшую хэш-таблицу.
#include <cstring> // memcmp, strlen
#include <cassert>
using namespace std;

/// Стандарт, увы, не определяет целочисленных вариантов pow.
size_t ipow(size_t x, size_t y)
{
  size_t p = 1;
  while (y != 0)
  {
    p *= y & 1? x: 1;
    x *= x;
    y >>= 1;
  }

  return p;
}


/// Основание системы счисления для rk_hash.
const size_t RK_HASH_RADIX = 39989; // Простое число.

/// Вычислить хэш строки.
size_t rk_hash(const char *text, size_t text_size)
{
  size_t s = 0; // сумма
  while (text_size--)
  {
    s *= RK_HASH_RADIX;
    s += (unsigned char)*text++;
  }
  return s;
}


/// Набор искомых строк с дополнительной информацией.
struct RK_subs
{
  /// Записи, соответствующие искомым подстрокам.
  struct Entry
  {
    const char *value;
    size_t len;
    size_t hash;
    size_t radix_pow;
    size_t window_hash;
  };

  Entry *entries = nullptr;
  size_t size = 0, active = 0;

  /// Инициализировать записи.
  RK_subs(const char *subs[], size_t subs_size)
    : entries(new Entry[subs_size]),
      size(subs_size),
      active(subs_size)
  {
    for (size_t i = 0; i < subs_size; ++i)
    {
      const auto len = strlen(subs[i]);
      entries[i] = Entry{ subs[i], len, rk_hash(subs[i], len) };
    }
  }

  ~RK_subs() { delete[] entries; }

  /// Инициализировать окна.
  void init_windows(const char *text, size_t text_len)
  {
    for (size_t i = 0; i < size; ++i)
    {
      auto &e = entries[i];
      if (text_len < e.len)
      {
        e.value = nullptr; // "деактивировать" запись.
        --active;
      }
      else
      {
        e.window_hash = rk_hash(text, e.len); // неэффективно...
        e.radix_pow = ipow(RK_HASH_RADIX, e.len-1);
      }
    }
  }

  /// Поправить окна: шаг вперёд на один символ.
  void shift_windows(const char *text, size_t text_len, size_t old_pos)
  {
    for (size_t i = 0; i < size; ++i)
    {
      auto &e = entries[i];
      if (e.value) // активна?
      {
        if (text_len <= old_pos + e.len)
        {
          e.value = nullptr; // деактивировать
          --active;
        }
        else
        {
          const unsigned char
            old_byte = text[old_pos],
            new_byte = text[old_pos + e.len];

          ((e.window_hash -= old_byte * e.radix_pow)
                          *= RK_HASH_RADIX)
                          += new_byte;
        }
      }
    }
  }

  /// Найти все совпадения.
  void find_matches(const char *text, size_t cur_pos,
    void (*report_match)(const char*, size_t)) const
  {
    const auto start = text + cur_pos;
    for (size_t i = 0; i < size; ++i)
    {
      const auto &e = entries[i];
      if (e.value && e.hash == e.window_hash
       && memcmp(start, e.value, e.len) == 0)
         report_match(e.value, cur_pos);
    }
  }
};


/// Выполнить поиск всех подстрок алгоритмом Рабина-Карпа.
void rk_search_all
  (const char *text,
   const char *what[], size_t what_size,
   void (*report_match)(const char*, size_t))
{
  const auto text_len = strlen(text);
  if (text_len == 0 || what_size == 0)
    return;

  // Создать окна поиска.
  RK_subs searcher(what, what_size);
  searcher.init_windows(text, text_len);

  // Искать, пока есть "помещающиеся" окна.
  for (size_t pos = 0; searcher.active != 0 && pos != text_len; ++pos)
  {
    searcher.find_matches(text, pos, report_match);
    searcher.shift_windows(text, text_len, pos);
  }
}



////////////////////////////////////////////////////////////////////////////////
// Тестирование.
#include <iostream>
#include <string>

const char *text = "alpha bet beta to send gamma into delta, that's all about alpha and beta";
int main()
{
  const char *subs[]
  {
    "alpha",
    "beta",
    "gamma",
    "delta",
  };

  rk_search_all(text, subs, 4,
    [](const char *sub, size_t pos)
    {
      cout << "found: " << sub << " at " << pos << '\n';
      cout << (text + pos) << '\n';
    });

  return 0;
}

// lcs.cpp
// longest common substring
// самя длинная общая подстрока
#include <cassert>
#include <cstring>
#include <string>
#include <algorithm>
#include <utility>
using namespace std;

/// Решение: позиции начала в первой и второй строках и длина.
struct Lcs_solution
{
  size_t pos1, pos2, len;
};

inline bool operator==(const Lcs_solution &a, const Lcs_solution &b)
{
  return a.pos1 == b.pos1 && a.pos2 == b.pos2 && a.len == b.len;
}

inline bool operator!=(const Lcs_solution &a, const Lcs_solution &b)
{
  return !(a == b);
}


/// Наивный переборный алгоритм: O(N3)-время, O(1)-память.
/// (На самом деле O(asz*bsz + sum(match(i, j), 0<=i<asz, 0<=j<bsz)),
///  что на практике близко O(N2).)
Lcs_solution lcs_naive(const string &a, const string &b)
{
  Lcs_solution result {};

  const size_t asz = a.size(), bsz = b.size();
  for (size_t i = 0; i < asz; ++i)
  {
    for (size_t j = 0; j < bsz; ++j)
    {
      const size_t max_match = min(asz - i, bsz - j);
      if (max_match <= result.len)
        break;

      size_t match = 0;
      while (match < max_match && a[i + match] == b[j + match])
        ++match;

      if (result.len < match)
        result = Lcs_solution { i, j, match };
    }

    if (asz - i <= result.len)
      break;
  }

  return result;
}


/// Переборный алгоритм + мемоизация: O(N2)-время, O(N2)-память.
/// "Словарь".
struct Lcs_memo
{
  const string &a, &b;
  size_t **m, asz, bsz;

  // Конструктор объектов Memo.
  Lcs_memo(const string &a, const string &b)
    : a(a), b(b), asz(a.size()), bsz(b.size())
  {
    // Создать таблицу в динамической памяти.
    assert(asz && bsz && asz < SIZE_MAX/bsz/sizeof(size_t));
    m = new size_t*[asz];
    m[0] = new size_t[asz * bsz];
    memset(m[0], -1, asz * bsz * sizeof(size_t)); // -1 == "not stored"
    for (size_t j = 1; j < asz; ++j)
      m[j] = m[j-1] + bsz;
  }

  // Деструктор объектов Memo. Освобождает память.
  ~Lcs_memo() { delete[] m[0]; delete[] m; }

  // Проверить, вычислена ли ячейка таблицы с индексами i, j.
  bool stored(size_t i, size_t j) const
  {
    return m[i][j] != size_t(-1);
  }

  // Вычислить совпадение, начиная с индексов i, j.
  size_t match(size_t i, size_t j)
  {
    if (i == asz || j == bsz)
      return 0;

    if (stored(i, j))
      return m[i][j];

    return m[i][j] = a[i] != b[j]? 0: match(i+1, j+1) + 1;
  }
};

Lcs_solution lcs_memo(const string &a, const string &b)
{
  Lcs_solution result {};
  Lcs_memo memo(a, b);

  const size_t asz = a.size(), bsz = b.size();
  for (size_t i = 0; i < asz; ++i)
  {
    for (size_t j = 0; j < bsz; ++j)
    {
      const size_t match = memo.match(i, j);
      if (result.len < match)
        result = Lcs_solution { i, j, match };
    }

    if (asz - i <= result.len)
      break;
  }

  return result;
}


/// Мемоизация с запоминанием только двух строк таблицы: O(N2)-время, O(N)-память.
Lcs_solution lcs_memo2(const string &a, const string &b)
{
  Lcs_solution result{};
  const size_t asz = a.size(), bsz = b.size();
  auto m_prev = new size_t[bsz + 1]{}; // в начале нули.
  auto m_next = new size_t[bsz + 1];
  m_next[bsz] = 0;

  // С конца к началу.
  for (size_t i = asz; i-- != 0;)
  {
    const auto ai = a[i];
    for (size_t j = bsz; j-- != 0;)
    {
      if (ai != b[j])
      {
        m_next[j] = 0;
      }
      else
      {
        const auto val = m_next[j] = m_prev[j + 1] + 1;
        if (result.len <= val)
          result = Lcs_solution{ i, j, val };
      }
    }

    swap(m_next, m_prev);
  }

  delete[] m_next;
  delete[] m_prev;
  return result;
}


/// Мемоизация с запоминанием только двух строк таблицы: O(N2)-время, O(N)-память.
Lcs_solution lcs_memo3(const string &a, const string &b)
{
  Lcs_solution result {};
  const size_t asz = a.size(), bsz = b.size();
  auto m_prev = new size_t[bsz+1]{}; // в начале нули.
  auto m_next = new size_t[bsz+1];
  m_next[bsz] = 0;

  // С конца к началу.
  for (size_t i = asz; i-- != 0;)
  {
    const auto ai = a[i];
    for (size_t j = 0; j < bsz; ++j)
    {
      const auto nz = m_prev[j+1] + 1;
      m_next[j] = ai != b[j]? 0: nz;
    }

    const auto max_match = max_element(m_next, m_next + bsz);
    if (result.len <= *max_match)
      result = Lcs_solution{ i, size_t(max_match - m_next), *max_match };

    swap(m_next, m_prev);
  }

  delete[] m_next;
  delete[] m_prev;
  return result;
}


////////////////////////////////////////////////////////////////////////////////
// Максимальная общая подпоследовательность.

/// Мемоизация с запоминанием только двух строк таблицы: O(N3)-время, O(N2)-память.
/// Из-за копирования строк получается O(N3)-время, можно улучшить до O(N2)-времени.
string longest_common_subsequence(const string &a, const string &b)
{
  const size_t asz = a.size(), bsz = b.size();
  auto m_prev = new string[bsz + 1]; // Строки будут изначально пустые,
  auto m_next = new string[bsz + 1]; // благодаря тому, что для string определён "конструктор по умолчанию".
  string result;

  // С конца к началу.
  for (size_t i = asz; i-- != 0;)
  {
    const auto ai = a[i];
    for (size_t j = bsz; j-- != 0;)
    {
      if (ai != b[j])
        m_next[j] = m_next[j+1].size() < m_prev[j].size()?
          m_prev[j]: m_next[j+1];
      else
        m_next[j] = ai + m_prev[j + 1];
    }

    const auto &candidate = *max_element(m_next, m_next+bsz,
      [](const string &a, const string &b) // Найти самую длинную строку.
      { return a.size() < b.size(); });
    if (result.size() < candidate.size())
      result = candidate;

    swap(m_next, m_prev);
  }

  delete[] m_next;
  delete[] m_prev;
  return result;
}


/// Мемоизация с запоминанием только двух строк таблицы: O(N3)-время, O(N2)-память.
/// Из-за копирования строк получается O(N3)-время, можно улучшить до O(N2)-времени.
/// Другой порядок: с начала в конец.
string longest_common_subsequence2(const string &a, const string &b)
{
  if (a.size() < b.size())
    return longest_common_subsequence2(b, a);
    
  const size_t asz = a.size(), bsz = b.size();
  const auto storage = new string[2*bsz + 2];
  auto m_prev = storage + 1, m_next = m_prev + bsz + 1;
  
  string result;

  for (size_t i = 0; i < asz; ++i)
  {
    const auto ai = a[i];
    for (size_t j = 0; j < bsz; ++j)
    {
      if (ai != b[j])
        m_next[j] = m_next[j-1].size() < m_prev[j].size()?
          m_prev[j]: m_next[j-1];
      else
        m_next[j] = m_prev[j-1] + ai;
    }

    const auto &candidate = *max_element(m_next, m_next+bsz,
      [](const string &a, const string &b) // Найти самую длинную строку.
      { return a.size() < b.size(); });
    if (result.size() < candidate.size())
      result = candidate;

    swap(m_next, m_prev);
  }

  delete[] storage;
  return result;
}


////////////////////////////////////////////////////////////////////////////////
// Тестирование.
#include <iostream>
#include <random>
#include <chrono>


/// Генерирует псевдослучайную строку.
string random_string(size_t len, unsigned seed = 421129)
{
  string result(len, '\0');
  mt19937_64 rng(seed);
  for (auto &ch: result)
  {
    unsigned code = rng() % 32;
    ch = code < 26? 'a'+code: ' '; // латинские буквы или пробелы.
  }
  return result;
}

/// Перемешивает строки (чтобы появилось много совпадений).
void crosscopy(string &a, string &b, unsigned seed = 2017)
{
  mt19937_64 rng(seed);
  const auto sz = min(a.size(), b.size());
  bool first = false;
  for (size_t i = 0, j = 0; i < sz; ++i)
  {
    if (j == 0)
    {
      j = (rng() % sz) / 16;
      first = rng() % 2 == 0;
    }
    else
    {
      (first? a: b)[i] = (first? b: a)[i];
      --j;
    }
  }

  rotate(a.begin(), a.begin() + (rng() % sz), a.end());
}


/// Тестирование на заранее заданных строках.
int test_lcs(Lcs_solution (*lcs)(const string&, const string&))
{
  if (lcs("", "") != Lcs_solution{})
    return 1;
  if (lcs("", "abc") != Lcs_solution{})
    return 2;
  if (lcs("abc", "") != Lcs_solution{})
    return 3;
  if (lcs("abcdef", "ghik") != Lcs_solution{})
    return 4;
  if (lcs("abcd", "efga") != Lcs_solution{0, 3, 1})
    return 5;
  if (lcs("abra cadabra", "mabra fabra") != Lcs_solution{0, 1, 5})
    return 6;
  if (lcs("this is a description of another situation",
          "that situation has been described already") != Lcs_solution{32, 4, 10})
    return 7;

  return 0;
}


/// Тестирование производительности поиска общей подстроки.
void time_test(Lcs_solution (*lcs)(const string&, const string&), size_t len = 1 << 14)
{
  using namespace chrono;
  auto
    a = random_string(len, 2011),
    b = random_string(len, 2017);
  crosscopy(a, b);

  const auto t0 = high_resolution_clock::now();
  const auto p = lcs(a, b);
  const auto t1 = high_resolution_clock::now();
  cout << (int)duration<double, std::milli>(t1 - t0).count() << " ms  "
       << p.pos1 << ' ' << p.pos2 << ' ' << p.len << endl;
}


int main()
{
  cout << test_lcs(lcs_naive);
  cout << test_lcs(lcs_memo);
  cout << test_lcs(lcs_memo2);
  cout << test_lcs(lcs_memo3) << '\n';

  cout << "LCsubsequence: ";
  cout << longest_common_subsequence("abcd bcd ab", "dbca cda") << '\n';
  cout << longest_common_subsequence("this is a description of another situation",
          "that situation has been described already") << '\n';
  cout << longest_common_subsequence2("this is a description of another situation",
          "that situation has been described already") << '\n';

  cout << "\nRandom string: " << random_string(20) << '\n';

  size_t len = 4'000;
  cout << "\n2 x " << len << " characters:";

  cout << "\nNaive: ";
  time_test(lcs_naive, len);

  cout << "Memo:  ";
  time_test(lcs_memo, len);

  cout << "Memo2: ";
  time_test(lcs_memo2, len);

  cout << "Memo3: ";
  time_test(lcs_memo3, len);

  len = 20'000;
  cout << "\n2 x " << len << " characters:";

  cout << "\nNaive: ";
  time_test(lcs_naive, len);

  cout << "Memo2: ";
  time_test(lcs_memo2, len);

  cout << "Memo3: ";
  time_test(lcs_memo3, len);

  cout << endl;
  return 0;
}

// vt_println.cpp
#include <iostream>

inline void print() {} // Нечего печатать.

template <class Ch, class Tr>
void print(std::basic_ostream<Ch, Tr>&) {} // Нечего печатать.

template <class Ch, class Tr, class First, class... Other>
void print(std::basic_ostream<Ch, Tr> &os,
  const First &first, const Other&... other)
{
  os << first;
  print(os, other...);
}

template <class... Other>
void print(const Other&... other)
{
  print(std::cout, other...);
}

template <class Ch, class Tr, class... Other>
void println(std::basic_ostream<Ch, Tr> &os, const Other&... other)
{
  print(os, other..., '\n');
}

template <class... Args>
void println(const Args&... args)
{
  println(std::cout, args...);
}


int main()
{
  println(42, " is the new ", 23);
  println(3.1415926535897932);
  return 0;
}

// vt_format.cpp
#include <string>
#include <cstddef>
#include <cctype>

namespace impl
{
  template <class T>
  std::string to_str(const T &val)
  {
    using std::to_string;
    return to_string(val);
  }

  inline std::string to_str(const char *p)
  {
    return p;
  }

  inline const std::string& to_str(const std::string &s)
  {
    return s;
  }
}


template <class... Params>
std::string format(const std::string &model, const Params&... params)
{
  using namespace std;
  // Конвертировать params в массив строк.
  string param[] { "", impl::to_str(params)... };

  // Подготовить место для результата.
  size_t result_possible_size = model.size();
  for (auto &p: param)
    result_possible_size += p.size();

  string result;
  result.reserve(result_possible_size);

  // Выполнить подстановку (конечный автомат).
  size_t param_index;
  auto p = model.c_str();
  auto stored_p = p;

outer:
  while (auto ch = *p++)
  {
    if (ch == '%')
      goto percent;
    result += ch;
  }

  return result;

percent:
  stored_p = p;
  switch (auto ch = *p++)
  {
  case '\0':
    result += '%';
    return result;

  case '%':
    result += '%';
    goto outer;

  default:
    if (!isdigit(ch))
    {
      (result += '%') += ch;
      goto outer;
    }

    param_index = ch - '0';
  }

  while (auto ch = *p++)
  {
    if (ch == '%' && param_index <= sizeof...(params))
    {
      result += param[param_index];
      goto outer;
    }

    if (!isdigit(ch))
      break;

    param_index = param_index * 10 + (ch - '0');
  }

  result += '%';
  p = stored_p;
  goto outer;
}


#include <iostream>
int main()
{
  std::cout << format(
      "%2% %3% of %1%\n"
      "%20%\n"
      "%2broken%\n"
      "%2%% of %1%, percent-sign: %%\n"
      "not closed %\n",
    "acid",
    37.5,
    "ounces",
    "not used?");
  return 0;
}

// vt_format_regex.cpp
#include <string>
#include <cstddef>
#include <regex>

namespace impl
{
  template <class T>
  std::string to_str(const T &val)
  {
    using std::to_string;
    return to_string(val);
  }

  inline std::string to_str(const char *p)
  {
    return p;
  }

  inline const std::string& to_str(const std::string &s)
  {
    return s;
  }
}

template <class... Params>
std::string format(std::string model, const Params&... params)
{
  using namespace std;
  // Конвертировать params в массив строк.
  string param[] { "", impl::to_str(params)... };

  // Выполнить подстановку (конечный автомат через std::regex).
  string result;

  regex preg("%(\\d*)%");
  smatch match;
  size_t i;
  while (regex_search(model, match, preg))
  {
    result += match.prefix();
    if (match.size() == 1 || match[1].length() == 0)
      result += '%';
    else if ((i = stoi(match[1])) <= sizeof...(params))
      result += param[i];
    else
      result += match[0];

    model = match.suffix();
  }

  result += match.suffix();
  return result;
}


#include <iostream>
int main()
{
  std::cout << format(
      "%2% %3% of %1%\n"
      "%20%\n"
      "%2broken%\n"
      "%2%% of %1%, percent-sign: %%\n"
      "not closed %\n",
    "acid",
    37.5,
    "ounces",
    "not used?");
  return 0;
}

// vt_type_list.cpp
#include <type_traits>
#include <cstddef>

template <class... Types>
struct Type_list
{
  static constexpr auto size = sizeof...(Types);
};


#define APPLY(meta, ...) \
  typename meta::template apply<__VA_ARGS__>::type


struct T_prepend
{
  template <class A, class B>
  struct apply {};

  template <class A, class... Other>
  struct apply<A, Type_list<Other...>>
  {
    using type = Type_list<A, Other...>;
  };
};


struct T_map
{
  template <class Map, class List>
  struct apply { using type = List; };

  template <class Map, class T, class... Tail>
  struct apply<Map, Type_list<T, Tail...>>
  {
    using type =
    APPLY(T_prepend, \
      APPLY(Map, T), \
      APPLY(T_map, Map, Type_list<Tail...>));
  };
};


struct T_fold
{
  template <class Fold, class List> struct apply {};

  template <class Fold, class T>
  struct apply<Fold, Type_list<T>>
  {
    using type = T;
  };

  template <class Fold, class T, class... Tail>
  struct apply<Fold, Type_list<T, Tail...>>
  {
    using type = APPLY(Fold, T, APPLY(T_fold, Fold, Type_list<Tail...>));
  };
};


struct T_sizeof
{
  template <class T>
  struct apply { using type = std::integral_constant<size_t, sizeof(T)>; };
};


struct T_add
{
  template <class A, class B>
  struct apply
  {
    using type = std::integral_constant<
      std::common_type_t<typename A::value_type, typename B::value_type>,
      A::value + B::value>;
  };
};

struct T_max
{
  template <class A, class B>
  struct apply
  {
    using type = std::integral_constant<
      std::common_type_t<typename A::value_type, typename B::value_type>,
      (A::value < B::value? B::value: A::value)>;
  };
};


#include <iostream>

template <class BinOp, class Types>
using F = APPLY(T_fold, \
  BinOp, \
  APPLY(T_map, T_sizeof, Types));

int main()
{
  using Ts = Type_list<int, char, bool, char>;
  std::cout
    << "Type list has " << Ts::size << " types in it\n"
    << "Type list max sizeof = "
    << F<T_max, Ts>{}
    << "\nType list sizeof sum = "
    << F<T_add, Ts>{} << '\n';

  return 0;
}

// vt_variant.cpp
#include <algorithm>
#include <cstdint>
#include <typeinfo>

template <class... Types>
struct Type_list
{
  static constexpr auto size = sizeof...(Types);
};

template <class T, class L>
constexpr size_t t_count = 0;

template <class T, class... Tail>
constexpr size_t t_count<T, Type_list<T, Tail...>> = 1 + t_count<T, Type_list<Tail...>>;

template <class T, class X, class... Tail>
constexpr size_t t_count<T, Type_list<X, Tail...>> = t_count<T, Type_list<Tail...>>;


template <class T, class L>
constexpr size_t t_index = -1;

template <class T, class... Tail>
constexpr size_t t_index<T, Type_list<T, Tail...>> = 0;

template <class T, class X, class... Tail>
constexpr size_t t_index<T, Type_list<X, Tail...>> = 1 + t_index<T, Type_list<Tail...>>;


template <size_t i, class L>
struct t_at {};

template <size_t i, class T, class... Tail>
struct t_at<i, Type_list<T, Tail...>>
{
  using type =
    typename t_at<i - 1, Type_list<Tail...>>
    ::type;
};

template <class T, class... Tail>
struct t_at<0, Type_list<T, Tail...>>
{
  using type = T;
};

template <size_t i, class L>
using t_at_t = typename t_at<i, L>::type;


template <class L>
constexpr size_t max_sizeof = 0;

template <class T, class... Tail>
constexpr size_t max_sizeof<Type_list<T, Tail...>> =
  std::max(sizeof(T), max_sizeof<Type_list<Tail...>>);

template <class L>
constexpr size_t max_alignof = 1;

template <class T, class... Tail>
constexpr size_t max_alignof<Type_list<T, Tail...>> =
  std::max(alignof(T), max_alignof<Type_list<Tail...>>);


template <class L1, class L2>
constexpr bool t_subset = true; // for empty list L1.

template <class T, class... Types1, class... Types2>
constexpr bool t_subset<
    Type_list<T, Types1...>,
    Type_list<Types2...>> =
  t_count<T, Type_list<Types2...>> != 0 &&
  t_subset<Type_list<Types1...>, Type_list<Types2...>>;


constexpr size_t ti_prop1(size_t a)
{
  return a != size_t(-1)? a + 1: -1;
}

template <class... Args>
constexpr size_t t_first_constructible = -1;

template <class T, class... Types, class... Args>
constexpr size_t t_first_constructible<
    Type_list<T, Types...>, Args...> =
  std::is_constructible<T, Args...>::value? 0
  : ti_prop1(t_first_constructible<Type_list<Types...>, Args...>);


/// Вариантный тип.
template <class... Types>
class Variant
{
public:
  static_assert(sizeof...(Types) > 0,
    "Variant: at least one type needed.");
  static_assert(sizeof...(Types) < 64,
    "Variant: too many possible types.");
  using types = Type_list<Types...>;

private:
  using Storage = std::aligned_storage_t<
    max_sizeof<types>, max_alignof<types>>;

  Storage _stor;
  uint8_t _type = -1;


  template <class T>
  static void dtor(void *p)
  {
    static_cast<T*>(p)->~T();
  }

  static auto fetch_dtor(uint8_t t)
  {
    static void (*d[])(void*)
    { &dtor<Types>... }; // Destructor table.
    return d[t];
  }


  template <class T>
  static void copy_ctor(void *p, const void *val)
  {
    new (p) T(*static_cast<const T*>(val));
  }

  static auto fetch_copy_ctor(uint8_t t)
  {
    static void (*cc[])(void*, const void*)
    { &copy_ctor<Types>... }; // Copy constructor table.
    return cc[t];
  }


  template <class T, class Visitor>
  static void do_visit(Visitor &&visitor, void *val)
  {
    visitor(*static_cast<T*>(val));
  }

  template <class Visitor>
  static auto fetch_do_visit(uint8_t t)
  {
    using DVT = void(Visitor&&, void*);
    static DVT* const vt[]
    { &do_visit<Types, Visitor>... };
    return vt[t];
  }

public:
  /// Извлечь индекс хранимого типа.
  uint8_t get_type() const noexcept { return _type; }

  /// Получить адрес хранилища.
  const void* data() const noexcept { return &_stor; }

  /// Проверить, "пуст" ли вариант (содержит ли он какое-либо значение).
  bool empty() const noexcept { return _type == -1; }

  /// Очистить содержимое, вернув к состоянию "пусто".
  void clear()
  {
    if (!empty())
    {
      // Dispatch destructor call:
      fetch_dtor(_type)(&_stor);
      _type = -1;
    }
  }

  // Деструктор, вызывает clear().
  ~Variant() { clear(); }

  // Конструктор из произвольного набора аргументов.
  // Выбирает первый подходящий тип из Types (слева направо) для инициализации.
  template <class... Args>
  Variant(Args&&... args)
  {
    static constexpr size_t i = t_first_constructible<types, Args&&...>;
    static_assert(i != -1, "Variant constructor: invalid arguments.");
    new (static_cast<void*>(&_stor))
      t_at_t<i, types>(std::forward<Args>(args)...);
    _type = uint8_t(i);
  }

  // Конструктор по умолчанию, создаёт пустой объект варианта.
  Variant() {}

  // Копирующий конструктор.
  Variant(const Variant &other) { *this = other; }

  // Копирующий конструктор (вариант для не const-ссылки).
  Variant(Variant &other) { *this = other; }
  
  // Перемещающий конструктор.
  // (По сути, не реализован -- вызывает копирование.)
  Variant(Variant&& other) { *this = other; }

  // Копирующий оператор присваивания.
  Variant& operator=(const Variant &other)
  {
    if (this != &other)
    {
      clear();
      // Dispatch copy constructor call:
      fetch_copy_ctor(other._type)(&_stor, &other._stor);
      _type = other._type;
    }
    return *this;
  }

  // Присваивание варианта с более узким множеством возможных типов значений.
  template <class... Subset>
  std::enable_if_t<
    t_subset<Type_list<Subset...>, types>,
  Variant<Types...>&> operator=(const Variant<Subset...> &other)
  {
    clear();
    if (other.empty())
      return *this;

    // Таблица перевода индексов other._type в this->_type.
    static const uint8_t tt[]
    { t_index<Subset, types>... };

    fetch_copy_ctor(tt[other.get_type()])(&_stor, other.data());
    _type = tt[other.get_type()];
    return *this;
  }

  // Оператор присваивания произвольного значения.
  // Доступен только если тип T есть в Types.
  template <class T>
  std::enable_if_t<
    t_count<T, types> != 0,
  Variant<Types...>&> operator=(const T &val)
  {
    return *this = Variant<Types...>(val); // possibly suboptimal.
  }

  // Получить значение type_info для хранимого типа.
  const std::type_info& get_typeid() const
  {
    if (empty())
      return typeid(void);

    static const std::type_info* const ti[]
    { &typeid(Types)... };

    return *ti[_type];
  }

  // Обработать значение варианта с помощью паттерна "посетитель".
  template <class Visitor>
  void visit(Visitor &&visitor)
  {
    fetch_do_visit<Visitor>(_type)
      (std::forward<Visitor>(visitor), &_stor);
  }
};


#include <iostream>
#include <string>

using namespace std;

struct Twice
{
  template <class T>
  void operator()(T &x) const
  {
    x += x;
  }
} const twice;

struct To_cout
{
  template <class T>
  ostream& operator()(const T &x) const
  {
    return cout << x;
  }
} const to_cout;


int main()
{
  Variant<int, string> v = 3;
  v.visit(twice);
  v.visit(to_cout);
  cout << ' ' << v.get_typeid().name();
  cout << '\n';

  v = "da";
  v.visit(twice);
  v.visit(to_cout);
  cout << ' ' << v.get_typeid().name();
  cout << '\n';

  return 0;
}