poly.c

这是一个C的源代码

#include <ctype.h>
#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include "pbc_field.h"
#include "pbc_darray.h"
#include "pbc_poly.h"
#include "pbc_utils.h"
#include "pbc_memory.h"

//implements R[x] for a given ring R
//also R[x]_{f(x)}

void poly_alloc(element_ptr e, int n)
{
    poly_field_data_ptr pdp = e->field->data;
    poly_element_ptr p = e->data;
    element_ptr e0;
    int k;
    k = p->coeff->count;
    while (k < n) {
	e0 = pbc_malloc(sizeof(element_t));
	element_init(e0, pdp->field);
	darray_append(p->coeff, e0);
	k++;
    }
    while (k > n) {
	k--;
	e0 = darray_at(p->coeff, k);
	element_clear(e0);
	pbc_free(e0);
	darray_remove_last(p->coeff);
    }
}

static void poly_init(element_ptr e)
{
    poly_element_ptr p;
    e->data = pbc_malloc(sizeof(poly_element_t));
    p = e->data;
    darray_init(p->coeff);
}

static void poly_clear(element_ptr e)
{
    poly_element_ptr p = e->data;

    poly_alloc(e, 0);
    darray_clear(p->coeff);
    pbc_free(e->data);
}

void poly_remove_leading_zeroes(element_ptr e)
{
    poly_element_ptr p = e->data;
    int n = p->coeff->count - 1;
    while (n >= 0) {
	element_ptr e0 = p->coeff->item[n];
	if (!element_is0(e0)) return;
	element_clear(e0);
	pbc_free(e0);
	darray_remove_last(p->coeff);
	n--;
    }
}

static void poly_set_si(element_ptr e, signed long int op)
{
    poly_element_ptr p = e->data;
    element_ptr e0;

    poly_alloc(e, 1);
    e0 = p->coeff->item[0];
    element_set_si(e0, op);
    poly_remove_leading_zeroes(e);
}

static void poly_set_mpz(element_ptr e, mpz_ptr op)
{
    poly_element_ptr p = e->data;
    element_ptr e0;

    poly_alloc(e, 1);
    e0 = p->coeff->item[0];
    element_set_mpz(e0, op);
    poly_remove_leading_zeroes(e);
}

void poly_set_coeff(element_ptr e, element_ptr a, int n)
{
    poly_element_ptr p = e->data;
    element_ptr e0;
    if (p->coeff->count < n + 1) {
	poly_alloc(e, n + 1);
    }
    e0 = p->coeff->item[n];
    element_set(e0, a);
}

void poly_setx(element_ptr f)
{
    poly_alloc(f, 2);
    element_set0(poly_coeff(f, 0));
    element_set1(poly_coeff(f, 1));
}

static void poly_set(element_ptr dst, element_ptr src)
{
    poly_element_ptr psrc = src->data;
    poly_element_ptr pdst = dst->data;
    int i;

    poly_alloc(dst, psrc->coeff->count);
    for (i=0; i<psrc->coeff->count; i++) {
	element_set(pdst->coeff->item[i], psrc->coeff->item[i]);
    }
}

static int poly_sgn(element_ptr f)
{
    int res = 0;
    int i;
    int n = poly_coeff_count(f);
    for (i=0; i<n; i++) {
	res = element_sgn(poly_coeff(f, i));
	if (res) break;
    }
    return res;
}

static void poly_add(element_ptr sum, element_ptr f, element_ptr g)
{
    int i, n, n1;
    element_ptr big;

    n = poly_coeff_count(f);
    n1 = poly_coeff_count(g);
    if (n > n1) {
	big = f;
	n = n1;
	n1 = poly_coeff_count(f);
    } else {
	big = g;
    }

    poly_alloc(sum, n1);
    for (i=0; i<n; i++) {
	element_add(poly_coeff(sum, i), poly_coeff(f, i), poly_coeff(g, i));
    }
    for (; i<n1; i++) {
	element_set(poly_coeff(sum, i), poly_coeff(big, i));
    }
    poly_remove_leading_zeroes(sum);
}

static void poly_sub(element_ptr diff, element_ptr f, element_ptr g)
{
    int i, n, n1;
    element_ptr big;

    n = poly_coeff_count(f);
    n1 = poly_coeff_count(g);
    if (n > n1) {
	big = f;
	n = n1;
	n1 = poly_coeff_count(f);
    } else {
	big = g;
    }

    poly_alloc(diff, n1);
    for (i=0; i<n; i++) {
	element_sub(poly_coeff(diff, i), poly_coeff(f, i), poly_coeff(g, i));
    }
    for (; i<n1; i++) {
	if (big == f) {
	    element_set(poly_coeff(diff, i), poly_coeff(big, i));
	} else {
	    element_neg(poly_coeff(diff, i), poly_coeff(big, i));
	}
    }
    poly_remove_leading_zeroes(diff);
}

