VEX/test/rounderr.c

   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