void Tate2_pairings()

in src/compression/pairing.c [156:287]
109 lines of code
11 McCabe index (conditional complexity)

void Tate2_pairings(const point_t P, const point_t Q, point_full_proj_t *Qj, f2elm_t* f)
{
    felm_t *x, *y, *x_, *y_, *l1;
    f2elm_t finv[2*t_points], one = {0};
    f2elm_t *x_first, *y_first, l1_first, t0, t1, g, h;
    
    fpcopy((digit_t*)&Montgomery_one, one[0]);

    for (int j = 0; j < t_points; j++) {
        fp2copy(one, f[j]);
        fp2copy(one, f[j+t_points]);
    }

    // Pairings with P
    x_first = (f2elm_t*)P->x;
    y_first = (f2elm_t*)P->y;

    x_ = (felm_t*)T_tate2_firststep_P + 0;
    y_ = (felm_t*)T_tate2_firststep_P + 1;
    fpcopy((digit_t*)T_tate2_firststep_P + 2*NWORDS_FIELD, l1_first[0]);         
    fpcopy((digit_t*)T_tate2_firststep_P + 3*NWORDS_FIELD, l1_first[1]);         
    
    for (int j = 0; j < t_points; j++) {
        fp2sub(Qj[j]->X, *x_first, t0);
        fp2sub(Qj[j]->Y, *y_first, t1);
        fp2mul_mont(l1_first, t0, t0);
        fp2sub(t0, t1, g);

        fpsub(Qj[j]->X[0], *x_, h[0]);
        fpcopy(Qj[j]->X[1], h[1]);
        fpneg(h[1]);
        fp2mul_mont(g, h, g);

        fp2sqr_mont(f[j], f[j]);
        fp2mul_mont(f[j], g, f[j]);
    }
    x = x_;
    y = y_;
    
    for (int k = 0; k < OALICE_BITS - 2; k++) {
        x_ = (felm_t*)T_tate2_P + 3 * k + 0;
        y_ = (felm_t*)T_tate2_P + 3 * k + 1;
        l1 = (felm_t*)T_tate2_P + 3 * k + 2;
        for (int j = 0; j < t_points; j++) {
            fpsub(*x, Qj[j]->X[0], t0[1]);
            fpmul_mont(*l1, t0[1], t0[1]);
            fpmul_mont(*l1, Qj[j]->X[1], t0[0]);
            fpsub(Qj[j]->Y[1], *y, t1[1]);
            fpsub(t0[1], t1[1], g[1]);
            fpsub(t0[0], Qj[j]->Y[0], g[0]);

            fpsub(Qj[j]->X[0], *x_, h[0]);
            fpcopy(Qj[j]->X[1], h[1]);
            fpneg(h[1]);
            fp2mul_mont(g, h, g);

            fp2sqr_mont(f[j], f[j]);
            fp2mul_mont(f[j], g, f[j]);
        }
        x = x_;
        y = y_;
    }
    for (int j = 0; j < t_points; j++) {
        fpsub(Qj[j]->X[0], *x, g[0]);
        fpcopy(Qj[j]->X[1], g[1]);
        fp2sqr_mont(f[j], f[j]);
        fp2mul_mont(f[j], g, f[j]);
    }

    // Pairings with Q
    x_first = (f2elm_t*)Q->x;
    y_first = (f2elm_t*)Q->y; 
    x_ = (felm_t*)T_tate2_firststep_Q + 0;
    y_ = (felm_t*)T_tate2_firststep_Q + 1;
    fpcopy(((felm_t*)T_tate2_firststep_Q)[2], l1_first[0]);
    fpcopy(((felm_t*)T_tate2_firststep_Q)[3], l1_first[1]);

    for (int j = 0; j < t_points; j++) {
        fp2sub(Qj[j]->X, *x_first, t0);
        fp2sub(Qj[j]->Y, *y_first, t1);
        fp2mul_mont(l1_first, t0, t0);
        fp2sub(t0, t1, g);

        fpsub(Qj[j]->X[0], *x_, h[0]);
        fpcopy(Qj[j]->X[1], h[1]);
        fpneg(h[1]);
        fp2mul_mont(g, h, g);

        fp2sqr_mont(f[j+t_points], f[j+t_points]);
        fp2mul_mont(f[j+t_points], g, f[j+t_points]);
    }
    x = x_;
    y = y_;
    
    for (int k = 0; k < OALICE_BITS - 2; k++) {
        x_ = (felm_t*)T_tate2_Q + 3*k + 0;
        y_ = (felm_t*)T_tate2_Q + 3*k + 1;
        l1 = (felm_t*)T_tate2_Q + 3*k + 2;
        for (int j = 0; j < t_points; j++) {
            fpsub(Qj[j]->X[0], *x, t0[0]);
            fpmul_mont(*l1, t0[0], t0[0]);
            fpmul_mont(*l1, Qj[j]->X[1], t0[1]);
            fpsub(Qj[j]->Y[0], *y, t1[0]);
            fpsub(t0[0], t1[0], g[0]);
            fpsub(t0[1], Qj[j]->Y[1], g[1]);

            fpsub(Qj[j]->X[0], *x_, h[0]);
            fpcopy(Qj[j]->X[1], h[1]);
            fpneg(h[1]);
            fp2mul_mont(g, h, g);

            fp2sqr_mont(f[j+t_points], f[j+t_points]);
            fp2mul_mont(f[j+t_points], g, f[j+t_points]);
        }
        x = x_;
        y = y_;
    }
    // Last iteration
    for (int j = 0; j < t_points; j++) {
        fpsub(Qj[j]->X[0], *x, g[0]);
        fpcopy(Qj[j]->X[1], g[1]);

        fp2sqr_mont(f[j+t_points], f[j+t_points]);
        fp2mul_mont(f[j+t_points], g, f[j+t_points]);
    }

    // Final exponentiation:
    mont_n_way_inv(f, 2*t_points, finv);
    for (int j = 0; j < 2*t_points; j++) {
        final_exponentiation_2_torsion(f[j], finv[j], f[j]);
    }
}