double        dotprod1(double a[], double b[], int n)
{
  double sum;
  int    i;

  sum = 0.0;
  for ( i = 0; i < n; i++ )
    sum = sum + a[i]*b[i];
  

  return sum;
}

double        dotprod4(double a[], double b[], int n)
{
  double sum0, sum1, sum2, sum3, sum;
  int    i, i4, n4;

  n4 = n / 4;
  /* do inner product in blocks of length 4 */
  sum0 = sum1 = sum2 = sum3 = 0.0;
  for ( i4 = 0; i4 < n4; i4++ )
  {
    sum0 = sum0 + a[4*i4  ]*b[4*i4  ];
    sum1 = sum1 + a[4*i4+1]*b[4*i4+1];
    sum2 = sum2 + a[4*i4+2]*b[4*i4+2];
    sum3 = sum3 + a[4*i4+3]*b[4*i4+3];
  }
  sum = sum0 + sum1 + sum2 + sum3;
  /* take care of fractional block */
  for ( i = 4*n4; i < n; i++ )
    sum = sum + a[i]*b[i];
  

  return sum;
}

int main(int argc, char *argv[])
{
#define        N 103
  double a[N], b[N], sum1, sum4;
  int    i;

  /* here is some test data */
  for ( i = 0; i < N; i++ )
    {
      a[i] = (double)i;
      b[i] = ( i & 1 ) ? -1.0 : 1.0;
    }

  sum1 = dotprod1(a,b,N);
  sum4 = dotprod4(a,b,N);
  printf("dotprod1 = %g, dotprod2 = %g, difference = %g\n",
         sum1, sum4, sum1-sum4);
}