Короткие примеры C++ кода-1 Кувшинов Д

using namespace std;

double left = a, f_left = f(left);

bool f_left_zero = abs(f_left) <= f_tol;



// Корень на левой границе исходного отрезка?

if (f_left_zero && !report(left))

return left;
while (left != b)

{

// Правая граница очередного участка.

const double right = fmin(b, left + step), f_right = f(right);

const bool f_right_zero = abs(f_right) <= f_tol;



// Корень на правой границе участка?

if (f_right_zero && !report(right))

return right;
// Есть корень внутри участка?

if (!(f_left_zero || f_right_zero) && signbit(f_left) != signbit(f_right))

{

const double root = solver(f, left, right, x_tol);

if (!report(root))

return root;

}
// Передвинуть левую границу.

left = right;

f_left = f_right;

f_left_zero = f_right_zero;

}
return b;

}

///////////////////////////////////////////////////////////////////////////////

// Тестирование.
#include

#include

using namespace std;
/// Модельная функция 1.

double poly5(double x)

{

return (x - 1.)*(x - 2.)*(x - 3.)*(x - 4.)*(x - 5.);

}
/// Модельная функция 2.

double loge(double x)

{

return log(x) - 1.;

}
/// Модельная функция 3.

double ipl(double x)

{

return x*x*x - pow(3., x);

}

/// Вывести очередной корень.

bool print_root(double root)

{

cout << root << endl;

return true;

}
int main()

{

struct { const char *msg; Unary_real_function f; double a, b, step; }

tasks[] =

{

{ "poly: ", poly5, -10.0, 10.0, 0.5 },

{ "poly (1): ", poly5, -10.5, 10.5, 0.3 },

{ "e: ", loge, 0.01, 100., 0.1 },

{ "ipl: ", ipl, -100., 100., 0.5 },

{ "sin: ", sin, 0., 30., 0.2 }

};
cout.precision(16);

for (auto &task : tasks)

{

cout << task.msg << endl;

repeated_nsolve(task.f, task.a, task.b, task.step, print_root);

cout << "---------------\n" << endl;

}
return EXIT_SUCCESS;

}

0610-global_solve.cpp

// global_solve.cpp

#include

#include

#include
/// Точность приближённого решения, используемая по умолчанию.

const double TOLERANCE = 1e-9;
/// Найти диапазон [left, right] такой, что x0 >= left и существует left <= x <= right, что f(x) = 0.

/// Предполагается, что f(x) -- непрерывная функция.

/// Минимальная ширина шага по x задаётся параметром tol.

typedef double (*Function)(double);

bool catch_root_to_the_right

(

Function f,

double &left,

double &right,

double x0 = 0.,

double tol = TOLERANCE,

double max_j = DBL_MAX

)

{

using namespace std;

max_j = fmin(max_j, 1. / DBL_EPSILON);

double f_left = f(x0), f_right;
// Найти левую границу отрезка.

if (isnan(f_left))

for (double j = 1.; j < max_j; j += 2.)

for (double i = tol * j; isfinite(i); i *= 2.)

if (isfinite(f_left = f(x0 + i)))

{

x0 += i;

goto left_found;

}
left_found:

// Уже нашли ноль?

if (f_left == 0.)

{

left = right = x0;

return true;

}
// Найти правую границу отрезка.

for (double j = 1.; j < max_j; j += 2.)

for (double i = tol * j; isfinite(i); i *= 2.)

if (isfinite(f_right = f(x0 + i)))

{

// Нашли ноль?

if (f_right == 0.)

{

left = right = x0 + i;

return true;

}
// Разные знаки на концах отрезка?

if (f_left * f_right < 0.)

{

left = x0;

right = x0 + i;

return true;

}

}
// Найти не удалось.

return false;

}

/// Алгоритм численного решения уравнения f(x) = 0 на отрезке [a, b] делением отрезка пополам.

/// Данный алгоритм является вариантом двоичного поиска.

double nsolve(Function f, double a, double b, double tol = TOLERANCE)

{

using namespace std;

assert(f != nullptr);

assert(a <= b);

assert(0. <= tol);

for (auto fa = f(a), fb = f(b);;)

{

// Проверим значения функции на концах отрезка.

if (fa == 0.)

return a;

if (fb == 0.)

return b;
// Делим отрезок пополам.

const auto mid = 0.5 * (a + b); // середина отрезка

if (mid <= a || b <= mid)

return abs(fa) < abs(fb) ? a : b;

if (b - a <= tol)

return mid;
// Выберем одну из половин в качестве уточнённого отрезка.

const auto fmid = f(mid);

if (signbit(fa) != signbit(fmid))

{

// Корень на левой половине.

b = mid;

fb = fmid;

}

else

{

assert(signbit(fb) != signbit(fmid));

// Корень на правой половине.

a = mid;

fa = fmid;

}

}

}

