fieldquadratic.c

这是一个C的源代码

    fq_data_ptr p = e->data;
    element_t e0, e1;
    element_ptr nqr = fq_nqr(e->field);
    int result;
    element_init(e0, p->x->field);
    element_init(e1, e0->field);
    element_square(e0, p->x);
    element_square(e1, p->y);
    element_mul(e1, e1, nqr);
    element_sub(e0, e0, e1);
    result = element_is_sqr(e0);
    element_clear(e0);
    element_clear(e1);
    return result;
}

static void fq_sqrt(element_ptr n, element_ptr e)
{
    fq_data_ptr p = e->data;
    fq_data_ptr r = n->data;
    element_ptr nqr = fq_nqr(n->field);
    element_t e0, e1, e2;

    //if (a+b sqrt(nqr))^2 = x+y sqrt(nqr) then
    //2a^2 = x +- sqrt(x^2 - nqr y^2)
    //(take the sign which allows a to exist)
    //and 2ab = y
    element_init(e0, p->x->field);
    element_init(e1, e0->field);
    element_init(e2, e0->field);
    element_square(e0, p->x);
    element_square(e1, p->y);
    element_mul(e1, e1, nqr);
    element_sub(e0, e0, e1);
    element_sqrt(e0, e0);
    //e0 = sqrt(x^2 - nqr y^2)
    element_add(e1, p->x, e0);
    element_set_si(e2, 2);
    element_invert(e2, e2);
    element_mul(e1, e1, e2);
    //e1 = (x + sqrt(x^2 - nqr y^2))/2
    if (!element_is_sqr(e1)) {
	element_sub(e1, e1, e0);
	//e1 should be a square
    }
    element_sqrt(e0, e1);
    element_add(e1, e0, e0);
    element_invert(e1, e1);
    element_mul(r->y, p->y, e1);
    element_set(r->x, e0);
    element_clear(e0);
    element_clear(e1);
    element_clear(e2);
}

static void field_clear_fq(field_ptr f)
{
    UNUSED_VAR(f);
    //f->order gets cleared automatically
}

static void fq_out_info(FILE *out, field_ptr f)
{
    field_ptr fbase = f->data;
    element_fprintf(out, "x^2 + %B quadratic extension, base field:\n", fq_nqr(f));
    field_out_info(out, fbase);
}

void field_init_quadratic(field_ptr f, field_ptr fbase)
{
    field_init(f);

    f->field_clear = field_clear_fq;
    f->data = fbase;

    f->init = fq_init;
    f->clear = fq_clear;
    f->set_si = fq_set_si;
    f->set_mpz = fq_set_mpz;
    f->to_mpz = fq_to_mpz;
    f->out_str = fq_out_str;
    f->snprint = fq_snprint;
    f->set_str = fq_set_str;
    f->sign = fq_sign;
    f->add = fq_add;
    f->sub = fq_sub;
    f->set = fq_set;
    f->mul = fq_mul;
    f->mul_mpz = fq_mul_mpz;
    f->mul_si = fq_mul_si;
    f->square = fq_square;
    f->doub = fq_double;
    f->neg = fq_neg;
    f->cmp = fq_cmp;
    f->invert = fq_invert;
    f->random = fq_random;
    f->from_hash = fq_from_hash;
    f->is1 = fq_is1;
    f->is0 = fq_is0;
    f->set0 = fq_set0;
    f->set1 = fq_set1;
    f->is_sqr = fq_is_sqr;
    f->sqrt = fq_sqrt;
    f->to_bytes = fq_to_bytes;
    f->from_bytes = fq_from_bytes;
    f->out_info = fq_out_info;

    mpz_mul(f->order, fbase->order, fbase->order);
    if (fbase->fixed_length_in_bytes < 0) {
	f->length_in_bytes = fq_length_in_bytes;
	f->fixed_length_in_bytes = -1;
    } else {
	f->fixed_length_in_bytes = 2 * fbase->fixed_length_in_bytes;
    }
}

void element_field_to_quadratic(element_ptr r, element_ptr a)
{
    fq_data_ptr p = r->data;
    element_set(p->x, a);
    element_set0(p->y);
}

static void fi_mul(element_ptr n, element_ptr a, element_ptr b)
{
    fq_data_ptr p = a->data;
    fq_data_ptr q = b->data;
    fq_data_ptr r = n->data;
    element_t e0, e1, e2;

    element_init(e0, p->x->field);
    element_init(e1, e0->field);
    element_init(e2, e0->field);
    /* Naive way
    element_mul(e0, p->x, q->x);
    element_mul(e1, p->y, q->y);
    element_sub(e0, e0, e1);
    element_mul(e1, p->x, q->y);
    element_mul(e2, p->y, q->x);
    element_add(e1, e1, e2);
    element_set(r->x, e0);
    element_set(r->y, e1);
    */
    //Karatsuba:
    element_add(e0, p->x, p->y);
    element_add(e1, q->x, q->y);
    element_mul(e2, e0, e1);
    element_mul(e0, p->x, q->x);
    element_sub(e2, e2, e0);
    element_mul(e1, p->y, q->y);
    element_sub(r->x, e0, e1);
    element_sub(r->y, e2, e1);

    element_clear(e0);
    element_clear(e1);
    element_clear(e2);
}

static void fi_square(element_ptr n, element_ptr a)
{
    fq_data_ptr p = a->data;
    fq_data_ptr r = n->data;
    element_t e0, e1;

    element_init(e0, p->x->field);
    element_init(e1, e0->field);
    //Re(n) = x^2 - y^2 = (x+y)(x-y)
    element_add(e0, p->x, p->y);
    element_sub(e1, p->x, p->y);
    element_mul(e0, e0, e1);
    //Im(n) = 2xy
    element_mul(e1, p->x, p->y);
    element_add(e1, e1, e1);
    element_set(r->x, e0);
    element_set(r->y, e1);
    element_clear(e0);
    element_clear(e1);
}

