<cfenv>

#include <cfenv>

// Opaque types
std::fenv_t    // complete floating-point environment state
std::fexcept_t // exception flag state (a subset of fenv_t)

// Exception flag macros (a platform may define any subset; unsupported ones == 0)
FE_DIVBYZERO   // x / 0.0 where x is finite and non-zero
FE_INEXACT     // result is not exactly representable
FE_INVALID     // e.g. 0.0/0.0, sqrt(-1.0), NaN comparisons
FE_OVERFLOW    // result magnitude exceeds the representable range
FE_UNDERFLOW   // result magnitude is too small (subnormal or flushed to zero)
FE_ALL_EXCEPT  // bitwise OR of all supported exception macros

// Rounding mode macros
FE_DOWNWARD    // toward -∞
FE_TONEAREST   // round-half-to-even (IEEE 754 default, banker's rounding)
FE_TOWARDZERO  // toward zero (truncation)
FE_UPWARD      // toward +∞

FE_DFL_ENV     // pointer to const fenv_t representing the default environment

// Exception flag functions
int std::feclearexcept(int excepts);                                  // clear flags
int std::feraiseexcept(int excepts);                                  // raise flags
int std::fetestexcept(int excepts);                                   // test flags (returns set subset)
int std::fegetexceptflag(std::fexcept_t* flagp, int excepts);        // snapshot flags
int std::fesetexceptflag(const std::fexcept_t* flagp, int excepts);  // restore snapshot

// Rounding mode
int std::fegetround();               // returns current rounding mode macro, or -1
int std::fesetround(int round);      // 0 on success

// Full environment save/restore
int std::fegetenv(std::fenv_t* envp);                // snapshot env (flags + rounding + more)
int std::fesetenv(const std::fenv_t* envp);          // restore env verbatim
int std::feholdexcept(std::fenv_t* envp);            // snapshot, clear all flags, non-stop mode
int std::feupdateenv(const std::fenv_t* envp);       // restore env; re-raise pre-existing exceptions

#pragma STDC FENV_ACCESS ON
#include <cfenv>
#include <cmath>
#include <cstdio>

// math_errhandling & MATH_ERREXCEPT must be non-zero for <cmath>
// functions to actually raise FP exception flags (true on most hosted platforms).
void probe(double x) {
    std::feclearexcept(FE_ALL_EXCEPT);  // flags are sticky — clear before the region

    double r = std::log(x);

    if (std::fetestexcept(FE_INVALID))   { std::puts("invalid: argument < 0"); return; }
    if (std::fetestexcept(FE_DIVBYZERO)) { std::puts("pole: log(0) = -inf");   return; }
    if (std::fetestexcept(FE_OVERFLOW))  { std::puts("overflow");               return; }

    std::printf("log(%g) = %.17g  (inexact=%d)\n",
                x, r, !!(std::fetestexcept(FE_INEXACT)));
}

int main() {
    probe(-1.0);   // FE_INVALID
    probe(0.0);    // FE_DIVBYZERO
    probe(1.0);    // exact; FE_INEXACT not set
    probe(1e300);  // FE_INEXACT
}

#pragma STDC FENV_ACCESS ON
#include <cfenv>
#include <cstdio>

struct Interval { double lo, hi; };

Interval iadd(Interval a, Interval b) {
    double lo, hi;

    std::fesetround(FE_DOWNWARD);   // C++11 — round toward -∞ for lower bound
    lo = a.lo + b.lo;

    std::fesetround(FE_UPWARD);     // C++11 — round toward +∞ for upper bound
    hi = a.hi + b.hi;

    std::fesetround(FE_TONEAREST);  // always restore default before returning
    return {lo, hi};
}

Interval imul(Interval a, Interval b) {
    // For positive intervals only (full case requires checking all four products)
    double lo, hi;
    std::fesetround(FE_DOWNWARD); lo = a.lo * b.lo;
    std::fesetround(FE_UPWARD);   hi = a.hi * b.hi;
    std::fesetround(FE_TONEAREST);
    return {lo, hi};
}

int main() {
    // 0.1 + 0.2 is not exactly representable; the interval must contain 0.3
    Interval x{0.1, 0.1}, y{0.2, 0.2};
    auto z = iadd(x, y);
    std::printf("[%.17g, %.17g]\n", z.lo, z.hi);
    // [0.29999999999999998, 0.30000000000000004] — contains 0.3
}

#pragma STDC FENV_ACCESS ON
#include <cfenv>
#include <cmath>
#include <stdexcept>

class FenvGuard {
    std::fenv_t saved_;
public:
    FenvGuard()  { std::feholdexcept(&saved_); }  // save, clear all flags, non-stop mode
    ~FenvGuard() { std::feupdateenv(&saved_); }   // restore; re-raise pre-existing exceptions
    FenvGuard(const FenvGuard&) = delete;
    FenvGuard& operator=(const FenvGuard&) = delete;
};

