poly.c

这是一个C的源代码

    element_add(c4, s2[1], s2[2]);
    element_mul(c12, c3, c4);

    element_mul(dst[1], s1[1], s2[1]);

    //constant term
    element_mul(dst[0], s1[0], s2[0]);

    //coefficient of x^4
    element_mul(c4, s1[2], s2[2]);

    //coefficient of x^3
    element_add(c3, dst[1], c4);
    element_sub(c3, c12, c3);

    //coefficient of x^2
    element_add(dst[2], c4, dst[0]);
    element_sub(c02, c02, dst[2]);
    element_add(dst[2], dst[1], c02);

    //coefficient of x
    element_sub(c01, c01, dst[0]);
    element_sub(dst[1], c01, dst[1]);

    polymod_const_mul(p0, c3, p->xpwr[0]);
    element_add(res, res, p0);
    polymod_const_mul(p0, c4, p->xpwr[1]);
    element_add(res, res, p0);

    element_clear(p0);
    element_clear(c3);
    element_clear(c4);
}
*/

static void polymod_mul(element_ptr res, element_ptr e, element_ptr f)
//NOTE: cases n = 3, 6 have dedicated routines
{
    polymod_field_data_ptr p = res->field->data;
    int n = p->n;
    element_t *dst;
    element_t *s1 = e->data, *s2 = f->data;
    element_t prod, p0, c0;
    int i, j;
    element_t *high; //coefficients of x^n,...,x^{2n-2}

    high = pbc_malloc(sizeof(element_t) * (n - 1));
    for (i=0; i<n-1; i++) {
	element_init(high[i], p->field);
	element_set0(high[i]);
    }
    element_init(prod, res->field);
    dst = prod->data;
    element_init(p0, res->field);
    element_init(c0, p->field);

    for (i=0; i<n; i++) {
	int ni = n - i;
	for (j=0; j<ni; j++) {
	    element_mul(c0, s1[i], s2[j]);
	    element_add(dst[i + j], dst[i + j], c0);
	}
	for (;j<n; j++) {
	    element_mul(c0, s1[i], s2[j]);
	    element_add(high[j - ni], high[j - ni], c0);
	}
    }

    for (i=0; i<n-1; i++) {
	polymod_const_mul(p0, high[i], p->xpwr[i]);
	element_add(prod, prod, p0);
	element_clear(high[i]);
    }
    pbc_free(high);

    element_set(res, prod);
    element_clear(prod);
    element_clear(p0);
    element_clear(c0);
}

static void polymod_square_degree3(element_ptr res, element_ptr e)
//TODO: investigate if squaring is significantly cheaper
//than multiplication. If so convert to Karatsuba.
{
    element_t *dst = res->data;
    element_t *src = e->data;
    polymod_field_data_ptr p = res->field->data;
    element_t p0;
    element_t c0, c2;
    element_ptr c1, c3;

    element_init(p0, res->field);
    element_init(c0, p->field);
    element_init(c2, p->field);

    c3 = p0->data;
    c1 = c3 + 1;

    element_mul(c3, src[0], src[1]);
    element_mul(c1, src[0], src[2]);
    element_square(dst[0], src[0]);

    element_mul(c2, src[1], src[2]);
    element_square(c0, src[2]);
    element_square(dst[2], src[1]);

    element_add(dst[1], c3, c3);

    element_add(c1, c1, c1);
    element_add(dst[2], dst[2], c1);

    polymod_const_mul(p0, c0, p->xpwr[1]);
    element_add(res, res, p0);

    element_add(c2, c2, c2);
    polymod_const_mul(p0, c2, p->xpwr[0]);
    element_add(res, res, p0);

    element_clear(p0);
    element_clear(c0);
    element_clear(c2);
}

