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

[Dan Sebald]

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

I've looked through the glibc library code, and as I suspected it is a big
collection of files concerning sine and cosine.  These are sort of the base
routines for which many of the other trig routines are derived (hyperbolic,
arc, etc.).  So, it is hard to pin down exactly which files and which routines
are pertinent.  To really dig in, one would have to put debug types of code in
spots and compile the library.  I really don't want to do such a thing.

My best guess of where to look in the library for 64-bit double (BTW, there is
128-bit long double code in the library, e.g., sinl, apparently not used by
Octave) would be in the library directory

glibc/sysdeps/ieee754/dbl-64

There is consideration for the various ranges of |x| for sin(x).  Ranges are x
small on the order of 1.0 (e.g., up to 2.4), various quadrants and so on.  It
uses a Taylor series polynomial approximation there.  Then there are ranges
for |x| moderately large, up to 105414350, for which a reduce_sincos_1() is
called to get some integer value that is passed into a slightly different
version of the sin/cos functions that eventually call the Taylor series
formula.  Then up to 281474976710656 is another range with a similar construct
but reduce_sincos_2().  Then there is one last range for huge numbers that
uses reduce_and_compute().

What's interesting is that those ranges roughly agree with the plot of
https://savannah.gnu.org/bugs/?48307#comment5.  For the case of 1, 10, ...,
10000000, the error shows a typical exponential-like growth.  But for 1e8 the
estimate jumps around and gets much worse in some cases.  So, I wonder if
there is a transition in there for the 105414350 boundary that the original
programmers recognized.

In any case, it is a lot of work to pick that code apart, understand what it
is doing, fix it, test it, and then (worst of all) convince some library
maintainer that it should be changed.  The best I think we could do as far as
an Octave bug goes, is to make the case the glibc sin/cos is not as accurate
as other estimates and submit a bug report to the GNU and by chance get some
feedback from someone familiar with the code.  Colin, if you could do some
exploration similar to what I did by plotting the sin(pi*x) for x=1, 10, 100,
etc. for a software program that you know doesn't use glibc, maybe we could
contrast the two.


    _______________________________________________________

Reply to this item at:

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

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