static void fi_invert(element_ptr n, element_ptr a)
{
    fq_data_ptr p = a->data;
    fq_data_ptr r = n->data;
    element_t e0, e1;

    element_init(e0, p->x->field);
    element_init(e1, e0->field);
    element_square(e0, p->x);
    element_square(e1, p->y);
    element_add(e0, e0, e1);
    element_invert(e0, e0);
    element_mul(r->x, p->x, e0);
    element_neg(e0, e0);
    element_mul(r->y, p->y, e0);

    element_clear(e0);
    element_clear(e1);
}

static int fi_is_sqr(element_ptr e)
{
    //x + yi is a square <=> x^2 + y^2 is (in the base field)

    // Proof: (=>) if x+yi = (a+bi)^2,
    // then a^2 - b^2 = x, 2ab = y,
    // thus (a^2 + b^2)^2 = (a^2 - b^2)^2 + (2ab)^2 =  x^2 + y^2
    // (<=) Suppose A^2 = x^2 + y^2
    // then if there exist a, b satisfying:
    //   a^2 = (+-A + x)/2, b^2 = (+-A - x)/2
    // then (a + bi)^2 = x + yi.
    // We show that exactly one of (A + x)/2, (-A + x)/2
    // is a quadratic residue (thus a, b do exist).
    // Suppose not. Then the product
    // (x^2 - A^2) / 4 is some quadratic residue, a contradiction
    // since this would imply x^2 - A^2 = -y^2 is also a quadratic residue,
    // but we know -1 is not a quadratic residue.
    fq_data_ptr p = e->data;
    element_t e0, e1;
    int result;
    element_init(e0, p->x->field);
    element_init(e1, e0->field);
    element_square(e0, p->x);
    element_square(e1, p->y);
    element_add(e0, e0, e1);
    result = element_is_sqr(e0);
    element_clear(e0);
    element_clear(e1);
    return result;
}

static void fi_sqrt(element_ptr n, element_ptr e)
{
    fq_data_ptr p = e->data;
    fq_data_ptr r = n->data;
    element_t e0, e1, e2;

    //if (a+bi)^2 = x+yi then
    //2a^2 = x +- sqrt(x^2 + y^2)
    //(take the sign such that a exists) and 2ab = y
    //[thus 2b^2 = - (x -+ sqrt(x^2 + y^2))]
    element_init(e0, p->x->field);
    element_init(e1, e0->field);
    element_init(e2, e0->field);
    element_square(e0, p->x);
    element_square(e1, p->y);
    element_add(e0, e0, e1);
    element_sqrt(e0, e0);
    //e0 = sqrt(x^2 + y^2)
    element_add(e1, p->x, e0);
    element_set_si(e2, 2);
    element_invert(e2, e2);
    element_mul(e1, e1, e2);
    //e1 = (x + sqrt(x^2 + y^2))/2
    if (!element_is_sqr(e1)) {
	element_sub(e1, e1, e0);
	//e1 should be a square
    }
    element_sqrt(e0, e1);
    element_add(e1, e0, e0);
    element_invert(e1, e1);
    element_mul(r->y, p->y, e1);
    element_set(r->x, e0);
    element_clear(e0);
    element_clear(e1);
    element_clear(e2);
}

void element_field_to_fi(element_ptr a, element_ptr b)
{
    fq_data_ptr p = a->data;

    element_set(p->x, b);
    element_set0(p->y);
}

static void fi_out_info(FILE *out, field_ptr f)
{
    field_ptr fbase = f->data;
    fprintf(out, "x^2 + 1 quadratic extension, base field:\n");
    field_out_info(out, fbase);
}

static void field_clear_fi(field_ptr f)
{
    UNUSED_VAR(f);
    //f->order gets cleared automatically
}

void field_init_fi(field_ptr f, field_ptr fbase)
{
    field_init(f);
    f->field_clear = field_clear_fi;
    f->data = fbase;
    f->init = fq_init;
    f->clear = fq_clear;
    f->set_si = fq_set_si;
    f->set_mpz = fq_set_mpz;
    f->to_mpz = fq_to_mpz;
    f->out_str = fq_out_str;
    f->snprint = fq_snprint;
    f->set_str = fq_set_str;
    f->sign = fq_sign;
    f->add = fq_add;
    f->sub = fq_sub;
    f->set = fq_set;
    f->mul = fi_mul;
    f->mul_mpz = fq_mul_mpz;
    f->mul_si = fq_mul_si;
    f->square = fi_square;
    f->doub = fq_double;
    f->neg = fq_neg;
    f->cmp = fq_cmp;
    f->invert = fi_invert;
    f->random = fq_random;
    f->from_hash = fq_from_hash;
    f->is1 = fq_is1;
    f->is0 = fq_is0;
    f->set0 = fq_set0;
    f->set1 = fq_set1;
    f->is_sqr = fi_is_sqr;
    f->sqrt = fi_sqrt;
    f->to_bytes = fq_to_bytes;
    f->from_bytes = fq_from_bytes;
    f->out_info = fi_out_info;

    mpz_mul(f->order, fbase->order, fbase->order);
    if (fbase->fixed_length_in_bytes < 0) {
	f->length_in_bytes = fq_length_in_bytes;
	f->fixed_length_in_bytes = -1;
    } else {
	f->fixed_length_in_bytes = 2 * fbase->fixed_length_in_bytes;
    }
}