indexcalculus.c

这是一个C的源代码

	pbc_free(matrix[faci]);
    }

    for (i=0; i<r; i++) {
	delcell(rel[i]);
	delcell(relm[i]);
    }

    pbc_free(prime);
    pbc_free(rel);
    pbc_free(relm);
    mpz_clear(k);
    mpz_clear(km);
    mpz_clear(z);
    mpz_clear(z0);
    mpz_clear(z1);
}

// brute-force: does not use the fact that matrices are sparse
void slow_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;
    //mpz_t rel[r + 1];
    //mpz_t relm[r + 1];
    mpz_t *rel = pbc_malloc(sizeof(mpz_t) * (r + 1));
    mpz_t *relm = pbc_malloc(sizeof(mpz_t) * (r + 1));
    unsigned int *prime = pbc_malloc(sizeof(unsigned int) * r);
    darray_t matrix;
    int faci;
    mpz_t k;
    int minfound[fac->count];

    for (i=0; i<r+1; i++) mpz_init(rel[i]);
    for (i=0; i<r+1; i++) mpz_init(relm[i]);

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

    darray_init(matrix);

    for (i=0; i<fac->count; i++) {
	darray_append(matrix, pbc_malloc(r * sizeof(mpz_t *)));
	minfound[i] = 0;
    }

    for (j=0; j<fac->count; j++) {
	mpz_t **row = matrix->item[j];
	for (i=0; i<r; i++) row[i] = NULL;
    }

    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);
    }
    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);
	for (i=0; i<r; i++) {
	    mpz_set_ui(rel[i], 0);
	    while (mpz_divisible_ui_p(z1, prime[i])) {
		mpz_add_ui(rel[i], rel[i], 1);
		mpz_divexact_ui(z1, z1, prime[i]);
	    }
	}
	if (mpz_cmp_ui(z1, 1)) {
	    continue;
	}
	mpz_set(rel[r], k);

