1
2/* peach (7400, altivec supported, 450MHz, gcc -O)
3   m1 = 1.20000000000000018,  exp = 1.19999999999999996
4   m2 = 1.19999999999998885,  exp = 1.19999999999999996
5*/
6
7/* peach (7400, altivec supported, 450MHz, gcc)
8   m1 = 1.20000000000000018,  exp = 1.19999999999999996
9   m2 = 1.19999999999998885,  exp = 1.19999999999999996
10*/
11
12/* phoenix, gcc -O
13   m1 = 1.19999999999999996,  exp = 1.19999999999999996
14   m2 = 1.19999999999999996,  exp = 1.19999999999999996
15*/
16
17/* phoenix, icc -O
18   m1 = 1.19999999999999996,  exp = 1.19999999999999996
19   m2 = 1.19999999999999996,  exp = 1.19999999999999996
20*/
21
22/* phoenix, gcc -O, iropt-level=2
23   m1 = 1.20000000000000040,  exp = 1.19999999999999996
24   m2 = 1.19999999999999440,  exp = 1.19999999999999996
25*/
26
27/* phoenix, gcc -O, iropt-level=1/0
28   m1 = 1.20000000000000018,  exp = 1.19999999999999996
29   m2 = 1.19999999999998885,  exp = 1.19999999999999996
30*/
31
32
33#include <stdlib.h>
34#include <stdio.h>
35#include <math.h>
36
37#define NNN 1000
38
39
40double
41my_mean1 (const double data[], size_t stride, const size_t size)
42{
43  /* Compute the arithmetic mean of a dataset using the recurrence relation
44     mean_(n) = mean(n-1) + (data[n] - mean(n-1))/(n+1)   */
45  long double mean = 0;
46  size_t i;
47
48  for (i = 0; i < size; i++)
49    {
50      mean += (data[i * stride] - mean) / (i + 1);
51    }
52  return mean;
53}
54
55double
56my_mean2 (const double data[], size_t stride, const size_t size)
57{
58  /* Compute the arithmetic mean of a dataset using the
59     obvious scheme. */
60  int i;
61  long double sum = 0;
62  for (i = 0; i < size; i++)
63    sum += data[i * stride];
64  return sum / (double)size;
65}
66
67
68int main (void)
69{
70  int i;
71  const size_t nacc2 = NNN+1;
72  double numacc2[NNN+1] ;
73
74  numacc2[0] = 1.2 ;
75
76  for (i = 1 ; i < NNN; i += 2)
77     numacc2[i] = 1.1 ;
78
79  for (i = 1 ; i < NNN; i += 2)
80      numacc2[i+1] = 1.3 ;
81
82#if 1
83  asm __volatile__("fninit");
84#endif
85
86  {
87    double m1 = my_mean1 (numacc2, 1, nacc2);
88    double m2 = my_mean2 (numacc2, 1, nacc2);
89    double expected_mean = 1.2;
90    printf("m1 = %19.17f,  exp = %19.17f\n", m1, expected_mean);
91    printf("m2 = %19.17f,  exp = %19.17f\n", m2, expected_mean);
92  }
93
94  return 0;
95}
96
97
98