efficient way to divide a very large number stored in 2 registers by a constant

Question

Let's say I want to calculate the following:

A/Z

Where A is of length 128 bit and Z is 64 bit long. A is stored in 2 64 bit registers since the registers of the system can store up to 64 bits. What would be an efficient way to calculate the result?

P.S: I've solved similar multiplication problems by using CSD representations. However, this would require calculating 1/Z first.

Answer 1

[Edit1] spotted bug repaired

I assume you want integer division so here is the math for 8bit analogy:

A = { a0 + (a1<<8) }
D = { d0 + (d1<<8) } ... division result
Z = { z0 }
D = (a0/z0) + ((a1*256)/z0) +                      (( (a0%z0) + ((a1*256)%z0) )/z0);
D = (a0/z0) + ((a1/z0)*256) + ((a1%z0)*(256/z0)) + (( (a0%z0) + ((a1%z0)*(256%z0)) )/z0);

Now the terms 256/z0 and 256%z0 can be computed like this (C++):

    i0=0xFF/z0; if ((z0&(z0-1))==0) i0++;   // i0 = 256/z0
    i1=i0*z0; i1^=0xFF; i1++;               // i1 = 256%z0

So the i0 is just incremented in case the z0 is power of 2, and i1 is just remainder computed from the division.

    a/b = d + r/b
    r = a - a*d

Here tested 8bit code:

//---------------------------------------------------------------------------
// unsigned 8 bit ALU in C++
//---------------------------------------------------------------------------
BYTE cy;                                                    // carry flag cy = { 0,1 }
void inc(BYTE &a);                                          // a++
void dec(BYTE &a);                                          // a--
void add(BYTE &c,BYTE a,BYTE b);                            // c = a+b
void adc(BYTE &c,BYTE a,BYTE b);                            // c = a+b+cy
void sub(BYTE &c,BYTE a,BYTE b);                            // c = a-b
void sbc(BYTE &c,BYTE a,BYTE b);                            // c = a-b-cy
void mul(BYTE &h,BYTE &l,BYTE a,BYTE b);                    // (h,l) = a/b
void div(BYTE &h,BYTE &l,BYTE &r,BYTE ah,BYTE al,BYTE b);   // (h,l) = (ah,al)/b ; r = (ah,al)%b
//---------------------------------------------------------------------------
void inc(BYTE &a) { if (a==0xFF) cy=1; else cy=0; a++; }
void dec(BYTE &a) { if (a==0x00) cy=1; else cy=0; a--; }
void add(BYTE &c,BYTE a,BYTE b)
    {
    c=a+b;
    cy=BYTE(((a &1)+(b &1)   )>>1);
    cy=BYTE(((a>>1)+(b>>1)+cy)>>7);
    }
void adc(BYTE &c,BYTE a,BYTE b)
    {
    c=a+b+cy;
    cy=BYTE(((a &1)+(b &1)+cy)>>1);
    cy=BYTE(((a>>1)+(b>>1)+cy)>>7);
    }
void sub(BYTE &c,BYTE a,BYTE b)
    {
    c=a-b;
    if (a<b) cy=1; else cy=0;
    }
void sbc(BYTE &c,BYTE a,BYTE b)
    {
    c=a-b-cy;
    if (cy) { if (a<=b) cy=1; else cy=0; }
    else    { if (a< b) cy=1; else cy=0; }
    }
void mul(BYTE &h,BYTE &l,BYTE a,BYTE b)
    {
    BYTE ah,al;
    h=0; l=0; ah=0; al=a;
    if ((a==0)||(b==0)) return;
    // long binary multiplication
    for (;b;b>>=1)
        {
        if (BYTE(b&1))
            {
            add(l,l,al);    // (h,l)+=(ah,al)
            adc(h,h,ah);
            }
        add(al,al,al);      // (ah,al)<<=1
        adc(ah,ah,ah);
        }
    }
void div(BYTE &d1,BYTE &d0,BYTE &r,BYTE a1,BYTE a0,BYTE z0)
    {
    //      D = (a0/z0) + ((a1*256)/z0) +                      (( (a0%z0) + ((a1*256)%z0) )/z0);
    //      D = (a0/z0) + ((a1/z0)*256) + ((a1%z0)*(256/z0)) + (( (a0%z0) + ((a1%z0)*(256%z0)) )/z0);
    // edge cases
    if (z0==0){ d0= 0; d1= 0; r=0; }
    if (z0==1){ d0=a0; d1=a1; r=0; }
    // normal division
    if (z0>=2)
        {
        BYTE i0,i1,e0,e1,f0,f1,t,dt;
        i0=0xFF/z0; if ((z0&(z0-1))==0) i0++;   // i0 = 256/z0
        i1=i0*z0; i1^=0xFF; i1++;               // i1 = 256%z0
        t=a1%z0;
        mul(e1,e0,t,i0);                        // e = (a1%z0)*(256/z0)
        mul(f1,f0,t,i1);                        // f = (a1%z0)*(256%z0)
        add(f0,f0,a0%z0);                       // f = (a0%z0) + (a1%z0)*(256%z0)
        adc(f1,f1,0);
        add(d0,a0/z0,e0);
        adc(d1,a1/z0,e1);
        // t = division of problematic term by z0
        t=0;
        for (;f1;)
            {
            dt=f1*i0;
            mul(e1,e0,dt,z0);
            sub(f0,f0,e0);
            sbc(f1,f1,e1);
            t+=dt;
            }
        if (f0>=z0) t+=f0/z0;
        // correct output
        add(d0,d0,t);
        adc(d1,d1,0);
        // remainder
        r=d0*z0;
        r=a0-r;
        }
    }
//---------------------------------------------------------------------------

The 8bit ALU is not optimized at all I just busted it to test it right now as original project is nowhere to found... I assume you are doing it in asm so you use can use CPU/ALU instructions carry instead. The only important function is the div.

Notes:

This is only 8 bit. To convert it to 64 bit just change all 0xFF to 0xFFFFFFFFFFFFFFFF and BYTE to your data type and <<8 to <<64.

Division result is in d0, d1 and remainder is in r Code does not handle negative values.

Sadly the term:

(( (a0%z0) + ((a1%z0)*(256%z0)) )/z0);

in its current state requires also 16 bit division (not full though as result is not arbitrary instead a composite of two mod z0 values). I managed to avoid long division by few (for 16bit:8bit is the worst case 7) iterations. However my guts are telling me it should be computed simpler using some modular math identity I do not know or cant think of right now. This makes this division relatively slow.

Answer 2

The right way to solve such a problem, is by returning to the basics:

• divide the most significant register by the denominator
• calculate the quotient Q and the rest R
• define a new temporary register preferrably with the same length as the other 2
• the rest should occupy the most significant bits in the temporary register
• shift the lesser significant register to the right by the same amount of bits contained iR and add to the result to the temporary register.
• go back to step 1

after the division, the resulting rest must be casted to double, divided by the denominator then added to the quotient.