static void polymod_square(element_ptr res, element_ptr e)
{
    element_t *dst;
    element_t *src = e->data;
    polymod_field_data_ptr p = res->field->data;
    int n = p->n;

    element_t prod, p0, c0;
    int i, j;

    element_t *high; //coefficients of x^n,...,x^{2n-2}

    high = pbc_malloc(sizeof(element_t) * (n - 1));
    for (i=0; i<n-1; i++) {
	element_init(high[i], p->field);
	element_set0(high[i]);
    }

    element_init(prod, res->field);
    dst = prod->data;
    element_init(p0, res->field);
    element_init(c0, p->field);

    for (i=0; i<n; i++) {
	int twicei = 2 * i;
	element_square(c0, src[i]);
	if (twicei < n) {
	    element_add(dst[twicei], dst[twicei], c0);
	} else {
	    element_add(high[twicei - n], high[twicei - n], c0);
	}

	for (j=i+1; j<n-i; j++) {
	    element_mul(c0, src[i], src[j]);
	    element_add(c0, c0, c0);
	    element_add(dst[i + j], dst[i + j], c0);
	}
	for (;j<n; j++) {
	    element_mul(c0, src[i], src[j]);
	    element_add(c0, c0, c0);
	    element_add(high[i + j - n], high[i + j - n], c0);
	}
    }

    for (i=0; i<n-1; i++) {
	polymod_const_mul(p0, high[i], p->xpwr[i]);
	element_add(prod, prod, p0);
	element_clear(high[i]);
    }
    pbc_free(high);

    element_set(res, prod);
    element_clear(prod);
    element_clear(p0);
    element_clear(c0);
}

static int polymod_is0(element_ptr e)
{
    polymod_field_data_ptr p = e->field->data;
    element_t *coeff = e->data;
    int i, n = p->n;

    for (i=0; i<n; i++) {
	if (!element_is0(coeff[i])) return 0;
    }
    return 1;
}

static int polymod_is1(element_ptr e)
{
    polymod_field_data_ptr p = e->field->data;
    element_t *coeff = e->data;
    int i, n = p->n;

    if (!element_is1(coeff[0])) return 0;
    for (i=1; i<n; i++) {
	if (!element_is0(coeff[i])) return 0;
    }
    return 1;
}

static void polymod_set0(element_ptr e)
{
    polymod_field_data_ptr p = e->field->data;
    element_t *coeff = e->data;
    int i, n = p->n;

    for (i=0; i<n; i++) {
	element_set0(coeff[i]);
    }
}

static void polymod_set1(element_ptr e)
{
    polymod_field_data_ptr p = e->field->data;
    element_t *coeff = e->data;
    int i, n = p->n;

    element_set1(coeff[0]);
    for (i=1; i<n; i++) {
	element_set0(coeff[i]);
    }
}

static int polymod_sgn(element_ptr e)
{
    polymod_field_data_ptr p = e->field->data;
    element_t *coeff = e->data;
    int res = 0;
    int i, n = p->n;
    for (i=0; i<n; i++) {
	res = element_sgn(coeff[i]);
	if (res) break;
    }
    return res;
}

static size_t polymod_out_str(FILE *stream, int base, element_ptr e)
{
    size_t result = 2, status;
    polymod_field_data_ptr p = e->field->data;
    element_t *coeff = e->data;
    int i, n = p->n;

    if (EOF == fputc('[', stream)) return 0;
    for (i=0; i<n; i++) {
	if (i) {
	    if (EOF == fputs(", ", stream)) return 0;
	    result += 2;
	}
	status = element_out_str(stream, base, coeff[i]);
	if (!status) return 0;
	result += status;
    }
    if (EOF == fputc(']', stream)) return 0;
    return result;
}

static int polymod_snprint(char *s, size_t size, element_ptr e)
{
    polymod_field_data_ptr p = e->field->data;
    element_t *coeff = e->data;
    int i, n = p->n;
    size_t result = 0, left;
    int status;

    void clip_sub(void)
    {
	result += status;
	left = result >= size ? 0 : size - result;
    }

    status = snprintf(s, size, "[");
    if (status < 0) return status;
    clip_sub();

    for (i=0; i<n; i++) {
	if (i) {
	    status = snprintf(s + result, left, ", ");
	    if (status < 0) return status;
	    clip_sub();
	}
	status = element_snprint(s + result, left, coeff[i]);
	if (status < 0) return status;
	clip_sub();
    }
    status = snprintf(s + result, left, "]");
    if (status < 0) return status;
    return result + status;
}