// Caller's FP environment is fully preserved regardless of what happens inside.
double hypot_safe(double a, double b) {
    FenvGuard guard;
    std::fesetround(FE_TONEAREST);
    double r = std::sqrt(a * a + b * b);
    if (std::fetestexcept(FE_OVERFLOW)) {
        // Handle internally before guard destructor restores caller's state
        r = std::numeric_limits<double>::infinity();
    }
    return r;
}

#pragma STDC FENV_ACCESS ON
#include <cfenv>
#include <cmath>
#include <cassert>
#include <cstdio>

// Returns true if the computation was free of all five exceptional conditions.
bool audit_computation(double* result, double x) {
    std::feclearexcept(FE_ALL_EXCEPT);  // C++11 — wipe slate before region

    *result = std::log1p(x) * std::sqrt(x) / (x + 1.0);

    int raised = std::fetestexcept(FE_ALL_EXCEPT);
    if (raised) {
        if (raised & FE_INVALID)   std::puts("  audit: FE_INVALID raised");
        if (raised & FE_DIVBYZERO) std::puts("  audit: FE_DIVBYZERO raised");
        if (raised & FE_OVERFLOW)  std::puts("  audit: FE_OVERFLOW raised");
        if (raised & FE_UNDERFLOW) std::puts("  audit: FE_UNDERFLOW raised");
        if (raised & FE_INEXACT)   std::puts("  audit: FE_INEXACT raised");
    }
    return raised == 0;
}

// In a test suite:
void test_well_behaved_input() {
    double result;
    bool clean = audit_computation(&result, 1.0);
    // FE_INEXACT is expected for most transcendental results; mask it if acceptable
    int raised = std::fetestexcept(FE_ALL_EXCEPT & ~FE_INEXACT);
    assert(raised == 0 && "computation raised unexpected FP exception");
}

#pragma STDC FENV_ACCESS ON
#include <cfenv>
#include <cmath>
#include <cstdio>
#include <vector>
#include <numeric>

// Sum 1000 doubles under FE_UPWARD, then FE_DOWNWARD.
// The difference is a lower-bound estimate of the accumulated rounding error.
double directed_sum(const std::vector<double>& v, int rounding_mode) {
    std::fesetround(rounding_mode);  // C++11
    double acc = 0.0;
    for (double x : v) {
        acc += x;  // each addition rounds in the directed direction
    }
    std::fesetround(FE_TONEAREST);  // restore default before returning
    return acc;
}

int main() {
    // Construct a vector where catastrophic cancellation is likely:
    // alternating large positive and negative values whose exact sum is small.
    std::vector<double> data;
    data.reserve(1000);
    for (int i = 0; i < 500; ++i) {
        data.push_back( 1e15 + i);   // large positive
        data.push_back(-1e15 - i);   // large negative (exact sum of each pair = 0)
    }
    // Add a small meaningful signal at the end
    data.push_back(42.0);

    double sum_up   = directed_sum(data, FE_UPWARD);    // C++11
    double sum_down = directed_sum(data, FE_DOWNWARD);  // C++11

    double spread = std::abs(sum_up - sum_down);

    std::printf("FE_UPWARD   result: %.17g\n", sum_up);
    std::printf("FE_DOWNWARD result: %.17g\n", sum_down);
    std::printf("Spread (error estimate): %.6g\n", spread);
    std::printf("True answer: 42.0\n");
    // If spread >> 1 ULP, the summation is numerically unstable for this input.
}

#pragma STDC FENV_ACCESS ON
#include <cfenv>
#include <cstdio>
#include <vector>

// Kahan compensated summation — O(1) rounding error regardless of n.  C++11
double kahan_sum(const std::vector<double>& v) {
    double sum  = 0.0;
    double comp = 0.0;  // running compensation for lost low-order bits

    for (double x : v) {
        double y = x - comp;       // compensated addend: correct for previous loss
        double t = sum + y;        // sum + y: low-order bits of y are lost here
        comp = (t - sum) - y;      // recover the lost bits (should be ~0 in exact arith)
        sum  = t;
    }
    return sum;
}

// Naive sum for comparison
double naive_sum(const std::vector<double>& v) {
    double acc = 0.0;
    for (double x : v) acc += x;
    return acc;
}

int main() {
    // 10^7 additions of 0.1: exact answer is 1,000,000.0
    std::vector<double> data(10'000'000, 0.1);

    double naive  = naive_sum(data);
    double kahan  = kahan_sum(data);

    std::printf("Naive:  %.17g  (error: %.6g)\n", naive,  naive  - 1e7 * 0.1);
    std::printf("Kahan:  %.17g  (error: %.6g)\n", kahan,  kahan  - 1e7 * 0.1);
    // Naive error grows proportionally with n; Kahan error is bounded.
}