/*
	printf("found r-smooth number after %d tries\n", try);
	gmp_printf("g^%Zd = %Zd:", rel[r], z);
	for (i=0; i<r; i++) {
	    if (mpz_sgn(rel[i])) {
		gmp_printf(" %u:%Zd", i, rel[i]);
	    }
	}
	printf("\n");
*/

	try = 0;

	for (faci=0; faci<fac->count; faci++) {
	    mpz_t **row = matrix->item[faci];
	    mpz_ptr order = fac->item[faci];
	    //gmp_printf("mod %Zd\n", order);
	    for (i=0; i<r+1; i++) {
		mpz_mod(relm[i], rel[i], order);
	    }

	    for (;;) {
		/*
		for (i=0; i<r && !mpz_sgn(relm[i]); i++);
		if (i == r) {
		    //printf("redundant relation\n");
		    break;
		}
		*/
		for (i=r-1; i>=0 && !mpz_sgn(relm[i]); i--);
		if (i < 0) {
		    //printf("redundant relation\n");
		    break;
		}
		if (i < minfound[faci]) {
		    break;
		}
		mpz_set(z0, relm[i]);
		if (!row[i]) {
		    row[i] = pbc_malloc(sizeof(mpz_t) * (r + 1));
		    mpz_invert(z1, z0, order);
		    for (j=0; j<r+1; j++) {
			mpz_init(row[i][j]);
			mpz_mul(row[i][j], z1, relm[j]);
			mpz_mod(row[i][j], row[i][j], order);
		    }
		    count++;
printf("%d / %d\n", count, r * fac->count);
/*
for (j=0; j<r; j++) {
    if (mpz_sgn(row[i][j])) {
	gmp_printf(" %d:%Zd", j, row[i][j]);
    }
}
gmp_printf(" %Zd\n", row[i][j]);
*/

		    if (i == minfound[faci]) {
			do {
			    if (!minfound[faci]) {
			    printf("found log p_%d\n", minfound[faci]);
			    gmp_printf("km = %Zd mod %Zd\n", row[i][r], order);
			    }
			    minfound[faci]++;
			    if (minfound[faci] >= r) break;
			} while (row[minfound[faci]]);
		    }
		    break;
		}

		/*
		printf("before:");
		for (j=0; j<r; j++) {
		    if (mpz_sgn(relm[j])) {
			gmp_printf(" %d:%Zd", j, relm[j]);
		    }
		}
		gmp_printf(" %Zd\n", relm[j]);

		printf("sub %d:", i);
		for (j=0; j<r; j++) {
		    if (mpz_sgn(row[i][j])) {
			gmp_printf(" %d:%Zd", j, row[i][j]);
		    }
		}
		gmp_printf(" %Zd\n", row[i][j]);
		*/

		for (j=0; j<r+1; j++) {
		    mpz_mul(z1, z0, row[i][j]);
		    mpz_sub(relm[j], relm[j], z1);
		    mpz_mod(relm[j], relm[j], order);
		}

		/*
		printf("after:");
		for (j=0; j<r; j++) {
		    if (mpz_sgn(relm[j])) {
			gmp_printf(" %d:%Zd", j, relm[j]);
		    }
		}
		gmp_printf(" %Zd\n", relm[j]);
		*/
	    }
	}

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

    for (faci=0; faci<fac->count; faci++) {
	mpz_t **row = matrix->item[faci];
	mpz_ptr order = fac->item[faci];
	/*
	gmp_printf("mod %Zd:\n", order);
	for (i=0; i<r; i++) {
	    for (j=0; j<r+1; j++) {
		gmp_printf(" %Zd", row[i][j]);
	    }
	    printf("\n");
	}
	printf("\n");
	*/

	for (i=1; i<r; i++) {
	    for (j=0; j<i; j++) {
		if (mpz_sgn(row[i][j])) {
		    mpz_mul(z0, row[i][j], row[j][r]);
		    mpz_sub(row[i][r], row[i][r], z0);
		    mpz_mod(row[i][r], row[i][r], order);
		}
	    }
	}
	/*
	for (i=r-2; i>=0; i--) {
	    for (j=i+1; j<r; j++) {
		if (mpz_sgn(row[i][j])) {
		    mpz_mul(z0, row[i][j], row[j][r]);
		    mpz_sub(row[i][r], row[i][r], z0);
		    mpz_mod(row[i][r], row[i][r], order);
		}
	    }
	}
	*/

	/*
	for (i=0; i<r; i++) {
	    mpz_set(rel[i], row[i][r]);
	    gmp_printf(" %Zd", row[i][r]);
	    printf("\n");
	}
	*/
    }

    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++) {
	    mpz_t **row = matrix->item[faci];
	    mpz_set(tmp[faci], row[i][r]);
	}
	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; faci<matrix->count; faci++) {
	mpz_t **row = matrix->item[faci];
	for (j=0; j<r; j++) {
	    for (i=0; i<r+1; i++) {
		mpz_clear(row[j][i]);
	    }
	    pbc_free(row[j]);
	}
	pbc_free(row);
    }
    darray_clear(matrix);
    for (i=0; i<r+1; i++) mpz_clear(rel[i]);
    for (i=0; i<r+1; i++) mpz_clear(relm[i]);
    pbc_free(prime);
    pbc_free(rel);
    pbc_free(relm);
    mpz_clear(k);
    mpz_clear(z);
    mpz_clear(z0);
    mpz_clear(z1);

    printf("step 1 completed\n");
    for (i=0; i<r; i++) element_printf(" %Zd", ind[i]);
    printf("\n");
}

