<cmath>

#include <cmath>

// All functions live in namespace std
double r = std::sqrt(2.0);
float  s = std::sqrtf(2.0f);   // C-style suffix form — still valid
float  t = std::sqrt(2.0f);    // C++11: preferred float overload

double x = std::sqrt(4);    // C++11: int promoted to double → 2.0

#include <cmath>
#include <numbers>  // C++20: std::numbers::pi
#include <print>    // C++23: std::println

struct Vec2 { double x, y; };

double magnitude(Vec2 v) {
    return std::hypot(v.x, v.y);  // avoids intermediate overflow
}

double angle_deg(Vec2 v) {
    double rad = std::atan2(v.y, v.x);
    return rad * (180.0 / std::numbers::pi);  // C++20
}

int main() {
    Vec2 v{3.0, 4.0};
    std::println("magnitude: {}", magnitude(v));   // 5.0
    std::println("angle: {}°", angle_deg(v));       // 53.13...
}

#include <cmath>
#include <limits>

bool safe_log(double x, double& result) {
    if (!std::isfinite(x) || x <= 0.0)
        return false;
    result = std::log(x);
    return true;
}

void inspect(double v) {
    switch (std::fpclassify(v)) {
        case FP_INFINITE:  /* handle inf  */ break;
        case FP_NAN:       /* handle NaN  */ break;
        case FP_SUBNORMAL: /* handle denormal */ break;
        case FP_ZERO:      /* handle zero */ break;
        case FP_NORMAL:    /* ordinary value  */ break;
    }
}

#include <cmath>

// C++11: all rounding variants exist
double v = 2.7;

std::floor(v);      // 2.0  — toward -∞
std::ceil(v);       // 3.0  — toward +∞
std::trunc(v);      // 2.0  — toward zero (C++11)
std::round(v);      // 3.0  — nearest, halfway away from zero
std::nearbyint(v);  // 3.0  — nearest, respects current rounding mode, no FE_INEXACT
std::rint(v);       // 3.0  — nearest, respects current rounding mode, raises FE_INEXACT

#include <cmath>

// C++11: std::fma computes (x * y) + z with a single rounding step
double dot(double ax, double bx, double ay, double by) {
    return std::fma(ax, bx, ay * by);
}

#include <cmath>
#include <cerrno>
#include <cstring>

double safe_sqrt(double x) {
    errno = 0;
    double r = std::sqrt(x);
    if (errno == EDOM)
        throw std::domain_error("sqrt of negative");
    return r;
}

#include <cmath>
#include <cfloat>

// naive:   two roundings, up to 2 ULP error
double naive = x * y + z;

// fma:     one rounding, at most 0.5 ULP error (correctly rounded)
double exact = std::fma(x, y, z);

#include <cmath>

#ifdef FP_FAST_FMA
    double r = std::fma(a, b, c);   // single instruction
#else
    double r = a * b + c;           // fast but two-rounded; or accept slower std::fma
#endif

// Evaluate p(x) = c0 + c1*x + c2*x^2 + c3*x^3
// Horner form: ((c3*x + c2)*x + c1)*x + c0
double horner(double x, double c0, double c1, double c2, double c3) {
    double r = c3;
    r = std::fma(r, x, c2);
    r = std::fma(r, x, c1);
    r = std::fma(r, x, c0);
    return r;
}

double dot(const double* a, const double* b, std::size_t n) {
    double sum = 0.0;
    for (std::size_t i = 0; i < n; ++i)
        sum = std::fma(a[i], b[i], sum);
    return sum;
}

// Exactly computes a*b - c*d without intermediate cancellation error
double diff_of_products(double a, double b, double c, double d) {
    double w  = c * d;
    double err = std::fma(-c, d, w);   // exact rounding error of c*d
    double r  = std::fma(a, b, -w);
    return r + err;
}

#include <cmath>
#include <limits>

double big = std::numeric_limits<double>::max() / 2.0;

// Naive overflows to inf even though the true answer is finite:
double bad = std::sqrt(big * big + big * big);   // +inf

// std::hypot handles it correctly:
double good = std::hypot(big, big);              // 1.27...e+307

// Two-argument form (C++11)
double      std::hypot(double x, double y);
float       std::hypot(float x, float y);
long double std::hypot(long double x, long double y);

