Во-первых, необходимо<#includes>.

#include <boost/math/tools/roots.hpp>
//using boost::math::policies::policy;
//using boost::math::tools::newton_raphson_iterate;
//using boost::math::tools::halley_iterate; //
//using boost::math::tools::eps_tolerance; // Binary functor for specified number of bits.
//using boost::math::tools::bracket_and_solve_root;
//using boost::math::tools::toms748_solve;
#include <boost/math/special_functions/next.hpp> // For float_distance.
#include <tuple> // for std::tuple and std::make_tuple.
#include <boost/math/special_functions/cbrt.hpp> // For boost::math::cbrt.

	Tip
	Для наглядности<using>приведен список функций, используемых в этом примере: вы можете, конечно, частично или полностью квалифицировать имена другими способами. (Для вашего приложения вы, возможно, захотите извлечь некоторые части в файлы заголовков, но вы никогда не должны использовать<using>заявления по всему миру в файлах заголовков).

Предположим, что мы хотим найти корень числаaи начать вычислять кубический корень.

Уравнение, которое мы хотим решить:

  f(x) = x³ -a

Сначала мы решим эту проблему, не используя никакой информации о наклоне или кривизне корневой функции куба.

К счастью, функция кубического корня «действительно хорошо себя ведет» в том смысле, что она монотонна и имеет только один корень (мы оставляем отрицательные значения «как упражнение для ученика»).

Затем мы покажем, как добавление того, что мы можем знать об этой функции, сначала только наклон или 1-я производнаяf'(x), ускорит наведение на решение.

Наконец, мы показываем, как добавление кривизныf''(x)тоже ускорит сходимость еще больше.

Cube root function without derivatives

Для начала определим объект функции (функтора):

template <class T>
struct cbrt_functor_noderiv
{
  //  cube root of x using only function - no derivatives.
  cbrt_functor_noderiv(T const& to_find_root_of) : a(to_find_root_of)
  { /* Constructor just stores value a to find root of. */ }
  T operator()(T const& x)
  {
    T fx = x*x*x - a; // Difference (estimate x^3 - a).
    return fx;
  }
private:
  T a; // to be 'cube_rooted'.
};

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

template <class T>
T cbrt_noderiv(T x)
{
  // return cube root of x using bracket_and_solve (no derivatives).
  using namespace std;                          // Help ADL of std functions.
  using namespace boost::math::tools;           // For bracket_and_solve_root.
  int exponent;
  frexp(x, &exponent);                          // Get exponent of z (ignore mantissa).
  T guess = ldexp(1., exponent/3);              // Rough guess is to divide the exponent by three.
  T factor = 2;                                 // How big steps to take when searching.
  const boost::uintmax_t maxit = 20;            // Limit to maximum iterations.
  boost::uintmax_t it = maxit;                  // Initally our chosen max iterations, but updated with actual.
  bool is_rising = true;                        // So if result if guess^3 is too low, then try increasing guess.
  int digits = std::numeric_limits<T>::digits;  // Maximum possible binary digits accuracy for type T.
  // Some fraction of digits is used to control how accurate to try to make the result.
  int get_digits = digits - 3;                  // We have to have a non-zero interval at each step, so
                                                // maximum accuracy is digits - 1.  But we also have to
                                                // allow for inaccuracy in f(x), otherwise the last few
                                                // iterations just thrash around.
  eps_tolerance<T> tol(get_digits);             // Set the tolerance.
  std::pair<T, T> r = bracket_and_solve_root(cbrt_functor_noderiv<T>(x), guess, factor, is_rising, tol, it);
  return r.first + (r.second - r.first)/2;      // Midway between brackets is our result, if necessary we could
                                                // return the result as an interval here.
}

Note

Последний параметр, определяющий максимальное количество итераций, является необязательным. Но это не то, что<boost::uintmax_t maxit= (std::numeric_limits<boost::uintmax_t>::max)();>, а<18446744073709551615>, и это больше, чем кто-либо хотел бы ждать! Таким образом, может быть разумно выбрать некоторую разумную оценку того, сколько итераций может потребоваться, в этом случае функция настолько хорошо ведет себя, что мы можем выбрать низкое значение 20. Внутреннее повышение. Математика использует эти функции, она устанавливает максимальные итерации к<policies::get_max_root_iterations<Policy>();>.

Если бы мы захотели, мы могли бы показать, сколько итераций было использовано в<bracket_and_solve_root>(эта информация потеряна за пределами<cbrt_noderiv>), например:

if (it >= maxit)
{
  std::cout << "Unable to locate solution in " << maxit << " iterations:"
    " Current best guess is between " << r.first << " and " << r.second << std::endl;
}
else
{
  std::cout << "Converged after " << it << " (from maximum of " << maxit << " iterations)." << std::endl;
}

для выхода как

Converged after 11 (from maximum of 20 iterations).

Этот фрагмент из<main()>вroot_finding_example.cppпоказывает, как его можно использовать.

