tab-co.c ver 0.7

#include <stdio.h>
#include <stdlib.h>
typedef unsigned long long U64;
typedef unsigned short gtype;

typedef struct filter_data
{
    unsigned int val;
    U64 amin;
    U64 bmin;
    gtype gmin;
} filter_data;

void init_filter_data(int bits, filter_data *f)
{
    U64 a = (1<<bits);
    U64 b = (f->val)&(a-1);
    gtype g = 0;

    while((a&1)==0)
    {
	if (b&1) //odd
	{
	    a = 3*a;
	    b = 3*b + 1;
	    g++;
	}
	a /= 2;
	b /= 2;
	g++;
	if(a< (1<<bits)/4) break;
    }

    f->amin = a;
    f->bmin = b;
    f->gmin = g;
}

// For binary search
typedef struct b_entry
{
    long long val;
    int hash;
    struct b_entry *next_gt, *next_lt;
} b_entry;

// Hash to make the binary tree evenly grow.
// Any mapping which randomize the integer order is OK
// 1-to-1 mapping is better (but not necessary)
#define BHASH(x) ((unsigned int)(x)*54321*54321*54321)

static b_entry *b_buffer;
static int n_b_buffer;

static b_entry *new_entry()
{
    return &b_buffer[n_b_buffer++];
}

static int find_b(long long b, b_entry *root)
{
    b_entry *p = root;
    b_entry **prev = NULL;
    int bhash = BHASH(b);

    //Binary search
    while(p!=NULL)
    {
	long long bcmp = p->hash - bhash;
	if(bcmp==0) bcmp = p->val - b;
        if(bcmp == 0) return 0;
	prev = (bcmp < 0)? &(p->next_lt):&(p->next_gt);
	p = *prev;
    }

    p = new_entry();

    p->val = b;
    p->hash = bhash;
    p->next_lt = p->next_gt = NULL;
    if (prev != NULL) *prev = p;
    return 1;
}

extern void init_cutoff_tab(int maxbits, int *tabsize, filter_data **co_tab)
{
    int pow3[maxbits+1];
    b_entry root[maxbits+1];
    int i;
    #define TABSIZE (1<<maxbits)
    unsigned char *tmp_buf = (unsigned char *)malloc(TABSIZE);
    int size = 0;
    int nskip = 0;
    b_buffer = (b_entry *)malloc(sizeof(b_entry)*TABSIZE/2);
    n_b_buffer = 0;

    printf("# Cutoff Init: Temporally Memory Usage: %f Mb\n",
	    (sizeof(b_entry)*TABSIZE/2 + TABSIZE)*1e-6);

    for (i=0;i<=maxbits;i++)
    {
	pow3[i] = (i==0)? 1:3*pow3[i-1];
	root[i].val = -1;
	root[i].next_gt = root[i].next_lt = NULL;
    }

    for(i=0;i<TABSIZE;i++)
    {
	int j;
	int a3 = 0;
	long long b = 0;
	int i0 = i;
        b_entry *p;

	// Transform 3^a3 * (2m+ 0 or 1) +b
	for(j=0;j<maxbits;j++)
	{
	    if (i0&1)
	        b += pow3[a3];
            if (b&1)
	    {
	        b = b*3+1;
	        a3++;
	    }
	    b >>=1;
	    i0>>=1;
	}
        p = &root[a3];

	if (p->val == -1)
	{
	    p->val = b;
	    tmp_buf[i]=1;
	}
	else
	{
	    tmp_buf[i] = find_b(b, p);
	}

	if (tmp_buf[i]==1)
	{
	    for(j=0;j<3;j++)
	    {
	        int m = (i+j*TABSIZE)%6;
#ifdef CHECK_EVEN
	        if (m==1 || m==3 || m==0) size ++;
#else
	        if (m==1 || m==3) size ++;
#endif
	    }
	}
    }
    printf("# %d Bit: Cutoff Ratio %6.3f (%d/%d)\n",
	    maxbits,size*1.0/(TABSIZE*3),size,TABSIZE*3);

    free(b_buffer);

    (*co_tab) = (filter_data *)malloc(size*sizeof(filter_data));

    printf("# Cutoff Table size  %6.3f Mb\n", size*sizeof(filter_data)*1e-6);

    *tabsize = size;
    size=0;

    for(i=0;i<TABSIZE*3;i++)
    {
        int m = i%6;
#ifdef CHECK_EVEN
	if (tmp_buf[i%TABSIZE]==1 && (m==1 || m==3 || m==0))
#else
	if (tmp_buf[i%TABSIZE]==1 && (m==1 || m==3))
#endif
	{
	    filter_data *f = &((*co_tab)[size]);
	    f->val = i;
	    init_filter_data(maxbits, f);
	    size++;
	    if (f->amin < TABSIZE/2) nskip++;
	}
    }
    free(tmp_buf);
    printf("# MAX a-posteriori skip ratio %f\n", nskip*1.0/size);
}