poly.c

这是一个C的源代码

	element_set0(coeff[i]);
    }
}

static int poly_cmp(element_ptr f, element_ptr g)
{
    int i;
    int n = poly_coeff_count(f);
    int n1 = poly_coeff_count(g);
    if (n != n1) return 1;
    for (i=0; i<n; i++) {
	if (element_cmp(poly_coeff(f, i), poly_coeff(g, i))) return 1;
    }
    return 0;
}

static void field_clear_poly(field_ptr f)
{
    poly_field_data_ptr p = f->data;
    pbc_free(p);
}

//2 bytes hold the number of terms
//then the terms follow
//bad for sparse polynomials

static int poly_length_in_bytes(element_t p)
{
    int count = poly_coeff_count(p);
    int result = 2;
    int i;
    for (i=0; i<count; i++) {
	result += element_length_in_bytes(poly_coeff(p, i));
    }
    return result;
}

static int poly_to_bytes(unsigned char *buf, element_t p)
{
    int count = poly_coeff_count(p);
    int result = 2;
    int i;
    buf[0] = (unsigned char) count;
    buf[1] = (unsigned char) (count >> 8);
    for (i=0; i<count; i++) {
	result += element_to_bytes(&buf[result], poly_coeff(p, i));
    }
    return result;
}

static int poly_from_bytes(element_t p, unsigned char *buf)
{
    int result = 2;
    int count = buf[0] + buf[1] * 256;
    int i;
    poly_alloc(p, count);
    for (i=0; i<count; i++) {
	result += element_from_bytes(poly_coeff(p, i), &buf[result]);
    }
    return result;
}

//Contrived? This returns to_mpz(constant term)
static void poly_to_mpz(mpz_t z, element_ptr e)
{
    if (!poly_coeff_count(e)) {
	mpz_set_ui(z, 0);
    } else {
	element_to_mpz(z, poly_coeff(e, 0));
    }
}

void poly_out_info(FILE *str, field_ptr f)
{
    poly_field_data_ptr p = f->data;
    fprintf(str, "Polynomial ring. Base field:\n");
    field_out_info(str, p->field);
}

void field_init_poly(field_ptr f, field_ptr base_field)
{
    poly_field_data_ptr p;

    field_init(f);
    f->data = pbc_malloc(sizeof(poly_field_data_t));
    p = f->data;
    p->field = base_field;
    p->mapbase = element_field_to_poly;
    f->field_clear = field_clear_poly;
    f->init = poly_init;
    f->clear = poly_clear;
    f->set_si = poly_set_si;
    f->set_mpz = poly_set_mpz;
    f->to_mpz = poly_to_mpz;
    f->out_str = poly_out_str;
    f->snprint = poly_snprint;
    f->set = poly_set;
    f->sign = poly_sgn;
    f->add = poly_add;
    f->doub = poly_double;
    f->is0 = poly_is0;
    f->is1 = poly_is1;
    f->set0 = poly_set0;
    f->set1 = poly_set1;
    f->sub = poly_sub;
    f->neg = poly_neg;
    f->mul = poly_mul;
    f->mul_mpz = poly_mul_mpz;
    f->mul_si = poly_mul_si;
    f->cmp = poly_cmp;
    f->out_info = poly_out_info;

    f->to_bytes = poly_to_bytes;
    f->from_bytes = poly_from_bytes;
    f->fixed_length_in_bytes = -1;
    f->length_in_bytes = poly_length_in_bytes;
}

static void field_clear_polymod(field_ptr f)
{
    polymod_field_data_ptr p = f->data;
    int i, n = p->n;

    for (i=0; i<n; i++) {
	element_clear(p->xpwr[i]);
    }
    pbc_free(p->xpwr);

    element_clear(p->poly);
    pbc_free(f->data);
}

static int polymod_is_sqr(element_ptr e)
{
    int res;
    mpz_t z;
    element_t e0;

    element_init(e0, e->field);
    mpz_init(z);
    mpz_sub_ui(z, e->field->order, 1);
    mpz_divexact_ui(z, z, 2);

    element_pow_mpz(e0, e, z);
    res = element_is1(e0);
    element_clear(e0);
    mpz_clear(z);
    return res;
}

