#include <string.h>
#include <stdio.h>
#include <errno.h>
#include <stdlib.h>
#include <math.h>

#include "nmc.h"

double nmc_mean (double *l, unsigned int len)
{
    /* given an array of doubles (l) of length len, cacluate their mean */

    double avg = 0;
    unsigned int i=len;
    if (len == 1) return (double) *l;
    while (i--) avg += *(l++);
    return avg / len;
}

double nmc_mean_c (unsigned char *l, unsigned int len)
{
    double avg = 0;
    unsigned int i = len;
    if (len == 1) return (double) *l;
    while (i--) avg += (double) *(l++);
    return avg / len;
}

double nmc (unsigned int n, unsigned char *s1, unsigned int s1_len, unsigned char *s2, unsigned int s2_len)
{
    unsigned int s1_chunk, s1_chunk_last;
    unsigned int s2_chunk, s2_chunk_last;
    unsigned int i;
    double *diffs, mean;

    if ( n > s1_len || n > s2_len )
    {
        errno = EINVAL;
        return -1.0;
    }

    if ( 0 == n )
		n = (s1_len < s2_len) ? s1_len : s2_len;

    /* allocate memory for the difference array. We need n+1, n pieces
     * plus one final slot for the difference in length */
    diffs = calloc(n+1, sizeof(double));

    /* breaking s1 and s2 apart into n-pieces gives chunk-sized pieces */
    s1_chunk = s1_len / n;
    s2_chunk = s2_len / n;

    /* to account for s1 and s2 possibly not being evenly divisible by n */
    s1_chunk_last = ( s1_len % n ) ? s1_len % n : s1_chunk;
    s2_chunk_last = ( s2_len % n ) ? s2_len % n : s2_chunk;

    for (i = 0; i < n-1; i++)
    {
        diffs[i] =
			fabs (
				nmc_mean_c ( s1 + i*s1_chunk, s1_chunk )
				-
				nmc_mean_c ( s2 + i*s2_chunk, s2_chunk )
			);
        diffs[i] = diffs[i] < 128 ? diffs[i] : 256 - diffs[i];
        diffs[i] *= diffs[i];
    }

    /* the final chunk */
    diffs[n-1] =
		fabs (
			nmc_mean_c ( s1 + (n-1)*s1_chunk, s1_chunk_last )
			-
			nmc_mean_c ( s2 + (n-1)*s2_chunk, s2_chunk_last )
		);
    diffs[n-1] = diffs[n-1] < 128 ? diffs[n-1] : 256 - diffs[n-1];
    diffs[n-1] *= diffs[n-1];

    /* final diff: size difference */
    diffs[n] = fabs ( (double) s1_len - s2_len );
	diffs[n] *= diffs[n];

    mean = sqrt(nmc_mean(diffs,n+1));

    free(diffs);

    return mean;
}