/*
  Native calls to do dot product and sliding dot product quickly

  @author Scott Fernandez
  @version $Id: dotprod.c,v 1.1 2011/10/14 18:52:41 ntn Exp $
 */

float inline dot_product (float* in1, float* in2, int toCalculate) {
  int unroll = 4;
  int left = toCalculate % unroll;
  float* end = &in1[toCalculate];
  float* unrolledEnd = end - left;
  float sum, sum1, sum2, sum3, sum4;
  sum1 = sum2 = sum3 = sum4 = 0;
  while(in1 != end) {
    sum1 += in1[0]*in2[0];
    sum2 += in1[1]*in2[1];
    sum3 += in1[2]*in2[2];
    sum4 += in1[3]*in2[3];
    in1 += 4;
    in2 += 4;
  }
  return sum1 + sum2 + sum3 + sum4;
}

/* Ugly but fast */
int inline sliding_dot_product (float* in, float* imp, float* out, int toCalculate, int impSize)
{
  /* Setup */
  int unrollOuter = 4;
  int unrollInner = 2;
  int outerLeft = toCalculate % unrollOuter;
  int innerLeft = impSize % unrollInner;
  int impMov = innerLeft - impSize;
  int inMov = impMov + unrollOuter;
  double sum, sum1, sum2, sum3, sum4;
  float impVal;
  float* outEnd = &out[toCalculate];
  float* outerUnrolledEnd = outEnd - outerLeft;
  float* impEnd = &imp[impSize];
  float* innerUnrolledEnd = impEnd - innerLeft;
  /* Unrolled outer looping--calculate 4 values at once */
  if (innerLeft) {
    while (out != outerUnrolledEnd) {
      sum1 = sum2 = sum3 = sum4 = 0.0;
      /* Unrolled inner looping--calculate 2 products before checking condition */
      while (imp != innerUnrolledEnd) {
        impVal = (imp++)[0];
        sum1 += (double)((in++)[0] * impVal);
        sum2 += (double)(in[0]     * impVal);
        sum3 += (double)(in[1]     * impVal);
        sum4 += (double)(in[2]     * impVal);

        impVal = (imp++)[0];
        sum1 += (double)((in++)[0] * impVal);
        sum2 += (double)(in[0]     * impVal);
        sum3 += (double)(in[1]     * impVal);
        sum4 += (double)(in[2]     * impVal);
      }
      impVal = imp[0];
      sum1 += (double)(in[0] * impVal);
      sum2 += (double)(in[1] * impVal);
      sum3 += (double)(in[2] * impVal);
      sum4 += (double)(in[3] * impVal);

      (out++)[0] = (float)sum1;
      (out++)[0] = (float)sum2;
      (out++)[0] = (float)sum3;
      (out++)[0] = (float)sum4;
      imp += impMov;
      in += inMov;
    }
  }
  else {
    while (out != outerUnrolledEnd) {
      sum1 = sum2 = sum3 = sum4 = 0.0;
      /* Unrolled inner looping--calculate 2 products before checking condition */
      while (imp != innerUnrolledEnd) {
        impVal = (imp++)[0];
        sum1 += (double)((in++)[0] * impVal);
        sum2 += (double)(in[0]     * impVal);
        sum3 += (double)(in[1]     * impVal);
        sum4 += (double)(in[2]     * impVal);

        impVal = (imp++)[0];
        sum1 += (double)((in++)[0] * impVal);
        sum2 += (double)(in[0]     * impVal);
        sum3 += (double)(in[1]     * impVal);
        sum4 += (double)(in[2]     * impVal);
      }
      (out++)[0] = (float)sum1;
      (out++)[0] = (float)sum2;
      (out++)[0] = (float)sum3;
      (out++)[0] = (float)sum4;
      imp += impMov;
      in += inMov;
    }
  }
  /* Finish up one at a time when close to the end */
  while (out != outEnd) {
    sum = 0.0;
    while (imp != impEnd) {
      impVal = (imp++)[0];
      sum += (double) ((in++)[0] * impVal);
    }
    imp -= impSize;
    in -= impSize - 1;

    (out++)[0] = (float)sum;
  }
  return toCalculate;
}

