Why that rotation does what it does
It's really probably better to think of multiplication as making rotations, although of course division and multiplication are really the same. In that spirit, look at this:
$$
\frac{z}{b-ai}=z(\frac{b+ai}{a^2+b^2})
$$
Here we can see that $z$ is being multiplied by a complex number of unit length (here you can see it's already normalized) on the right hand side, which we know effects a rotation on the complex plane. It isn't normalized , so it will contribute dilation, but that's OK: two steps later the inverse transformation contributes a dilation factor that cancels this one's dilation out. All that matters is that this transformation carries our line onto the $x$-axis.
So on the step you're at, the line has been translated so that it goes through the origin, so its points all fall on $ax+by=0$.
Let $(x',y')$ be on $ax+by=0$. After rotation by the above mapping, $(x',y')\mapsto (bx'-ay',ax'+by')/(a^2+b^2)$. But since $ax'+by'=0$, this is just $((bx'-ay')/(a^2+b^2),0)$, so the line has been laid on the $x$-axis, as advertised.
From the last paragraph you can see how this can be reverse engineered. If one had applied an arbitrary rotation $(x,y)\mapsto((cx-dy)/(c^2+d^2),cx+dy)$, and they were armed with the equation $ax+by=0$, and you asked them "What choice of $c$ and $d$ makes this rotation put $x,y$ on the $x$ axis? They might say "how about $c=a$ and $d=b$?"
Moved this here since it turned out to be a digression.
Looking at the map as an $\Bbb R$ linear transformation on $\Bbb C$, and viewing the inputs as column vectors $[x,y]^\top$ that correspond to $x+yi$, you can confirm that this matrix accomplishes the transformation by multiplication on the left of the column vectors:
$$
\frac{1}{a^2+b^2}\begin{bmatrix}b&-a\\a&b\end{bmatrix}
$$
It's a matrix with determinant $1$ (of course it should be, since it represents a complex number of unit modulus) and so it's a rotation of $\Bbb R^2$ (which is being identified with $\Bbb C$.)
Best Answer
This solutions depends on whether $a ≠ 0$ or $b ≠ 0$.
Translate the problem to the origin. $z ↦ z-\tfrac{c}{a}$ or $z ↦ z - i\tfrac{c}{b}$
Spin the line onto the real line. $z ↦ \tfrac{z}{b-ai}$
Reflect the plane at the real line. $z ↦ \overline z$
Spin back to the former line . $z ↦ (b-ai) z$
Translate to your former point. $z ↦ z + \tfrac{c}{a}$ or $z ↦ z + i\tfrac{c}{b}$
All in all I get: $(b-ai)\overline{\left(\tfrac{z-\tfrac{c}{a}}{b-ai}\right)} + \tfrac{c}{a} = \tfrac{2ic + (b-ai)\overline z}{b+ai}$ (And the same with the other variant.)
Note: I think something's wrong with this as it should be $2ic$ instead of $2iac$. I can't figure it out, but I really want to know what my error is.
This is now corrected.