static void poly_neg(element_ptr f, element_ptr g)
{
    poly_element_ptr pf = f->data;
    poly_element_ptr pg = g->data;
    int i, n;

    n = pg->coeff->count;
    poly_alloc(f, n);
    for (i=0; i<n; i++) {
	element_neg(pf->coeff->item[i], pg->coeff->item[i]);
    }
}

static void poly_double(element_ptr f, element_ptr g)
{
    poly_element_ptr pf = f->data;
    poly_element_ptr pg = g->data;
    int i, n;

    n = pg->coeff->count;
    poly_alloc(f, n);
    for (i=0; i<n; i++) {
	element_double(pf->coeff->item[i], pg->coeff->item[i]);
    }
}

static void poly_mul_mpz(element_ptr f, element_ptr g, mpz_ptr z)
{
    poly_element_ptr pf = f->data;
    poly_element_ptr pg = g->data;
    int i, n;

    n = pg->coeff->count;
    poly_alloc(f, n);
    for (i=0; i<n; i++) {
	element_mul_mpz(pf->coeff->item[i], pg->coeff->item[i], z);
    }
}

static void poly_mul_si(element_ptr f, element_ptr g, signed long int z)
{
    poly_element_ptr pf = f->data;
    poly_element_ptr pg = g->data;
    int i, n;

    n = pg->coeff->count;
    poly_alloc(f, n);
    for (i=0; i<n; i++) {
	element_mul_si(pf->coeff->item[i], pg->coeff->item[i], z);
    }
}

static void poly_mul(element_ptr r, element_ptr f, element_ptr g)
{
    poly_element_ptr pprod;
    poly_element_ptr pf = f->data;
    poly_element_ptr pg = g->data;
    poly_field_data_ptr pdp = r->field->data;
    int fcount = pf->coeff->count;
    int gcount = pg->coeff->count;
    int i, j, n;
    element_t prod;
    element_t e0;

    if (!fcount || !gcount) {
	element_set0(r);
	return;
    }
    element_init(prod, r->field);
    pprod = prod->data;
    n = fcount + gcount - 1;
    poly_alloc(prod, n);
    element_init(e0, pdp->field);
    for (i=0; i<n; i++) {
	element_ptr x = pprod->coeff->item[i];
	element_set0(x);
	for (j=0; j<=i; j++) {
	    if (j < fcount && i - j < gcount) {
		element_mul(e0, pf->coeff->item[j], pg->coeff->item[i - j]);
		element_add(x, x, e0);
	    }
	}
    }
    poly_remove_leading_zeroes(prod);
    element_set(r, prod);
    element_clear(e0);
    element_clear(prod);
}

void poly_random_monic(element_ptr f, int deg)
{
    int i;
    poly_alloc(f, deg + 1);
    for (i=0; i<deg; i++) {
	element_random(poly_coeff(f, i));
    }
    element_set1(poly_coeff(f, i));
}

element_ptr polymod_coeff(element_ptr e, int i)
{
    element_t *coeff = e->data;

    return coeff[i];
}

inline int polymod_field_degree(field_t f)
{
    polymod_field_data_ptr p = f->data;
    return p->n;
}

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

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

static void polymod_from_hash(element_ptr e, void *data, int len)
{
    //TODO: improve this
    element_t *coeff = e->data;
    int i, n = polymod_field_degree(e->field);
    for (i=0; i<n; i++) {
	element_from_hash(coeff[i], data, len);
    }
}

void polymod_const_mul(element_ptr res, element_ptr a, element_ptr e)
    //a lies in R, e in R[x]
{
    element_t *coeff = e->data, *dst = res->data;
    int i, n = polymod_field_degree(e->field);

    for (i=0; i<n; i++) {
	element_mul(dst[i], coeff[i], a);
    }
}

static void poly_set0(element_ptr e)
{
    poly_alloc(e, 0);
}

static void poly_set1(element_ptr e)
{
    poly_element_ptr p = e->data;
    element_ptr e0;

    poly_alloc(e, 1);
    e0 = p->coeff->item[0];
    element_set1(e0);
}

static int poly_is0(element_ptr e)
{
    poly_element_ptr p = e->data;
    return !p->coeff->count;
}

static int poly_is1(element_ptr e)
{
    poly_element_ptr p = e->data;
    if (p->coeff->count == 1) {
	return element_is1(p->coeff->item[0]);
    }
    return 0;
}

static size_t poly_out_str(FILE *stream, int base, element_ptr e)
{
    int i;
    int n = poly_coeff_count(e);
    size_t result = 2, status;

    /*
    if (!n) {
	if (EOF == fputs("[0]", stream)) return 0;
	return 3;
    }
    */
    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, poly_coeff(e, i));
	if (!status) return 0;
	result += status;
    }
    if (EOF == fputc(']', stream)) return 0;
    return result;
}

