// Native methods implementation of
// /home/rahul/NetBeansProjects/MotifKmer/src/motifkmer/FilterSequence.java

#include "FilterSequence.h"
#include<jni.h>

int gcd(a,b)
        int a,b;
{
        int c;

        if (b<0) b= -b;
        if (b>a) { c=a; a=b; b=c; }
        for (;b;b=c) { c=a%b; a=b; }
        return a;
}





jdoubleArray JNICALL Java_motifkmer_FilterSequence_CalculateVar
  (JNIEnv * env, jobject object, jint low, jint high, jdoubleArray r) 
{
    jdoubleArray ret;
        ret =(*env)->NewDoubleArray(env,2);
        jdouble *pr = (*env)->GetDoubleArrayElements(env, r, 0);
        jdouble  K;                     /* Pointer to parameter K             */
        jdouble  H;                     /* Pointer to parameter H             */
        jint     i,j,range,lo,hi,first,last;
        jdouble  y,slope,delta,sum,Sum,beta,temp,temp2;
        jdouble  *P,*ptrP,*ptr1,*ptr2,exp();
        jchar    *calloc();
        
        
         jdouble p[high-low+1];
         for (i=range=high-low;i> -low && !pr[i];--i);
         sum=i=0;
         for (;i<=range;i++)
         {
             sum+=pr[i];
         }
           for (Sum=low,i=0;i<=range;++i) 
           {
               p[i]=pr[i]/sum;
               Sum+=i*p[i];
           } 
          jdouble lambda = -log(p[range])/high;
             do {
        beta=exp(lambda);
        temp=exp(lambda*(low-1));
        y=slope=0;
        ptr1=p;
        for (i=low;i<=high;++i) {
            y+=temp2= *ptr1++ * (temp*=beta);
            slope+= temp2*i;
            }
        lambda-=delta=(y-1)/slope;
      } while (delta / lambda > 1e-6);
     H= lambda*slope/log(2.0);
     Sum=lo=hi=0;
     P= (jdouble *) calloc(500*range+1,sizeof(jdouble));
     
     for (*P=sum=j=1;j<=500 && sum>0.0001;Sum+=sum/=j++) {
                first=last=range;
                for (ptrP=P+(hi+=high)-(lo+=low);ptrP>=P;*ptrP-- =sum) {
                        ptr1=ptrP-first;
                        ptr2=p+first;
                        for (sum=0,i=first;i<=last;++i) sum+= *ptr1-- * *ptr2++;
                        if (first) --first;
                        if (ptrP-P<=range) --last;
                }
                temp=exp(lambda*(lo-1));
                for (sum=0,i=lo;i;++i) sum+= *++ptrP * (temp*=beta);
                for (;i<=hi;++i) sum+= *++ptrP;
        }
        free(P);
        for (i=low;!p[i-low];++i);
        for (j= -i;i<high && j>1;) if (p[++i-low]) j=gcd(j,i);
        K = (j*exp(-2*Sum))/(slope*(1.0-exp(- lambda*j)));
       
        jdouble arr[]={lambda,K};
        (*env)->SetDoubleArrayRegion(env,ret,0,2,arr) ;
        
        return (ret);
}