try
{
  double threecubed = 27.;   // Value that has an *exactly representable* integer cube root.
  double threecubedp1 = 28.; // Value whose cube root is *not* exactly representable.
  std::cout << "cbrt(28) " << boost::math::cbrt(28.) << std::endl; // boost::math:: version of cbrt.
  std::cout << "std::cbrt(28) " << std::cbrt(28.) << std::endl;    // std:: version of cbrt.
  std::cout <<" cast double " << static_cast<double>(3.0365889718756625194208095785056696355814539772481111) << std::endl;
  // Cube root using bracketing:
  double r = cbrt_noderiv(threecubed);
  std::cout << "cbrt_noderiv(" << threecubed << ") = " << r << std::endl;
  r = cbrt_noderiv(threecubedp1);
  std::cout << "cbrt_noderiv(" << threecubedp1 << ") = " << r << std::endl;

  cbrt_noderiv(27) = 3
  cbrt_noderiv(28) = 3.0365889718756618

Результатом<bracket_and_solve_root>являетсяпаразначений, которые могут быть отображены.

Количество разделяющих их битов можно найти с помощью<float_distance(r.first,r.second)>. Расстояние равно нулю для 3³= 27, но<float_distance(r.first,r.second)=3>для кубического корня 28 с этой функцией. Результат (избегание переполнения) находится на полпути между этими двумя значениями.

Cube root function with 1st derivative (slope)

Теперь мы решаем ту же проблему, но используя больше информации о функции, чтобы показать, как это может ускорить поиск лучшей оценки корня.

Для корневой функции известен 1-й дифференциал (наклон касательной к кривой в любой точке).

Этот алгоритм похож на этотn-й корневой алгоритм.

Если вам нужны какие-то напоминания, топроизводные элементарных функциймогут помочь.

Используя правило, что производнаяxⁿдля положительного n (фактически все ненулевые n) являетсяn x^n-1, позволяет получить 1-й дифференциал как3x².

Чтобы увидеть, как эта дополнительная информация используется для поиска корня, просмотритеитерации Ньютона-Рафсонаианимацию.

Мы определяем лучший функтор<cbrt_functor_deriv>, который возвращает как оценку функции для решения, так и ее первую производную:

Длявозвратадвух значений мы используемstd:: паразначений с плавающей точкой.

template <class T>
struct cbrt_functor_deriv
{ // Functor also returning 1st derivative.
  cbrt_functor_deriv(T const& to_find_root_of) : a(to_find_root_of)
  { // Constructor stores value a to find root of,
    // for example: calling cbrt_functor_deriv<T>(a) to use to get cube root of a.
  }
  std::pair<T, T> operator()(T const& x)
  {
    // Return both f(x) and f'(x).
    T fx = x*x*x - a;                // Difference (estimate x^3 - value).
    T dx =  3 * x*x;                 // 1st derivative = 3x^2.
    return std::make_pair(fx, dx);   // 'return' both fx and dx.
  }
private:
  T a;                               // Store value to be 'cube_rooted'.
};

Корневая функция куба теперь:

template <class T>
T cbrt_deriv(T x)
{
  // return cube root of x using 1st derivative and Newton_Raphson.
  using namespace boost::math::tools;
  int exponent;
  frexp(x, &exponent);                                // Get exponent of z (ignore mantissa).
  T guess = ldexp(1., exponent/3);                    // Rough guess is to divide the exponent by three.
  T min = ldexp(0.5, exponent/3);                     // Minimum possible value is half our guess.
  T max = ldexp(2., exponent/3);                      // Maximum possible value is twice our guess.
  const int digits = std::numeric_limits<T>::digits;  // Maximum possible binary digits accuracy for type T.
  int get_digits = static_cast<int>(digits * 0.6);    // Accuracy doubles with each step, so stop when we have
                                                      // just over half the digits correct.
  const boost::uintmax_t maxit = 20;
  boost::uintmax_t it = maxit;
  T result = newton_raphson_iterate(cbrt_functor_deriv<T>(x), guess, min, max, get_digits, it);
  return result;
}

Результатом функции<newton_raphson_iterate>является единственное значение.

Tip

Существует компромисс между точностью и скоростью при выборе значения<digits>. Заманчиво просто выбрать<std::numeric_limits<T>::digits>, но это может означать некоторые неэффективные и ненужные итерации, поскольку функция трэширует, пытаясь найти последний бит. Теоретически, поскольку точность удваивается с каждым шагом, достаточно остановиться, когда половина битов верна: так как последний шаг удвоит это до полной точности. Конечно, функция не может сказать, действительно ли это так, если она не сделает еще один шаг, чтобы быть уверенным. На практике установка точности чуть больше<std::numeric_limits<T>::digits/ 2>является хорошим выбором.

Обратите внимание, что абонент функции должен проверить количество итераций после вызова, чтобы увидеть, прекратилась ли итерация в результате исчерпания итераций, а не соответствия требуемой точности.

Используя тестовые данные в/test/test_cbrt.cpp, мы обнаружили, что корень куба точен до последней цифры в каждом случае и не более чем в 6 итерациях с двойной точностью. Тем не менее, вы заметите, что в этом примере использовалась высокая точность, именно то, что было предупреждено ранее в этих документах! В этом конкретном случае можно вычислитьf(x)точно и без неоправданной ошибки отмены, поэтому высокий предел не является слишком большой проблемой.

Однако снижение предела до<std::numeric_limits<T>::digits *2/3>дало полную точность во всех, кроме одного из тестовых случаев (и что один был выключен всего на один бит). Максимальное количество итераций оставалось 6, но в большинстве случаев сокращалось на одну.

Обратите внимание, что вышеупомянутый код опускает вероятную оптимизацию путем вычисления z² и при повторном использовании он опускает обработку ошибок и не обрабатывает отрицательные значения z правильно. (Это обычное упражнение для читателя!)

Функция<boost::math::cbrt>также включает в себя эти и другие улучшения: самое главное, она использует гораздо лучшее первоначальное предположение, которое уменьшает количество итераций до 1 почти во всех случаях.

Cube root with 1st & 2nd derivative (slope & curvature)

Далее мы определяем еще один, еще лучший функтор<cbrt_functor_2deriv>, который возвращает как оценку функции для решения, так и ее первуюи вторуюпроизводную:

  f''(x) = 6x

использование информации о наклоне и кривизне для ускорения конвергенции.

Длявозврататрех значений мы используем<tuple>трех значений с плавающей точкой:

template <class T>
struct cbrt_functor_2deriv
{
  // Functor returning both 1st and 2nd derivatives.
  cbrt_functor_2deriv(T const& to_find_root_of) : a(to_find_root_of)
  { // Constructor stores value a to find root of, for example:
    // calling cbrt_functor_2deriv<T>(x) to get cube root of x,
  }
  std::tuple<T, T, T> operator()(T const& x)
  {
    // Return both f(x) and f'(x) and f''(x).
    T fx = x*x*x - a;                     // Difference (estimate x^3 - value).
    T dx = 3 * x*x;                       // 1st derivative = 3x^2.
    T d2x = 6 * x;                        // 2nd derivative = 6x.
    return std::make_tuple(fx, dx, d2x);  // 'return' fx, dx and d2x.
  }
private:
  T a; // to be 'cube_rooted'.
};

Корневая функция куба теперь:

template <class T>
T cbrt_2deriv(T x)
{
  // return cube root of x using 1st and 2nd derivatives and Halley.
  //using namespace std;  // Help ADL of std functions.
  using namespace boost::math::tools;
  int exponent;
  frexp(x, &exponent);                                // Get exponent of z (ignore mantissa).
  T guess = ldexp(1., exponent/3);                    // Rough guess is to divide the exponent by three.
  T min = ldexp(0.5, exponent/3);                     // Minimum possible value is half our guess.
  T max = ldexp(2., exponent/3);                      // Maximum possible value is twice our guess.
  const int digits = std::numeric_limits<T>::digits;  // Maximum possible binary digits accuracy for type T.
  // digits used to control how accurate to try to make the result.
  int get_digits = static_cast<int>(digits * 0.4);    // Accuracy triples with each step, so stop when just
                                                      // over one third of the digits are correct.
  boost::uintmax_t maxit = 20;
  T result = halley_iterate(cbrt_functor_2deriv<T>(x), guess, min, max, get_digits, maxit);
  return result;
}

Функция<halley_iterate>также возвращает одно значение, и число итераций покажет, соответствует ли она критерию конвергенции, установленному<get_digits>.

Непроизводный метод дает результат

cbrt_noderiv(28) = 3.0365889718756618

с расстоянием в 3 бита между скобками, в то время как производные методы сходятся к одному значению

cbrt_2deriv(28) = 3.0365889718756627

Мы можем сравнить сповышением::матема::cbrt

cbrt(28)              = 3.0365889718756627

Обратите внимание, что итерации должны останавливаться только на половине полной точности, и все же, даже в этом случае, ни один из тестовых случаев не ошибся. Более того, максимальное количество итераций теперь составляло всего 4.

Чтобы завершить картину, мы могли бы назвать<schroder_iterate>в последнем примере, и на самом деле это не имеет никакого значения для точности или количества итераций в этом конкретном случае. Однако относительная производительность этих двух методов может варьироваться в зависимости от характераf(x)и точности, с которой может быть вычислено первоначальное предположение. Нет никаких обобщений, которые могут быть сделаны, кроме как «попробуй и посмотри».

Наконец, если бы мы позвонили<cbrt>сNTL::RR, установленной с точностью 1000 бит (около 300 десятичных цифр), то полная точность может быть получена всего с 7 итерациями. Чтобы представить это в перспективе, увеличение точности в 20 раз менее чем удвоило количество итераций. Это просто подчеркивает, что большинство итераций используются, чтобы исправить первые несколько цифр: после этого эти методы могут создавать дополнительные цифры с замечательной эффективностью.

Или, другими словами:ничто не превосходит действительно хорошую первоначальную догадку!

Полный код этого примера находится по адресуroot_finding_example.cpp,