static int polymod_set_str(element_ptr e, char *s, int base)
{
    polymod_field_data_ptr p = e->field->data;
    element_t *coeff = e->data;
    int i, n = p->n;
    char *cp = s;
    element_set0(e);
    while (*cp && isspace(*cp)) cp++;
    if (*cp++ != '[') return 0;
    for (i=0; i<n; i++) {
	cp += element_set_str(coeff[i], cp, base);
	while (*cp && isspace(*cp)) cp++;
	if (i<n-1 && *cp++ != ',') return 0;
    }
    if (*cp++ != ']') return 0;
    return cp - s;
}

//Contrived? This returns to_mpz(constant term)
static void polymod_to_mpz(mpz_t z, element_ptr e)
{
    element_to_mpz(z, polymod_coeff(e, 0));
}

static void compute_x_powers(field_ptr field, element_ptr poly)
//compute x^n,...,x^{2n-2} mod poly
{
    polymod_field_data_ptr p = field->data;
    element_t p0;
    element_ptr pwrn;
    element_t *coeff, *coeff1;
    int i, j;
    int n = p->n;
    element_t *xpwr;
    
    xpwr = p->xpwr;

    element_init(p0, field);
    for (i=0; i<n; i++) {
	element_init(xpwr[i], field);
    }
    pwrn = xpwr[0];
    element_poly_to_polymod_truncate(pwrn, poly);
    element_neg(pwrn, pwrn);

    for (i=1; i<n; i++) {
	coeff = xpwr[i-1]->data;
	coeff1 = xpwr[i]->data;

	element_set0(coeff1[0]);
	for (j=1; j<n; j++) {
	    element_set(coeff1[j], coeff[j - 1]);
	}
	polymod_const_mul(p0, coeff[n - 1], pwrn);
	element_add(xpwr[i], xpwr[i], p0);
    }
    element_clear(p0);
}

void polymod_out_info(FILE *str, field_ptr f)
{
    polymod_field_data_ptr p = f->data;
    element_fprintf(str, "Field extension. Min poly = %B\n", p->poly);
    fprintf(str, "Base field:\n");
    field_out_info(str, p->field);
}

void field_init_polymod(field_ptr f, element_ptr poly)
    //assumes poly is monic
{
    poly_field_data_ptr pdp = poly->field->data;
    polymod_field_data_ptr p;
    int n;

    field_init(f);
    f->data = pbc_malloc(sizeof(polymod_field_data_t));
    p = f->data;
    p->field = pdp->field;
    p->mapbase = element_field_to_poly;
    element_init(p->poly, poly->field);
    element_set(p->poly, poly);
    n = p->n = poly_degree(p->poly);
    f->field_clear = field_clear_polymod;
    f->init = polymod_init;
    f->clear = polymod_clear;
    f->set_si = polymod_set_si;
    f->set_mpz = polymod_set_mpz;
    f->out_str = polymod_out_str;
    f->snprint = polymod_snprint;
    f->set_str = polymod_set_str;
    f->set = polymod_set;
    f->sign = polymod_sgn;
    f->add = polymod_add;
    f->doub = polymod_double;
    f->sub = polymod_sub;
    f->neg = polymod_neg;
    f->is0 = polymod_is0;
    f->is1 = polymod_is1;
    f->set0 = polymod_set0;
    f->set1 = polymod_set1;
    f->cmp = polymod_cmp;
    f->to_mpz = polymod_to_mpz;
    switch(n) {
	case 3:
	    f->mul = polymod_mul_degree3;
	    f->square = polymod_square_degree3;
	    break;
	case 6:
	    f->mul = polymod_mul_degree6;
	    f->square = polymod_square;
	    break;
	default:
	    f->mul = polymod_mul;
	    f->square = polymod_square;
	    break;
    }

    f->mul_mpz = polymod_mul_mpz;
    f->mul_si = polymod_mul_si;
    f->random = polymod_random;
    f->from_hash = polymod_from_hash;
    f->invert = polymod_invert;
    f->is_sqr = polymod_is_sqr;
    f->sqrt = polymod_sqrt;
    f->to_bytes = polymod_to_bytes;
    f->from_bytes = polymod_from_bytes;
    f->out_info = polymod_out_info;

    if (pdp->field->fixed_length_in_bytes < 0) {
	f->fixed_length_in_bytes = -1;
	f->length_in_bytes = polymod_length_in_bytes;
    } else {
	f->fixed_length_in_bytes = pdp->field->fixed_length_in_bytes * poly_degree(poly);
    }
    mpz_pow_ui(f->order, p->field->order, n);

    p->xpwr = pbc_malloc(sizeof(element_t) * n);
    compute_x_powers(f, poly);
}

