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
40 double
my_mean1(const double data[],size_t stride,const size_t size)41 my_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
55 double
my_mean2(const double data[],size_t stride,const size_t size)56 my_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
main(void)68 int 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