∑MathIntermediate

Numerical Stability

Key Points

•
Numerical stability measures how much rounding and tiny input changes can distort an algorithm’s output on real computers using floating-point arithmetic.
•
Forward error asks “How far is my computed answer from the true answer?”, while backward error asks “For which nearby input would my computed answer be exact?”.
•
Backward-stable algorithms are ideal: they solve a slightly perturbed problem exactly, and the perturbation is on the order of machine precision.
•
The condition number links problem sensitivity to errors: forward error ≈ condition number × backward error.
•
Subtractive cancellation, catastrophic reordering of operations, and lack of pivoting are common causes of instability.
•
Simple fixes like Kahan summation, Horner’s method for polynomials, and partial pivoting in Gaussian elimination greatly improve stability.
•
Stability and conditioning are different: even a perfectly stable algorithm will struggle on an ill-conditioned problem.
•
Stable methods often have the same big-O time as naive ones, with small constant-factor overhead for accuracy.

Prerequisites

→IEEE floating-point representation (binary32/binary64) — Understanding rounding, subnormal numbers, overflow/underflow, and unit roundoff is essential for stability reasoning.
→Vector and matrix norms — Condition numbers and error bounds are expressed in norms; you need norms to interpret and compute these bounds.
→Big-O notation — To compare algorithmic cost among stable and naive methods and to understand complexity statements.
→Basic calculus (derivatives) — Scalar condition number κ(f,x) = |x f'(x)/f(x)| relies on derivatives.
→Linear algebra (LU/QR, Gaussian elimination) — Stability of solving Ax=b depends on factorization choices and pivoting.
→C++ programming basics — Implementing and benchmarking stable algorithms requires facility with loops, arrays, and numeric types.
→Error analysis concepts — Absolute/relative error, residuals, and perturbation analysis are the language of numerical stability.
→Algorithmic reformulation — Recognizing and applying algebraic transformations (e.g., Horner, log-sum-exp) improves stability.

Detailed Explanation

Tap terms for definitions

01Overview

Computers use floating-point numbers, which are approximations to real numbers with limited precision. Every arithmetic operation introduces a tiny rounding error. Numerical stability studies how these small errors grow—or remain controlled—when an algorithm runs. Two complementary lenses help: forward error (how wrong the final answer is) and backward error (how much we would need to nudge the input so that the computed answer becomes exactly right). A backward-stable algorithm is usually considered ideal because it guarantees the computed answer is the exact solution to a nearby problem, with the “nearness” on the order of machine precision. The gap between the true and computed answers is also influenced by the problem’s sensitivity, captured by the condition number. Well-conditioned problems don’t magnify input errors much; ill-conditioned ones do. The art of numerical stability is to design and choose algorithms that keep rounding errors from exploding, especially for sensitive problems. Hook → Concept → Example: Imagine adding dollar amounts in a huge ledger. If you add pennies to billions in the wrong order, you can lose track of the pennies due to rounding. That’s instability. The concept is to control such error growth, often by rearranging computations. For example, Kahan’s compensated summation tracks lost pennies (rounding) so they aren’t forgotten.

02Intuition & Analogies

Think of carrying water with a cup that spills a drop every time you pour. If you pour once, no big deal. If you pour millions of times, those drops add up. Floating-point arithmetic is like that cup: each operation "spills" a tiny amount of information. Numerical stability asks: after all these pours, how wet is the floor and how much water did we lose from the bucket we care about? Another image: adding a small coin to a giant pile of cash. If your counting machine can’t register tiny additions on top of a huge total, the small coins may get ignored—this is what happens when you add small numbers to large ones in floating point. A stable summation algorithm ensures the coins aren’t lost by keeping a separate side-pocket for the small change and combining it later (Kahan summation). For solving linear equations, imagine using a teeter-totter balance with nearly equal weights on both sides. A tiny nudge can flip it drastically; that system is ill-conditioned. Even if your method of weighing (algorithm) is precise (stable), the setup itself is so sensitive that small disturbances create large changes in the result. Thus, two forces interact: the problem’s condition number (how twitchy the setup is) and the algorithm’s stability (how well your procedure avoids injecting extra disturbances). Backward stability is like saying, “My weighing method is so reliable that the reading equals what you’d get if the weights were adjusted by an imperceptibly tiny amount.”

