Stable Reformulations Near Removable Singularities

Why Equivalent Formulas Can Behave Differently

Two algebraically equivalent formulas can give very different numerical results when evaluated on a computer. The issue is especially common near a removable singularity: a formula is undefined at an input, but its limit there exists and is finite. Direct substitution may involve subtracting nearly equal quantities, so rounding error can overwhelm the small difference we want to compute.

A reliable strategy is to replace the original expression with an equivalent one whose calculation does not discard most of its significant digits. Taylor expansions are particularly useful because they expose the small leading terms directly. Algebraic identities can help as well, and numerical libraries often provide specialized functions such as expm1 for computing \(e^x-1\) accurately when \(x\) is small.

Definition: A formula \(F(x)\) has a removable singularity at \(x=a\) if it is undefined at \(a\), but \(\lim_{x\to a}F(x)=L\) exists and is finite. Defining \(F(a)=L\) gives a continuous extension at \(a\). A stable reformulation is an equivalent expression or method of evaluation that avoids unnecessary loss of numerical accuracy.

The limiting value and the numerical calculation are related, but they are not the same question. The limit identifies the value of the continuous extension. A stable reformulation addresses how to compute values near the singularity without losing accuracy.

What Cancellation Does to a Calculation

Suppose a computation forms \(a-b\), where \(a\) and \(b\) are close. Each value may already have a small rounding error, but their difference can be much smaller than either value. Subtracting them can therefore leave very few reliable digits. The relative effect is especially large when \(|a-b|\) is tiny compared with \(|a|+|b|\).

For example, \(1-\cos x\) is of order \(x^2\) near zero, even though both \(1\) and \(\cos x\) are close to \(1\). A computer may represent \(\cos x\) as exactly \(1\) after rounding, making the computed difference zero. The exact expression is positive for small nonzero \(x\), so the rounded calculation has lost the quantity of interest.

In floating-point arithmetic, one often describes cancellation by saying that subtraction can amplify relative errors in its inputs. This does not mean every subtraction of close numbers is inaccurate: the concern is that the result is small while the input errors are not comparably small. In these examples, the exact result tends to zero faster than the individual terms being subtracted.

A Stable Form for \(1-\cos x\)

The identity \(1-\cos x=2\sin^2(x/2)\) replaces subtraction of two values near \(1\) with a square of a small sine. For nonzero \(x\), dividing by \(x^2\) gives an equivalent expression. The Taylor expansion of cosine gives another useful form.

\frac{1-\cos x}{x^2} = 2\left(\frac{\sin(x/2)}{x}\right)^2 = \frac12-\frac{x^2}{24}+\frac{x^4}{720}-\frac{x^6}{40320}+\cdots .

The half-angle version is useful when a reliable sine function is available. The series is useful when a polynomial evaluation is preferred and the input is small enough that a suitable remainder bound is available. In either case, the value at \(x=0\) for the continuous extension is \(1/2\), not the value of the original quotient, which is undefined there.

Worked Example: Evaluating a Cosine Quotient at \(x=0.1\)

Consider \(F(x)=(1-\cos x)/x^2\). Direct evaluation uses the difference of \(1\) and \(\cos(0.1)\), which are close. Instead use the Taylor expansion through the \(x^6\) term:

P_6(0.1) = \frac12-\frac{0.1^2}{24}+\frac{0.1^4}{720}-\frac{0.1^6}{40320}.

The individual terms are \(0.5\), \(0.0004166666666666667\), \(0.0000001388888888889\), and \(0.0000000000248015873\). Thus

P_6(0.1)\approx 0.4995834721974206349.

The alternating cosine series has decreasing term magnitudes here. Its first omitted term in the quotient is \(0.1^8/10!=2.7557319\times10^{-15}\), so the absolute truncation error is at most \(2.76\times10^{-15}\). That bound does not certify every digit displayed in \(P_6\). Including the next term gives approximately \(0.4995834721974233907\); the following term is less than \(2.1\times10^{-19}\). Therefore the correctly rounded value to 15 decimal places is \(0.499583472197423\).

