The zeta function provides incorrect numerical output for numbers that aren't even integers greater than some number between 33.44 and 33.45. The output for 33.44 is 1.000000000128721, and for 33.45 is 0.0. It should be clear from the definition ( http://dlmf.nist.gov/25.2 ) that for finite positive real values zeta should be > 1 (>= 1.0 given the accuracy limitations of floating point). For values this large, direct calculation using the definition should suffice.
Another potentially useful thing to do would be to provide a function zetam1(s) that is formally equal to zeta(s) - 1 that can work for large values of s where the difference between zeta(s) and 1 would be lost in floating point cancellation.