#include <iostream>
#include <iomanip>
#include <cmath>
#include <algorithm>
#include <limits>
#include <cstdlib>
#include <ctime>
#include <string>

template<typename T>
bool equal_naive (T a, T b) {
  return a == b;
}

template<typename T>
bool equal_abs (T a, T b, T eps) {
  return std::fabs (a - b) <= eps;
}

template<typename T>
bool equal_rel (T a, T b, T eps) {
  const T scale = std::max (std::fabs (a), std::fabs (b));
  return std::fabs (a - b) <= eps * scale;
}

template<typename T>
bool equal_combined (T a, T b, T abs_eps, T rel_eps) {
  const T scale = std::max (std::fabs (a), std::fabs (b));
  return std::fabs (a - b) <= std::max (abs_eps, rel_eps * scale);
}

template<typename T>
void print_bool_result (const std::string &label, bool value) {
  std::cout << "  " << std::left << std::setw (28) << label
            << ": " << (value ? "true" : "false") << "\n";
}

template<typename T>
void print_value_line (const std::string &label, T value) {
  std::cout << "  " << std::left << std::setw (28) << label
            << ": " << std::setprecision (std::numeric_limits<T>::digits10 + 2)
            << value << "\n";
}

void print_section (const std::string &title) {
  std::cout << "\n============================================================\n";
  std::cout << title << "\n";
  std::cout << "============================================================\n";
}

void example_basic_01_plus_02() {
  print_section ("1. Basic classic: 0.1 + 0.2 != 0.3");

  const double a = 0.1 + 0.2;
  const double b = 0.3;

  print_value_line ("a", a);
  print_value_line ("b", b);
  print_value_line ("abs(a - b)", std::fabs (a - b));

  print_bool_result<double> ("a == b", equal_naive (a, b));
  print_bool_result<double> ("abs eps = 1e-12", equal_abs (a, b, 1e-12));
  print_bool_result<double> ("abs eps = 1e-9", equal_abs (a, b, 1e-9));
  print_bool_result<double> ("combined", equal_combined (a, b, 1e-12, 1e-12));
}

void example_large_numbers() {
  print_section ("2. Large numbers: fixed absolute epsilon becomes useless");

  const double a = 1e9;
  const double b = a + 0.1;

  print_value_line ("a", a);
  print_value_line ("b", b);
  print_value_line ("abs(a - b)", std::fabs (a - b));

  print_bool_result<double> ("a == b", equal_naive (a, b));
  print_bool_result<double> ("abs eps = 1e-12", equal_abs (a, b, 1e-12));
  print_bool_result<double> ("abs eps = 1e-9", equal_abs (a, b, 1e-9));
  print_bool_result<double> ("abs eps = 1e-3", equal_abs (a, b, 1e-3));

  print_bool_result<double> ("rel eps = 1e-12", equal_rel (a, b, 1e-12));
  print_bool_result<double> ("rel eps = 1e-9", equal_rel (a, b, 1e-9));
  print_bool_result<double> ("combined", equal_combined (a, b, 1e-9, 1e-9));
}

void example_small_numbers() {
  print_section ("3. Very small numbers: fixed absolute epsilon can be too large");

  const double a = 1e-12;
  const double b = 2e-12;

  print_value_line ("a", a);
  print_value_line ("b", b);
  print_value_line ("abs(a - b)", std::fabs (a - b));

  print_bool_result<double> ("a == b", equal_naive (a, b));
  print_bool_result<double> ("abs eps = 1e-9", equal_abs (a, b, 1e-9));
  print_bool_result<double> ("abs eps = 1e-12", equal_abs (a, b, 1e-12));
  print_bool_result<double> ("rel eps = 1e-9", equal_rel (a, b, 1e-9));
  print_bool_result<double> ("rel eps = 0.5", equal_rel (a, b, 0.5));
  print_bool_result<double> ("combined", equal_combined (a, b, 1e-15, 1e-6));
}

void example_near_zero_relative_problem() {
  print_section ("4. Near zero: relative comparison alone is weak");

  const double a = 0.0;
  const double b = 1e-15;

  print_value_line ("a", a);
  print_value_line ("b", b);
  print_value_line ("abs(a - b)", std::fabs (a - b));

  print_bool_result<double> ("a == b", equal_naive (a, b));
  print_bool_result<double> ("rel eps = 1e-6", equal_rel (a, b, 1e-6));
  print_bool_result<double> ("abs eps = 1e-12", equal_abs (a, b, 1e-12));
  print_bool_result<double> ("combined", equal_combined (a, b, 1e-12, 1e-6));
}