///////////////////////////////////////////////////////////////////////////////

// Тестирование
#include

#include
struct Test_result

{

bool root_found; // корень был найден?

double a, b, root, f_root; // диапазон, корень, значение функции в корне

};
// Для заданной функции f найти диапазон, содержащий корень f(x) = 0.

// Затем найти на этом диапазоне сам корень (приближённо) и вычислить значение функции в нём.

void test(Function f, Test_result &result)

{

result.root_found = false;

if (catch_root_to_the_right(f, result.a, result.b))

{

result.root_found = true;

result.root = nsolve(f, result.a, result.b);

result.f_root = f(result.root);

}

}
int main()

{

using namespace std;

cout.precision(9);

// Набор заданий: функции и текстовые обозначения.

struct Task {

Function f;

const char *name;

Test_result result;

} task[] =

{

{ sin, "sin" },

{ cos, "cos" },

#define F(expr, name) \

{ [](double x) { return expr; }, name }

F(tan(x - 1.), "tan1"),

F(exp(x) - 1000., "expm"),

F(20. / (x*x) - 0.1, "rep"),

F(pow(3., x) - x*x*x, "W3"),

F(exp(sin(0.1*x)) - exp(cos(sqrt(2.)*x)), "cexp"),

F(log((x + 100.) / x) - 1 / sqrt(x), "logr"),

F(log(sin(x) / x) / tan(x - 2.5)*cos(x), "f1"),

F(10. / cbrt(x) - 1 / (1. + log(x + 1.)) + 100. / sqrt(x), "f2"),

F(atan(log(abs(tan(x) - exp(x)))) - 0.5, "f3"),

F(atan(log(abs(tan(x + 0.5)*tan(sqrt(2.)*x + 0.3)))) - 1.4, "f4"),

F(tanh(-sin(1. / cosh(x - 100.)*sin(sinh(x - 100.)))) + 0.15, "chirp"),

F(exp(-1. / pow(x - 100123., 8)) - 0.7, "b100k"),

F(log(x) / exp(1. / x) + 0.09726, "touch"),

F(0.1 + cbrt(x) - pow(x, 1./5.) - pow(x, 1./7.), "357"),

F((x < 1e+6? abs(sin(x+0.1)*sin(sqrt(30.)*x+0.1)) - 0.001: 1.), "s30")

#undef F

};
// Параллельное вычисление результатов (использует OpenMP).

static const int tasks = sizeof(task) / sizeof(Task);

#pragma omp parallel for

for (int i = 0; i < tasks; ++i)

test(task[i].f, task[i].result);
// Последовательный вывод результатов.

for (auto &t: task)

{

cout << t.name << ": ";

auto &r = t.result;

if (r.root_found)

{

cout << '[' << r.a << ", " << r.b << "], root = ";

cout << r.root << ", f(root) = " << r.f_root << '\n';

}

else

{

cout << "no roots found\n";

}

}
cout << "done" << endl;

cin.ignore();

return EXIT_SUCCESS;

}

0620-isqrt.cpp

// isqrt.cpp

// Квадратный корень целого числа, округлённый вниз.

#include

#include

// Через стандартную функцию для вычисления квадратного корня в плавающей точке.

// Должно работать быстро на современных процессорах.

// "Внутри" обычно реализовано на основе метода Ньютона.

inline uint32_t isqrt_fpu(uint32_t n)

{

return uint32_t(sqrt(double(n)));

}

// Метод деления отрезка пополам.

uint32_t isqrt_dt1(uint32_t n)

{

// Границы диапазона: [0, 65535]

// (65535 в квадрате представимо в 32-битным числом, 65536 -- уже нет).

uint32_t a = 0, b = 65536;

// Цикл выполняет 16 итераций.

for (int i = 0; i < 16; ++i)

{

// Середина отрезка [a, b], с округлением вниз.

const auto m = (a + b) >> 1;

if (m * m > n)

b = m;

else

a = m;

}

return a;

}

// Метод деления отрезка пополам: на одну переменную меньше.

uint32_t isqrt_dt2(uint32_t n)

{

register uint32_t a = 0, d = 32768;

do

{

// Середина отрезка [a, b], с округлением вниз.

const auto m = a + d;

if (m * m <= n)

a = m;

} while (d >>= 1);

return a;

}

/////////////////////////////////////////////////////////////////////

// Тестирование

#include

#include