indexcalculus.c

这是一个C的源代码

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gmp.h>
#include "pbc.h"
#include "pbc_utils.h"

static int is_gen(mpz_t x, mpz_t q, darray_ptr fac, darray_ptr mul)
{
    int result = 1;
    mpz_t z;
    mpz_t q1;
    int i;
    UNUSED_VAR(mul);

    mpz_init(z);
    mpz_init(q1);

    mpz_sub_ui(q1, q, 1);
    for (i=0; i<fac->count; i++) {
	mpz_divexact(z, q1, fac->item[i]);
	mpz_powm(z, x, z, q);
	if (!mpz_cmp_ui(z, 1)) {
	    result = 0;
	    break;
	}
    }

    mpz_clear(q1);
    mpz_clear(z);
    return result;
}

static void CRT(mpz_t x, mpz_ptr *v, mpz_ptr *m, int t)
//Garner's Algorithm
//see Algorithm 14.71, Handbook of Cryptography
{
    mpz_t u;
    mpz_t C[t];
    int i, j;

    mpz_init(u);
    for (i=1; i<t; i++) {
	mpz_init(C[i]);
	mpz_set_ui(C[i], 1);
	for (j=0; j<i; j++) {
	    mpz_invert(u, m[j], m[i]);
	    mpz_mul(C[i], C[i], u);
	    mpz_mod(C[i], C[i], m[i]);
	}
    }
    mpz_set(u, v[0]);
    mpz_set(x, u);
    for (i=1; i<t; i++) {
	mpz_sub(u, v[i], x);
	mpz_mul(u, u, C[i]);
	mpz_mod(u, u, m[i]);
	for (j=0; j<i; j++) {
	    mpz_mul(u, u, m[j]);
	}
	mpz_add(x, x, u);
    }

    for (i=1; i<t; i++) mpz_clear(C[i]);
    mpz_clear(u);
}

