/ ulldiv.s (emx+gcc) -- Copyright (c) 1992-1993 by Eberhard Mattes

/ There's a better solution: set x':=x, y':=y. Shift right both x' and
/ y' until x' < 2^32. Then divide x' by y'. The result can be off by -1
/ due to the approximation. This can be checked by multiplying the
/ result by y and comparing to x-y. This verification step is not
/ required for x < 2^32.
/
        .globl  __ulldiv

        .text

        .align  2, 0x90

/ _ulldiv_t _ulldiv (unsigned long long x, unsigned long long y)

/define saved_ebp     0(%ebp)
/define ret_addr      4(%ebp)
#define dst           8(%ebp)
#define x_lo         12(%ebp)
#define x_hi         16(%ebp)
#define y_lo         20(%ebp)
#define y_hi         24(%ebp)

#define quot_lo       0(%edi)
#define quot_hi       4(%edi)
#define rem_lo        8(%edi)
#define rem_hi       12(%edi)

#define tmp_lo       %ebx
#define tmp_hi       %esi

__ulldiv:
        pushl   %ebp
        movl    %esp, %ebp
        pushl   %ebx
        pushl   %esi
        pushl   %edi
        movl    dst, %edi
        cmpl    $0, y_hi
        jz      Ldiv32
        movl    x_lo, %eax
        movl    x_hi, %edx
        movl    %eax, tmp_lo
        movl    %edx, tmp_hi
        xorl    %eax, %eax
        xorl    %edx, %edx
        movl    $64, %ecx
        .align  2, 0x90
1:      shll    $1, tmp_lo                      / Shift 128 bits
        rcll    $1, tmp_hi
        rcll    $1, %eax
        rcll    $1, %edx
        subl    y_lo, %eax                      / Subtract divisor
        sbbl    y_hi, %edx
        jb      2f                              / Negative ->
        incl    tmp_lo                          / Set bit 0 of quotient
        loop    1b
        jmp     3f

        .align  2, 0x90
2:      addl    y_lo, %eax                      / Restore divisor
        adcl    y_hi, %edx
        loop    1b
3:      movl    %eax, rem_lo
        movl    %edx, rem_hi
        movl    tmp_lo, quot_lo
        movl    tmp_hi, quot_hi
4:      movl    %edi, %eax
        popl    %edi
        popl    %esi
        popl    %ebx
        leave
        ret

/
/ Divide 64-bit number by 32-bit number
/
/ 2^32 a + b         a    2^32 (a mod c) + b
/ ---------- = 2^32 --- + ------------------   (integer division)
/      c             c           c
/
/
/ (2^32 a + b) mod c = (2^32 (a mod c) + b) mod c
/
/ Proof:

/ Divide by 2^32 c by successive subtraction (k steps) of 2^32 c.
/
/ 2^32 a + b            2^32 a + b - 2^32 k c            2^32 (a-kc) + b
/ ---------- = 2^32 k + --------------------- = 2^32 k + ---------------
/     c                            c                            c
/
/ For k = [a/c], the equality 0 <= a-kc < c is true. Then, a-kc = a mod c.
/ For this k, the final division cannot overflow.
/
/
        .align  2, 0x90
Ldiv32:
        movl    x_hi, %eax
        xorl    %edx, %edx
        divl    y_lo
        movl    %eax, quot_hi           / high-order 32 bits of result
        movl    x_lo, %eax              / remainder of prev. division in %edx
        divl    y_lo
        movl    %eax, quot_lo
        movl    %edx, rem_lo
        movl    $0, rem_hi
        jmp     4b