void trial_divide(darray_ptr factor, darray_ptr mult, mpz_t n, mpz_ptr limit)
{
    mpz_t p, m;
    mpz_ptr fac;
    unsigned int mul;

    mpz_init(p);
    mpz_init(m);
    mpz_set(m ,n);
    mpz_set_ui(p, 2);

    while (mpz_cmp_ui(m, 1)) {
	if (mpz_probab_prime_p(m, 10)) {
	    mpz_set(p, m);
	}
	if (limit && mpz_cmp(p, limit) > 0) {
	    mpz_set(p, m);
	}
	if (mpz_divisible_p(m, p)) {
	    fac = pbc_malloc(sizeof(mpz_t));
	    mul = 0;
	    mpz_init(fac);
	    mpz_set(fac, p);
	    darray_append(factor, fac);
	    do {
		mpz_divexact(m, m, p);
		mul++;
	    } while (mpz_divisible_p(m, p));
	    darray_append(mult, int_to_voidp(mul));
	}
	mpz_nextprime(p, p);
    }
   
    mpz_clear(m);
    mpz_clear(p);
}

void poly_gcd(element_ptr d, element_ptr f, element_ptr g)
{
    element_t a, b, q, r;
    element_init(a, d->field);
    element_init(b, d->field);
    element_init(q, d->field);
    element_init(r, d->field);

    element_set(a, f);
    element_set(b, g);
    for(;;) {
	//TODO: don't care about q
	poly_div(q, r, a, b);
	if (element_is0(r)) break;
	element_set(a, b);
	element_set(b, r);
    }
    element_set(d, b);
    element_clear(a);
    element_clear(b);
    element_clear(q);
    element_clear(r);
}

int poly_is_irred_degfac(element_ptr f, darray_t factor)
    //called by poly_is_irred
    //needs to be passed a list of the factors of deg f
{
    int res;
    element_t xpow, x, g;
    mpz_t deg, z;
    int i, n;
    field_ptr basef = poly_base_field(f);
    field_t rxmod;

    mpz_init(deg);
    mpz_init(z);

    field_init_polymod(rxmod, f);

    mpz_set_ui(deg, poly_degree(f));
    element_init(xpow, rxmod);
    element_init(x, rxmod);
    element_init(g, f->field);
    n = factor->count;
    //element_set1(((element_t *) x->data)[1]);
    element_set1(polymod_coeff(x, 1));

    res = 0;
    for (i=0; i<n; i++) {
	mpz_divexact(z, deg, factor->item[i]);
	mpz_pow_ui(z, basef->order, mpz_get_ui(z));
	element_pow_mpz(xpow, x, z);
	element_sub(xpow, xpow, x);
	if (element_is0(xpow)) {
	    goto done;
	}
	element_polymod_to_poly(g, xpow);
	poly_gcd(g, f, g);
	if (poly_degree(g) != 0) goto done;
    }

    mpz_pow_ui(z, basef->order, poly_degree(f));
    element_pow_mpz(xpow, x, z);
    element_sub(xpow, xpow, x);
    if (element_is0(xpow)) res = 1;

done:
    element_clear(g);
    element_clear(xpow);
    element_clear(x);
    mpz_clear(deg);
    mpz_clear(z);
    field_clear(rxmod);
    return res;
}

int poly_is_irred(element_ptr f)
{
    darray_t fac, mul;
    mpz_t deg;
    int res;
    void clear(void *p) {
	mpz_clear(p);
	pbc_free(p);
    }

    if (poly_degree(f) <= 1) return 1;

    darray_init(fac);
    darray_init(mul);
    mpz_init(deg);

    mpz_set_ui(deg, poly_degree(f));

    trial_divide(fac, mul, deg, NULL);
    res = poly_is_irred_degfac(f, fac);

    darray_forall(fac, clear);

    darray_clear(fac);
    darray_clear(mul);
    mpz_clear(deg);
    return res;
}