static void polymod_sqrt(element_ptr res, element_ptr a)
    //use Cantor-Zassenhaus aka Legendre's method
    //TODO: use another, faster method?
{
    field_t kx;
    element_t f;
    element_t r, s;
    element_t e0;
    mpz_t z;

    field_init_poly(kx, a->field);
    mpz_init(z);
    element_init(f, kx);
    element_init(r, kx);
    element_init(s, kx);
    element_init(e0, a->field);

    poly_alloc(f, 3);
    element_set1(poly_coeff(f, 2));
    element_neg(poly_coeff(f, 0), a);

    mpz_sub_ui(z, a->field->order, 1);
    mpz_divexact_ui(z, z, 2);
    for (;;) {
	int i;
	element_ptr x;
	element_ptr e1, e2;

	poly_alloc(r, 2);
	element_set1(poly_coeff(r, 1));
	x = poly_coeff(r, 0);
	element_random(x);
	element_mul(e0, x, x);
	if (!element_cmp(e0, a)) {
	    element_set(res, x);
	    break;
	}
	element_set1(s);
	//TODO: this can be optimized greatly
	//since we know r has the form ax + b
	for (i = mpz_sizeinbase(z, 2) - 1; i >=0; i--) {
	    element_mul(s, s, s);
	    if (poly_degree(s) == 2) {
		e1 = poly_coeff(s, 0);
		e2 = poly_coeff(s, 2);
		element_mul(e0, e2, a);
		element_add(e1, e1, e0);
		poly_alloc(s, 2);
		poly_remove_leading_zeroes(s);
	    }
	    if (mpz_tstbit(z, i)) {
		element_mul(s, s, r);
		if (poly_degree(s) == 2) {
		    e1 = poly_coeff(s, 0);
		    e2 = poly_coeff(s, 2);
		    element_mul(e0, e2, a);
		    element_add(e1, e1, e0);
		    poly_alloc(s, 2);
		    poly_remove_leading_zeroes(s);
		}
	    }
	}
	if (poly_degree(s) < 1) continue;
	element_set1(e0);
	e1 = poly_coeff(s, 0);
	e2 = poly_coeff(s, 1);
	element_add(e1, e1, e0);
	element_invert(e0, e2);
	element_mul(e0, e0, e1);
	element_mul(e2, e0, e0);
	if (!element_cmp(e2, a)) {
	    element_set(res, e0);
	    break;
	}
    }

    mpz_clear(z);
    element_clear(f);
    element_clear(r);
    element_clear(s);
    element_clear(e0);
    field_clear(kx);
}

void poly_make_monic(element_t f, element_t g)
{
    int n = poly_coeff_count(g);
    int i;
    element_ptr e0;
    poly_alloc(f, n);
    if (!n) return;

    e0 = poly_coeff(f, n - 1);
    element_invert(e0, poly_coeff(g, n - 1));
    for (i=0; i<n-1; i++) {
	element_mul(poly_coeff(f, i), poly_coeff(g, i), e0);
    }
    element_set1(e0);
}

static int polymod_to_bytes(unsigned char *data, element_t f)
{
    polymod_field_data_ptr p = f->field->data;
    element_t *coeff = f->data;
    int i, n = p->n;
    int len = 0;
    for (i=0; i<n; i++) {
	len += element_to_bytes(data + len, coeff[i]);
    }
    return len;
}

static int polymod_length_in_bytes(element_t f)
{
    polymod_field_data_ptr p = f->field->data;
    element_t *coeff = f->data;
    int i, n = p->n;
    int res = 0;

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

    return res;
}

static int polymod_from_bytes(element_t f, unsigned char *data)
{
    polymod_field_data_ptr p = f->field->data;
    element_t *coeff = f->data;
    int i, n = p->n;
    int len = 0;

    for (i=0; i<n; i++) {
	len += element_from_bytes(coeff[i], data + len);
    }
    return len;
}

