gistfile1.c

#include <stdio.h>
#include <gmp.h>
#include <math.h>
#include <assert.h>
#include <float.h>
#include "mpunrank.h"

/*
int same(double a, double b)
{
	return (a+0.01 > b) && (b + 0.01 > a);
}
*/
	const int xsiz = 537;
	const int ysiz = 40;

	const double lc_upps = +1.000000000000001;
	const double lc_down = -0.0000000000000000557;
	const double lc_safe = +1.0000000000000010557;

double lc(unsigned long int xx, unsigned long int yy)
{
	if (xx == 0 ||  yy  == 0)
		return 1.;

	double x = (xx);
	double y = (yy);

	double w = (lgammal(x+y-1.)-lgammal(x)-lgammal(y))/log(2.0);

	return w;
}

int gmpccmp(mpz_t dst, mpz_t bi[xsiz][ysiz], unsigned int x, unsigned int y)
{
	if (x == 0 ||  y == 0)
		return mpz_cmp_ui(dst, 0);
	if (x == 1 ||  y == 1)
		return mpz_cmp_ui(dst, 1);
	if (x == 2)
		return mpz_cmp_ui(dst, y);
	if (y == 2)
		return mpz_cmp_ui(dst, x);

	if (x < y)
		return mpz_cmp(dst, bi[y-3][x-3]);
	else
		return mpz_cmp(dst, bi[x-3][y-3]);
}

void gmpcsub(mpz_t dst, mpz_t bi[xsiz][ysiz], unsigned int x, unsigned int y)
{
	if (x == 0 ||  y  == 0)
		return;
	if (x == 1 ||  y  == 1)
		return mpz_sub_ui(dst, dst, 1);
	if (x == 2)
		return mpz_sub_ui(dst, dst, y);
	if (y == 2)
		return mpz_sub_ui(dst, dst, x);

	if (x < y)
		return mpz_sub(dst, dst, bi[y-3][x-3]);
	else
		return mpz_sub(dst, dst, bi[x-3][y-3]);
}


void gmpcadd(mpz_t dst, mpz_t bi[xsiz][ysiz], unsigned int x, unsigned int y)
{
	if (x == 0 ||  y  == 0)
		return;
	if (x == 1 ||  y  == 1)
		return mpz_add_ui(dst, dst, 1);
	if (x == 2)
		return mpz_add_ui(dst, dst, y);
	if (y == 2)
		return mpz_add_ui(dst, dst, x);

	if (x < y)
		return mpz_add(dst, dst, bi[y-3][x-3]);
	else
		return mpz_add(dst, dst, bi[x-3][y-3]);
}

void gmpgenc(mpz_t bi[xsiz][ysiz])
{
	unsigned int x, y;

	if (0 < xsiz && 0 < ysiz)
		mpz_init_set_ui(bi[0][0], 6);

	for (x = 1; x < xsiz; x++) {
		mpz_init(bi[x][0]);
		mpz_add_ui(bi[x][0], bi[x-1][0], x+3);
	}
	for (y = 1; y < ysiz; y++) {
		mpz_init(bi[0][y]);
		mpz_add_ui(bi[0][y], bi[0][y-1], y+3);
	}

	for (y = 1; y < ysiz; y++) {
		for (x = 1; x < xsiz; x++) {
			mpz_init(bi[x][y]);
			mpz_add(bi[x][y], bi[x][y-1], bi[x-1][y]);
		}
	}
}

void gmpfrec(mpz_t bi[xsiz][ysiz])
{
	unsigned int x, y;

	for (y = 0; y < ysiz; y++)
		for (x = 0; x < xsiz; x++)
			mpz_clear(bi[x][y]);
}

void unrank(mpz_t bi[xsiz][ysiz], ucallback_f f, void *p, mpz_t tr, int z1, int z0)
{
	unsigned long int n = 0;

	while (z0 + z1) {

		/* test 0 - logarithmic test */
		double tl = lc(z0, z1 + 1);
		double utl = tl + lc_safe;
		unsigned long l = mpz_sizeinbase(tr, 2);
		double ltl = tl - lc_safe;
#if 0
		int reality;

			/* test 1 - treshold comparation test*/
			if (gmpccmp(tr, bi, z0, z1 + 1) < 0) {
				reality = 0;
			} else {
				reality = 1;
			}
#endif
		if (utl < l) {
#if 0
			assert(reality);
#endif
			f(p, n);
			gmpcsub(tr, bi, z0, z1 + 1);
			assert(z1);
			z1--;
		} else if (ltl > l) {
#if 0
			if (reality) {
				printf("%u %u\n", z0, z1 + 1);
				printf("%lf > %lu\n", ltl, l);
			}
			assert(!reality);
#endif
			assert(z0);
			z0--;
		} else {
			/* test 1 - treshold comparation test*/
			if (gmpccmp(tr, bi, z0, z1 + 1) < 0) {
				assert(z0);
				z0--;
			} else {
				f(p, n);
				gmpcsub(tr, bi, z0, z1 + 1);
				assert(z1);
				z1--;
			}
		}
		/**/

		n++;
	}

}

void ftunrankz(mpz_t bi[xsiz][ysiz], ucallback_f f, unsigned char *ra, unsigned char *rb, mpz_t rc, int z1, int z0)
{
	struct nested n;
	n.bi = n.bj = 0;
	n.p_ra = ra;
	n.p_rb[0] = &rb[0];
	n.p_rb[1] = &rb[(z1+7)/8];
	n.p_rb[2] = &rb[(z1+3)/4];
	n.bin = &bi[0][0];

	unrank(bi, f, &n, rc, z1, z0);
}

void funrankz(ucallback_f f, unsigned char *ra, unsigned char *rb, mpz_t rc, int z1, int z0)
{
	mpz_t bi[xsiz][ysiz];
	gmpgenc(bi);

	ftunrankz(bi, f, ra, rb, rc, z1, z0);
	gmpfrec(bi);
}

void funrank(ucallback_f f, unsigned char *ra, unsigned char *rb, mpz_t rc)
{
	int z0 = 536;	/*560 = 4 14 or 5 15 */
	int z1 = 24;	/* 24 clausules */

	funrankz(f, ra, rb, rc, z1, z0);
}

int cgbit(unsigned char *c, unsigned long int n)
{
	return c[n / 8] & (1<<(n & 7));
}

void csbit(unsigned char *c, unsigned long int n)
{
	c[n / 8] |= (1<<(n & 7));
}