It's different from Eq. (b). Is $f(q_k)$ equal to $f\left(\frac{q_{k+1}+q_k}{2}\right)$?
Notice in the very first line that $f(q_k)$ multiplies a delta function $\delta(q_k-q_{k+1})$. This delta forces, in anything multipliying it, that $q_k = q_{k+1}$. This is the well-known property:
$$f(x)\delta(x-a)=f(a)\delta(x-a).$$
Now, this means that $$f(q_k)\delta(q_k-q_{k+1})=f\left(\frac{q_k+q_k}{2}\right)\delta(q_k-q_{k+1})=f\left(\frac{q_k+q_{k+1}}{2}\right)\delta(q_k-q_{k+1}).$$
This is, in fact, a convention. Here it may seem weird, but it appears because when you evaluate matrix elements $\langle q'|H|q\rangle$ of the Hamiltonian you need assume one operator ordering convention. The point is that while this seems superflous for a function of just the coordinates $f(q)$, in general $H$ is a function $H(q,p)$ and $q$ and $p$ do not commute.
With the operator ordering that Peskin & Schroeder choose you end up evaluating things at the average of the endpoints of the subintervals.
But again this is just a convention. Check Weinberg's The Quantum Theory of Fields, Vol. 1, Chapter 9. He does the same derivation but he employs another convention. Specifically he orders $H$ such that all position are to the left of all momenta.
In the Eq. (c), he said "we introduce a complete set of momentum eigenstates", but how to introduce it?
Well, now you have a function $f(p)$. What are its matrix elements in coordinate basis $\langle q'|f(P)|q\rangle$? The answer is that by definition you know how to compute $f(P)$ acting on the momentum eigenstates $$f(P)|p\rangle = f(p)|p\rangle,$$
therefore you should use this. Introducing a complete set means using the identity $$\mathbf{1}=\int |p\rangle \langle p| dp,$$
well-known from quantum mechanics.
This will lead you to just evaluate $\langle q|p\rangle$ which is just $\frac{1}{\sqrt{2\pi}}e^{ipq}$.