Project Nayuki


Numerically stable law of cosines

Introduction

Triangle diagram

The law of cosines: \(c^2 = a^2 + b^2 - 2ab \cos C\).

This trigonometric equation relates four values in a triangle – the three side lengths and one angle. Depending on which variables you know, you can solve the equation for a side or for an angle.

However when this formula is implemented straightforwardly in computer floating-point arithmetic, some situations will yield highly inaccurate results. When the given triangle is very long and skinny, the numerical result is inaccurate (large relative error) due to intermediate calculations that subtract very similar numbers.

This article explains how to modify the law of cosines to calculate accurately in these situations, giving mathematical derivations to alternative formulas that can be implemented. (In fact they are implemented on my JavaScript triangle solver page already.)

Solving for side

From the law of cosines formula, this should be clear: \(c = \sqrt{a^2 + b^2 - 2ab \cos C}\).

If \(a \approx b\) (they have similar magnitudes) and angle \(C \approx 0\) (close to zero), then \(\cos C \approx 1\).
Hence \(a^2 + b^2 - 2ab \cos C \:\approx\: a^2 + a^2 - 2a^2 \cos 0 \:\approx\: 2a^2 - 2a^2 \:\approx\: 0\).
Therefore \(c \approx \sqrt{0} \approx 0\). But \(c\) could be a small number like \(10^{-20}\), which is still representable in floating-point.

We can do better than this. The Taylor series expansion of the cosine function (for \(C\) in radians) yields: \(\cos C = 1 - \frac{C^2}{2!} + \frac{C^4}{4!} - \cdots\).

Then \(a^2 + b^2 - 2ab \cos C\)
\(= a^2 + b^2 - 2ab (1 + \cos C - 1)\)
\(= (a^2 + b^2 - 2ab) - 2ab (\cos C - 1)\)
\(= (a^2 - 2ab + b^2) - 2ab \left(-\frac{C^2}{2!} + \frac{C^4}{4!} - \cdots\right)\)
\(\approx (a - b)^2 + 2ab \left(\frac{C^2}{2!} - \frac{C^4}{4!}\right)\)
\(=(a - b)^2 + ab C^2 \left(1 - \frac{C^2}{12}\right)\).

Therefore our numerically stable formula is \(c \approx \sqrt{(a - b)^2 + ab C^2 \left(1 - \frac{C^2}{12}\right)}\), empirically suitable for \(|C| < 10^{-3}\).

Notes:

Solving for angle

Rearranging the law of cosines formula, we get: \(C = \cos^{-1}\! \left(\frac{a^2 + b^2 - c^2}{2ab}\right)\).

When side \(c \approx 0\) (which also implies \(a \approx b\)), we should have angle \(C \approx 0\). But \(\cos C = \frac{a^2 + b^2 - c^2}{2ab} \approx 1\), and numbers very close to \(1\) will be rounded to \(1\). This means \(\cos^{-1}\! \left(\frac{a^2 + b^2 - c^2}{2ab}\right)\) will be exactly \(0\) instead of a small number like \(10^{-20}\) which is still representable in floating-point.

By using the Taylor expansion again, we can derive the following approximation:
Assume that \(x \approx 0\) and \(y \approx 1\). Then \(y = \cos x \approx 1 - \frac{x^2}{2}\). Thus \(x \approx \sqrt{2 - 2y}\).

Now let’s avoid calculating the subformula \(a^2 + b^2 - c^2\) because it subtracts a tiny value \(c^2\) from a huge value \(a^2 + b^2\) (which is a different type of cancellation):
\(C \approx \sqrt{2 - 2\left(\frac{a^2 + b^2 - c^2}{2ab}\right)}\)
\(= \sqrt{2 - \frac{a^2 + b^2 - c^2}{ab}}\)
\(= \sqrt{\frac{2ab - a^2 - b^2 + c^2}{ab}}\)
\(= \sqrt{\frac{c^2 - (a-b)^2}{ab}}\).

Therefore our numerically stable formula is \(C \approx \sqrt{\frac{c^2 - (a-b)^2}{ab}}\), which is empirically suitable for when the arccosine argument \(\frac{a^2 + b^2 - c^2}{2ab} > 1 - 10^{-8}\) or equivalently \(\frac{c^2 - (a-b)^2}{2ab} < 10^{-8}\).


