Mass-calculation of 1/k mod p with a fixed k and prime p

[Ilya Zakharevich <nospam-abuse@ilyaz.org>]

I think about a cheap way to implement forprimestep() — without current
pessimizations.  To do this, it seems that given k, I need to find 1/k
mod p for all “small” p mutually prime with k.  (Say, the list of p ≤
pMAX<10¹¹ is already pre-calculated.  Already covering the case when k
pMAX < 2⁶⁴ may be quite useful.)

  Is there a way to do it quickly (possibly pre-calculating suitable
  data depending on k)?

⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜

So far I can only see the following (silly) way: if k is “small”,
pre-calculate 1/l mod k for every l mod k; store in array INV[l].
Then looping through p, it is cheap to find p mod k, hence r ≔ 1/p mod k.
Hence rp ≡ 1 mod k, say rp - 1 = sk.  So 1/k ≡ -s mod p, and all one
needs to find is p-(rp-1)/k.

Likewise if k=k₁k₂ with small k₁ and k₂ (twice slower).  Etc.

Do I miss something obvious?!

Thanks in advance,
Ilya