// Example: when FMA contraction is desirable — better accuracy, one rounding.  C++11
#include <cmath>

double dot_product_fast(const double* a, const double* b, int n) {
    // With FP_CONTRACT ON (the default), the compiler may emit FMA instructions.
    // This is almost always what you want: faster and more accurate.
#pragma STDC FP_CONTRACT ON
    double sum = 0.0;
    for (int i = 0; i < n; ++i)
        sum += a[i] * b[i];  // may become: sum = fma(a[i], b[i], sum)
    return sum;
}

// Example: when FMA contraction must be disabled — exact IEEE 754 reproducibility.
double dot_product_reproducible(const double* a, const double* b, int n) {
    // With FP_CONTRACT OFF, every multiply and every add is a distinct
    // IEEE 754 operation with its own rounding.  Results are identical
    // across compilers, architectures, and optimization levels.
#pragma STDC FP_CONTRACT OFF
    double sum = 0.0;
    for (int i = 0; i < n; ++i)
        sum += a[i] * b[i];
    return sum;
}

// feenableexcept is POSIX/Linux only — guard all uses.  C++11
#ifdef __linux__
#include <cfenv>
#include <csignal>
#include <cstdio>
#include <cstdlib>
#include <execinfo.h>   // GNU extension for backtrace()

// SIGFPE handler: print a stack trace, then abort.
static void fp_exception_handler(int /*sig*/, siginfo_t* info, void* /*ctx*/) {
    std::puts("\n*** Floating-point exception detected ***");

    // Decode which exception triggered the signal (si_code)
    switch (info->si_code) {
        case FPE_FLTDIV: std::puts("  cause: division by zero");    break;
        case FPE_FLTINV: std::puts("  cause: invalid operation");   break;
        case FPE_FLTOVF: std::puts("  cause: overflow");            break;
        case FPE_FLTUND: std::puts("  cause: underflow");           break;
        case FPE_FLTRES: std::puts("  cause: inexact result");      break;
        default:         std::puts("  cause: unknown FP exception"); break;
    }

    // Print a raw stack trace (requires linking with -rdynamic for symbol names)
    void* frames[64];
    int   depth = backtrace(frames, 64);
    backtrace_symbols_fd(frames, depth, STDERR_FILENO);

    std::abort();
}

// Call once at program startup (in main, before any numeric work).
void install_fp_traps() {
    struct sigaction sa{};
    sa.sa_flags     = SA_SIGINFO;
    sa.sa_sigaction = fp_exception_handler;
    sigemptyset(&sa.sa_mask);
    sigaction(SIGFPE, &sa, nullptr);

    // Enable traps for the three most dangerous silent exceptions.
    // FE_INEXACT fires on almost every FP operation — do not enable it here.
    feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
}

int main() {
    install_fp_traps();

    double x = 0.0;
    double y = 1.0 / x;   // <- SIGFPE delivered here; handler prints backtrace
    (void)y;
    return 0;
}
#endif  // __linux__

#pragma STDC FENV_ACCESS ON
#include <cfenv>
#include <cstdio>

struct FpCapabilities {
    bool rounding_modes;    // fesetround actually works
    bool exception_flags;   // fetestexcept actually works
    bool math_errexcept;    // <cmath> raises exception flags
};

// Probe FP capabilities at startup.  C++11
FpCapabilities probe_fp_capabilities() {
    FpCapabilities caps{};

    // Test rounding mode support: try to set FE_UPWARD and read it back.
    if (std::fesetround(FE_UPWARD) == 0) {
        caps.rounding_modes = (std::fegetround() == FE_UPWARD);
    }
    std::fesetround(FE_TONEAREST);  // restore default regardless

    // Test exception flag support: raise FE_INEXACT manually and test it.
    std::feclearexcept(FE_ALL_EXCEPT);
    if (std::feraiseexcept(FE_INEXACT) == 0) {
        caps.exception_flags = !!(std::fetestexcept(FE_INEXACT));
    }
    std::feclearexcept(FE_ALL_EXCEPT);

    // Test whether <cmath> raises flags (MATH_ERREXCEPT bit).
    caps.math_errexcept = !!(math_errhandling & MATH_ERREXCEPT);

    return caps;
}

int main() {
    FpCapabilities caps = probe_fp_capabilities();
    std::printf("Rounding modes:    %s\n", caps.rounding_modes  ? "yes" : "no (fesetround is a no-op)");
    std::printf("Exception flags:   %s\n", caps.exception_flags ? "yes" : "no (flags always zero)");
    std::printf("math_errexcept:    %s\n", caps.math_errexcept  ? "yes" : "no (use errno instead)");

    if (!caps.rounding_modes) {
        // Fall back: use software interval arithmetic or avoid directed rounding.
    }
    if (!caps.exception_flags) {
        // Fall back: check errno after every math call, or use isnan/isinf.
    }
}