static void polymod_init(element_t e)
{
    int i;
    polymod_field_data_ptr p = e->field->data;
    int n = p->n;
    element_t *coeff;
    coeff = e->data = pbc_malloc(sizeof(element_t) * n);

    for (i=0; i<n; i++) {
	element_init(coeff[i], p->field);
    }
}

static void polymod_clear(element_t 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_clear(coeff[i]);
    }
    pbc_free(e->data);
}

static void polymod_set_si(element_t e, signed long int x)
{
    polymod_field_data_ptr p = e->field->data;
    element_t *coeff = e->data;
    int i, n = p->n;
    element_set_si(coeff[0], x);
    for (i=1; i<n; i++) {
	element_set0(coeff[i]);
    }
}

static void polymod_set_mpz(element_t e, mpz_t z)
{
    polymod_field_data_ptr p = e->field->data;
    element_t *coeff = e->data;
    int i, n = p->n;
    element_set_mpz(coeff[0], z);
    for (i=1; i<n; i++) {
	element_set0(coeff[i]);
    }
}

static void polymod_set(element_t e, element_t f)
{
    polymod_field_data_ptr p = e->field->data;
    element_t *dst = e->data, *src = f->data;
    int i, n = p->n;
    for (i=0; i<n; i++) {
	element_set(dst[i], src[i]);
    }
}

static void polymod_neg(element_t e, element_t f)
{
    polymod_field_data_ptr p = e->field->data;
    element_t *dst = e->data, *src = f->data;
    int i, n = p->n;
    for (i=0; i<n; i++) {
	element_neg(dst[i], src[i]);
    }
}

static int polymod_cmp(element_ptr f, element_ptr g)
{
    polymod_field_data_ptr p = f->field->data;
    element_t *c1 = f->data, *c2 = g->data;
    int i, n = p->n;
    for (i=0; i<n; i++) {
	if (element_cmp(c1[i], c2[i])) return 1;
    }
    return 0;
}

static void polymod_add(element_t r, element_t e, element_t f)
{
    polymod_field_data_ptr p = r->field->data;
    element_t *dst = r->data, *s1 = e->data, *s2 = f->data;
    int i, n = p->n;
    for (i=0; i<n; i++) {
	element_add(dst[i], s1[i], s2[i]);
    }
}

static void polymod_double(element_t r, element_t f)
{
    polymod_field_data_ptr p = r->field->data;
    element_t *dst = r->data, *s1 = f->data;
    int i, n = p->n;
    for (i=0; i<n; i++) {
	element_double(dst[i], s1[i]);
    }
}

static void polymod_sub(element_t r, element_t e, element_t f)
{
    polymod_field_data_ptr p = r->field->data;
    element_t *dst = r->data, *s1 = e->data, *s2 = f->data;
    int i, n = p->n;
    for (i=0; i<n; i++) {
	element_sub(dst[i], s1[i], s2[i]);
    }
}

static void polymod_mul_mpz(element_t e, element_t f, mpz_ptr z)
{
    polymod_field_data_ptr p = e->field->data;
    element_t *dst = e->data, *src = f->data;
    int i, n = p->n;
    for (i=0; i<n; i++) {
	element_mul_mpz(dst[i], src[i], z);
    }
}

static void polymod_mul_si(element_t e, element_t f, signed long int z)
{
    polymod_field_data_ptr p = e->field->data;
    element_t *dst = e->data, *src = f->data;
    int i, n = p->n;
    for (i=0; i<n; i++) {
	element_mul_si(dst[i], src[i], z);
    }
}

//Karatsuba for degree 2 polynomials
static void kar_poly_2(element_t *dst, element_t c3, element_t c4, element_t *s1, element_t *s2, element_t *scratch)
{
    element_ptr c01, c02, c12;

    c12 = scratch[0];
    c02 = scratch[1];
    c01 = scratch[2];

    element_add(c3, s1[0], s1[1]);
    element_add(c4, s2[0], s2[1]);
    element_mul(c01, c3, c4);
    element_add(c3, s1[0], s1[2]);
    element_add(c4, s2[0], s2[2]);
    element_mul(c02, c3, c4);
    element_add(c3, s1[1], s1[2]);
    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]);
}