03Formal Definition

Let f map inputs x to outputs

y = f (x

). Given a floating-point algorithm that returns \(

\overset{y}{^}

\), the absolute forward error is \(

∥

\overset{y}{^}

- y

∥

\) and the relative forward error is \(

\frac{∥ y ^ - y ∥}{∥ y ∥}

\). Backward error measures how much we perturb the input so that the computed output becomes exact for a nearby problem: find \(

Δ

x\) such that \(

\overset{y}{^}

= f(x +

Δ

x)\); the relative backward error is \(

\frac{∥ Δ x ∥}{∥ x ∥}

\). The condition number \(

κ

(f,x)\) quantifies problem sensitivity. In the scalar case with differentiable f, \(

κ

(f,x) =

\frac{x f ^{'} ( x )}{f ( x )}

\). In vector/matrix settings with norms, \(

κ

(f,x) =

lim_{ε \to 0}

sup_{∥ Δ x ∥ \leq ε ∥ x ∥}

\frac{∥ f ( x + Δ x ) - f ( x )∥ / ∥ f ( x )∥}{ε}

\). A key relation is that, to first order, relative forward error is bounded by condition number times relative backward error. An algorithm is backward stable if its backward error is bounded by a modest multiple of machine precision u (often times a low-degree polynomial in problem size). For linear systems \(

Ax=b\

), the matrix condition number is \(

κ

(A) =

∥

∥

∥

A^{- 1}

∥

\). Gaussian elimination with partial pivoting is typically backward stable, with bounds that may involve a growth factor but are practically small for most matrices.

04When to Use

Use numerical stability analysis whenever you design or select algorithms that run in floating point—practically always in scientific computing, data analysis, graphics, and machine learning. Choose Kahan or pairwise summation for long sums of mixed-magnitude numbers, such as aggregating gradients, probabilities, or financial ledgers. Prefer Horner’s method when evaluating polynomials, especially near roots where terms nearly cancel. In linear algebra, solve (Ax=b) with Gaussian elimination with partial pivoting (or use robust libraries like LAPACK) rather than naive elimination; pivoting is essential when diagonal entries are small. Apply stable transforms like log-sum-exp for softmax or partition functions to prevent overflow/underflow. Scale inputs (e.g., normalize vectors, rescale units) to avoid catastrophic cancellation and to keep values within comfortable dynamic ranges. When problems are ill-conditioned (large (\kappa)), expect limited accuracy regardless of algorithm; in such cases, consider reformulations (regularization, better parameterizations, orthogonal factorizations like QR) that improve conditioning.

⚠️Common Mistakes

Ignoring conditioning: Blaming the algorithm when the problem is inherently ill-conditioned. Always estimate or reason about (\kappa); if it’s huge, even perfect stability won’t rescue accuracy.
No pivoting in Gaussian elimination: Dividing by tiny pivots explodes rounding errors. Always use at least partial pivoting or better factorizations (QR with Householder) when solving linear systems.
Subtractive cancellation: Directly computing differences of nearly equal numbers (e.g., (\sqrt{x+1}-\sqrt{x})) causes loss of significant digits. Use algebraic reformulations (rationalize) or compensated techniques.
Bad summation order: Adding small numbers to very large ones loses the small terms. Sum from small to large, use pairwise/Kahan summation, or both.
Overflow/underflow: Computing (e^{x}) for large (x) or products of many terms without scaling leads to inf/0. Use log-domain tricks (log-sum-exp, softmax shifting) and normalization.
Using float instead of double by default: Single precision halves the number of significant digits versus double; choose precision intentionally.
Misinterpreting backward stability: A backward-stable algorithm can still yield a result with large forward error if the problem is ill-conditioned. Stability is not a silver bullet.
Skipping validation: Not checking residuals (e.g., (\lVert A\hat{x}-b \rVert)) or sensitivity diagnostics can hide serious numerical issues.

Key Formulas

Floating-point rounding model

fl (a \circ b) = (a \circ b) (1 + δ), ∣ δ ∣ \leq u

Explanation: Each basic operation (add, subtract, multiply, divide) is rounded to the nearest representable number, introducing a small relative error bounded by the unit roundoff u. This model underpins most stability analyses.

Forward error

abs. forward error = ∥ \overset{y}{^} - y ∥, rel. forward error = \frac{∥ y ^ - y ∥}{∥ y ∥}

Explanation: Measures how far the computed solution is from the true one, either absolutely or relative to the true magnitude. Use this to assess final accuracy.

Backward error

rel. backward error = Δ x : \overset{y}{^} = f (x + Δ x) min \frac{∥ Δ x ∥}{∥ x ∥}

Explanation: Quantifies the smallest relative input change that makes the computed output exact. Small backward error indicates a backward-stable algorithm.

Condition number (general)

κ (f, x) = ε \to 0 lim ∥ Δ x ∥ \leq ε ∥ x ∥ sup \frac{∥ f ( x + Δ x ) - f ( x )∥ / ∥ f ( x )∥}{ε}

Explanation: This is the worst-case relative output change per unit relative input change near x. It measures the intrinsic sensitivity of the problem itself.

Scalar condition number

κ (f, x) = \frac{x f ^{'} ( x )}{f ( x )}

Explanation: For scalar differentiable functions, this simple formula measures sensitivity. If κ is large, small input errors will be amplified in the output.

Matrix condition number

κ (A) = ∥ A ∥ ∥ A^{- 1} ∥

Explanation: For solving Ax=b, κ(A) bounds how much input perturbations affect x. In 2-norm, κ is the ratio of largest to smallest singular values.

Forward vs backward error

\frac{∥ y ^ - y ∥}{∥ y ∥} ≲ κ (f, x) \cdot \frac{∥ Δ x ∥}{∥ x ∥} + O (u^{2})

Explanation: Forward error is approximately the product of condition number and backward error. If the algorithm is backward stable (small backward error), the forward error is controlled unless κ is large.

Unit roundoff

u = \frac{1}{2} β^{1 - p}

Explanation: For base-β floating point with p digits and round-to-nearest, unit roundoff u bounds one-step relative error. In IEEE binary64, $β = 2$ and $p = 53$ .

Error accumulation factor

γ_{n} = \frac{n u}{1 - n u}

Explanation: In naive summation or dot products, the relative error can grow like $γ_{n} .$ For small nu, $γ_{n}$ ≈ n u, showing linear growth with the number of operations.

Summation error bound

fl (i = 1 \sum n a_{i}) - i = 1 \sum n a_{i} \leq γ_{n} i = 1 \sum n ∣ a_{i} ∣

Explanation: Naive left-to-right summation accumulates rounding error approximately proportional to n u times the 1-norm of inputs. This motivates compensated or pairwise summation.

Kahan summation bound

∣ s_{Kahan} - s ∣ \leq (u + O (n u^{2})) i = 1 \sum n ∣ a_{i} ∣

Explanation: Kahan’s method reduces the leading error growth from O(nu) to essentially O(u) for well-behaved data. It tracks and compensates lost low-order bits.

Log-sum-exp stabilization

lo g i = 1 \sum n e^{x_{i}} = m + lo g i = 1 \sum n e^{x_{i} - m}, m = i max x_{i}

Explanation: Shifting by the maximum prevents overflow and underflow when exponentiating large or very negative numbers. This is widely used in softmax and partition functions.

Growth factor in GEPP

ρ = i, j, k max \frac{∣ a _{ij}^{(k)} ∣}{max _{i, j} ∣ a _{ij}^{(0)} ∣}

Explanation: The growth factor ρ measures how much entries grow during Gaussian elimination with partial pivoting. Error bounds often scale with ρ, which is typically modest in practice.

Complexity Analysis

Stable algorithms often have the same asymptotic cost as naive ones, with modest constant-factor overhead. For summation of n numbers, both naive and Kahan summation take O(n) time and O(1) extra space. Kahan adds a couple of floating-point operations per term to maintain a compensation variable, usually a negligible overhead for substantial accuracy gains. Polynomial evaluation via Horner’s method costs O(n) time and O(1) space, compared with O(

n^{2}

) if one computes each power from scratch; Horner is both faster and generally more stable due to fewer operations and less cancellation. Solving linear systems with Gaussian elimination has O(

n^{3}

) time and O(

n^{2}

) space with or without pivoting. Partial pivoting adds an O(

n^{2}

) search for each column (comparable to the O(

n^{2}

) updates in that step), so the overall complexity remains O(

n^{3}

). The memory footprint is also essentially unchanged: storing the LU factors overwrites A in-place, using O(

n^{2}

) memory. In practice, the extra cost of stability (pivoting, compensation) is small relative to the cost of the core arithmetic, and modern BLAS/LAPACK implementations are highly optimized. It’s important to note that while stability affects error growth, it doesn’t change the big-O scalability. However, instability can lead to meaningless results that might force recomputation or algorithm changes, which is far more expensive than a small constant-factor overhead for a stable method. Thus, prefer stable methods by default, and benchmark only if the overhead matters for your workload.

Code Examples

Accurate summation with Kahan compensation vs naive sum

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 double naive_sum(const vector<double>& a) {
5     double s = 0.0;
6     for (double v : a) s += v; // susceptible to rounding, esp. mixing magnitudes
7     return s;
8 }
9 
10 double kahan_sum(const vector<double>& a) {
11     double s = 0.0;     // running sum
12     double c = 0.0;     // compensation for lost low-order bits
13     for (double v : a) {
14         double y = v - c;        // first remove the small error we carried
15         double t = s + y;        // add corrected value
16         c = (t - s) - y;         // new rounding error stored as compensation
17         s = t;
18     }
19     return s;
20 }
21 
22 int main() {
23     ios::sync_with_stdio(false);
24     cout.setf(std::ios::fixed); cout << setprecision(17);
25 
26     // Construct a challenging sum: one huge number, then many ones.
27     // Adding small numbers to a huge accumulator loses precision in naive summation.
28     const int m = 500000; // half a million ones
29     vector<double> a;
30     a.reserve(m + 1);
31     a.push_back(1e8);           // a large term
32     for (int i = 0; i < m; ++i) a.push_back(1.0); // many small terms
33 
34     // Reference using long double (higher precision) and exact math
35     long double ref = 1e8L + static_cast<long double>(m);
36 
37     double s_naive = naive_sum(a);
38     double s_kahan = kahan_sum(a);
39 
40     cout << "Reference (long double): " << (double)ref << "\n";
41     cout << "Naive sum:                " << s_naive << "\n";
42     cout << "Kahan sum:                " << s_kahan << "\n";
43 
44     cout << "Abs error naive:  " << fabs(s_naive - (double)ref) << "\n";
45     cout << "Abs error Kahan:  " << fabs(s_kahan - (double)ref) << "\n";
46 
47     // Tip: Pairwise or sorted-by-magnitude summation can also help when Kahan is unavailable.
48     return 0;
49 }
50

This program compares naive left-to-right summation with Kahan’s compensated summation on data where many small terms are added to a very large term. Naive summation loses low-order bits when adding small to large; Kahan keeps a running compensation to recover those bits. The long double reference approximates the true value closely, showing how much each method deviates.

Time: O(n)Space: O(1)

Stable polynomial evaluation: Horner’s method vs expanded form

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 // Evaluate polynomial given coefficients [a0, a1, ..., ak] representing:
5 // P(x) = a0*x^k + a1*x^(k-1) + ... + ak
6 // Naive: compute powers explicitly (more operations, more rounding)
7 double poly_naive(const vector<double>& coef, double x) {
8     int k = (int)coef.size() - 1;
9     double y = 0.0;
10     for (int i = 0; i <= k; ++i) {
11         int power = k - i;
12         y += coef[i] * pow(x, power);
13     }
14     return y;
15 }
16 
17 // Horner's method: nested multiplication
18 // P(x) = (...((a0*x + a1)*x + a2)*x + ... + ak)
19 double poly_horner(const vector<double>& coef, double x) {
20     double y = 0.0;
21     for (double a : coef) {
22         y = y * x + a;
23     }
24     return y;
25 }
26 
27 int main() {
28     cout.setf(std::ios::fixed); cout << setprecision(20);
29 
30     // Consider P(x) = (x - 1)^5 = x^5 - 5x^4 + 10x^3 - 10x^2 + 5x - 1
31     vector<double> coef = {1.0, -5.0, 10.0, -10.0, 5.0, -1.0};
32     double x = 1.0001; // near a root; naive evaluation suffers cancellation
33 
34     // True value using high precision
35     long double xL = (long double)x;
36     long double ref = pow(xL - 1.0L, 5.0L);
37 
38     double y_naive = poly_naive(coef, x);
39     double y_horner = poly_horner(coef, x);
40 
41     cout << "Reference (long double): " << (double)ref << "\n";
42     cout << "Naive evaluation:        " << y_naive << "\n";
43     cout << "Horner evaluation:       " << y_horner << "\n";
44     cout << "|Naive - ref|:           " << fabs(y_naive - (double)ref) << "\n";
45     cout << "|Horner - ref|:          " << fabs(y_horner - (double)ref) << "\n";
46 
47     return 0;
48 }
49

Evaluating (x−1)^5 near x=1 leads to severe cancellation if you use the expanded form, because large terms with alternating signs nearly cancel. Horner’s method reduces the number of operations and the exposure to cancellation, yielding a more accurate result with the same O(n) time and O(1) space.

Time: O(n)Space: O(1)

Gaussian elimination: unstable no-pivot vs backward-stable partial pivoting

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 using Matrix = vector<vector<double>>;
5 using Vec = vector<double>;
6 
7 Vec back_substitute(const Matrix& U, const Vec& b) {
8     int n = (int)U.size();
9     Vec x(n, 0.0);
10     for (int i = n - 1; i >= 0; --i) {
11         double s = b[i];
12         for (int j = i + 1; j < n; ++j) s -= U[i][j] * x[j];
13         x[i] = s / U[i][i];
14     }
15     return x;
16 }
17 
18 optional<Vec> solve_no_pivot(Matrix A, Vec b) {
19     int n = (int)A.size();
20     // Forward elimination without pivoting (UNSTABLE)
21     for (int k = 0; k < n - 1; ++k) {
22         if (fabs(A[k][k]) == 0.0) return nullopt; // breakdown
23         for (int i = k + 1; i < n; ++i) {
24             double m = A[i][k] / A[k][k];
25             for (int j = k; j < n; ++j) A[i][j] -= m * A[k][j];
26             b[i] -= m * b[k];
27         }
28     }
29     // Back substitution
30     return back_substitute(A, b);
31 }
32 
33 optional<Vec> solve_gepp(Matrix A, Vec b) {
34     int n = (int)A.size();
35     // Gaussian Elimination with Partial Pivoting (GEPP)
36     for (int k = 0; k < n - 1; ++k) {
37         // Find pivot row with largest |A[i][k]|
38         int p = k;
39         double best = fabs(A[k][k]);
40         for (int i = k + 1; i < n; ++i) {
41             double val = fabs(A[i][k]);
42             if (val > best) { best = val; p = i; }
43         }
44         if (best == 0.0) return nullopt; // singular
45         if (p != k) { swap(A[p], A[k]); swap(b[p], b[k]); }
46         for (int i = k + 1; i < n; ++i) {
47             double m = A[i][k] / A[k][k];
48             A[i][k] = 0.0; // store L's entry implicitly if desired
49             for (int j = k + 1; j < n; ++j) A[i][j] -= m * A[k][j];
50             b[i] -= m * b[k];
51         }
52     }
53     return back_substitute(A, b);
54 }
55 
56 int main() {
57     cout.setf(std::ios::fixed); cout << setprecision(17);
58 
59     // A tiny leading pivot causes trouble without pivoting
60     Matrix A = {{1e-20, 1.0},
61                 {1.0,    1.0}};
62     Vec b = {1.0, 2.0};
63 
64     auto x_np = solve_no_pivot(A, b);
65     auto x_pp = solve_gepp(A, b);
66 
67     // Reference solution via long double (solve analytically for 2x2)
68     long double a11 = 1e-20L, a12 = 1.0L, a21 = 1.0L, a22 = 1.0L;
69     long double det = a11*a22 - a12*a21; // ~= -1 + 1e-20
70     long double x1 = (1.0L*a22 - 1.0L*a12) / det; // Cramer's rule
71     long double x2 = (a11*2.0L - a21*1.0L) / det;
72 
73     cout << "No pivoting solution:  ";
74     if (x_np) cout << "[" << (*x_np)[0] << ", " << (*x_np)[1] << "]\n"; else cout << "failed\n";
75 
76     cout << "Partial pivoting sol.: [" << (*x_pp)[0] << ", " << (*x_pp)[1] << "]\n";
77     cout << "Reference (long dbl):  [" << (double)x1 << ", " << (double)x2 << "]\n";
78 
79     return 0;
80 }
81

This example contrasts naive Gaussian elimination (no pivoting) with partial pivoting on a tiny 2×2 system where the (1,1) entry is extremely small. Without pivoting, we divide by a tiny number, amplifying rounding errors and risking catastrophic growth. Partial pivoting swaps rows to use a better pivot, producing a backward-stable solution that matches a high-precision reference much more closely.

Time: O(n^3)Space: O(n^2)

1	#include <bits/stdc++.h>
2	using namespace std;
3
4	double naive_sum(const vector<double>& a) {
5	double s = 0.0;
6	for (double v : a) s += v; // susceptible to rounding, esp. mixing magnitudes
7	return s;
8	}
9
10	double kahan_sum(const vector<double>& a) {
11	double s = 0.0; // running sum
12	double c = 0.0; // compensation for lost low-order bits
13	for (double v : a) {
14	double y = v - c; // first remove the small error we carried
15	double t = s + y; // add corrected value
16	c = (t - s) - y; // new rounding error stored as compensation
17	s = t;
18	}
19	return s;
20	}
21
22	int main() {
23	ios::sync_with_stdio(false);
24	cout.setf(std::ios::fixed); cout << setprecision(17);
25
26	// Construct a challenging sum: one huge number, then many ones.
27	// Adding small numbers to a huge accumulator loses precision in naive summation.
28	const int m = 500000; // half a million ones
29	vector<double> a;
30	a.reserve(m + 1);
31	a.push_back(1e8); // a large term
32	for (int i = 0; i < m; ++i) a.push_back(1.0); // many small terms
33
34	// Reference using long double (higher precision) and exact math
35	long double ref = 1e8L + static_cast<long double>(m);
36
37	double s_naive = naive_sum(a);
38	double s_kahan = kahan_sum(a);
39
40	cout << "Reference (long double): " << (double)ref << "\n";
41	cout << "Naive sum: " << s_naive << "\n";
42	cout << "Kahan sum: " << s_kahan << "\n";
43
44	cout << "Abs error naive: " << fabs(s_naive - (double)ref) << "\n";
45	cout << "Abs error Kahan: " << fabs(s_kahan - (double)ref) << "\n";
46
47	// Tip: Pairwise or sorted-by-magnitude summation can also help when Kahan is unavailable.
48	return 0;
49	}
50

1	#include <bits/stdc++.h>
2	using namespace std;
3
4	// Evaluate polynomial given coefficients [a0, a1, ..., ak] representing:
5	// P(x) = a0x^k + a1x^(k-1) + ... + ak
6	// Naive: compute powers explicitly (more operations, more rounding)
7	double poly_naive(const vector<double>& coef, double x) {
8	int k = (int)coef.size() - 1;
9	double y = 0.0;
10	for (int i = 0; i <= k; ++i) {
11	int power = k - i;
12	y += coef[i] * pow(x, power);
13	}
14	return y;
15	}
16
17	// Horner's method: nested multiplication
18	// P(x) = (...((a0x + a1)x + a2)*x + ... + ak)
19	double poly_horner(const vector<double>& coef, double x) {
20	double y = 0.0;
21	for (double a : coef) {
22	y = y * x + a;
23	}
24	return y;
25	}
26
27	int main() {
28	cout.setf(std::ios::fixed); cout << setprecision(20);
29
30	// Consider P(x) = (x - 1)^5 = x^5 - 5x^4 + 10x^3 - 10x^2 + 5x - 1
31	vector<double> coef = {1.0, -5.0, 10.0, -10.0, 5.0, -1.0};
32	double x = 1.0001; // near a root; naive evaluation suffers cancellation
33
34	// True value using high precision
35	long double xL = (long double)x;
36	long double ref = pow(xL - 1.0L, 5.0L);
37
38	double y_naive = poly_naive(coef, x);
39	double y_horner = poly_horner(coef, x);
40
41	cout << "Reference (long double): " << (double)ref << "\n";
42	cout << "Naive evaluation: " << y_naive << "\n";
43	cout << "Horner evaluation: " << y_horner << "\n";
44	cout << "\|Naive - ref\|: " << fabs(y_naive - (double)ref) << "\n";
45	cout << "\|Horner - ref\|: " << fabs(y_horner - (double)ref) << "\n";
46
47	return 0;
48	}
49

1	#include <bits/stdc++.h>
2	using namespace std;
3
4	using Matrix = vector<vector<double>>;
5	using Vec = vector<double>;
6
7	Vec back_substitute(const Matrix& U, const Vec& b) {
8	int n = (int)U.size();
9	Vec x(n, 0.0);
10	for (int i = n - 1; i >= 0; --i) {
11	double s = b[i];
12	for (int j = i + 1; j < n; ++j) s -= U[i][j] * x[j];
13	x[i] = s / U[i][i];
14	}
15	return x;
16	}
17
18	optional<Vec> solve_no_pivot(Matrix A, Vec b) {
19	int n = (int)A.size();
20	// Forward elimination without pivoting (UNSTABLE)
21	for (int k = 0; k < n - 1; ++k) {
22	if (fabs(A[k][k]) == 0.0) return nullopt; // breakdown
23	for (int i = k + 1; i < n; ++i) {
24	double m = A[i][k] / A[k][k];
25	for (int j = k; j < n; ++j) A[i][j] -= m * A[k][j];
26	b[i] -= m * b[k];
27	}
28	}
29	// Back substitution
30	return back_substitute(A, b);
31	}
32
33	optional<Vec> solve_gepp(Matrix A, Vec b) {
34	int n = (int)A.size();
35	// Gaussian Elimination with Partial Pivoting (GEPP)
36	for (int k = 0; k < n - 1; ++k) {
37	// Find pivot row with largest \|A[i][k]\|
38	int p = k;
39	double best = fabs(A[k][k]);
40	for (int i = k + 1; i < n; ++i) {
41	double val = fabs(A[i][k]);
42	if (val > best) { best = val; p = i; }
43	}
44	if (best == 0.0) return nullopt; // singular
45	if (p != k) { swap(A[p], A[k]); swap(b[p], b[k]); }
46	for (int i = k + 1; i < n; ++i) {
47	double m = A[i][k] / A[k][k];
48	A[i][k] = 0.0; // store L's entry implicitly if desired
49	for (int j = k + 1; j < n; ++j) A[i][j] -= m * A[k][j];
50	b[i] -= m * b[k];
51	}
52	}
53	return back_substitute(A, b);
54	}
55
56	int main() {
57	cout.setf(std::ios::fixed); cout << setprecision(17);
58
59	// A tiny leading pivot causes trouble without pivoting
60	Matrix A = {{1e-20, 1.0},
61	{1.0, 1.0}};
62	Vec b = {1.0, 2.0};
63
64	auto x_np = solve_no_pivot(A, b);
65	auto x_pp = solve_gepp(A, b);
66
67	// Reference solution via long double (solve analytically for 2x2)
68	long double a11 = 1e-20L, a12 = 1.0L, a21 = 1.0L, a22 = 1.0L;
69	long double det = a11a22 - a12a21; // ~= -1 + 1e-20
70	long double x1 = (1.0La22 - 1.0La12) / det; // Cramer's rule
71	long double x2 = (a112.0L - a211.0L) / det;
72
73	cout << "No pivoting solution: ";
74	if (x_np) cout << "[" << (x_np)[0] << ", " << (x_np)[1] << "]\n"; else cout << "failed\n";
75
76	cout << "Partial pivoting sol.: [" << (x_pp)[0] << ", " << (x_pp)[1] << "]\n";
77	cout << "Reference (long dbl): [" << (double)x1 << ", " << (double)x2 << "]\n";
78
79	return 0;
80	}
81