For a check, multiplying the series for \(F(x)\) by \(x^2\) gives \(x^2/2-x^4/24+x^6/720-\cdots\), which is the Taylor series for \(1-\cos x\). The half-angle identity provides a second equivalent check: at \(x=0.1\), \(2(\sin(0.05)/0.1)^2\) represents the same quotient without the subtraction of near-equal values.

Using a Dedicated Exponential Function

The quotient \((e^x-1)/x\) has a removable singularity at zero. Directly computing \(e^x-1\) can lose precision when \(x\) is small, since \(e^x\) is then close to \(1\). Many numerical libraries provide expm1(x), which is designed to compute \(e^x-1\) accurately in this regime. When such a function is available, use expm1(x)/x rather than forming \(e^x\) and then subtracting \(1\).

The Taylor expansion also gives an explicit stable method. For \(x\ne0\), the exponential series yields

\frac{e^x-1}{x} = 1+\frac{x}{2!}+\frac{x^2}{3!}+\frac{x^3}{4!}+\cdots .

Its continuous extension at zero has value \(1\). Because the constant term is visible in the series, this method does not require subtracting two numbers close to \(1\).

Worked Example: Computing \((e^x-1)/x\) at \(x=0.1\)

Use the series through the term \(x^8/9!\):

S_8 = 1+\frac{0.1}{2!}+\frac{0.1^2}{3!} +\frac{0.1^3}{4!}+\cdots+\frac{0.1^8}{9!}.

The successive terms from \(x^0\) through \(x^8\) are approximately \(1\), \(0.05\), \(0.0016666666666667\), \(0.0000416666666667\), \(0.0000008333333333\), \(0.0000000138888889\), \(0.0000000001984127\), \(0.0000000000024802\), and \(0.0000000000000276\). Their sum is approximately \(1.0517091807564761\).

The first omitted term is \(0.1^9/10!=2.7557319\times10^{-16}\). Each subsequent term is at most \(0.1/11\) times the preceding term, so the full omitted tail is less than \(2.79\times10^{-16}\). The series therefore supports rounding the result to 13 decimal places: \(1.0517091807565\).

To check the formula, multiplying the infinite series by \(x\) gives \(x+x^2/2!+x^3/3!+\cdots=e^x-1\), so division by \(x\) recovers the original quotient for \(x\ne0\). The limiting value at zero also agrees with the series’ constant term, \(1\).

When the Numerator Has More Than One Cancellation

Some expressions have successive cancellations. For example, in \(e^x-1-x\), the constant and linear terms of the exponential expansion are removed. The remaining expression is of order \(x^2\). Computing it as written can therefore lose accuracy twice: first when forming \(e^x-1\), and again when subtracting \(x\).

Instead, start from the exponential series and remove the terms that cancel exactly. Dividing the remaining series by \(x^2\) gives a formula with a finite limit:

\frac{e^x-1-x}{x^2} = \frac12+\frac{x}{3!}+\frac{x^2}{4!}+\frac{x^3}{5!}+\cdots .

This expansion is appropriate near zero. For values farther from zero, the number of terms needed for a given accuracy can increase; a Taylor polynomial should not be used without considering its error at the input in question.

Worked Example: Evaluating a Doubly Cancellation-Prone Quotient

Evaluate \(G(x)=(e^x-1-x)/x^2\) at \(x=0.1\) using the series through \(x^4/6!\):

T_4 = \frac12+\frac{0.1}{3!}+\frac{0.1^2}{4!} +\frac{0.1^3}{5!}+\frac{0.1^4}{6!}.

The terms are \(0.5\), \(0.0166666666666667\), \(0.0004166666666667\), \(0.0000083333333333\), and \(0.0000001388888889\). Their sum is approximately \(0.5170918055555556\). The next term is \(0.1^5/7!=1.984126984\times10^{-9}\). The ratio of each following term to the preceding one is at most \(0.1/8\), so the remaining tail is less than \(1.984126984\times10^{-9}/(1-0.1/8)<2.01\times10^{-9}\).

