Numerical issues with sin()/cos()

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
8 messages Options
Reply | Threaded
Open this post in threaded view
|

Numerical issues with sin()/cos()

vlad ionescu
Hello everyone

I am trying to modify fir1.m to create FIRs using the correct sum of
sines or cosines. I ended up making a new .m file, but when I got to
the "%! assert()" lines, I stumbled across some numerical issues. The
.m file was built after a wxMaxima file, which was used to generate
the default results, against which assert() could check, and I noticed
that both sin(), and cos(), seem to not produce (anti-)symmetrical
results. The differences can be more than 1e-15. "format long" is
enabled, but it's the same. Example, lowpass, type I, 33rd order,
wc=0.37, in wxMaxima and Octave (in this order):

0.006248855546642409" "
-0.0151891541716538" "
-0.02000749470712845" "
3.703551396872747*10^-4" "
0.02352639910729612" "
0.01987710435709456" "
-0.0107156702183056" "
-0.03346910762203006" "
-0.01647493965701502" "
0.02756343124616879" "
0.0468059024239567" "
0.006350821104335874" "
-0.06144311077699433" "
-0.07272786646455694" "
0.02972318687964276" "
0.2090466916579446" "
0.3495187814185787" "
0.3495187814185787" "
0.2090466916579446" "
0.02972318687964276" "
-0.07272786646455694" "
-0.06144311077699433" "
0.006350821104335874" "
0.0468059024239567" "
0.02756343124616879" "
-0.01647493965701502" "
-0.03346910762203006" "
-0.0107156702183056" "
0.01987710435709456" "
0.02352639910729612" "
3.703551396872747*10^-4" "
-0.02000749470712845" "
-0.0151891541716538" "
0.006248855546642409

6.24885554664247e-03
  -1.51891541716537e-02
  -2.00074947071285e-02
   3.70355139687233e-04
   2.35263991072961e-02
   1.98771043570946e-02
  -1.07156702183055e-02
  -3.34691076220301e-02
  -1.64749396570151e-02
   2.75634312461687e-02
   4.68059024239567e-02
   6.35082110433597e-03
  -6.14431107769943e-02
  -7.27278664645570e-02
   2.97231868796425e-02
   2.09046691657944e-01
   3.49518781418579e-01
   3.49518781418577e-01
   2.09046691657945e-01
   2.97231868796429e-02
  -7.27278664645570e-02
  -6.14431107769943e-02
   6.35082110433577e-03
   4.68059024239567e-02
   2.75634312461688e-02
  -1.64749396570150e-02
  -3.34691076220301e-02
  -1.07156702183056e-02
   1.98771043570945e-02
   2.35263991072961e-02
   3.70355139687233e-04
  -2.00074947071284e-02
  -1.51891541716538e-02
   6.24885554664247e-03

The simplest generator for this would be (sinc() works the same):