Numerical comparison

Setup notes:

Solving for side:

\(C\)\(c\) (naive)\(c\) (stable)\(c\) (precise)Naive rel errStable rel err
8.254×10−018.0217203041092455×10+78.0163093884354800×10+78.0217203041092457453923973×10+7−2.58×10−17−6.75×10−4
3.831×10−013.8077989132683828×10+73.8076838697627425×10+73.8077989132683806640465265×10+7+5.73×10−16−3.02×10−5
1.778×10−011.7759372471968371×10+71.7759347755013548×10+71.7759372471968379146724865×10+7−4.35×10−16−1.39×10−6
8.254×10−028.2516989633582728×10+68.2516984311624831×10+68.2516989633581144029010656×10+6+1.92×10−14−6.45×10−8
3.831×10−023.8309525449705063×10+63.8309525335063864×10+63.8309525449707342923629298×10+6−5.95×10−14−2.99×10−9
1.778×10−021.7782559792425837×10+61.7782559789960068×10+61.7782559792429917482831550×10+6−2.29×10−13−1.39×10−10
8.254×10−038.2540184218112822×10+58.2540184217583749×10+58.2540184218115887490007133×10+5−3.71×10−14−6.45×10−12
3.831×10−033.8311845064418390×10+53.8311845064677182×10+53.8311845064688626372997494×10+5−7.05×10−12−2.99×10−13
1.778×10−031.7782791757201680×10+51.7782791757300220×10+51.7782791757300465671183906×10+5−5.56×10−12−1.38×10−14
8.254×10−048.2540416185042341×10+48.2540416183712921×10+48.2540416183713007562273147×10+4+1.61×10−11−1.04×10−15
3.831×10−043.8311868239489442×10+43.8311868261264004×10+43.8311868261263991925039066×10+4−5.68×10−10+3.18×10−16
1.778×10−041.7782794156149928×10+41.7782794076958340×10+41.7782794076958339471918206×10+4+4.45×10−9+3.85×10−17
8.254×10−058.2540419189630975×10+38.2540418503371002×10+38.2540418503370954020360420×10+3+8.31×10−9+5.80×10−16
3.831×10−053.8311867613051704×10+33.8311868493229758×10+33.8311868493229788136203039×10+3−2.30×10−8−7.78×10−16
1.778×10−051.7782789432482184×10+31.7782794100154920×10+31.7782794100154919126759159×10+3−2.62×10−7+3.57×10−17
8.254×10−068.2540414343520229×10+28.2540418526567589×10+28.2540418526567533682467109×10+2−5.07×10−8+6.73×10−16
3.831×10−063.8311878053679385×10+23.8311868495549419×10+23.8311868495549446102570243×10+2+2.49×10−7−7.11×10−16
1.778×10−061.7783138080777533×10+21.7782794100386886×10+21.7782794100386884923399252×10+2+1.93×10−5+5.02×10−17
8.254×10−078.2534841127853397×10+18.2540418526799556×10+18.2540418526799499479107929×10+1−6.76×10−5+6.83×10−16
3.831×10−073.8314488121336034×10+13.8311868495572618×10+13.8311868495572642682234341×10+1+6.84×10−5−6.42×10−16
1.778×10−071.7776388834631177×10+11.7782794100389204×10+11.7782794100389204581365662×10+1−3.60×10−4−5.48×10−17
8.254×10−088.2462112512353212×10+08.2540418526801869×10+08.2540418526801819137074339×10+0−9.49×10−4+6.00×10−16
3.831×10−084.0000000000000000×10+03.8311868495572847×10+03.8311868495572874648030982×10+0+4.41×10−2−7.15×10−16
1.778×10−082.0000000000000000×10+01.7782794100389228×10+01.7782794100389227777945326×10+0+1.25×10−1−1.07×10−17
8.254×10−090.0000000000000000×10−18.2540418526801729×10−18.2540418526801842333654003×10−1−1.00×10+0−1.38×10−15
3.831×10−090.0000000000000000×10−13.8311868495572932×10−13.8311868495572876967688948×10−1−1.00×10+0+1.43×10−15
1.778×10−090.0000000000000000×10−11.7782794100389229×10−11.7782794100389228009911123×10−1−1.00×10+0+6.99×10−17
8.254×10−100.0000000000000000×10−28.2540418526801732×10−28.2540418526801842565619800×10−2−1.00×10+0−1.35×10−15
3.831×10−100.0000000000000000×10−23.8311868495572929×10−23.8311868495572876990885528×10−2−1.00×10+0+1.35×10−15
1.778×10−100.0000000000000000×10−21.7782794100389229×10−21.7782794100389228012230781×10−2−1.00×10+0+6.98×10−17
8.254×10−110.0000000000000000×10−38.2540418526801732×10−38.2540418526801842567939458×10−3−1.00×10+0−1.35×10−15
3.831×10−110.0000000000000000×10−33.8311868495572929×10−33.8311868495572876991117494×10−3−1.00×10+0+1.35×10−15
1.778×10−110.0000000000000000×10−31.7782794100389225×10−31.7782794100389228012253978×10−3−1.00×10+0−1.50×10−16
8.254×10−120.0000000000000000×10−48.2540418526801725×10−48.2540418526801842567962654×10−4−1.00×10+0−1.42×10−15
3.831×10−120.0000000000000000×10−43.8311868495572923×10−43.8311868495572876991119813×10−4−1.00×10+0+1.21×10−15
1.778×10−120.0000000000000000×10−41.7782794100389227×10−41.7782794100389228012254210×10−4−1.00×10+0−5.83×10−17
8.254×10−130.0000000000000000×10−58.2540418526801736×10−58.2540418526801842567962886×10−5−1.00×10+0−1.29×10−15

