Function Divide

src/cryptlib/integer.cpp:922–974

// B[NB] -------- divisor

void Divide(word *R, word *Q, word *T, const word *A, unsigned int NA, const word *B, unsigned int NB)
{
    assert(NA && NB && NA%2==0 && NB%2==0);
    assert(B[NB-1] || B[NB-2]);
    assert(NB <= NA);

    // set up temporary work space
    word *const TA=T;
    word *const TB=T+NA+2;
    word *const TP=T+NA+2+NB;

    // copy B into TB and normalize it so that TB has highest bit set to 1
    unsigned shiftWords = (B[NB-1]==0);
    TB[0] = TB[NB-1] = 0;
    CopyWords(TB+shiftWords, B, NB-shiftWords);
    unsigned shiftBits = WORD_BITS - BitPrecision(TB[NB-1]);
    assert(shiftBits < WORD_BITS);
    ShiftWordsLeftByBits(TB, NB, shiftBits);

    // copy A into TA and normalize it
    TA[0] = TA[NA] = TA[NA+1] = 0;
    CopyWords(TA+shiftWords, A, NA);
    ShiftWordsLeftByBits(TA, NA+2, shiftBits);

    if (TA[NA+1]==0 && TA[NA] <= 1)
    {
        Q[NA-NB+1] = Q[NA-NB] = 0;
        while (TA[NA] || Compare(TA+NA-NB, TB, NB) >= 0)
        {
            TA[NA] -= Subtract(TA+NA-NB, TA+NA-NB, TB, NB);
            ++Q[NA-NB];
        }
    }
    else
    {
        NA+=2;
        assert(Compare(TA+NA-NB, TB, NB) < 0);
    }

    word B0 = TB[NB-2] + 1;
    word B1 = TB[NB-1] + (B0==0);

    // start reducing TA mod TB, 2 words at a time
    for (unsigned i=NA-2; i>=NB; i-=2)
    {
        AtomicDivide(Q[i-NB], Q[i-NB+1], TA+i-2, B0, B1);
        CorrectQuotientEstimate(TA+i-NB, TP, Q[i-NB], Q[i-NB+1], TB, NB);
    }

    // copy TA into R, and denormalize it
    CopyWords(R, TA+shiftWords, NB);
    ShiftWordsRightByBits(R, NB, shiftBits);
}

static inline unsigned int EvenWordCount(const word *X, unsigned int N)
{