N=33; M=N/2; n=-M:1:M; h=sin(0.37*pi*n)./(pi*n); h(ceil((N+1)/2)=0.37;

Note that the code in Octave directly follows the one in wxMaxima, and
both use the range -M:1:M. There is no cheating in wxMaxima, so
sin(-1)=-sin(1), while in Octave it doesn't seem to be so. What's
more, I am not using big floats in wxMaxima, only default numbers,
float() when necessary.

The assert() that I was aiming for was 1e-15. In the case above, it
fails at h(18):

Location  |  Observed  |  Expected  |  Reason
    (18)       0.34952      0.34952      Abs err 1.6653e-15 exceeds tol 1e-15

Is this a known behaviour? Should I just relax the tolerance to 1e-14,
or less? Or am I doing something wrong? I'll attach the .m script as
it is, so far.

Vlad



fir.m (11K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Numerical issues with sin()/cos()

vlad ionescu
I have just seen a previous question in the mailing lists, "If pi is
so accurate why it's not producing that accurate result". There
probably are more of the same. Just to clarufy (if there's a
misunderstanding), I am not questioning the inherent floating point
errors, but their consistency. In Octave, sin(-x) doesn't seem to be
equal to -sin(x). This is my main point. As you can see, the the order
of the FIR is odd, so it has two midpoints of equal values coming from
sin(x)/x and sin(-x)/(-x), thus, for all intents and purposes,
sin(x)/x, yet they are not equal. That they differ in Octave and
wxMaxima, so be it, to each their own (numerical evaluation), but why
the two are different whithin Octave is my concern.

Vlad


Reply | Threaded
Open this post in threaded view
|

Re: Numerical issues with sin()/cos()

Mike Miller-4
In reply to this post by vlad ionescu
On Mon, Mar 26, 2018 at 16:01:04 +0000, vlad ionescu wrote:
> The simplest generator for this would be (sinc() works the same):
>
> N=33; M=N/2; n=-M:1:M; h=sin(0.37*pi*n)./(pi*n); h(ceil((N+1)/2)=0.37;
>
> Note that the code in Octave directly follows the one in wxMaxima, and
> both use the range -M:1:M. There is no cheating in wxMaxima, so
> sin(-1)=-sin(1), while in Octave it doesn't seem to be so.

Setting the fir function aside, I don't see any problems with this
particular example

    >> N = 33; M = N/2; n = (-M:M)'; h = sin (0.37 * pi * n) ./ (pi * n);
    >> assert (h, flipud (h))
    >> assert (h, x, 1e-15)

It doesn't look to me like the errors are in the sin function.

FWIW the h I obtain above is not equal to the h I get from the attached
fir function. The error is worst at index 18, the same as the error you
are seeing.

--
mike



signature.asc (849 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Numerical issues with sin()/cos()

Mike Miller-4
On Tue, Mar 27, 2018 at 11:29:54 -0700, Mike Miller wrote:
> It doesn't look to me like the errors are in the sin function.

In fact, I think what you want to look at are the differences between
the following v1 and v2

    N = 33; M = N/2; n = -M:M;
    v1 = sin (0.37 * pi * n) ./ (pi * n);
    v2 = sin (0.37 * pi * [n]) ./ (pi * [n]);

I think what is hurting your precision here is the lazy evaluation of
ranges in Octave.

--
mike



signature.asc (849 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Numerical issues with sin()/cos()

Gordon Haverland
I ran across a fairly recent (2008) paper on some of the "dark arts" of
finding good approximations to mathematical functions.  It does in part
look at sin(x).

The paper has a preferred citation:

Florent De Dinechin, Christoph Lauter. Optimizing polynomials for
floating-point implementation. 12 pages. 2008. <ensl-00260563>

I believe the URL is:

https://hal-ens-lyon.archives-ouvertes.fr/ensl-00260563

Gord



Reply | Threaded
Open this post in threaded view
|

Re: Numerical issues with sin()/cos()

Sergei Steshenko





________________________________
From: Gordon Haverland <[hidden email]>
To: [hidden email]
Sent: Wednesday, March 28, 2018 1:53 AM
Subject: Re: Numerical issues with sin()/cos()



I ran across a fairly recent (2008) paper on some of the "dark arts" of
finding good approximations to mathematical functions.  It does in part
look at sin(x).

The paper has a preferred citation:

Florent De Dinechin, Christoph Lauter. Optimizing polynomials for
floating-point implementation. 12 pages. 2008. <ensl-00260563>

I believe the URL is:

https://hal-ens-lyon.archives-ouvertes.fr/ensl-00260563


Gord

-----------------------------------------------------

Since Intel 80387 (see https://en.wikipedia.org/wiki/X86_instruction_listings ->
https://en.wikipedia.org/wiki/X86_instruction_listings#Added_with_80387 ) sine and cosine functions are part of the instruction set, so if one wants fast implementation, he's limited by what his x86 (if it's the case) provides.

--Sergei.


Reply | Threaded
Open this post in threaded view
|

Re: Numerical issues with sin()/cos()

vlad ionescu
In reply to this post by vlad ionescu
Hello

I am sorry, despite being able to see all the (current) 4 replies in
the gnu.org list, Gmail still didn't bother to let me know, so I am
replying to my own sent mail.

>    v1 = sin (0.37 * pi * n) ./ (pi * n);
>    v2 = sin (0.37 * pi * [n]) ./ (pi * [n]);
>
> I think what is hurting your precision here is the lazy evaluation of
> ranges in Octave.

This is the first time I see this, so I did some search. Correct me if
I am wrong, but, apparently, it's about the way Octave stores the
range. If it's without square braces, it uses less storage, which
probably means it uses floats instead of doubles, which might explain
the precision (in fact, I am pleasantly surprised by the accuracy,
now), but that still leaves the question of why sin(-x) != -sin(x),
though I guess it's the floating point representation of the numbers,
as they go up, or down, different fractions (just guessing), but
then...

> FWIW the h I obtain above is not equal to the h I get from the attached
> fir function. The error is worst at index 18, the same as the error you
> are seeing.

Could it be that different CPU architectures calculate differently the
floats? I can't imagine why you would get different results. Did you
use fir(33,[0 0.37])? Still, the eror is, again, at index 18, so maybe
it really is about the floating point representation.

At any rate, as soon as I made my range as n=[...], everything went
spot on. Thank you.

> I ran across a fairly recent (2008) paper on some of the "dark arts" of
> finding good approximations to mathematical functions.  It does in part
> look at sin(x).

Well, for my part, I'm fine with the way Octave does it, thank you for
the suggestion.

Vlad


Reply | Threaded
Open this post in threaded view
|

Re: Numerical issues with sin()/cos()

Mike Miller-4
On Wed, Mar 28, 2018 at 06:46:47 +0000, vlad ionescu wrote:
> > I think what is hurting your precision here is the lazy evaluation of
> > ranges in Octave.
>
> This is the first time I see this, so I did some search. Correct me if
> I am wrong, but, apparently, it's about the way Octave stores the
> range. If it's without square braces, it uses less storage, which
> probably means it uses floats instead of doubles

Actually it's that Octave stores ranges as a (start, step, stop)
triplet. Let's say the value n is the range -16.5:1:16.5 as in your
example. Octave stores that as (start=-16.5, step=1, stop=16.5).

The expression 0.37*pi*n evaluates to another range, with the original
triplet values each multiplied by 0.37*pi. Each of the 34 values of the
range has not been evaluated yet.

It's only when the range is passed in to the sin function that it is
evaluated into an array of discrete doubles. So I think you can see how
this lazy evaluation can lead to differences when adding subsequent
increments of 0.37*pi to the start value of -0.37*pi*16.5, compared to
what you actually wanted.

> > FWIW the h I obtain above is not equal to the h I get from the attached
> > fir function. The error is worst at index 18, the same as the error you
> > are seeing.
>
> Could it be that different CPU architectures calculate differently the
> floats? I can't imagine why you would get different results. Did you
> use fir(33,[0 0.37])? Still, the eror is, again, at index 18, so maybe
> it really is about the floating point representation.

Turns out this was only because I had tested with n=(-16.5:16.5)', where
the transpose operator also expands the range into a column vector.

--
mike



signature.asc (849 bytes) Download Attachment