// Three-argument form (C++17): sqrt(x^2 + y^2 + z^2)
double      std::hypot(double x, double y, double z);
float       std::hypot(float x, float y, float z);
long double std::hypot(long double x, long double y, long double z);

// Conceptual implementation (simplified)
// The standard library version handles all IEEE 754 edge cases.
double safe_hypot2(double x, double y) {
    x = std::fabs(x);
    y = std::fabs(y);
    if (x < y) std::swap(x, y);
    if (x == 0.0) return 0.0;
    double t = y / x;
    return x * std::sqrt(std::fma(t, t, 1.0));
}

#include <cmath>
#include <array>

// 3-D distance between two points (C++17)
struct Vec3 { double x, y, z; };

double distance(Vec3 a, Vec3 b) {
    return std::hypot(a.x - b.x, a.y - b.y, a.z - b.z);
}

// Magnitude of a complex number without <complex>
double complex_abs(double re, double im) {
    return std::hypot(re, im);
}

// math_errhandling is a macro, NOT in namespace std
// It holds a bitwise OR of:
//   MATH_ERRNO    (1) — errors set errno
//   MATH_ERREXCEPT (2) — errors raise IEEE 754 floating-point exceptions

if (math_errhandling & MATH_ERRNO) {
    // errors are reported via errno
}
if (math_errhandling & MATH_ERREXCEPT) {
    // errors raise FE_INVALID, FE_DIVBYZERO, or FE_OVERFLOW
}

#include <cmath>
#include <cfenv>

#pragma STDC FENV_ACCESS ON    // required pragma: tell the compiler not to
                               // reorder or optimise away fenv accesses

double checked_log(double x) {
    std::feclearexcept(FE_ALL_EXCEPT);         // clear all flags first

    double result = std::log(x);

    if (std::fetestexcept(FE_INVALID)) {
        // x was NaN — result is NaN
    }
    if (std::fetestexcept(FE_DIVBYZERO)) {
        // x was exactly 0.0 — result is -inf (pole error)
    }
    if (std::fetestexcept(FE_OVERFLOW)) {
        // result too large for double
    }
    return result;
}

#include <cmath>
#include <limits>

double x = std::numeric_limits<double>::quiet_NaN();

// All of these produce NaN — the "poison" propagates:
double a = x + 1.0;         // NaN
double b = std::sqrt(x);    // NaN
double c = std::fma(x, 2.0, 3.0); // NaN
double d = std::sin(x);     // NaN

// But comparison operators always return false with NaN:
bool eq = (x == x);  // false  — the classic NaN self-inequality test
bool ne = (x != x);  // true   — reliable NaN check without <cmath>

// Prefer the standard predicate:
bool is_bad = std::isnan(x);  // true

#include <limits>

float trap = std::numeric_limits<float>::signaling_NaN();
// Reading 'trap' in arithmetic raises FE_INVALID on platforms
// that support signaling NaN (most x86 and ARM systems do).

#include <cmath>

void describe(double v) {
    switch (std::fpclassify(v)) {
        case FP_NAN:       // either quiet or signaling NaN
        case FP_INFINITE:  // +inf or -inf; use std::signbit to distinguish
        case FP_ZERO:      // +0.0 or -0.0; std::signbit(v) for signed zero
        case FP_SUBNORMAL: // |v| < DBL_MIN (gradual underflow)
        case FP_NORMAL:    // ordinary representable value
    }
}

#include <cmath>  // special functions are in <cmath> from C++17

// All signatures follow the pattern:
//   double   func_name(args...);
//   float    func_namef(args...);
//   long double func_namel(args...);

// Laguerre polynomials L_n(x) — quantum harmonic oscillator,
// hydrogen atom radial wave functions
double std::laguerre(unsigned n, double x);

// Associated Laguerre polynomials L_n^m(x) — hydrogen atom
// wave functions with orbital angular momentum quantum number m
double std::assoc_laguerre(unsigned n, unsigned m, double x);

// Legendre polynomials P_l(x) — solutions on [-1, 1],
// gravitational/electrostatic multipole expansions, x in [-1,1]
double std::legendre(unsigned l, double x);

