Understanding Intel's algorithm for reducing a polynomial modulo an irreducible polynomial
The Galois field $GF(2^{128})$ has many different "concrete" representations. One popular representation is using polynomials in $GF(2)[x]$ (i.e. with coefficients in $GF(2)$) modulo some irreducible polynomial of degree $128$, say $x^{128}+x^7+x^2+x^1+x^0$. The magic constants $7,2,1,0$ come from this particular irreducible polynomial. You don't see $0$ since there is no need to shift by $0$. Similarly, the magic constants $57,62,63,64$ are the complements of the constants above with respect to $64$. Again, you don't need to shift by $64$ since the registers are $64$ bits wide. If we had used some other irreducible polynomial, the constants would have been different.

Regarding step 4, $[H_1:H_0]$ results from XORing $[E_1:E_0],[F_1:F_0],[G_1:G_0],[X_3:D]$. This implements the operation of MODing by $x^{128} + x^7 + x^2 + x^1 + x^0$. The idea is that modulo this polynomial, $x^{128} = x^7+x^2+x^1+^0$; addition here is XOR. The different summands corresponds to these different monomials – see if you can find the correspondence.