Thus \(T_4\) is a controlled approximation, but it does not support many decimal places. To check the reformulation, multiply its full series by \(x^2\): the result is \(x^2/2+x^3/3!+x^4/4!+\cdots=e^x-1-x\). The first omitted term also confirms the positive direction of the error at \(x=0.1\).

Why the Reformulations Work

The central reason is that Taylor expansions separate the small leading behavior from terms that cancel. For \((1-\cos x)/x^2\), the numerator starts with \(x^2/2\); dividing by \(x^2\) leaves a nonzero constant. For \((e^x-1)/x\), the numerator starts with \(x\); dividing by \(x\) leaves a constant. These leading terms explain both the finite limits and the appropriate forms for numerical evaluation.

The half-angle identity follows directly from the cosine double-angle formula \(\cos x=1-2\sin^2(x/2)\). Rearranging gives \(1-\cos x=2\sin^2(x/2)\). This derivation is useful computationally because the sine is small near zero, while its square directly represents the small difference we need. For \((e^x-1)/x\) and \((e^x-1-x)/x^2\), the corresponding stable series follow by subtracting the terms that cancel in the exponential Taylor expansion, then dividing each remaining term by the appropriate power of \(x\).

Practical rule: If a formula near a point subtracts quantities that approach the same value, look for an algebraically equivalent identity or a Taylor expansion that exposes the first nonzero term. Use a remainder bound to decide how many digits the chosen approximation supports.

A Common Mistake: Trusting a Small Direct Difference

Consider evaluating \((1-\cos(10^{-8}))/(10^{-8})^2\) using ordinary double-precision arithmetic. The exact value is close to \(1/2\). But \(\cos(10^{-8})\) differs from \(1\) by only about \(5\times10^{-17}\), which is smaller than the spacing between representable numbers near \(1\) in double precision. A typical computation rounds \(\cos(10^{-8})\) to \(1\), then calculates \(1-1=0\). It returns zero for the quotient instead of a value close to \(1/2\).

The correction is to avoid that subtraction. The half-angle formula gives \(2(\sin(5\times10^{-9})/10^{-8})^2\), which is close to \(2(1/2)^2=1/2\). A Taylor evaluation, \(1/2-x^2/24+x^4/720-\cdots\), also begins with the desired value directly. In exact arithmetic the original and revised formulas agree for \(x\ne0\); in floating-point arithmetic their accuracy can differ sharply.

A related mistake is to keep adding Taylor terms without checking the error. A stable formula does not automatically mean a sufficiently accurate approximation. The input size, number of retained terms, and a valid remainder estimate all matter. Round only after the error bound is small enough for the requested precision.

Choosing a Method

Situation	Useful approach
The expression is \(e^x-1\) with small \(x\)	Use a library’s `expm1(x)`, if available.
The expression is \(1-\cos x\) near zero	Use \(2\sin^2(x/2)\) or a cosine Taylor expansion.
Several low-order terms cancel	Expand first, remove the cancelling terms, and then divide.
A specified number of digits is needed	Bound the omitted terms before rounding the result.

Stable reformulations matter because numerical calculations are not judged only by algebraic equivalence. The arrangement of operations affects how rounding errors enter the result. Near a removable singularity, a good representation preserves the finite limiting behavior instead of allowing cancellation to hide it.

Check Your Understanding

Use the ideas in this tutorial to answer the following questions.

Why can direct evaluation of \(1-\cos x\) lose accuracy when \(x\) is small?
Use the half-angle identity to rewrite \((1-\cos x)/x^2\) for \(x\ne0\). What is the value of its continuous extension at zero?
Write the first four terms of the Taylor expansion for \((e^x-1)/x\).
Why is subtracting \(1\) and then \(x\) from \(e^x\) especially risky when \(x\) is small?
What does a remainder bound tell you before you round a Taylor approximation?