How are these AES MixColumn multiplication tables calculated? I'm using the mul2, mul3, mul9, mul11, mul13 and mul14 tables for the MixColumn and InvMixColumn steps in AES-128. However, I got these off some Github repository, and now I'm looking for an actual source explaining how these were calculated and so on.

I cannot find anything proper on Google, can anyone help me out?
It looks like those are just precalculated tables for multiplication by the constants 2, 3, 9, 11, 13 and 14 in the AES finite field.  We already have several questions about how multiplication in finite fields in general, and in the AES representation of GF(28) specifically, is done, such as:


Galois fields in cryptography
How is multiplication in field $GF(2^{8})$ done?
How to optimise a finite field multiplication?
How to do Hexadecimal multiplication in GF(2^8)


and even:


Efficient pen-and-paper calculation of the Galois Field multiplication?


That said, I would not recommend using lookup tables for these multiplications at all, since they can be easily calculated at runtime.  Not only do lookup tables introduce the risk of timing attacks and other side channel attacks (since memory lookups don't always run in constant time), but they are also often less efficient than calculating the multiplication directly, both on high-end CPUs (where memory access is relatively slow) and on low-end microcontrollers (where memory itself is often a scarce resource).

The usual way to do multiplication in GF(2n) is using the double-and-add algorithm, also known as "Russian peasant multiplication", which uses only two primitive arithmetic operations:


addition, which in GF(2n) is simply bitwise XOR, and
doubling, i.e. multiplication by 2 (or by $x$, in the polynomial representation), which requires a left shift followed by a polynomial modular reduction, and can be efficiently implemented e.g. like this in C:



static inline uint8_t dbl(uint8_t a) {
    return (a &lt;&lt; 1) ^ (0x11b &amp; -(a &gt;&gt; 7));
}



Both of these operations can be fairly easily implemented without branches or table lookups, which means they can be usually expected to run in constant time (although one must, of course, always watch out for possibly compiler and/or CPU optimizations that could reintroduce data-dependent timing variations).  Using these primitive operations, we can write a general AES finite field multiplication routine e.g. like this:



uint8_t gmul(uint8_t a, uint8_t b) {
    uint8_t p = 0;
    while (b &gt; 0) {
        p ^= a &amp; -(b &amp; 1);         /* add a to p if the lowest bit of b is set */
        a = dbl(a);                /* multiply a by 2 in the finite field */
        b &gt;&gt;= 1;                   /* shift the next bit of b to the lowest position */
    }
    return p;
}


Note that the running time of the code above depends on b, and specifically on the position of its highest set bit, but (hopefully, assuming a reasonable compiler and CPU) not on a.  If we wanted, we could also make the timing fully independent of b by replacing while (b &gt; 0) with a fixed 8 iteration loop such as for (int i = 0; i &lt; 8; i++), but this isn't really needed if b is a public constant.

However, for multiplication by a constant b, we don't really need this general routine at all; knowing the value of b, we can manually unroll the loop and eliminate any steps we know are not needed, leaving just e.g.:



uint8_t mul2(uint8_t a) {
    return dbl(a);
}
uint8_t mul3(uint8_t a) {  /* 3 = 2 + 1 */
    return dbl(a) ^ a;
}
uint8_t mul9(uint8_t a) {  /* 9 = 8 + 1 */
    return dbl(dbl(dbl(a))) ^ a;
}
uint8_t mul11(uint8_t a) {  /* 11 = 8 + 2 + 1 */
    uint8_t a2 = dbl(a), a4 = dbl(a2), a8 = dbl(a4);
    return a8 ^ a2 ^ a;
}
uint8_t mul13(uint8_t a) {  /* 13 = 8 + 4 + 1 */
    uint8_t a2 = dbl(a), a4 = dbl(a2), a8 = dbl(a4);
    return a8 ^ a4 ^ a;
}
uint8_t mul14(uint8_t a) {  /* 14 = 8 + 4 + 2 */
    uint8_t a2 = dbl(a), a4 = dbl(a2), a8 = dbl(a4);
    return a8 ^ a4 ^ a2;
}


Also note that, for AES (Inv)MixColumns, you're going to need to multiply each input byte by four different constants in a specific pattern and sum (i.e. XOR) the results.  It's possible to do this more efficiently by splitting the finite field multiplications into doublings and XORs and applying the distributive law dbl(x) ^ dbl(y) == dbl(x ^ y), e.g. like this:



void mix_columns(uint8_t cols[4]) {
    uint8_t a = cols[0], b = cols[1], c = cols[2], d = cols[3];
    cols[0] = dbl(a ^ b) ^ b ^ c ^ d;  /* 2a + 3b + c + d */
    cols[1] = dbl(b ^ c) ^ c ^ d ^ a;  /* 2b + 3c + d + a */
    cols[2] = dbl(c ^ d) ^ d ^ a ^ b;  /* 2c + 3d + a + b */
    cols[3] = dbl(d ^ a) ^ a ^ b ^ c;  /* 2d + 3a + b + c */
}

void inv_mix_columns(uint8_t cols[4]) {
    uint8_t a = cols[0], b = cols[1], c = cols[2], d = cols[3];
    uint8_t x = dbl(a ^ b ^ c ^ d), y = dbl(x ^ a ^ c), z = dbl(x ^ b ^ d);
    cols[0] = dbl(y ^ a ^ b) ^ b ^ c ^ d;  /* 14a + 11b + 13c + 9d */
    cols[1] = dbl(z ^ b ^ c) ^ c ^ d ^ a;  /* 14b + 11c + 13d + 9a */
    cols[2] = dbl(y ^ c ^ d) ^ d ^ a ^ b;  /* 14c + 11d + 13a + 9b */
    cols[3] = dbl(z ^ d ^ a) ^ a ^ b ^ c;  /* 14d + 11a + 13b + 9c */
}


(Online demo with self-test.)

Further performance improvements may be possible e.g. through the use of bitslicing techniques.