static void polymod_mul_degree3(element_ptr res, element_ptr e, element_ptr f)
{
    polymod_field_data_ptr p = res->field->data;
    element_t *dst = res->data, *s1 = e->data, *s2 = f->data;
    element_t c3, c4;
    element_t p0;

    element_init(p0, res->field);
    element_init(c3, p->field);
    element_init(c4, p->field);

    kar_poly_2(dst, c3, c4, s1, s2, p0->data);

    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_degree6(element_ptr res, element_ptr e, element_ptr f)
{
    polymod_field_data_ptr p = res->field->data;
    element_t *dst = res->data, *s0, *s1 = e->data, *s2 = f->data;
    element_t *a0, *a1, *b0, *b1;
    element_t p0, p1, p2, p3;

    a0 = s1;
    a1 = &s1[3];
    b0 = s2;
    b1 = &s2[3];

    element_init(p0, res->field);
    element_init(p1, res->field);
    element_init(p2, res->field);
    element_init(p3, res->field);

    s0 = p0->data;
    s1 = p1->data;
    s2 = p2->data;
    element_add(s0[0], a0[0], a1[0]);
    element_add(s0[1], a0[1], a1[1]);
    element_add(s0[2], a0[2], a1[2]);

    element_add(s1[0], b0[0], b1[0]);
    element_add(s1[1], b0[1], b1[1]);
    element_add(s1[2], b0[2], b1[2]);

    kar_poly_2(s2, s2[3], s2[4], s0, s1, p3->data);
    kar_poly_2(s0, s0[3], s0[4], a0, b0, p3->data);
    kar_poly_2(s1, s1[3], s1[4], a1, b1, p3->data);

    element_set(dst[0], s0[0]);
    element_set(dst[1], s0[1]);
    element_set(dst[2], s0[2]);

    element_sub(dst[3], s0[3], s0[0]);
    element_sub(dst[3], dst[3], s1[0]);
    element_add(dst[3], dst[3], s2[0]);

    element_sub(dst[4], s0[4], s0[1]);
    element_sub(dst[4], dst[4], s1[1]);
    element_add(dst[4], dst[4], s2[1]);

    element_sub(dst[5], s2[2], s0[2]);
    element_sub(dst[5], dst[5], s1[2]);

    //start reusing part of s0 as scratch(!)
    element_sub(s0[0], s2[3], s0[3]);
    element_sub(s0[0], s0[0], s1[3]);
    element_add(s0[0], s0[0], s1[0]);

    element_sub(s0[1], s2[4], s0[4]);
    element_sub(s0[1], s0[1], s1[4]);
    element_add(s0[1], s0[1], s1[1]);

    polymod_const_mul(p3, s0[0], p->xpwr[0]);
    element_add(res, res, p3);
    polymod_const_mul(p3, s0[1], p->xpwr[1]);
    element_add(res, res, p3);
    polymod_const_mul(p3, s1[2], p->xpwr[2]);
    element_add(res, res, p3);
    polymod_const_mul(p3, s1[3], p->xpwr[3]);
    element_add(res, res, p3);
    polymod_const_mul(p3, s1[4], p->xpwr[4]);
    element_add(res, res, p3);

    element_clear(p0);
    element_clear(p1);
    element_clear(p2);
    element_clear(p3);
}

/*
void polymod_mul_degree3(element_ptr res, element_ptr e, element_ptr f)
{
    polymod_field_data_ptr p = res->field->data;
    element_t *dst = res->data, *s1 = e->data, *s2 = f->data;
    element_t *p0data;
    element_ptr c01, c02, c12;
    element_t c3, c4;
    element_t p0;

    element_init(p0, res->field);
    element_init(c3, p->field);
    element_init(c4, p->field);
    p0data = p0->data;

    c12 = p0data[0];
    c02 = p0data[1];
    c01 = p0data[2];

    element_add(c3, s1[0], s1[1]);
    element_add(c4, s2[0], s2[1]);
    element_mul(c01, c3, c4);
    element_add(c3, s1[0], s1[2]);
    element_add(c4, s2[0], s2[2]);
    element_mul(c02, c3, c4);
    element_add(c3, s1[1], s1[2]);