mpDiv  多倍長整数同士の除算

int mpDiv(int *q, int *r, int *a, int *b);

q     （出力）多倍長整数同士の除算の商
r     （出力）多倍長整数同士の除算の余り
a, b  （入力）多倍長整数（a が被除数、b が除数）

正常に計算できたときは0。失敗したときは-1。除数が0のときは-2。

#define N 10000

int mpDiv(int *q, int *r, int *za, int *zb)
{
    int  i;
    int  m, n;
    int  *aa, *bb, *qq, *rr, *t;
    int  *a, *b;
    int  ca;
    long x;
    int  k, Q;

    *q = 0, *r = 0;
    if (*zb == 0) return -2;
    if (*za == 0) return 0;

    if (*za < *zb) {
        for (aa = za, rr = r, i = *za; i >= 0; i--) *rr++ = *aa++; 
        return 0;
    }
    
    if (*zb == 1) {
        *q = *za;
        zb++;
        for (ca = 0, aa = za + *za, qq = q + *q, i = *za; i > 0; i--) {
            x = (long)N * ca + *aa--;
            *qq-- = x / *zb, ca = x % *zb;
        }
        if (*(q + *q) == 0) (*q)--;
        if (ca > 0) {
            *r++ = 1;
            *r = ca;
        } else *r = 0;
        return 0;
    }

    if ((a = (int *)malloc((sizeof(*za)+1) * sizeof(int))) == NULL) return -1;
    if ((b = (int *)malloc((sizeof(*zb)+2) * sizeof(int))) == NULL) {
        free(a);
        return -1;
    }
    for (aa = a, i = *za; i >= 0; i--) *aa++ = *za++; 
    for (bb = b, i = *zb; i >= 0; i--) *bb++ = *zb++; 

    if ((k = (N/2-1) / *(b + *b) + 1) > 1) {
        for (ca = 0, aa = a, qq = q, i = 0; i < *a; i++) {
            x = (long)k * *++aa + ca;      /* a = a * k */
            *++qq = x % N, ca = x / N;
        }
        *++qq = ca;
        *q = i + (ca > 0);
        for (qq = q, aa = a, i = *q; i >= 0; i--) *aa++ = *qq++;

        for (ca = 0, bb = b, qq = q, i = 0; i < *b; i++) {
            x = (long)k * *++bb + ca;      /* b = b * k */
            *++qq = x % N, ca = x / N;
        }
        *++qq = ca;
        *q = i + (ca > 0);
        for (qq = q, bb = b, i = *q; i >= 0; i--) *bb++ = *qq++;
    }

    *q = *a - *b + 1;
    for (qq = q, i = *q; i > 0; i--) *++qq = 0;
    n = *b;
    while ((m = *a) >= n) {
        if (*(a + *a) >= *(b + *b)) {
            for (aa = a + *a, bb = b + *b; bb != b; aa--, bb--) 
                if (*aa != *bb) break;
            if (bb == b) {
                *a -= *b;
                *(q + m - n + 1) = 1;
                continue;
            } else if (*aa > *bb) {
                for (t = bb, ca = 0, aa = a + m - n, bb = b; bb < t; ) {
                    *++aa -= *++bb + ca;
                    ca = 0;
                    if (*aa < 0) *aa += N, ca = 1;
                }
                while (*aa == 0) aa--;
                *a = aa - a;
                *(q + m - n + 1) = 1;
                continue;
            }
            Q = N - 1;
        } else Q = ((long)N * *(a + *a) + *(a + *a - 1)) / *(b + *b);
        if (m == n) break;
        
        for ( ; ; ) {
            if (Q == 1) {
                *(b + *b + 1) = 0;
                for (ca=0, aa=a+(*a - *b-1), bb=b, i= *b; i >= 0; i--) {
                    *++aa -= *++bb + ca;
                    ca = 0;
                    if (*aa < 0) *aa += N, ca = 1;
                }
                while (*aa == 0) aa--;
                *a = aa - a;
                break;
            }
            for (ca = 0, rr = r, bb = b, i = 0; i < *b; i++) {
                x = (long)Q * *++bb + ca;
                *++rr = x % N, ca = x / N;
            }
            *++rr = ca;
            *r = i + 1;
            for (aa = a + *a, rr = r + *r; rr != r; aa--, rr--)
                if (*aa != *rr) break;
            if (rr == r) {
                *a -= *r;
                break;
            } else if (*aa > *rr) {
                for (t = rr, ca = 0, aa = a + (*a - *r), rr = r; rr < t; ) {
                    *++aa -= *++rr + ca;
                    ca = 0;
                    if (*aa < 0) *aa += N, ca = 1;
                }
                while (*aa == 0) aa--;
                *a = aa - a;
                break;
            } else Q--;
        }
        *(q + m - n) = Q;
    }
    if (*(q + *q) == 0) (*q)--;
    
    if (k > 1) {
        for (ca = 0, aa = a + *a, rr = r + *a, i = *a; i > 0; i--) {
            x = (long)N * ca + *aa--;
            *rr-- = x / k, ca = x % k;
        }
        *r = *a - (*(r + *a) == 0);
    } else for (aa = a, rr = r, i = *a; i >= 0; i--) *rr++ = *aa++; 

    free(a);
    free(b);
    return 0;
}