Here's how to do Bezier patches on the CPU, where cpp_dec_float_100 is a high-precision real number:
vertex_3 bezier(const double u, const double v, vector<vector<vector<float> > > &control_points, size_t num_wide, size_t num_tall)
{
vertex_3 ret;
cpp_dec_float_100 ret_x;
cpp_dec_float_100 ret_y;
cpp_dec_float_100 ret_z;
cpp_dec_float_100 u_long = u;
cpp_dec_float_100 v_long = v;
const size_t Ni = num_wide - 1;
const size_t Nj = num_tall - 1;
for (size_t i = 0; i <= Ni; i++)
{
for (size_t j = 0; j <= Nj; j++)
{
cpp_dec_float_100 binpow_u = binomial(Ni, i) * cached_pow(u_long, cpp_dec_float_100(i)) * cached_pow(cpp_dec_float_100(1 - u_long), cpp_dec_float_100(Ni - i));
cpp_dec_float_100 binpow_v = binomial(Nj, j) * cached_pow(v_long, cpp_dec_float_100(j)) * cached_pow(cpp_dec_float_100(1 - v_long), cpp_dec_float_100(Nj - j));
ret_x += binpow_u * binpow_v * control_points[i][j][0];
ret_y += binpow_u * binpow_v * control_points[i][j][1];
ret_z += binpow_u * binpow_v * control_points[i][j][2];
}
}
ret.x = mp_to_double(ret_x);
ret.y = mp_to_double(ret_y);
ret.z = mp_to_double(ret_z);
return ret;
}