double bernoulli[] = { 1/6, 1/30, ... , 236464091/2730 };

double coth(double x, double epsilon)
{
   double last_n = x/3;
   double sum = last_n;
   double n = 0;

   for(int i = 2; i<=12, i++)
   {
       n = (pow(-1,i+1)*pow(2, 2*i))/faculty(2*i) * bernoulli[i] * pow(x, 2*i-1)
       sum += n;
       if((n < epsilon) && (last_n < epsilon))
           return sum;

   }
   return sum;
}