Invariant Properties

  • rss
  • Home

Calculating sin(x) / x

Bear Giles | April 20, 2014

I’ve been busy with my new job so I haven’t been able to spend time on the major themes – time to catch up on some one-off topics.

Every so often I see a blog discussion of calculating sin(x) / x as an example of the horrors of numerical methods. It’s actually a good example of why most if not all computer science programs require a year of calculus. 🙂

(I may be biased since I my undergraduate degrees are in math and physics, not computer science, so this is all very basic material for me.)

L’Hospital’s Rule

There’s an obvious place where sin(x) / x (also known as sinc(x)) is problematic: x = 0. It appears to be division-by-zero and in fact it’s usually explicitly defined as equal to 1 at that point.

Can we do better (for varying definitions of “better”)?

We can start with L’Hospital’s Rule:

limx→0 f(x)/g(x) = limx→0 f'(x)/g'(x)

In this case

limx→0 sin(x)/x = limx→0 cos(x)/1 = cos(0) = 1

This does not mean that sinc(x) is equal to 1 at x = 0. It just means that the closer you get to zero, from either side, the closer the value gets to 1.

The existence of the limit also does not mean that sinc(x) doesn’t do funny stuff as you approach zero. In this case we know it’s well-behaved but in general you’ll want to double-check that.

Taylor Series Expansion

We have an additional tool available: the Taylor series around x = 0. Specifically we know that the Taylor series of sin(x) around x = 0 is:

sin(x) = x – x3/3! + x5/5! – x7/7! + …

We can immediately use term-wise division by x to get

sinc(x) = 1 – x2/3! + x4/5! – x6/7! + …

With sufficiently small x we can use the first-order approximation:

sinc(x) = 1 – x2/6 (first-order approximation)

Sinc(x) is a close approximation of a simple parabola near the origin.

Computing sinc(x)

Finally we have the question of what happens when we try to compute sin(x)/x using a computer.

Floating point numbers are almost always represented using the IEEE 754-2008 representation. For a double that’s 52(+1) bits. (The +1 is implicit.)

I don’t know if sin(x) is computed in hardware, using the Taylor series, using CORDIC functions, or by consulting the Oracle at Delphi. It doesn’t matter since we know, from the Taylor series, that the limited precision of floating point numbers means that at some value |x| < ε the higher terms of the Taylor series will drop off and

sin(x) = x(1 – x2/3!) (computed)

That second term will drop out when x2/3! < 2-53, or roughly |x| < 2-25. Beyond that point the computed sin(x) is equal to x and sinc(x) will be equal to 1 until you hit x = 0.

We can avoid special handling of x = 0 by noting that we’ve already factored the calculation above:

sinc(x) = sin(x) / x = 1 – x2/6 (computed)

If we cut over to this calculation when the next term in the expansion drops out (x4/5! < 2-53, or roughly |x| < 2-10) then we don’t have to handle x = 0 separately.

It’s debatable whether it’s cleaner to treat x = 0 or |x| < 2-10 as the special case but it shows how we can use basic analysis to turn an apparent problem into “no big deal”.

It’s important to note that the fact that we’re dividing by a small number is irrelevant. What’s critical is the precision of numbers involved. A small number is not a problem a priori, a small number due to subtraction of two nearly identical values is.

Calculating (x – sin(x)) / x

This brings up a variant of the earlier problem. x and sin(x) are nearly identical for sufficiently small x. This can cause problems due to loss of precision.

Our first step is to verify that there is a well-defined limit at 0. We could use L’Hospital’s Rule or simply observe:

(x – sin(x)) / x = (x/x) – (sin(x)/x) = 1 – sinc(x)

We can now use the Taylor series to see

(x – sin(x)) / x = 1 – sinc(x) = 1 – (1 – x2/3! + …) = x2/3! – …

This calculation will return 0 when |x| < 2-511. It is 2-511, not 2-25, since the maximum precision of an IEEE double is only 52(+1) bits but the minimum exponent is -1022. That means there’s a broad region (in terms of orders of magnitude) where the direct subtraction will return zero but the Taylor series will return a small but nonzero value.

