There are three important points here to consider.
1. We work in $\mathbb{F}_2[X]$. This means that we do additions and multiplications of binary polynomials, i.e. polynomials whose coefficients are 0 or 1. The addition of two polynomials is then a bitwise XOR; there is no carry propagation. Similarly, the multiplication is called a "carry-less" multiplication. There are some explanations about that here.
In particular, the "modular reduction" is about divisions of polynomials, not integers. The modulus is the binary polynomial $X^{128} + X^7 + X^2 + X + 1$ (it is an irreducible polynomial over $\mathbb{F}_2$, which implies that when computing stuff modulo that polynomial, we work in a finite field).
2. There is a fair amount of confusion about encoding. This is a problem from the GCM specification. The input bytes are to be interpreted as bits, i.e. the coefficients of a binary polynomial. "Endianness" is about how to map chunks of an element (e.g. bits in a byte) to their mathematical value. For instance, if you consider a byte of numerical value 0xC9 (201 in decimal) and write it down as a sequence of bits, you normally use big-endian notation, i.e. you write it as "11001001": the leftmost bit is the one with the highest value ($2^7 = 128$), the rightmost is the one with the lowest value ($2^0 = 1$). This is what happens with integers; with binary polynomials, the big-endian interpretation of a byte is to map the leftmost bit to the coefficient of $X^7$, and the rightmost to the coefficient of $X^0$; hence, "11001001" would be decoded as $X^7 + X^6 + X^3 + 1$.
The concept of endianness works at all levels; i.e. not only bits within bytes, but also for bytes within a sequence of bytes. The "pure big-endian" interpretation of a sequence of 16 bytes would be to decode each byte in big-endian (as above), then assemble all bytes with the most significant byte coming first, i.e. the first byte in the sequence corresponding to binary coefficients $X^{127}$ to $X^{120}$, the second byte being for coefficients $X^{119}$ to $X^{112}$, and so on. BearSSL decodes input values in pure big-endian interpretation; internally, it stores them not into bytes, but into 32-bit words:
/*
* Decode the block. The GHASH standard mandates
* big-endian encoding.
*/
yw[3] ^= br_dec32be(src);
yw[2] ^= br_dec32be(src + 4);
yw[1] ^= br_dec32be(src + 8);
yw[0] ^= br_dec32be(src + 12);
Here, yw[0] will contain coefficients $X^0$ to $X^{31}$, yw[1] coefficients $X^{32}$ to $X^{63}$, and so on.
Now there is a lie here. The GHASH standard does not really mandate big-endian encoding. It, in fact, mandates the complete opposite, which is pure little-endian encoding, i.e. not only the first byte should be for coefficients $X^0$ to $X^7$, but the bits within that byte are interpreted in "unnatural" order, i.e. the most significant bit (of numerical value 128) is for coefficient $X^0$, not $X^7$. This little-endian convention of bits within bytes is very inconvenient for computations, especially since BearSSL implements the carryless multiplication by way of multiplications on plain integers (i.e. multiplications with carries, but using masks and holes to remove these carries where necessary). If BearSSL were to use pure little-endian convention, it would then have to swap bits around within bytes, a cumbersome operation which would degrade performance and increase code size.
Instead, BearSSL uses a trick, which works on binary polynomials. Namely, if you define bit reversal:
$$ \mathrm{rev}_n \Big(\sum_{i=0}^{n-1} a_i X^i\Big) = \sum_{i=0}^{n-1} a_{n-1-i} X^i $$
then, when multiplying two polynomials $A$ and $B$ of degree less than $n$ (i.e. that fit on $n$ bits):
$$ \mathrm{rev}_n(A) \times \mathrm{rev}_n(B) = \mathrm{rev}_{2n-1}(A \times B) $$
i.e. if you reverse the two operands, you get the reverse of the product. This does not work at all on plain integers, because additions on plain integers have carries, and carries propagate in a definite direction (right-to-left, in the usual notations), and that does not work "in a mirror". But with carryless multiplication, there are no carries, and everything can be mirrored.
What that means here is that, by decoding 16 bytes with pure big-endian notations, BearSSL does not get the expected polynomial $A$, but $\mathrm{rev}_{128}(A)$. The product is then computed over reversed operands, and yields the reverse of the expected value. If that result is to be reencoded back with pure big-endian notation, then everything works fine! The compatibility of operations with bit reversal ("mirroring") means that you can really choose the encoding convention you want, and you get the right result, provided that you stick to your convention throughout.
... Except for one small detail. If you look at the equation above, the product of two bit-reversed 128-bit polynomials yields the bit-reversed result over 255 bits, not 256. The BearSSL code ends up with a 256-bit result in zw[], and that value is shifted by one bit, because of that reversed convention issue. Thus, the code must include a shifting step to put it back where it should, and this is the first half of the code you reproduce in your question:
/*
* GHASH specification has the bits "reversed" (most
* significant is in fact least significant), which does
* not matter for a carryless multiplication, except that
* the 255-bit result must be shifted by 1 bit.
*/
zw[0] = c0 << 1;
zw[1] = (c1 << 1) | (c0 >> 31);
zw[2] = (c2 << 1) | (c1 >> 31);
zw[3] = (c3 << 1) | (c2 >> 31);
zw[4] = (d0 << 1) | (c3 >> 31);
zw[5] = (d1 << 1) | (d0 >> 31);
zw[6] = (d2 << 1) | (d1 >> 31);
zw[7] = (d3 << 1) | (d2 >> 31);
Thus, these lines are really about making it as if everything had been done in pure little-endian convention, whereas it was done in pure big-endian convention.
3. Modular reduction is over binary polynomials. We work modulo polynomial $M = X^{128} + X^7 + X^2 + X + 1$. This means that for any polynomial $P$, we can add (or subtract, but in $\mathbb{F}_2[X]$, addition and subtraction are the same thing, i.e. a XOR) arbitrary multiples of $M$. The "modular reduction" of the 255-bit result $P$ really is an addition of multiples of $M$ until we get something that fits on 128 bits.
Write the polynomial to reduce $P = \sum_{i=0}^{254} p_i X^i$. Take one coefficient $p_i$ for some $i \ge 128$; we want to "cancel that bit", in case it as value $1$; the trick is to add $M\times X^{i-128}$. In other words, we add:
$$ p_i M \times X^{i-128} = p_i X^i + p_i X^{i-121} + p_i X^{i-126} + p_i X^{i-127} + p_i X^{i-128} $$
The first part ($p_i X^i$) merely clears the bit at position $i$ (it's a XOR of $p_i$ with itself), but the other parts "inject" (XOR) bit $p_i$ at other places down the sequence of coefficients. Crucially, these are relatively far (121 or more bits down the sequence), which allows us to process bits by large chunks (this is the reason the modulus was chosen with that value). Consider the first iteration of that loop:
/*
* We now do the reduction modulo the field polynomial
* to get back to 128 bits.
*/
for (i = 0; i < 4; i ++) {
uint32_t lw;
lw = zw[i];
zw[i + 4] ^= lw ^ (lw >> 1) ^ (lw >> 2) ^ (lw >> 7);
zw[i + 3] ^= (lw << 31) ^ (lw << 30) ^ (lw << 25);
}
The word zw[0] contains the coefficients $X^{255}$ to $X^{224}$. Since the coefficients for $X^i$ must be injected at positions $X^{i-121}$, $X^{i-126}$, $X^{i-127}$ and $X^{i-128}$, these 32 bits end up being injected in various positions in zw[3] and zw[4]. In the line that updates zw[i+4] (i.e. zw[4] during the first iteration), we see the XOR with lw, which is the injection of bits at positions $X^{i-128}$; the lw >> 1 is the injection of bits at positions $X^{i-127}$; and so on.
Nominally, we should clear the high bits (zw[0] to zw[3]) but this is optimized out, since we don't look at these bits again; the useful result, after reduction, is in words zw[4] to zw[7]:
memcpy(yw, zw + 4, sizeof yw);
And that's the complete modular reduction.
