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
Description: Text Data