Conclusion

This is a simple example but it illustrates a common problem in numeric methods. Some equations look like they’ll be a problem, e.g., division by zero or subtracting nearly identical values, but rephrasing the problem can produce much better results.

A fair question is whether this matters. Why not make use the most direct approach since the numbers involved are so small? We need to ask two questions:

  • Is this value used in subsequent calculations or an interation? If so the loss of few bits of precision can result in significant errors downstream.
  • Is the most precise value required? Are there contractual or legal requirements?

The bottom line is that this is not our call to make. We can and should ask for clarification before investing a significant amount of time into an analysis but these examples should only take a few hours for someone comfortable writing numerical code.

Footnote: Quadratic Equations

I should also mention quadratic equation. We all know the form taught in high school but there are two problematic scenarios:

First, when b2 ≅ 4ac then the equation involves subtraction of two nearly equal values and there will be a significant loss of precision.

Second, when b2 ≫ 4ac then the contribution by the second will be diminished or lost due to the limited precision of a floating point number. In this case where will be full precision but quantization noise.

In both cases there are alternative formulations that can give better results.

Comments
No Comments »
Categories
math
Tags
L'Hospital, sin, sinc, Taylor expansion
Comments rss Comments rss
Trackback Trackback

Archives

  • May 2020 (1)
  • March 2019 (1)
  • August 2018 (1)
  • May 2018 (1)
  • February 2018 (1)
  • November 2017 (4)
  • January 2017 (3)
  • June 2016 (1)
  • May 2016 (1)
  • April 2016 (2)
  • March 2016 (1)
  • February 2016 (3)
  • January 2016 (6)
  • December 2015 (2)
  • November 2015 (3)
  • October 2015 (2)
  • August 2015 (4)
  • July 2015 (2)
  • June 2015 (2)
  • January 2015 (1)
  • December 2014 (6)
  • October 2014 (1)
  • September 2014 (2)
  • August 2014 (1)
  • July 2014 (1)
  • June 2014 (2)
  • May 2014 (2)
  • April 2014 (1)
  • March 2014 (1)
  • February 2014 (3)
  • January 2014 (6)
  • December 2013 (13)
  • November 2013 (6)
  • October 2013 (3)
  • September 2013 (2)
  • August 2013 (5)
  • June 2013 (1)
  • May 2013 (2)
  • March 2013 (1)
  • November 2012 (1)
  • October 2012 (3)
  • September 2012 (2)
  • May 2012 (6)
  • January 2012 (2)
  • December 2011 (12)
  • July 2011 (1)
  • June 2011 (2)
  • May 2011 (5)
  • April 2011 (6)
  • March 2011 (4)
  • February 2011 (3)
  • October 2010 (6)
  • September 2010 (8)

Recent Posts

  • 8-bit Breadboard Computer: Good Encapsulation!
  • Where are all the posts?
  • Better Ad Blocking Through Pi-Hole and Local Caching
  • The difference between APIs and SPIs
  • Hadoop: User Impersonation with Kerberos Authentication

Meta

  • Log in
  • Entries RSS
  • Comments RSS
  • WordPress.org

Pages

  • About Me
  • Notebook: Common XML Tasks
  • Notebook: Database/Webapp Security
  • Notebook: Development Tips

Syndication

Java Code Geeks

Know Your Rights

Support Bloggers' Rights
Demand Your dotRIGHTS

Security

  • Dark Reading
  • Krebs On Security Krebs On Security
  • Naked Security Naked Security
  • Schneier on Security Schneier on Security
  • TaoSecurity TaoSecurity

Politics

  • ACLU ACLU
  • EFF EFF

News

  • Ars technica Ars technica
  • Kevin Drum at Mother Jones Kevin Drum at Mother Jones
  • Raw Story Raw Story
  • Tech Dirt Tech Dirt
  • Vice Vice

Spam Blocked

53,313 spam blocked by Akismet
rss Comments rss valid xhtml 1.1 design by jide powered by Wordpress get firefox