Let me guess about the next optimizations you have made after this?

Firstly, bitwise xor is associative, so there's no point in computing the ^ HASH_LEN in the innermost loop, you can put that to m8..i8 instead of m9..i9. Eg.

m8 = (m7 ^ q8) * H_PRIME ^ HASH_LEN;
m9 = m8 ^ q9;
[download]

Secondly, I wouldn't bother computing c8..i8 unless m9 is correct.

Thirdly, there is no point looping over q9. Instead, notice that as q9 takes all possible values from 0..127, m9 will take the values from (m8&~127)..(m8&~127)+127 in a different order. Among those 127 possible values of m9, there can only be one that has the right value modulo 1001. You can compute that value directly, and try the corresponding single value for q9, checking if it lies in the range.

However, I wonder, does the magic formula really have to use exactly 1001 as the modulus? I think it could use any modulus between 1001 and 9999.

To search for the magic formula in that case, first compute t = gcd(m9 - 1000, d9 - 500) (be careful with underflows). Then compute the small prime factors of t, and using this, all divisors of t between 1001 and 9999, and check all those divisors as the modulus against the rest of the numbers. In the rare lucky case that t is 0, you will have to try all moduluses (or all divisors of c9 - 100).

Most of the time, t should be between 1 and 1000, in case you don't have to compute a factorization because no modulus can work. Because of this, you should probably not bother to compute c9..i9 and c8..i8 unless t is greater than 1000 or zero.

This would certainly be more difficult to implement, and might take more time for each individual magic string, but it also has a higher chance of finding a formula, so it might be worth.

Update: fixed a typo from (m8&!127)+127 to (m8&~127)+127 that eyepopslikeamosquito noticed.