#include <iostream>
#include <vector>
#include <cmath>

double eps, a, b, tab_step;

std::vector<double> root_segs;

double f(double x) {
  return sqrt(4.0*x+7.0) - 3.0 * cos(x);
}

double fprime(double x) {
  return 3 * sin(x) + 2.0/sqrt(4*x+7);
}

void in_par() {
  //std::cin >> tab_step >> eps;
  //std::cin >> a >> b;
  a = -1.5; b = 2;
  eps = 1e-8; tab_step = 0.01;
}

void print_seg(int i) {
  std::cout << "[" << root_segs[i] << ", " << (root_segs[i]+tab_step) << "]";
}

void sep() {
  for (double x = a; x <= b; x += tab_step)
    if (f(x) * f(x + tab_step) <= 0)
      root_segs.push_back(x);
  if (!root_segs.size()) {
    std::cout << "No roots present in requested segment, terminating\n";
    exit(0);
  }
  std::cout << "Root separation yielded following segments:" << std::endl;
  print_seg(0);
  for (int i = 1; i < root_segs.size(); i++)
    std::cout << ", ", print_seg(i);
  std::cout << std::endl;
}

void bisect(double l, double r) {
  std::cout << "Bisection method; [" << l << ", " << r << "]: ";
  int step = 0;
  while (r - l > 2 * eps) {
    ++step;
    double mid = (l + r) / 2.0;
    if (f(mid) * f(r) < 0)
      l = mid;
    else
      r = mid;
  }
  double ans = (r+l)/2.0;
  std::cout
    << step << " steps, ans " << ans
    << ", |r| = " << fabs(f(ans))
    << ", fsl = " << (r-l)
    << std::endl;
}

void newton(double xs) {
  std::cout << "Newton's method; x_0 = " << xs << "; ";
  int steps = 0;
  double nx = xs, px = xs;
  do {
    ++steps;
    px = nx;
    nx = nx - f(nx) / fprime(nx);
  } while (nx - px > eps);
  std::cout
    << steps << " steps, ans " << nx
    << ", |r| = " << fabs(f(nx))
    << std::endl;
}

void mod_newton(double xs) {
  std::cout << "Modified Newton's method; x_0 = " << xs << "; ";
  int steps = 0;
  double nx = xs, px = xs;
  double fpr0 = fprime(xs);
  do {
    ++steps;
    px = nx;
    nx = nx - f(nx) / fpr0;
  } while (nx - px > eps);
  std::cout
    << steps << " steps, ans " << nx
    << ", |r| = " << fabs(f(nx))
    << std::endl;
}

void secant(double k, double l) {
  std::cout << "Secant method; x_0 = " << k << ", x_1 = " << l;
  int steps = 0;
  while (fabs(k - l) > eps) {
    ++steps;
    double fs = f(l);
    double n = l - fs * ((l - k) / (fs - f(k)));
    k = l; l = n;
  }
  std::cout
    << ", " << steps
    << " steps, ans " << l
    << ", |r| = " << fabs(f(l))
    << std::endl;
}

int main(int argc, char** argv) {
  std::cout << "LW 1: numerical methods of solving nonlinear equations\n";
  std::cout << "Variant 4, f(x) := sqrt(4*x + 7) - 3*cos(x)\n";

  std::cout.precision(10);
  in_par();

  std::cout
    << "A = " << a << "; B = " << b
    << "; h = " << tab_step
    << "; eps = " << eps << std::endl;

  std::cout << std::endl;
  sep();
  std::cout << std::endl;
  for (int i = 0; i < root_segs.size(); i++)
    bisect(root_segs[i], root_segs[i] + tab_step);
  for (int i = 0; i < root_segs.size(); i++)
    newton(root_segs[i]);
  for (int i = 0; i < root_segs.size(); i++)
    mod_newton(root_segs[i]);
  for (int i = 0; i < root_segs.size(); i++)
    secant(root_segs[i], root_segs[i] + tab_step);

  return 0;
}