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

[Dan Sebald]

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

Attached is a second version of the sinc tests in which I placed the double
value of the input, i.e.,


    double x = 10000001;
    x /= 3;


in addition to the long double value version of the above number into the
long-double version of the tgamma-based routine.  So here is the contrast:


DOUBLE
sin(pi*x)/(pi*x)   = -8.26993260956150473721458331322065e-08
tgamma_sinc(x)     = -8.26993260666192699379488145923489e-08
LONG DOUBLE
sin(pi*xl)/(pi*xl) = -8.26993260433311404283465912825143e-08
tgamma_sincl(xl)   = -8.26993260433248344169467932760349e-08
tgamma_sincl(x)    = -8.26993260666192689621341857211990e-08


Octave produces the first result in the list above.  The second result is the
custom double version of sinc that improves on the estimate in the
library--note that it only cuts the "error" in half, roughly, so the
improvement isn't the orders of magnitude needed to achieve 3 times eps.  The
fifth result is what happens when the double floating point value is passed
into the higher resolution long-double routine--note how close it is to the
custom double routine result.

So, yes, it seems that most of the loss of accuracy is a result of the fact
the input (i.e., Octave's storage class) is double and not long double.  I
don't know what can be done without a total overhaul of Octave's floating
point class.  Maybe symbolic can be made to work better for division problems
where the dividend and divisor are integers...or maybe that is already
happening (see below).

Let's go back to the original discrepancy that raised the issue.  It involved
the difference between the non-symbolic and symbolic result:


>> x = sym(10000001)/3
x = (sym) 10000001/3
>> d = double(x)              % large non-integer
d =    3.3333e+06
>> P1 = double(sin(pi*x) / (pi*x))
P1 =   -8.2699e-08
>> P2 = double(sinc(x))
P2 =   -8.2699e-08

>> Q = sinc(d)             % this is the value we're checking
Q =   -8.2699e-08

>> (P1-Q)/P1               % relative errors
ans =   -6.3216e-10
>> (P2-Q)/P2
ans =   -6.3222e-10


The relative errors only suggest that symbolic and non-symbolic values are
noticeably different.  Why?  In order to find out, I run the symbolic package
and get a "sympy version too old" error.  I went to upgrade, but then SymPy
wants a version 19 mpmath where I only have version 18.  I stop there.  So, if
someone with symbolic package can repeat the above test, but print out the
full length of the result with


sprintf("%0.100g", <both results>)


then we can determine how the symbolic result compares to the
higher-resolution results.  Maybe it is the case that symbolic IS using long
double because it is actually being done with the SymPi library.  If the
symbolic library casts 10000001 and 3 to long-doubles then does the
computation, it would compare to the high-resolution results of the C program.
 All conjecture, so someone should explore that to answer one of the original
issues.


    _______________________________________________________

Reply to this item at:

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

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