/*
Compilation:

 gcc -fPIC -O3 -Wall -shared mulp.c -o mulp.so

*/

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

void
gmp_isqrt(mp_limb_t *rp, mp_size_t rn, /* const */ mp_limb_t *sp, mp_size_t sn)
{
    if (sp[sn - 1] == 0) {
        sn--;
        rp[rn - 1] = 0;
    }
    mpn_sqrtrem(rp, 0, sp, sn);
}

int
count_zeros(mp_limb_t x)
{
    int res = 0;
    if (!(x & ((1l<<32) - 1))) {
        res = 32;
        x >>= 32;
    }
    if (!(x & ((1l<<16) - 1))) {
        res += 16;
        x >>= 16;
    }
    if (!(x & ((1l<<8) - 1))) {
        res += 8;
        x >>= 8;
    }
    if (!(x & ((1l<<4) - 1))) {
        res += 4 ;
        x >>= 4;
    }
    if (!(x & ((1l<<2) - 1))) {
        res += 2;
        x >>= 2;
    }
    if (!(x & 1)) {
        res += 1;
    }
    return res;
}

mp_size_t
gmp_gcd(mp_limb_t * rp, mp_limb_t *s1p,
       mp_size_t s1n, mp_limb_t *s2p, mp_size_t s2n)
{
    mp_limb_t rc;
    mp_size_t res;
    mp_size_t k1;
    mp_size_t k;
    mp_size_t z1l;
    mp_size_t z2l;
    mp_size_t zp;
    int z1;
    int z2;
    k1 = 1;
    while(k1 <= s1n && s1p[s1n - k1] == 0) k1++;
    s1n -= (k1 - 1);
    k1 = 0;
    while(k1 < s1n && s1p[k1] == 0) k1++;
    s1p += k1;
    s1n -= k1;
    z1l = k1*64;
    if (s1n == 0) {
        *rp = 0;
        return 0;
    }
    k1 = 1;
    while(k1 <= s2n && s2p[s2n - k1] == 0) k1++;
    s2n -= (k1 - 1);
    k1 = 0;
    while(k1 < s2n && s2p[k1] == 0) k1++;
    s2p += k1;
    s2n -= k1;
    z2l = k1*64;
    if (s2n == 0) {
        *rp = 0;
        return 0;
    }
    z2 = count_zeros(*s2p);
    mpn_rshift(s2p, s2p, s2n, z2);
    if (s2p[s2n - 1] == 0) {
        s2n--;
    }
    z1 = count_zeros(*s1p);
    mpn_rshift(s1p, s1p, s1n, z1);
    if (s1p[s1n - 1] == 0) {
        s1n--;
    }
    z1l += z1;
    z2l += z2;
    if (s2n > s1n || ((s2n == s1n) && s2p[s2n - 1] > s1p[s1n -1])) {
        mp_size_t tmp = s2n;
        mp_limb_t * ptmp = s2p;
        s2n = s1n;
        s2p = s1p;
        s1n = tmp;
        s1p = ptmp;
    }
    zp = (z1l < z2l)? z1l : z2l;
    k = zp / 64;
    zp = zp % 64;
    for(k1 = 0; k1 < k; k1++) {
        rp[k1] = 0;
    }
    rp += k;
    res = mpn_gcd(rp, s1p, s1n, s2p, s2n);
    rc = mpn_lshift(rp, rp, res, zp);
    if (rc) {
        rp[res] = rc;
        res++;
    }
    if (rp[res - 1] & (1ul << 63)) {
        rp[res] = 0;
        res++;
    }
    return res+k;
}

#if 0
int
main(void)
{
    int i;
    mp_limb_t num1[3] = {0};
    mp_limb_t num2[3] = {0};
    mp_limb_t res[3] = {0};
    num2[1] = 15;
    num1[0] = 20;
    gmp_gcd(res, num1, 3, num2, 3);
    printf("%lu\n", res[0]);
    return 0;
}
#endif