[Octave-bug-tracker] [bug #48307] sinc loses precision for large arguments

[Dan Sebald]

Follow-up Comment #19, bug #48307 (project octave):

I've attached some C code that explores various Taylor series solutions for
sin() with double and long double precision.  (Then to get sinc() do division
explicitely.)  I implemented the series with the first ten terms, because the
1/x! coefficient over/underflows for any higher terms.  Is ten terms way more
than enough?  Don't know, would have to test.  To extend the series to twenty
terms I used the Gamma function (tgamma) to compute factorial in floating
point form.  (One could pre-compute those Gamma values to speed up such an
approach.)  Then, for the ultimate test, I created a long double version of
the Gamma coefficient series.

Now, how about the Taylor (power) series for sin(x)/x?  Well, on a piece of
scratch paper, I get:


f    =  sin(x)/x
f'   =  cos(x)/x -   sin(x)/x^2
f''  = -sin(x)/x - 2 cos(x)/x^2 + 2 sin(x)/x^3
f''' = -cos(x)/x + 3 sin(x)/x^2 + 4 cos(x)/x^3 - 6 sin(x)/x^4


which just doesn't seem to lead anywhere convenient.  Expansion about x - x_0
where x_0 is some factor of 2 times pi gets rid of the sine and cosine of
course, but the formulas still aren't convenient, and I expect it can't be any
more accurate or converge any faster because that 1/x term hangs around at the
front.

Anyway, the results I am seeing with the "custom" sin() function are:


sebald@ ~/octave/bug_48307 $ gcc -o sinc sinc.c -lm
sebald@ ~/octave/bug_48307 $ ./sinc 
DOUBLE
M_PI = 3.141592653589793115997963468544185161590576171875
10000001/3 = 3333333.6666666665114462375640869140625
sin(pi*x)/(pi*x) =
-8.26993260956150473721458331322065049562297645024955272674560546875e-08
table_sinc(x) =
-8.26993260666192699379488145923489117450344565440900623798370361328125e-08
tgamma_sinc(x) =
-8.26993260666192699379488145923489117450344565440900623798370361328125e-08
LONG DOUBLE
M_PIl = 3.14159265358979323851280895940618620443274267017841339111328125
10000001/3 = 3333333.666666666666742457891814410686492919921875
sin(pi*x)/(pi*x) =
-8.2699326043331140428346591282514346722860854033143596097943373024463653564453125e-08
tgamma_sincl(x) =
-8.2699326043324834416946793276034964538696858671329437129315920174121856689453125e-08


>From that result I take away that we can in fact improve on the double
implementation:


sin(pi*x)/(pi*x) = -8.269932609561[snip]e-08
table_sinc(x)    = -8.269932606661[snip]e-08


Is it an adequate improvement in exchange for speed?  Or do we need to go to a
long double implementation


tgamma_sincl(x) =
-8.2699326043324834416946793276034964538696858671329437129315920174121856689453125e-08


which looks to be along the lines of Maple.

So I'm assuming Maple is using long double math and Octave is using trig
functions that toss out a bit of accuracy for increased speed.

Someone feel free to explore further to see just exactly how many terms are
needed, whether the algorithm is stable over various regions, so on.  Should
we have a fix for all the trig functions?  Or just the sin() computation
associated with sinc()?  Can we be content with just using the math libraries
sinl() routine?



(file #37708)
    _______________________________________________________

Additional Item Attachment:

File name: sinc.c                         Size:2 KB


    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?48307>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/