double precision and round-off errors

[Micko]

Hi, I'm writing program that calculates coefficients of a FIR filter. Anyway I have a problem regarding precision of numbers.
I need to use so called Blackman window:

w[i]=0.42-0.5*cos(2*pi*i/(Len-1))+0.08*cos(4*pi*i/(Len-1));

;

As you can notice coefficients are symetrical and for Len=11;
there is symetry around element w[5] which is equal 1.00

I wrote function that performs calculations:

#define pi 3.14159265
...
double* blackman(int Len)
{
	double* array=malloc(sizeof(double)*Len);
	int i;
	for(i=0;i<Len;i++)
		array[i]=0.42-0.5*cos(2*pi*i/(Len-1))+0.08*cos(4*pi*i/(Len-1));
	return array;
}

When I execute this code I get following results:
-1.387779e-017
4.021286e-002
2.007701e-001
5.097871e-001
8.492299e-001
1.000000e+000
8.492299e-001
5.097871e-001
2.007701e-001
4.021286e-002
-2.775558e-017

As you can see the first and the second element should be 0.00

I have idea to use if(){} to check if argumet to cos function is 2*pi*k (k=0,1,2,3...) and to set appropriate element to 0.

Is there any other way to increase precision? I assume problem is coming from cos function because cosine function is usually calculated through Taylor's series.
All your ideas are welcome.
Thanks

[Salem]

http://docs.sun.com/source/806-3568/ncg_goldberg.html

[Dave Evans]

Before you ask, "How can I increase the precision of the results?", I think you should as, "What is the best precision that I can hope for?"

Standard C compilers use "double" for all floating point calculations, and IEEE standard double precision is about 16 or 17 significant decimal digits. If your compiler has "long double" it may add about three more significant decimal digits.

Note that your first term is off about 1 unit it the 17th significant digit, and the last term has a little more.

Your approximation for pi is only nine significant digits; here's a way to get a better one without having to memorize (or store somwhere) more digits:

#include <math.h>
  .
  .
  .
  double pi;
  .
  .
  .
  pi = 4.0 * atan(1.0);
  .
  .
  .

If your compiler has "long double", try the same thing with the long atan function:

  long double pi_long;

  pi_long = 4.0L * atanl(1.0L);

Here's a program that I compiled and ran with Borland bcc32 version 5.1.1.

#include <math.h>
#include <stdio.h>

int main()
{
  int i;
  int Len;

  double pi1;
  double pi2;
  long double pi3;
  long double value3;

  double value1;
  double value2;

  pi1 = 3.14159265; /* only nine significant digits! */
  pi2 = 4.0 * atan(1.0);
  pi3 = 4.0L * atanl(1.0L);

  printf("\n\nCheck cos() with zero:\n");
  printf("  cos(0) = %21.18le\n\n", cos(0));

  printf("Look at approximations to pi:\n");
  printf("  pi1  = %21.18lf --- defined with 9 significant digits\n", pi1);
  printf("  pi2  = %21.18lf --- from 4.0 * atan(1.0)\n", pi2);
  printf("  pi3  = %21.18Lf --- long double from 4.0 * atanl(1.0)\n\n",pi3);

  printf("Look at roundoff with cos(4*pi):\n");
  printf("  cos(4.0 * pi1)  = %21.18lf\n", cos(4.0 * pi1));
  printf("  cos(4.0 * pi2)  = %21.18lf\n", cos(4.0 * pi2));
  printf("  cos(4.0 * pi3)  = %21.18Lf\n\n", cosl(4.0 * pi3));

  printf("Now, just a simple addition of constants:\n");
  printf("  with double:      0.42  - 0.5  + 0.08  = %21.18le\n", 
         0.42 - 0.5 + 0.08);
  printf("  with long double: 0.42L - 0.5L + 0.08L = %21.18Le\n\n", 
         0.42L - 0.5L + 0.08L);

  Len = 11;
  for (i = 0; i < Len; i++) {
    printf("Calculate term %d of transform:\n", i);

    value1 = 0.42 - 0.5 * cos(2.0 * pi1 * (double)i/(Len-1.0)) + 
             0.08 * cos(4.0 * pi1 * (double)i/(Len-1));

    printf("  value with pi1:  %21.18le\n", value1);

    value2 = 0.42 - 0.5 * cos(2.0 * pi2 * (double)i/(Len-1.0)) + 
             0.08 * cos(4.0 * pi2 * (double)i/(Len-1));

    printf("  value with pi2:  %21.18le\n", value2);


    value3 = 0.42L - 0.5L * cosl(2.0L * pi3 * (long double)i/(Len-1.0L)) + 
             0.08L * cosl(4.0L * pi3 * (long double)i/(Len-1.0L));

    printf("  value with pi3:  %21.18Le\n\n", value3);
  }

  return 0;
}

Regards,

Dave