Evaluate the Riemann Zeta Function at a Complex Number Introduction
I found this question that was closed because it was unclear, yet it was a nice idea. I'll do my best to make this into a clear challenge.
The Riemann Zeta function is a special function that is defined as the analytical continuation of

to the complex plane. There are many equivalent formulas for it which makes it interesting for code golf.
Challenge
Write a program that takes 2 floats as input (the real and imaginary part of a complex number) and evaluates the Riemann Zeta function at that point.
Rules

Input and output via console OR function input and return value
Built in complex numbers are not allowed, use floats (number, double, ...)
No mathematical functions except + - * / pow log and real valued trig functions (if you want to integrate, use the gamma function, ... you must include this functions definition in the code)
Input: 2 floats
Output: 2 floats
Your code must contain value that gives theoretically arbitrary precision when made arbitrary large/small
The behaviour at input 1 is not important (this is the only pole of this function)

Shortest code in bytes wins!
Example Input and Output
Input:

2, 0

Output:

1.6449340668482266, 0

Input:

1, 1

Output:

0.5821580597520037, -0.9268485643308071

Input:

-1, 0

Output:

-0.08333333333333559, 0
Python - 385

This is a straightforward implementation of Equation 21 from http://mathworld.wolfram.com/RiemannZetaFunction.html This uses Python's convention for optional arguments; if you want to specify a precision, you can pass a third argument to the function, otherwise it uses 1e-24 by default.

import numpy as N
def z(r,i,E=1e-24):
 R=0;I=0;n=0;
 while(True):
  a=0;b=0;m=2**(-n-1)
  for k in range(0,n+1):
   M=(-1)**k*N.product([x/(x-(n-k))for x in range(n-k+1,n+1)]);A=(k+1)**-r;t=-i*N.log(k+1);a+=M*A*N.cos(t);b+=M*A*N.sin(t)
  a*=m;b*=m;R+=a;I+=b;n+=1
  if a*a+b*b&lt;E:break
 A=2**(1-r);t=-i*N.log(2);a=1-A*N.cos(t);b=-A*N.sin(t);d=a*a+b*b;a=a/d;b=-b/d
 print(R*a-I*b,R*b+I*a)