binary_split.c

#include <stdio.h>
#include <gmp.h>

void binary_split(mpz_t a, mpz_t b, mpz_t Pab, mpz_t Qab, mpz_t Rab) {
    mpz_t c;
    mpz_init(c);

    mpz_add_ui(c, a, 1);
    if (mpz_cmp(b, c) == 0) {
        mpz_mul_ui(Pab, a, 6);
        mpz_sub_ui(Pab, Pab, 5);
        mpz_neg(Pab, Pab);

        mpz_mul_ui(c, a, 2);
        mpz_sub_ui(c, c, 1);
        mpz_mul(Pab, Pab, c);

        mpz_mul_ui(c, a, 6);
        mpz_sub_ui(c, c, 1);
        mpz_mul(Pab, Pab, c);

        mpz_pow_ui(Qab, a, 3);
        mpz_mul_ui(Qab, Qab, 10939058860032000);

        mpz_mul_ui(Rab, a, 545140134);
        mpz_add_ui(Rab, Rab, 13591409);
        mpz_mul(Rab, Rab, Pab);
    } else {
        mpz_t m;
        mpz_init(m);

        mpz_add(m, a, b);
        mpz_tdiv_q_2exp(m, m, 1);

        mpz_t Pam, Qam, Ram, Pmb, Qmb, Rmb;
        mpz_init(Pam);
        mpz_init(Qam);
        mpz_init(Ram);
        mpz_init(Pmb);
        mpz_init(Qmb);
        mpz_init(Rmb);

        binary_split(a, m, Pam, Qam, Ram);
        binary_split(m, b, Pmb, Qmb, Rmb);

        mpz_mul(Pab, Pam, Pmb);
        mpz_mul(Qab, Qam, Qmb);
        mpz_mul(Rab, Qmb, Ram);
        mpz_mul(c, Pam, Rmb);
        mpz_add(Rab, Rab, c);
    }
}

int main(int argc, char *argv[]) {
    mpz_t a, b, P1n, Q1n, R1n;
    mpz_init(a);
    mpz_init(b);
    mpz_init(P1n);
    mpz_init(Q1n);
    mpz_init(R1n);

    mpz_set_ui(a, 1);
    mpz_set_ui(b, 100);

    binary_split(a, b, P1n, Q1n, R1n);
    gmp_printf("P1n = %Zd\nQ1n = %Zd\nR1n = %Zd\n", P1n, Q1n, R1n);
    return 0;
}