; 32-bit x86 assembly
        BITS 32


        ; function entry
        global gcd
gcd:
        ; load arguments to eax (a) and edx (b).
        mov  eax, [esp + 4]
        mov  edx, [esp + 8]
        ; cdecl, must save ebx
        push ebx
        xor  ecx, ecx
        xor  ebx, ebx
        inc  ebx        ; 1

        ; check that eax and edx are non-zero
        test eax, eax
        jz   bail_out

        test edx, edx
        jz   bail_out

        ; main loop
do_loop:
        test edx, edx
        jz   end

        test edx, ebx
        jnz  edx_odd

        ; edx (b) ==  even
        shr  edx, 1

        test eax, ebx

        ; eax odd - shr eax (gcd(a, b) = gcd(a, b/2)
        jnz  do_loop

        ; eax is even as well - gcd(a, b) = 2gcd(a/2, b/2)
        shr  eax, 1
        inc  ecx        ; shift ++
        jmp  do_loop

        ; edx odd...
edx_odd:
        test eax, ebx
        jnz  eax_odd

        ; eax even - shr eax (gcd(a, b) = gcd(a/2, b)
        shr  eax, 1
        jmp  do_loop

eax_odd:
        ;  both are odd... swap a, b (eax, edx) if b (edx) < a (eax)!
        cmp  edx, eax
        jae  do_sub
        xchg edx, eax

do_sub:
        ; gcd(a, b) = gcd(a, (b - a) / 2)
        sub  edx, eax
        shr  edx, 1
        jmp do_loop

bail_out:
        ; returns 0
        xor  eax, eax

end:
        ; edx = 0, eax contains gcd / 2^(ecx)
        ; now we have to shift left by ecx bits...
        shl  eax, cl

        pop  ebx
        ; result in eax, return
        ret