Observations for the data above: The naive algorithm has less relative error when approximately \(C > 10^{-3}\) radians, because the stable algorithm truncates the cosine series. But for \(C\) below this, the stable algorithm is more accurate. Moreover, the naive algorithm completely falls apart when \(C < 10^{-8}\), always giving a result of 0 instead of a tiny positive number.

Solving for angle:

\(c\)\(C\) (naive)\(C\) (stable)\(C\) (precise)Naive rel errStable rel err
8.254×10+078.5083716972718348×10−018.2540418526801895×10−018.5083716972718281127650820×10−01+7.91×10−16−2.99×10−2
3.831×10+073.8550133141801396×10−013.8311868495572848×10−013.8550133141801425715491534×10−01−7.72×10−16−6.18×10−3
1.778×10+071.7806308740167129×10−011.7782794100389226×10−011.7806308740167164643291790×10−01−2.01×10−15−1.32×10−3
8.254×10+068.2563867392267828×10−028.2540418526801898×10−028.2563867392267940604361858×10−02−1.36×10−15−2.84×10−4
3.831×10+063.8314211971419307×10−023.8311868495572853×10−023.8314211971420579196363537×10−02−3.32×10−14−6.12×10−5
1.778×10+061.7783028417606830×10−021.7782794100389229×10−021.7783028417610801004879759×10−02−2.23×10−13−1.32×10−5
8.254×10+058.2540652837349952×10−038.2540418526801732×10−038.2540652837483225589699241×10−03−1.61×10−12−2.84×10−6
3.831×10+053.8311891926739367×10−033.8311868495572929×10−033.8311891926500117494692299×10−03+6.24×10−12−6.12×10−7
1.778×10+051.7782796442869352×10−031.7782794100389228×10−031.7782796443478916540907761×10−03−3.43×10−11−1.32×10−7
8.254×10+048.2540420878769543×10−048.2540418526801736×10−048.2540420869890877114652504×10−04+1.08×10−10−2.84×10−8
3.831×10+043.8311868745189868×10−043.8311868495572923×10−043.8311868729881766356180823×10−04+4.00×10−10−6.12×10−9
1.778×10+041.7782794191261900×10−041.7782794100389227×10−041.7782794123820116645208938×10−04+3.79×10−9−1.32×10−9
8.254×10+038.2540419783121356×10−058.2540418526801898×10−058.2540418550232731135519455×10−05+1.49×10−8−2.84×10−10
3.831×10+033.8311867817738467×10−053.8311868495572856×10−053.8311868497915965846466533×10−05−1.78×10−8−6.12×10−11
1.778×10+031.7782792544822105×10−051.7782794100389229×10−051.7782794100623536897758526×10−05−8.75×10−8−1.32×10−11
8.254×10+028.2540462802021769×10−068.2540418526801902×10−068.2540418527036151453460663×10−06+5.36×10−7−2.84×10−12
3.831×10+023.8311925772508779×10−063.8311868495572849×10−063.8311868495596307879669473×10−06+1.50×10−6−6.12×10−13
1.778×10+021.7783023543096177×10−061.7782794100389229×10−061.7782794100391571101109173×10−06+1.29×10−5−1.32×10−13
8.254×10+018.2536831045905480×10−078.2540418526801902×10−078.2540418526804185656817849×10−07−4.35×10−5−2.77×10−14
3.831×10+013.8310766614027379×10−073.8311868495572850×10−073.8311868495573111300005333×10−07−2.88×10−5−6.82×10−15
1.778×10+011.7756782901008428×10−071.7782794100389230×10−071.7782794100389251443142762×10−07−1.46×10−3−1.23×10−15
8.254×10+008.2966154259890665×10−088.2540418526801902×10−088.2540418526801865998851438×10−08+5.16×10−3+4.35×10−16
3.831×10+003.9424766765007238×10−083.8311868495572849×10−083.8311868495572879334208692×10−08+2.90×10−2−8.02×10−16
1.778×10+002.1073424255447017×10−081.7782794100389228×10−081.7782794100389228246563097×10−08+1.85×10−1−1.80×10−18
8.254×10−010.0000000000000000×10−098.2540418526801727×10−098.2540418526801842802271774×10−09−1.00×10+0−1.41×10−15
3.831×10−010.0000000000000000×10−093.8311868495572935×10−093.8311868495572877014550725×10−09−1.00×10+0+1.50×10−15
1.778×10−010.0000000000000000×10−091.7782794100389230×10−091.7782794100389228014597301×10−09−1.00×10+0+1.28×10−16
8.254×10−020.0000000000000000×10−108.2540418526801727×10−108.2540418526801842570305977×10−10−1.00×10+0−1.41×10−15
3.831×10−020.0000000000000000×10−103.8311868495572929×10−103.8311868495572876991354146×10−10−1.00×10+0+1.37×10−15
1.778×10−020.0000000000000000×10−101.7782794100389229×10−101.7782794100389228012277643×10−10−1.00×10+0+4.04×10−17
8.254×10−030.0000000000000000×10−118.2540418526801737×10−118.2540418526801842567986319×10−11−1.00×10+0−1.28×10−15
3.831×10−030.0000000000000000×10−113.8311868495572927×10−113.8311868495572876991122180×10−11−1.00×10+0+1.30×10−15
1.778×10−030.0000000000000000×10−111.7782794100389227×10−111.7782794100389228012254446×10−11−1.00×10+0−6.86×10−17
8.254×10−040.0000000000000000×10−128.2540418526801727×10−128.2540418526801842567963123×10−12−1.00×10+0−1.40×10−15
3.831×10−040.0000000000000000×10−123.8311868495572933×10−123.8311868495572876991119860×10−12−1.00×10+0+1.47×10−15
1.778×10−040.0000000000000000×10−121.7782794100389226×10−121.7782794100389228012254214×10−12−1.00×10+0−9.13×10−17
8.254×10−050.0000000000000000×10−138.2540418526801739×10−138.2540418526801842567962891×10−13−1.00×10+0−1.25×10−15

Observations for the data above: Let \(z = \frac{a^2 + b^2 - c^2}{2ab}\), which is the argument to the arccosine function. The naive algorithm has less relative error when \(z < 1 - 10^{-8}\), because the stable algorithm truncates the cosine series. But for \(z\) above this (i.e. very close to 1), the stable algorithm is more accurate. Moreover, the naive algorithm completely falls apart when \(z > 1 - 10^{-16.5}\), always giving a result of 0 instead of a tiny positive number because the argument \(z\) is closer to 1 than the machine epsilon of double precision.