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

[Dan Sebald]

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

>> Colin was initially assuming that the relative error should be of the order
of eps.
>> However, that is only true of approximately linear functions. 
>
>
>I'm not sure that is true.

I agree that the error might have a tighter bound than O(x * eps).  Take sin(x
+- eps) and use sine formulas, then use |sin(x)| < 1 and from there find the
error limit.  (One must take into account that this is floating point, so eps
may be something that varies with the magnitude of the x itself.)  But
regardless of the case, the issue is that there isn't much we can do.

The reason for looking at sin(x) is that all the implementations we've seen in
what we've compared are calculated as 1) compute sin(x*pi) 2) divide that
quantity by x*pi.  So, it is either the sin() function or the division that is
introducing the inaccuracy (and note that I compared against the C result so
it isn't Octave doing something strange).  I'm presuming that it is not the
division because the division is a fundamental math operation and has most
likely gotten a lot of scrutiny, not least of which could be the CPU designers
if the division is done with a coprocessor.  Plus, I looked at the C code for
the library computing the sin() and it's a complex set of tweaks revolving
around Taylor series expansion based upon the magnitude of the argument.  The
designers of sin() did use a fairly low order Taylor series expansion, maybe
because they wanted to keep the function call efficient.

So, the alternatives are:

1) Improve the sin(x) in the glibc library.  This could be a lot of work.  One
has to become comfortable with that particular code in the library, figure out
a way to improve accuracy, test it for the various regions.  After several
weeks of work then comes the task of convincing the glibc maintainers that
they should modify the code.  Because of how tricky the code is, and the fact
it appears to originate from a reputable author, I think they would be
inclined to leave it as is.  I would first go to the maintainers with a good
example of there being inaccuracy in sin() and then ask if someone were to
improve the accuracy would the patch be accepted.  I just found a post with a
very good description of the sin() function

http://stackoverflow.com/a/2285277

and it sounds like I was looking at the correct code.

2) Change the Octave code to use "long double" trig functions.  That might
improve the accuracy some.  It would slow things a bit since long double
(128-bit) requires software for basic arithmetic.

3) Find a different library.  (See link above.)

4) As Lachlan suggested, a new routine for sinc() that is more accurate.  Is
there one out there in some library Octave already has linked in?  Somewhere
else?  Who knows?  Maybe if one writes out the Taylor series for sinc() it
comes out to something that converges more quickly than sin() since sinc()
does decay.  (I would attempt such a thing starting from the simple C code I
gave below because it would take just seconds to recompile the small C file
whereas in Octave it would take minutes.  Once satisfied, then the code could
be moved into Octave.)  The drawback is that this might be a slower routine
than glibc's sin() routine, and it doesn't fix the accuracy of Octave's sin()
if there is a problem.


    _______________________________________________________

Reply to this item at:

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

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