My initial idea would be to exploit a function like
$$ g(z) = \frac{f(z)}{\sin z}.$$
What is nice about $g(z)$ is that it is holomorphic, and it is holomorphic everywhere. It is true that $1 / \sin z$ has poles at $z = k\pi$ for $k \in \mathbb N$, but these poles are cancelled out by the zeroes of $f(z)$ in the numerator.
Our ultimate goal is to prove that $g(z) = 1$. Having established that $g(z)$ is holomorphic, we can try to show that $g(z) = 1$ using familiar theorems about holomorphic functions such as the maximum modulus principle and Liouville's theorem.
First, let us obtain a bound on $|g(x + iy)|$. Assuming I did my algebra correctly, we have an inequality:
$$|g(x + iy)| \leq \frac{e^{|y|}}{|\sin(x + iy)|}=\frac{\sqrt{2}e^{|y|}}{\sqrt{\cosh 2y - \cos 2x}}$$
Notice that this inequality is not uniformly useful! For example, at $x + iy = 0$, the expression on the right-hand side is infinite, so the inequality tells us nothing useful at $x + iy = 0$. However, there are some special contours within the complex plane where the expression on the RHS is very useful for bounding $|g(x + iy)|$, and also easy to analyse:
- On vertical lines of the form $x = (N + \frac 1 2) \pi$ where $N \in \mathbb N$, the $\cos 2x$ term vanishes. It is easy to see that the expression on the RHS is bounded by $2$.
- On horizontal lines of the form $y = M$ where $|M| \gg 0$, we can approximate $\cosh 2y$ by $e^{2|y|}$. The expression on the RHS is bounded by $2 / \sqrt{1 - 2e^{-2|M|}}$, which tends to $2$ as $M \to \infty$.
The maximum modulus theorem allows us to use our knowledge about $|g(x + iy)|$ on these special contours to bound $|g(x + iy)|$ in other places. Let us define the rectangular domain
$$D_{N, M} =\{x + iy \in \mathbb C \ : -(N+\tfrac 1 2)\pi \leq x \leq (N + \tfrac 1 2)\pi, \ -M \leq y \leq M \},$$ for arbitrary $ N \in \mathbb Z,$ $M \gg 0.$ Using our above observations about the behaviour of $|g(x + iy)|$ on the boundary of this rectangular domain, the maximum modulus principle tells us that
$$ z \in D_{N,M} \implies |g(z)| \leq \frac{2}{\sqrt{1 - 2 e^{-2M}}} \ .$$
And since any $z \in \mathbb C$ is contained in $D_{N,M}$ for sufficiently large $N$ and $M$, we have
$$ z \in \mathbb C \implies |g(z)| \leq 2 \ .$$
Having shown that $|g(z)| \leq 2$, Liouville tells us that $g(z)$ is constant, i.e. that $f(z) = C \sin z$ for some $C \in \mathbb C$. We can fix the value of $C$ to be $1$ using the fact that $f'(0) = 1$, and this completes the proof.