//TODO: http://www.cecm.sfu.ca/CAG/abstracts/aaron27Jan06.pdf
//TODO: don't need to store last element of list in row[i]
//TODO: linked lists might be better than dynamic arrays (avoids memmove())
//TODO: allow holes in the table
//(if drought lasts too long)
void index_calculus_step1(mpz_t *ind, int r, mpz_t g, mpz_t q,
	darray_ptr fac, darray_ptr mul)
{
    int count = 0;
    int i, j;
    mpz_t z, z0, z1;
    int relcount;
    unsigned int *prime = pbc_malloc(sizeof(unsigned int) * r);
    int bundlecount = (r - 10 + 19) / 20;
    mpz_t *bundle = pbc_malloc(sizeof(mpz_t) * bundlecount);
    int faci;
    mpz_t k, km;

    struct cell_s {
	int ind;
	mpz_t data;
    };
    typedef struct cell_s *cell_ptr;

    cell_ptr newcell(void)
    {
	cell_ptr res;
	res = pbc_malloc(sizeof(struct cell_s));
	//res->data = pbc_malloc(sizeof(mpz_t));
	//mpz_init(res->data);
	mpz_init(res->data);
	return res;
    }
    void delcell(void *p)
    {
	cell_ptr cp = p;
	mpz_clear(cp->data);
	pbc_free(p);
    }

    cell_ptr *rel = pbc_malloc(sizeof(cell_ptr) * r);
    cell_ptr *relm = pbc_malloc(sizeof(cell_ptr) * r);
    //''matrix'' is actually a list of matrices
    //(we solve over different moduli and combine using CRT)
    //darray_t **matrix = pbc_malloc(sizeof(darray_t *) * fac->count);
    darray_t *matrix[fac->count];
    int minfound[fac->count];

    for (i=0; i<r; i++) {
	rel[i] = newcell();
	relm[i] = newcell();
    }
    for (i=0; i<fac->count; i++) {
	//similarly ''row'' refers to a list of rows
	darray_t *row = pbc_malloc(sizeof(darray_t) * r);
	for (j=0; j<r; j++) {
	    darray_init(row[j]);
	}
	matrix[i] = row;
	minfound[i] = 0;
    }

    mpz_init(k);
    mpz_init(km);
    mpz_init(z);
    mpz_init(z1);
    mpz_init(z0);

    printf("building prime table...\n");
    prime[0] = 2;
    mpz_set_ui(z, 2);
    for (i=1; i<r; i++) {
	mpz_nextprime(z, z);
	prime[i] = mpz_get_ui(z);
    }

    for (i=0; i<bundlecount; i++) {
	mpz_init(bundle[i]);
	mpz_set_ui(bundle[i], 1);
	for (j=0; j<20; j++) {
	    int jj = 10 + 20 * i + j;
	    if (jj >= r) break;
	    mpz_mul_ui(bundle[i], bundle[i], prime[jj]);
	}
	element_printf("bundle %d: %Zd\n", i, bundle[i]);
    }
    printf("searching for r-smooth numbers\n");

    mpz_set_ui(z, 1);
    mpz_init(k);
    int try = 0;
    do {
	mpz_mul(z, z, g);
	mpz_mod(z, z, q);
	mpz_add_ui(k, k, 1);

	/*
	pbc_mpz_random(k, q);
	mpz_powm(z, g, k, q);
	*/

	try++;

	mpz_set(z1, z);
	relcount = 0;
	for (i=0; i<10; i++) {
	    if (i >= r) break;
	/*
	for (i=0; i<r; i++) {
	*/
	    j = 0;
	    while (mpz_divisible_ui_p(z1, prime[i])) {
		mpz_divexact_ui(z1, z1, prime[i]);
		j++;
	    }
	    if (j) {
		rel[relcount]->ind = i;
		mpz_set_ui(rel[relcount]->data, j);
		relcount++;
		if (!mpz_cmp_ui(z1, 1)) goto found;
	    }
	}
	//continue;
	for (i=0; i<bundlecount; i++) {
	    mpz_gcd(z0, bundle[i], z1);
	    if (mpz_cmp_ui(z0, 1)) {
		int ii;
		for (ii = 0; ii < 20; ii++) {
		    int jj = 10 + i * 20 + ii;
		    if (jj >= r) break;
		    j = 0;
		    while (mpz_divisible_ui_p(z1, prime[jj])) {
			mpz_divexact_ui(z1, z1, prime[jj]);
			j++;
		    }
		    if (j) {
			rel[relcount]->ind = jj;
			mpz_set_ui(rel[relcount]->data, j);
			relcount++;
			if (!mpz_cmp_ui(z1, 1)) goto found;
		    }
		}
	    }
	}
	continue;
found:

/*
	printf("found r-smooth number after %d tries\n", try);

	gmp_printf("g^%Zd = %Zd:", k, z);
	for (i=0; i<relcount; i++) {
	    gmp_printf(" %u:%Zd", rel[i]->ind, rel[i]->data);
	}
	printf("\n");
*/
	try = 0;

	for (faci=0; faci<fac->count; faci++) {
	    darray_t *row = matrix[faci];
	    mpz_ptr order = fac->item[faci];
	    int relmcount = 0;
	    mpz_mod(km, k, order);

	    //gmp_printf("mod %Zd\n", order);
	    for (i=0; i<relcount; i++) {
		mpz_mod(z0, rel[i]->data, order);
		if (mpz_sgn(z0)) {
		    mpz_set(relm[relmcount]->data, z0);
		    relm[relmcount]->ind = rel[i]->ind;
		    relmcount++;
		}
	    }

	    while (relmcount) {
		//start from the sparse end
		int rind = relm[relmcount - 1]->ind;
		darray_ptr rp = row[rind];

		if (rind < minfound[faci]) break;

		mpz_set(z0, relm[relmcount - 1]->data);
		if (!rp->count) {
		    mpz_invert(z0, z0, order);
		    cell_ptr cnew = newcell();
		    cnew->ind = -1;
		    mpz_mul(z1, km, z0);
		    mpz_mod(cnew->data, z1, order);
		    darray_append(rp, cnew);
		    for (j=0; j<relmcount; j++) {
			cnew = newcell();
			cnew->ind = relm[j]->ind;
			mpz_mul(z1, relm[j]->data, z0);
			mpz_mod(cnew->data, z1, order);
			darray_append(rp, cnew);
		    }
		    count++;
printf("%d / %d\n", count, r * fac->count);
/*
for (i=1; i<rp->count; i++) {
    cnew = rp->item[i];
    gmp_printf(" %u:%Zd", cnew->ind, cnew->data);
}
cnew = rp->item[0];
gmp_printf(" %Zd\n", cnew->data);
*/

		    if (rind == minfound[faci]) {
			do {
			    if (!minfound[faci]) {
			    printf("found log p_%d\n", minfound[faci]);
			    cnew = rp->item[0];
			    gmp_printf("km = %Zd mod %Zd\n", cnew->data, order);
			    }
			    minfound[faci]++;
			    if (minfound[faci] >= r) break;
			    rp = row[minfound[faci]];
			} while (rp->count);
		    }
		    break;

		}

/*
{
//gmp_printf("mod = %Zd\n", order);
printf("before:");
for (i=0; i<relmcount; i++) {
    gmp_printf(" %u:%Zd", relm[i]->ind, relm[i]->data);
}
gmp_printf(" %Zd\n", km);
cell_ptr cp;
printf("sub %d:", rind);
for (i=1; i<rp->count; i++) {
    cp = rp->item[i];
    gmp_printf(" %u:%Zd", cp->ind, cp->data);
}
cp = rp->item[0];
gmp_printf(" %Zd\n", cp->data);
}
*/
		cell_ptr cpi, cpj;
		relmcount--;
		i=0; j=1;
		while (i<relmcount && j<rp->count - 1) {
		    cpi = relm[i];
		    cpj = rp->item[j];
		    if (cpi->ind == cpj->ind) {
			mpz_mul(z1, z0, cpj->data);
			mpz_mod(z1, z1, order);
			int res = mpz_cmp(z1, cpi->data);
			if (!res) {
			    memmove(&relm[i], &relm[i + 1], (relmcount - i - 1) * sizeof(cell_ptr));
			    relm[relmcount - 1] = cpi;
			    relmcount--;
			    j++;
			} else if (res > 0) {
			    mpz_sub(z1, order, z1);
			    mpz_add(cpi->data, cpi->data, z1);
			    i++;
			    j++;
			} else {
			    mpz_sub(cpi->data, cpi->data, z1);
			    i++;
			    j++;
			}
		    } else if (cpi->ind > cpj->ind) {
			cpi = relm[relmcount];
			memmove(&relm[i + 1], &relm[i], (relmcount - i) * sizeof(cell_ptr));
			relm[i] = cpi;
			relmcount++;

			cpi->ind = cpj->ind;
			mpz_mul(z1, z0, cpj->data);
			mpz_mod(z1, z1, order);
			mpz_sub(cpi->data, order, z1);
			//cpi->data = order - ((u0 * cpj->data) % order);
			i++;
			j++;
		    } else {
			i++;
		    }
		}

		if (i == relmcount) {
		    while (j < rp->count - 1) {
			cpi = relm[relmcount];
			cpj = rp->item[j];
			cpi->ind = cpj->ind;
			mpz_mul(z1, z0, cpj->data);
			mpz_mod(z1, z1, order);
			mpz_sub(cpi->data, order, z1);
			//cpi->data = order - ((u0 * cpj->data) % order);
			relmcount++;
			j++;
		    }
		}

		cpj = rp->item[0];
		mpz_mul(z1, z0, cpj->data);
		mpz_mod(z1, z1, order);
		//u1 = (u0 * cpj->data) % order;
		if (mpz_cmp(km, z1) >= 0) {
		    mpz_sub(km, km, z1);
		} else {
		    mpz_sub(z1, order, z1);
		    mpz_add(km, km, z1);
		}

/*
printf("after:");
for (i=0; i<relmcount; i++) {
    gmp_printf(" %u:%Zd", relm[i]->ind, relm[i]->data);
}
gmp_printf(" %Zd\n", km);
*/
	    }
	}

    } while (count < r * fac->count);

    for (faci=0; faci<fac->count; faci++) {
	darray_t *row = matrix[faci];
	mpz_ptr order = fac->item[faci];
	for (i=1; i<r; i++) {
	    darray_ptr rp = row[i];
	    cell_ptr c0 = rp->item[0];
	    for (j=1; j<rp->count-1; j++) {
		cell_ptr cp = rp->item[j];
		darray_ptr r2 = row[cp->ind];
		cell_ptr c2 = r2->item[0];
		mpz_mul(z0, cp->data, c2->data);
		mpz_sub(c0->data, c0->data, z0);
		mpz_mod(c0->data, c0->data, order);
	    }
	}
    }

    mpz_ptr *tmp = pbc_malloc(sizeof(mpz_ptr) * fac->count);
    for (i=0; i<fac->count; i++) {
	tmp[i] = pbc_malloc(sizeof(mpz_t));
	mpz_init(tmp[i]);
	mpz_pow_ui(fac->item[i], fac->item[i], (unsigned int) mul->item[i]);
    }

    for (i=0; i<r; i++) {
	for (faci=0; faci<fac->count; faci++) {
	    darray_t *row = matrix[faci];
	    cell_ptr cp = row[i]->item[0];
	    mpz_set(tmp[faci], cp->data);
	}
	CRT(ind[i], tmp, (mpz_ptr *) fac->item, fac->count);
    }

    for (i=0; i<fac->count; i++) {
	mpz_clear(tmp[i]);
    }
    pbc_free(tmp);

    for (faci=0; i<fac->count; faci++) {
	//similarly ''row'' refers to a list of rows
	darray_t *row = matrix[faci];
	for (j=0; j<r; j++) {
	    darray_forall(row[j], delcell);
	    darray_clear(row[j]);
	}