void example_associativity() {
  print_section ("5. Same math on paper, different result in floating point");

  const double a = 1e16;
  const double b = -1e16;
  const double c = 1.0;

  const double x = (a + b) + c;
  const double y = a + (b + c);

  print_value_line ("x = (a + b) + c", x);
  print_value_line ("y = a + (b + c)", y);
  print_value_line ("abs(x - y)", std::fabs (x - y));

  print_bool_result<double> ("x == y", equal_naive (x, y));
  print_bool_result<double> ("combined", equal_combined (x, y, 1e-12, 1e-12));
}

void example_accumulation() {
  print_section ("6. Accumulation of error: repeated addition");

  double x = 0.0;
  for (int i = 0; i < 1000000; ++i)
    x += 0.1;

  const double y = 100000.0;

  print_value_line ("x", x);
  print_value_line ("y", y);
  print_value_line ("abs(x - y)", std::fabs (x - y));

  print_bool_result<double> ("x == y", equal_naive (x, y));
  print_bool_result<double> ("abs eps = 1e-9", equal_abs (x, y, 1e-9));
  print_bool_result<double> ("abs eps = 1e-6", equal_abs (x, y, 1e-6));
  print_bool_result<double> ("combined", equal_combined (x, y, 1e-6, 1e-12));
}

void example_float_vs_double() {
  print_section ("7. float vs double: same idea, different precision");

  float af = 0.1f + 0.2f;
  float bf = 0.3f;

  double ad = 0.1 + 0.2;
  double bd = 0.3;

  print_value_line ("float a", af);
  print_value_line ("float b", bf);
  print_value_line ("float abs diff", std::fabs (af - bf));
  print_bool_result<float> ("float a == b", equal_naive (af, bf));

  std::cout << "\n";

  print_value_line ("double a", ad);
  print_value_line ("double b", bd);
  print_value_line ("double abs diff", std::fabs (ad - bd));
  print_bool_result<double> ("double a == b", equal_naive (ad, bd));
}

double random_noise (double amplitude) {
  const double unit = static_cast<double> (std::rand()) / static_cast<double> (RAND_MAX);
  return (unit * 2.0 - 1.0) * amplitude;
}

void run_sensor_simulation (double noise_amplitude, double abs_eps) {
  const double target = 25.0;
  const double rel_eps = 1e-6;

  std::cout << "\nNoise amplitude: +/-" << noise_amplitude
            << ", abs_eps: " << abs_eps << "\n";

  std::cout << "  " << std::left
            << std::setw (10) << "sample"
            << std::setw (18) << "value"
            << std::setw (10) << "=="
            << std::setw (10) << "abs"
            << std::setw (10) << "comb"
            << "\n";

  for (int i = 0; i < 12; ++i) {
    const double value = target + random_noise (noise_amplitude);

    const bool naive = equal_naive (value, target);
    const bool abs_ok = equal_abs (value, target, abs_eps);
    const bool comb_ok = equal_combined (value, target, abs_eps, rel_eps);

    std::cout << "  " << std::left
              << std::setw (10) << i
              << std::setw (18) << std::setprecision (10) << value
              << std::setw (10) << (naive ? "true" : "false")
              << std::setw (10) << (abs_ok ? "true" : "false")
              << std::setw (10) << (comb_ok ? "true" : "false")
              << "\n";
  }
}

void example_sensor_noise() {
  print_section ("8. Sensor-like values: 'equal' is often the wrong question");

  run_sensor_simulation (0.05, 0.1);
  run_sensor_simulation (0.005, 0.01);
  run_sensor_simulation (0.0005, 0.0001);
}

