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 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 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 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