// Associated Legendre functions P_l^m(x) — building block
// for spherical harmonics
double std::assoc_legendre(unsigned l, unsigned m, double x);

// Hermite polynomials H_n(x) — quantum harmonic oscillator
// energy eigenstates, probabilists' form
double std::hermite(unsigned n, double x);

// Spherical harmonics angular part Y_l^m(theta, phi=0):
// sph_legendre gives the theta-dependent factor
double std::sph_legendre(unsigned l, unsigned m, double theta);

// Bessel functions of the first kind J_nu(x):
// bounded at the origin; oscillatory for large x
double std::cyl_bessel_j(double nu, double x);

// Neumann functions (Bessel functions of the second kind) Y_nu(x):
// diverge at the origin; oscillatory for large x
double std::cyl_neumann(double nu, double x);

// Modified Bessel functions of the first kind I_nu(x):
// exponentially growing; arise in probability (Rice distribution)
double std::cyl_bessel_i(double nu, double x);

// Modified Bessel functions of the second kind K_nu(x):
// exponentially decaying; arise in heat transfer, elasticity
double std::cyl_bessel_k(double nu, double x);

// Spherical Bessel functions of the first kind j_n(x):
// appear in quantum scattering, acoustic resonators
double std::sph_bessel(unsigned n, double x);

// Spherical Neumann functions n_n(x)
double std::sph_neumann(unsigned n, double x);

// Complete elliptic integral of the first kind K(k):
// period of a nonlinear pendulum, k is the modulus in [0,1)
double std::comp_ellint_1(double k);

// Complete elliptic integral of the second kind E(k):
// arc length of an ellipse with eccentricity k
double std::comp_ellint_2(double k);

// Complete elliptic integral of the third kind Pi(nu, k)
double std::comp_ellint_3(double k, double nu);

// Incomplete forms with upper limit phi:
double std::ellint_1(double k, double phi);
double std::ellint_2(double k, double phi);
double std::ellint_3(double k, double nu, double phi);

// Beta function B(x, y) = Gamma(x)*Gamma(y)/Gamma(x+y):
// normalisation constant for the beta distribution,
// combinatorics, Bayesian statistics
double std::beta(double x, double y);

// Exponential integral Ei(x) = -P.V. integral from -x to inf of e^{-t}/t dt:
// radiative transfer, plasma physics
double std::expint(double x);

// Riemann zeta function zeta(s):
// analytic number theory, quantum field theory, statistical mechanics
double std::riemann_zeta(double s);

#include <cmath>
#include <numbers>  // C++20

// Radial probability density |R_10(r)|^2 for hydrogen 1s orbital
// using the associated Laguerre polynomial (n=1, l=0, so L_0^1)
double hydrogen_1s_radial_density(double r_over_a0) {
    double rho = 2.0 * r_over_a0;          // dimensionless radius
    double L   = std::assoc_laguerre(0, 1, rho);  // L_0^1(rho) = 1
    double R   = 2.0 * std::exp(-r_over_a0) * L;
    return R * R;  // |R|^2 is the probability density
}

#include <cmath>   // C++20: std::lerp

float blend_colour(float from, float to, float t) {
    return std::lerp(from, to, t);  // guaranteed exact at t=0 and t=1
}

// Why the naive formula breaks:
float a = 0.1f, b = 0.2f;
float t = 1.0f;

float naive = a + t * (b - a);   // may not equal b exactly due to
                                  // floating-point subtraction error
float safe  = std::lerp(a, b, t); // == b exactly

#include <numeric>  // C++20: std::midpoint

// Classic bug in binary search:
int mid_bad  = (lo + hi) / 2;     // overflows for large lo and hi
int mid_good = std::midpoint(lo, hi);  // never overflows

// Works for floating-point too:
double fmid = std::midpoint(1e308, -1e308);  // 0.0, not NaN

// Binary search using std::midpoint (C++20):
#include <algorithm>
#include <numeric>
#include <span>

int binary_search_safe(std::span<const int> arr, int target) {
    int lo = 0, hi = static_cast<int>(arr.size()) - 1;
    while (lo <= hi) {
        int mid = std::midpoint(lo, hi);
        if (arr[mid] == target)  return mid;
        if (arr[mid] < target)   lo = mid + 1;
        else                     hi = mid - 1;
    }
    return -1;
}