The magic numbers are best explained here:
/**
* Lots of magic numbers :)
* To understand what's going on below, note that
*
* (1) q = 2^252 + q0 where q0 = 27742317777372353535851937790883648493.
* (2) s11 is the coefficient of 2^(11*21), s23 is the coefficient of 2^(^23*21) and 2^252 = 2^((23-11) * 21)).
* (3) 2^252 congruent -q0 modulo q.
* (4) -q0 = 666643 * 2^0 + 470296 * 2^21 + 654183 * 2^(2*21) - 997805 * 2^(3*21) + 136657 * 2^(4*21) - 683901 * 2^(5*21)
*
* Thus
* s23 * 2^(23*11) = s23 * 2^(12*21) * 2^(11*21) = s3 * 2^252 * 2^(11*21) congruent
* s23 * (666643 * 2^0 + 470296 * 2^21 + 654183 * 2^(2*21) - 997805 * 2^(3*21) + 136657 * 2^(4*21) - 683901 * 2^(5*21)) * 2^(11*21) modulo q =
* s23 * (666643 * 2^(11*21) + 470296 * 2^(12*21) + 654183 * 2^(13*21) - 997805 * 2^(14*21) + 136657 * 2^(15*21) - 683901 * 2^(16*21)).
*
* The same procedure is then applied for s22,...,s18.
*/
If it's not yet clear, the constants 666643, 470296, 654183, -997805, 136657, -683901 are the 21-bit limbs of the 125-bit number -q0.
Proof (in python):
>>> limbs = [666643<<(0*21), 470296<<(1*21), 654183<<(2*21), -997805<<(3*21), 136657<<(4*21), -683901<<(5*21)]
>>> sum(limbs)
-27742317777372353535851937790883648493L