PHDay CTF mach-o solution

from struct import pack, unpack

# the crackme is a huge mach-o executable with a nice and original obfuscation

# I decided to study it under ollydbg by loading it with a custom loader (there is only few imports and all are standards (except the MD5 call))
# it's ugly but it works well :) (and olly is definitively the best debugger to study obfuscated codes :) )

'''
#include <stdio.h>
#include <stdlib.h>
#include <openssl/md5.h>
#include <windows.h>

int main()
{
    PBYTE maclol = VirtualAlloc(NULL, 0xB8744,  MEM_COMMIT | MEM_RESERVE, PAGE_EXECUTE_READWRITE);
    HANDLE hFile = CreateFileA("a.out", GENERIC_READ, 0, NULL, OPEN_EXISTING, 0, NULL);
    DWORD numberOfBytesRead;
    int (*macmain)(int argc, char** argv);
    char *argv[2] = {"a.out"};

    ReadFile(hFile, maclol, 0xB8744, &numberOfBytesRead, NULL);
    maclol -= 0x1000;
    maclol[0x056C3E] = 0x68;
    *(PDWORD)(&maclol[0x056C3E+1]) = (DWORD)fread;
    maclol[0x056C3E+5] = 0xC3;
    maclol[0x056C44] = 0x68;
    *(PDWORD)(&maclol[0x056C44+1]) = (DWORD)fseek;
    maclol[0x056C44+5] = 0xC3;
    maclol[0x056C4A] = 0x68;
    *(PDWORD)(&maclol[0x056C4A+1]) = (DWORD)ftell;
    maclol[0x056C4A+5] = 0xC3;
    maclol[0x056C56] = 0x68;
    *(PDWORD)(&maclol[0x056C56+1]) = (DWORD)memcmp;
    maclol[0x056C56+5] = 0xC3;
    maclol[0x056C5C] = 0x68;
    *(PDWORD)(&maclol[0x056C5C+1]) = (DWORD)memcpy;
    maclol[0x056C5C+5] = 0xC3;
    maclol[0x056C62] = 0x68;
    *(PDWORD)(&maclol[0x056C62+1]) = (DWORD)memset;
    maclol[0x056C62+5] = 0xC3;
    maclol[0x056C68] = 0x68;
    *(PDWORD)(&maclol[0x056C68+1]) = (DWORD)printf;
    maclol[0x056C68+5] = 0xC3;
    maclol[0x056C38] = 0x68;
    *(PDWORD)(&maclol[0x056C38+1]) = (DWORD)fopen;
    maclol[0x056C38+5] = 0xC3;
    maclol[0x056BEA] = 0x68;
    *(PDWORD)(&maclol[0x056BEA+1]) = (DWORD)MD5;
    maclol[0x056BEA+5] = 0xC3;
    macmain = (int (*)(int argc, char** argv))(&maclol[0x40E40]);
    macmain(2, argv);
    return 0;
}
'''

# the verification involved 4 modular polynomials equations with 4 variables :
# the key consist of 0x2000 4-uples of polynoms (P1, P2, P3, P4)

# C1, C2, C3 and C4 are calculated by using the following equations and must be equals to those stored in the binary :
# C1 = P1*(x^1+x^2+x^3 [0xE]) + P2*(x^0+x^1+x^3 [0xB]) + P3*(x^0+x^2+x^3 [0xD]) + P4*(x^0+x^3     [0x9])
# C2 = P1*(x^0+x^3     [0x9]) + P2*(x^1+x^2+x^3 [0xE]) + P3*(x^0+x^1+x^3 [0xB]) + P4*(x^0+x^2+x^3 [0xD])
# C3 = P1*(x^0+x^2+x^3 [0xD]) + P2*(x^0+x^3     [0x9]) + P3*(x^1+x^2+x^3 [0xE]) + P4*(x^0+x^1+x^3 [0xB])
# C4 = P1*(x^0+x^1+x^3 [0xB]) + P2*(x^0+x^2+x^3 [0xD]) + P3*(x^0+x^3     [0x9]) + P4*(x^1+x^2+x^3 [0xE])

# I coded the basic operation on this field and
# a little Gaussian elimination to solve those equations and so recover the 4uples P1, P2, P3 and P4.

def mul_mod(P1, P2) :
    tmp1 = P1
    tmp2 = P2
    r = 0
    while tmp2 :
        if tmp2 & 1 :
            r ^= tmp1
        tmp2 >>= 1
        if tmp1 & 0x80000000 :
            tmp1 ^= 0x8000001B
        tmp1 <<= 1
    return r

def mul(P1, P2) :
    tmp1 = P1
    tmp2 = P2
    r = 0
    while tmp2 :
        if tmp2 & 1 :
            r ^= tmp1
        tmp1 <<= 1
        tmp2 >>= 1
    return r

def div(P1, P2) :
    assert P1 > P2
    i = 0
    r = 0
    m = P1
    log2P2 = 0
    while P2 >> log2P2 :
        log2P2 += 1
    log2P2 -= 1
    while P1 >= (P2 << i) :
        i += 1
    while i >= 0 :
        if m & (1 << (i + log2P2)) :
            m ^= P2 << i
            r |= 1 << i
        i -= 1
    return r, m

def eea(u, v):
    u1 = 1
    u2 = 0
    u3 = u
    v1 = 0
    v2 = 1
    v3 = v
    while v3 != 0:
        q, t3 = div(u3, v3)
        t1 = u1 ^ mul(q, v1)
        t2 = u2 ^ mul(q, v2)
        u1 = v1
        u2 = v2
        u3 = v3
        v1 = t1
        v2 = t2
        v3 = t3
    return u1, u2, u3

def inv(P) :
    u1, u2, u3 = eea(0x100000036, P)
    assert(u3 == 1)
    return u2

gauss = [
 ([0, 0, 1, 0], [0xD, 0x9, 0xE, 0xB]),
 ([0, 1, 0, 0], [0x9, 0xE, 0xB, 0xD]),
 ([0, 0, 0, 1], [0xB, 0xD, 0x9, 0xE]),
 ([1, 0, 0, 0], [0xE, 0xB, 0xD, 0x9])
]

for i in xrange(4) :
    invx = inv(gauss[i][1][i])

    for j in xrange(4) :
        gauss[i][0][j] = mul_mod(gauss[i][0][j], invx)
        gauss[i][1][j] = mul_mod(gauss[i][1][j], invx)

    for j in xrange(4) :
        if j != i :
            coeff = gauss[j][1][i]
            for k in xrange(4) :
                gauss[j][0][k] ^= mul_mod(coeff, gauss[i][0][k])
                gauss[j][1][k] ^= mul_mod(coeff, gauss[i][1][k])

# table of the reference (C1, C2, C3, C4) 4-uple values checked at the end of the execution
table = file("table", "rb").read()
key = ""
while table :
    out = unpack("<IIII", table[:0x10])
    table = table[0x10:]
    for i in xrange(4) :
        r = 0
        for j in xrange(4) :
            r ^= mul_mod(gauss[i][0][j], out[j])
        key += pack("<I", r)
file("ctf-pass-file.txt", "wb").write(key)

# by the way, there is a bug in the crackme, only the first 0x11 bytes of the 0x20 000 bytes table are checked :)