static int poly_snprint(char *s, size_t size, element_ptr e)
{
    int i;
    int n = poly_coeff_count(e);
    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, poly_coeff(e, i));
	if (status < 0) return status;
	clip_sub();
    }
    status = snprintf(s + result, left, "]");
    if (status < 0) return status;
    return result + status;
}

void poly_div(element_ptr quot, element_ptr rem,
	element_ptr a, element_ptr b)
{
    poly_element_ptr pq, pr;
    poly_field_data_ptr pdp = a->field->data;
    element_t q, r;
    element_t binv, e0;
    element_ptr qe;
    int m, n;
    int i, k;

    if (element_is0(b)) {
	fprintf(stderr, "BUG! division by zero!\n");
	exit(1);
    }
    n = poly_degree(b);
    m = poly_degree(a);
    if (n > m) {
	element_set(rem, a);
	element_set0(quot);
	return;
    }
    element_init(r, a->field);
    element_init(q, a->field);
    element_init(binv, pdp->field);
    element_init(e0, pdp->field);
    pq = q->data;
    pr = r->data;
    element_set(r, a);
    k = m - n;
    poly_alloc(q, k + 1);
    element_invert(binv, poly_coeff(b, n));
    while (k >= 0) {
	qe = pq->coeff->item[k];
	element_mul(qe, binv, pr->coeff->item[m]);
	for (i=0; i<=n; i++) {
	    element_mul(e0, qe, poly_coeff(b, i));
	    element_sub(pr->coeff->item[i + k], pr->coeff->item[i + k], e0);
	}
	k--;
	m--;
    }
    poly_remove_leading_zeroes(r);
    element_set(quot, q);
    element_set(rem, r);

    element_clear(q);
    element_clear(r);
    element_clear(e0);
    element_clear(binv);
}

void poly_const_mul(element_ptr res, element_ptr a, element_ptr poly)
{
    int i, n = poly_coeff_count(poly);
    poly_alloc(res, n);
    for (i=0; i<n; i++) {
	element_mul(poly_coeff(res, i), a, poly_coeff(poly, i));
    }
    poly_remove_leading_zeroes(res);
}

void poly_invert(element_ptr res, element_ptr f, element_ptr m)
{
    element_t q, r0, r1, r2;
    element_t b0, b1, b2;
    element_t inv;

    element_init(b0, res->field);
    element_init(b1, res->field);
    element_init(b2, res->field);
    element_init(q, res->field);
    element_init(r0, res->field);
    element_init(r1, res->field);
    element_init(r2, res->field);
    element_init(inv, poly_base_field(res));
    element_set0(b0);
    element_set1(b1);
    element_set(r0, m);
    element_set(r1, f);

    for (;;) {
	poly_div(q, r2, r0, r1);
	if (element_is0(r2)) break;
	element_mul(b2, b1, q);
	element_sub(b2, b0, b2);
	element_set(b0, b1);
	element_set(b1, b2);
	element_set(r0, r1);
	element_set(r1, r2);
    }
    element_invert(inv, poly_coeff(r1, 0));
    poly_const_mul(res, inv, b1);
    element_clear(inv);
    element_clear(q);
    element_clear(r0);
    element_clear(r1);
    element_clear(r2);
    element_clear(b0);
    element_clear(b1);
    element_clear(b2);
}

void element_polymod_to_poly(element_ptr f, element_ptr e)
{
    polymod_field_data_ptr p = e->field->data;
    element_t *coeff = e->data;
    int i, n = p->n;
    poly_alloc(f, n);
    for (i=0; i<n; i++) {
	element_set(poly_coeff(f, i), coeff[i]);
    }
    poly_remove_leading_zeroes(f);
}

void element_poly_to_polymod_truncate(element_ptr e, element_ptr f)
{
    polymod_field_data_ptr p = e->field->data;
    element_t *coeff = e->data;
    int i;
    int n;
    n = poly_coeff_count(f);
    if (n > p->n) n = p->n;

    for (i=0; i<n; i++) {
	element_set(coeff[i], poly_coeff(f, i));
    }
    for (; i<p->n; i++) {
	element_set0(coeff[i]);
    }
}

static void polymod_invert(element_ptr r, element_ptr e)
{
    polymod_field_data_ptr p = r->field->data;
    element_ptr minpoly = p->poly;
    element_t f, r1;

    element_init(f, minpoly->field);
    element_init(r1, minpoly->field);
    element_polymod_to_poly(f, e);

    poly_invert(r1, f, p->poly);

    element_poly_to_polymod_truncate(r, r1);

    element_clear(f);
    element_clear(r1);
}

void element_field_to_poly(element_ptr f, element_ptr g)
{
    poly_alloc(f, 1);
    element_set(poly_coeff(f, 0), g);
    poly_remove_leading_zeroes(f);
}

void element_field_to_polymod(element_ptr f, element_ptr g)
{
    polymod_field_data_ptr p = f->field->data;
    element_t *coeff = f->data;
    int i, n = p->n;
    element_set(coeff[0], g);
    for (i=1; i<n; i++) {