void example_scaled_integer_vs_double() {
  print_section ("9. Scaled integer values: often better to compare raw domain");

  const int raw_value = 1234;     // imagine 123.4 in deci-units
  const int raw_target = 1234;

  const double scaled_value = raw_value * 0.1;
  const double scaled_target = 123.4;

  print_value_line ("raw_value", raw_value);
  print_value_line ("raw_target", raw_target);
  print_bool_result<int> ("raw_value == raw_target", raw_value == raw_target);

  std::cout << "\n";

  print_value_line ("scaled_value", scaled_value);
  print_value_line ("scaled_target", scaled_target);
  print_value_line ("abs diff", std::fabs (scaled_value - scaled_target));
  print_bool_result<double> ("scaled_value == scaled_target",
                             equal_naive (scaled_value, scaled_target));
  print_bool_result<double> ("combined",
                             equal_combined (scaled_value, scaled_target, 1e-12, 1e-12));

  std::cout << "\n  Note: if your system is naturally discrete, comparing raw units can be simpler\n";
  std::cout << "  and more honest than converting everything to floating point too early.\n";
}

void example_rounding_is_not_magic() {
  print_section ("10. Rounding is not a universal fix");

  const double a = 1.0049;
  const double b = 1.0051;

  const double rounded_a_2 = std::round (a * 100.0) / 100.0;
  const double rounded_b_2 = std::round (b * 100.0) / 100.0;

  const double rounded_a_3 = std::round (a * 1000.0) / 1000.0;
  const double rounded_b_3 = std::round (b * 1000.0) / 1000.0;

  print_value_line ("a", a);
  print_value_line ("b", b);

  std::cout << "\n";
  print_value_line ("round(a, 2)", rounded_a_2);
  print_value_line ("round(b, 2)", rounded_b_2);
  print_bool_result<double> ("equal after round(2)", rounded_a_2 == rounded_b_2);

  std::cout << "\n";
  print_value_line ("round(a, 3)", rounded_a_3);
  print_value_line ("round(b, 3)", rounded_b_3);
  print_bool_result<double> ("equal after round(3)", rounded_a_3 == rounded_b_3);

  std::cout << "\n  Rounding may hide differences or invent equality depending on chosen precision.\n";
}

void example_hysteresis() {
  print_section ("11. Hysteresis: in real systems you often want state logic, not equality");

  const double target = 25.0;
  const double on_threshold = target - 0.1;
  const double off_threshold = target + 0.1;

  const double signal[] = {
    24.92, 24.97, 25.01, 25.05, 24.99, 25.08, 25.11, 25.06,
    24.98, 24.91, 24.89, 24.95, 25.02, 25.12, 25.07, 24.93
  };

  bool heater_simple = true;
  bool heater_hysteresis = true;

  std::cout << "  " << std::left
            << std::setw (8)  << "step"
            << std::setw (12) << "value"
            << std::setw (16) << "simple_ctrl"
            << std::setw (16) << "hyst_ctrl"
            << "\n";

  for (size_t i = 0; i < (sizeof (signal) / sizeof (signal[0])); ++i) {
    const double value = signal[i];

    /*
     * Simple control:
     *   heater ON below target
     *   heater OFF at or above target
     * This tends to chatter near the threshold.
     */
    heater_simple = value < target;

    /*
     * Hysteresis control:
     *   heater turns ON only below lower threshold
     *   heater turns OFF only above upper threshold
     * This reduces chatter around the target.
     */
    if (heater_hysteresis && value > off_threshold)
      heater_hysteresis = false;
    else if (!heater_hysteresis && value < on_threshold)
      heater_hysteresis = true;

    std::cout << "  " << std::left
              << std::setw (8)  << i
              << std::setw (12) << std::setprecision (6) << value
              << std::setw (16) << (heater_simple ? "HEATER ON" : "HEATER OFF")
              << std::setw (16) << (heater_hysteresis ? "HEATER ON" : "HEATER OFF")
              << "\n";
  }

  std::cout << "\n  This is the key engineering point:\n";
  std::cout << "  sometimes the answer is not a better equality test,\n";
  std::cout << "  but a better control model.\n";
}

int main() {
  std::srand (static_cast<unsigned int> (std::time (NULL)));

  std::cout << std::boolalpha;

  example_basic_01_plus_02();
  example_large_numbers();
  example_small_numbers();
  example_near_zero_relative_problem();
  example_associativity();
  example_accumulation();
  example_float_vs_double();
  example_sensor_noise();
  example_scaled_integer_vs_double();
  example_rounding_is_not_magic();
  example_hysteresis();

  std::cout << "\nDone.\n";
  return 0;
}