Without doing it, I can point out a few things:
1. The situation is symmetric (it's cos(z1)cos(z2) = 1, symmetric under switching z1 and z2), so you only need to derive one of those equations.
2. If you just expand out by the definition, you'll get something LIKE a(z1, z2) + ib(z1, z2) = 1, and so two equations, a(z1,z2)=1 and b(z1,z2)=0.
Then substitue one into the other, trying to eliminate the parameter that doesn't show up in your ultimately desired equation. This is the direct way, and seems like it should work in principle.
3. There are all kinds of factoring and algebraic relationships going on (cosh^2 - sinh^2 = 1, sin^2 + cos^2 = 1, cosh(a+b) + cosh(a-b) = 2cosh(a)sinh(b), etc etc) for you to exploit.
4. There's so much symmetry going here that there might be a really easy clever way of deriving it. I don't see it at the moment, but I bet it's there.
5. In addition to getting two equations from a(z1, z2)+b(z1, z2)i = 1, you could look at it in polar form to get two equations from something like
R(z1, z2)exp(theta(z1, z2)i) = 1 exp(2 pi n i) (n an integer). In particular, look at the norms, right off the bat.