static void index_calculus_step2(mpz_t x, mpz_t *ind, int r,
	mpz_t g, mpz_t h, mpz_t q)
{
    mpz_t prime;
    mpz_t s;
    mpz_t z, z1;
    mpz_t rel[r];
    int i;

    mpz_init(z);
    mpz_init(z1);
    mpz_init(s);
    mpz_init(prime);
    for (i=0; i<r; i++) mpz_init(rel[i]);

    mpz_set(z, h);

    for (;;) {
	mpz_mul(z, z, g);
	mpz_mod(z, z, q);
	mpz_add_ui(s, s, 1);

	mpz_set(z1, z);
	mpz_set_ui(prime, 1);
	for (i=0; i<r; i++) {
	    mpz_set_ui(rel[i], 0);
	    mpz_nextprime(prime, prime);
	    while (mpz_divisible_p(z1, prime)) {
		mpz_add_ui(rel[i], rel[i], 1);
		mpz_divexact(z1, z1, prime);
	    }
	}
	if (mpz_cmp_ui(z1, 1)) continue;
	element_printf("found r-smooth number on try #%Zd\n", s);
	mpz_set_ui(x, 0);
	for (i=0; i<r; i++) {
	    mpz_mul(z, rel[i], ind[i]);
	    mpz_add(x, x, z);
	}
	mpz_sub(x, x, s);
	mpz_sub_ui(z, q, 1);
	mpz_mod(x, x, z);
	break;
    }
}

void index_calculus_dlog(mpz_t x, mpz_t g, mpz_t h, mpz_t q)
{
    darray_t fac, mul;
    int i, r;
    mpz_t q1, z0;

    void mpzclear(void *p)
    {
	mpz_clear(p);
	pbc_free(p);
    }

    darray_init(fac);
    darray_init(mul);

    mpz_init(q1);
    mpz_init(z0);

    mpz_sub_ui(q1, q, 1);
    mpz_setbit(z0, 6);
    trial_divide(fac, mul, q1, z0);

    for (i=0; i<mul->count; i++) {
	unsigned int m = (unsigned int) mul->item[i];
	if (m != 1) {
	    //TODO
	    printf("p-adics not implemented yet\n");
	    return;
	}
    }

    {
	double dq = mpz_get_d(q);
	//r = exp(sqrt(log(dq)*log(log(dq))));
	//printf("r = %d\n", r);
	r = exp(1.2 * sqrt(log(dq)));
	printf("r = %d\n", r);
    }
    mpz_t *ind = pbc_malloc(sizeof(mpz_t) * r);
    for (i=0; i<r; i++) mpz_init(ind[i]);

    if (is_gen(g, q, fac, mul)) {

	index_calculus_step1(ind, r, g, q, fac, mul);

	index_calculus_step2(x, ind, r, g, h, q);
    } else {
	mpz_t y, z;
	mpz_t d;

	mpz_init(d);
	mpz_init(y);
	mpz_init(z);
	do {
	    pbc_mpz_random(z, q);
	} while (!is_gen(z, q, fac, mul));

	gmp_printf("new gen: %Zd\n", z);

	index_calculus_step1(ind, r, z, q, fac, mul);
	//slow_index_calculus_step1(ind, r, z, q, fac, mul);

	index_calculus_step2(x, ind, r, z, g, q);
	index_calculus_step2(y, ind, r, z, h, q);
	//want y / x mod q-1
	mpz_gcd(d, x, q1);
	mpz_divexact(q1, q1, d);
	mpz_divexact(x, x, d);
	//if y not divisible by d there is no solution
	mpz_divexact(y, y, d);
	mpz_invert(x, x, q1);
	mpz_mul(x, y, x);
	mpz_mod(x, x, q1);

	do {
	    mpz_powm(z0, g, x, q);
	    if (!mpz_cmp(z0, h)) {
		break;
	    }
	    mpz_add(x, x, q1);
	    mpz_sub_ui(d, d, 1);
	} while (mpz_sgn(d));

	mpz_clear(d);
	mpz_clear(y);
	mpz_clear(z);
    }

    for (i=0; i<r; i++) mpz_clear(ind[i]);
    pbc_free(ind);

    darray_forall(fac, mpzclear);
    darray_clear(mul);
    darray_clear(fac);
    mpz_clear(q1);
    mpz_clear(z0);
}