#include <assert.h>
#include <math.h>


long factorial(int n)
{
	int res=1 ;
	for (int i=2; i<=n; i++)
		res*=i ;
	return res ;
}

long nchoosek(int n, int k)
{
	long res=1 ;
	
	for (int i=n-k+1; i<=n; i++)
		res*=i ;
	
	return res/factorial(k) ;
}

double compute(char *avec, int alen, char* bvec, int blen, int degree, int max_mismatch, double* weights)
{
  double sum=0;
  assert(alen==blen) ;
  int *match_vector=new int[alen] ;

  for (int i=0; i<alen; i++)
	  match_vector[i]=(avec[i]!=bvec[i]) ;
  
  for (int i=0; i<alen-degree; i++)
  {
	  int mismatches=0;
	  
	  for (int j=0; j<degree; j++)
	  {
		  if (match_vector[i+j])
		  {
			  mismatches++ ;
			  if (mismatches>max_mismatch)
				  break ;
		  } ;
		  sum += weights[j+degree*mismatches];
	  }
  }
  delete[] match_vector ;
  
  return (double) sum ;
}

double compute_normalized(char *avec, int alen, char* bvec, int blen, 
						  int degree, int max_mismatch)
{
	/* you may want to precompute this ...*/
	double* weights=new double[degree*(1+max_mismatch)];
	double sum=0;

	for (int i=0; i<degree; i++)
	{
		weights[i]=degree-i;
		sum+=weights[i];
	}
	for (int i=0; i<degree; i++)
		weights[i]/=sum;
	
	for (int i=0; i<degree; i++)
	{
		for (int j=1; j<=max_mismatch; j++)
		{
			if (j<i+1)
			{
				int nk=nchoosek(i+1, j) ;
				weights[i+j*degree]=weights[i]/(nk*pow(3,j)) ;
			}
			else
				weights[i+j*degree]= 0;
		} ;
	} ;
	
	/* the diagonal may also be precomputed */
	double x2=compute(avec, alen, avec, alen, degree, max_mismatch, weights) ;
	double y2=compute(bvec, blen, bvec, blen, degree, max_mismatch, weights) ;
	double xy=compute(avec, alen, bvec, blen, degree, max_mismatch, weights) ;

	return xy/sqrt(x2*y2) ;
}