/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  COMPUTE NUMBER OF FRACTIONS						      C
C  12/19/14 (DKC)							      C
C                                                                             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
#include <math.h>
#include <stdlib.h>
void add64(unsigned int *A, unsigned int *B);
void sub64(unsigned int *A, unsigned int *B);
void haros9(unsigned int *N, unsigned int *H, unsigned int *K,
		    unsigned int *HP, unsigned int *KP, unsigned int *D,
		    unsigned int *ND);
void mertenau(unsigned int *N, unsigned int *count, unsigned int *I) {
unsigned int J[2],ND[6],H[2],K[2],HP[2],KP[2],T[2],U[2];
J[1]=(I[1]>>1)|(I[0]<<31);
J[0]=I[0]>>1;
//J=I/2;
T[0]=J[0];
T[1]=J[1];
T[1]=T[1]+1;
if (T[1]==0)
   T[0]=T[0]+1;
if ((T[0]&0x80000000)!=0) {	 // N<=J
   count[0]=0;
   count[1]=0;
   count[2]=0;
   count[3]=0;
   return;
   }
U[0]=I[0];
U[1]=I[1];
sub64(N,U);
if (((T[0]&0x80000000)==0)&&((U[0]&0x80000000)!=0)) {  // N>=J+1 and N<I
   count[0]=0;
   count[1]=0;
   count[2]=J[0];
   count[3]=J[1];
   sub64(N,&count[2]);	     // N-J
   return;
   }
if ((U[0]==0)&&(U[1]==0)) {    // N=I
   count[0]=0;
   count[1]=0;
   if ((I[1]&1)==0) {
      count[2]=0;
      count[3]=1;
      sub64(J,&count[2]);     // J-1
      }
   else {
      count[2]=J[0];
      count[3]=J[1];
      }
   return;
   }
H[0]=0;
H[1]=0;
K[0]=0;
K[1]=1;
HP[0]=0;
HP[1]=1;
KP[0]=N[0];
KP[1]=N[1];
haros9(N,H,K,HP,KP,I,ND);
count[0]=ND[0];
count[1]=ND[1];
H[0]=0;
H[1]=1;
K[0]=I[0];
K[1]=I[1];
HP[0]=ND[2];
HP[1]=ND[3];
KP[0]=ND[4];
KP[1]=ND[5];
if ((I[1]&1)==0)
   haros9(N,H,K,HP,KP,J,ND);
else
   haros9(N,H,K,HP,KP,I,ND);
count[2]=ND[0];
count[3]=ND[1];
return;
}