/* WARNING: imp is assumed to be semi-normalized--see correlate.java */
int inline normalized_sliding_dot_product (float* in, float* imp, float* out,
    int toCalculate, int impSize, float* in_norm_factors)
{
  /* Setup */
  int unrollOuter = 4;
  int unrollInner = 2;
  int outerLeft = toCalculate % unrollOuter;
  int innerLeft = impSize % unrollInner;
  int impMov = innerLeft - impSize;
  int inMov = impMov + unrollOuter;
  double sum, sum1, sum2, sum3, sum4;
  float impVal;
  float* outEnd = &out[toCalculate];
  float* outerUnrolledEnd = outEnd - outerLeft;
  float* impEnd = &imp[impSize];
  float* innerUnrolledEnd = impEnd - innerLeft;

  /* Unrolled outer looping--calculate 4 values at once */
  if (innerLeft) {
    while (out != outerUnrolledEnd) {
      sum1 = sum2 = sum3 = sum4 = 0.0;
      /* Unrolled inner looping--calculate 2 products before checking condition */
      while (imp != innerUnrolledEnd) {
        impVal = (imp++)[0];
        sum1 += (double)((in++)[0] * impVal);
        sum2 += (double)(in[0]     * impVal);
        sum3 += (double)(in[1]     * impVal);
        sum4 += (double)(in[2]     * impVal);

        impVal = (imp++)[0];
        sum1 += (double)((in++)[0] * impVal);
        sum2 += (double)(in[0]     * impVal);
        sum3 += (double)(in[1]     * impVal);
        sum4 += (double)(in[2]     * impVal);
      }
      impVal = imp[0];
      sum1 += (double)(in[0]     * impVal);
      sum2 += (double)(in[1]     * impVal);
      sum3 += (double)(in[2]     * impVal);
      sum4 += (double)(in[3]     * impVal);

      (out++)[0] = (float) (sum1/(in_norm_factors++)[0]);
      (out++)[0] = (float) (sum2/(in_norm_factors++)[0]);
      (out++)[0] = (float) (sum3/(in_norm_factors++)[0]);
      (out++)[0] = (float) (sum4/(in_norm_factors++)[0]);
      imp += impMov;
      in += inMov;
    }
  }
  else {
    while (out != outerUnrolledEnd) {
      sum1 = sum2 = sum3 = sum4 = 0.0;
      /* Unrolled inner looping--calculate 2 products before checking condition */
      while (imp != innerUnrolledEnd) {
        impVal = (imp++)[0];
        sum1 += (double)((in++)[0] * impVal);
        sum2 += (double)(in[0]     * impVal);
        sum3 += (double)(in[1]     * impVal);
        sum4 += (double)(in[2]     * impVal);

        impVal = (imp++)[0];
        sum1 += (double)((in++)[0] * impVal);
        sum2 += (double)(in[0]     * impVal);
        sum3 += (double)(in[1]     * impVal);
        sum4 += (double)(in[2]     * impVal);
      }
      (out++)[0] = (float) (sum1/(in_norm_factors++)[0]);
      (out++)[0] = (float) (sum2/(in_norm_factors++)[0]);
      (out++)[0] = (float) (sum3/(in_norm_factors++)[0]);
      (out++)[0] = (float) (sum4/(in_norm_factors++)[0]);
      imp += impMov;
      in += inMov;
    }
  }
  /* Finish up one at a time when close to the end */
  while (out != outEnd) {
    sum = 0.0;
    while (imp != impEnd) {
      impVal = (imp++)[0];
      sum += (double)((in++)[0] * impVal);
    }
    imp -= impSize;
    in -= impSize - 1;

    (out++)[0] = (float) (sum/(in_norm_factors++)[0]);
  }

  return toCalculate;
}
