I discovered a way to evaluate infinite products of even analytic functions with high accuracy.
$$ \prod_{k=1}^{\infty} f(k^2) \approx \prod_{k=1}^{\infty} \left(1-\frac{A_1}{k^2}+\frac{A_2}{k^4}-\frac{A_3}{k^6}+\cdots+\frac{A_n}{k^{2n}} \right) = \frac{\sin (\pi a_1) \sin (\pi a_2) \cdots \sin (\pi a_n)}{\pi^n a_1 a_2 \cdots a_n} $$
Here:
$$A_1=\sum^{n}_{j=1} a_j^2,~~~~~A_2=\sum^{n}_{j\neq l} a_j^2 a_l^2,~~~~~A_3=\sum^{n}_{j\neq l\neq m} a_j^2 a_l^2 a_m^2,~~~\cdots~~~A_n=\prod^{n}_{j=1} a_j^2$$
The idea is as follows:
$$ \prod_{k=1}^{\infty} \left(1-\frac{a^2}{k^2} \right)=\frac{\sin (\pi a)}{\pi a}$$
Since the product converges, we can multiply two or more products by term:
$$ \prod_{k=1}^{\infty} \left(1-\frac{a^2}{k^2} \right)\left(1-\frac{b^2}{k^2} \right)=\prod_{k=1}^{\infty} \left(1-\frac{a^2+b^2}{k^2}+\frac{a^2 b^2}{k^4} \right)=\frac{\sin (\pi a) \sin (\pi b)}{\pi^2 a b}$$
We can easily recognize the Vieta formulas:
$$a^2+b^2=A_1,~~~~~~a^2b^2=A_2$$
$$x^2-A_1x+A_2=0$$
The connection to the coefficients of the Taylor series in general is obvious and we can always find $a_j$ by solving the equation:
$$x^n-A_1 x^{n-1}+A_2 x^{n-2}- \cdots + A_{n}=0$$
$$a^2_j=x_j,~~~~~~j=1,\dots,n$$
Of course $a_j$ can be complex, the product will still be real.
Now some examples to see how it works. Consider the following two products:
$$\prod_{k=1}^{\infty} k \tanh \frac{1}{k} = \prod_{k=1}^{\infty} \left(1-\frac{1}{3k^2}+\frac{2}{15k^4}-\frac{17}{35k^6}+\frac{62}{2835k^8}-\cdots \right)$$
$$\prod_{k=1}^{\infty} \frac{\sqrt{\pi}}{2} k ~\mathbb{erf} \frac{1}{k} = \prod_{k=1}^{\infty} \left(1-\frac{1}{3k^2}+\frac{1}{10k^4}-\frac{1}{42k^6}+\frac{1}{216k^8}-\cdots \right)$$
Coincidentally, the second term is the same for both functions, so we have the first approximation:
$$1>\prod_{k=1}^{\infty} k \tanh \frac{1}{k}>\frac{\sqrt{3}}{\pi} \sin \frac{\pi}{\sqrt{3}}=0.535131$$
$$1>\prod_{k=1}^{\infty} \frac{\sqrt{\pi}}{2} k ~\mathbb{erf} \frac{1}{k}>\frac{\sqrt{3}}{\pi} \sin \frac{\pi}{\sqrt{3}}=0.535131$$
To get the following approximations we solve quadratic, cubic and fourth-order equations and obtain the roots.
For example, the first product gives us:
$$a,b=\sqrt{\frac{5\pm i \sqrt{65}}{30}}$$
$$\prod_{k=1}^{\infty} k \tanh \frac{1}{k}<\frac{\sin \pi a \sin \pi b}{\pi^2 a b}=0.620862$$
The second product:
$$a,b=\sqrt{\frac{5\pm i \sqrt{95}}{30}}$$
$$\prod_{k=1}^{\infty} \frac{\sqrt{\pi}}{2} k ~\mathbb{erf} \frac{1}{k}<\frac{\sin \pi a \sin \pi b}{\pi^2 a b}=0.649760$$
The sequence of partial products up to $50$ is compared to the numerical value of the infinite product and three approximants in the graphs below.
Another interesting product is connected to the number $e$ (see here):
$$\prod_{k=2}^{\infty} e \left(1-\frac{1}{k^2} \right)^{k^2}=\prod_{k=2}^{\infty} \left(1-\frac{1}{2k^2}-\frac{5}{24k^4}-\frac{5}{48k^6}-\cdots \right)=\frac{\pi}{e^{3/2}}=0.70098$$
It converges extremely slow. But the second approximant already does a very good job of approximating the exact value.
Finally, let's consider much simpler case:
$$\prod_{k=2}^{\infty} \cos \frac{1}{k}=\prod_{k=2}^{\infty} \left(1-\frac{1}{2k^2}+\frac{1}{24k^4}-\cdots \right)=0.719109$$
The first approximant is the same as in the previous product:
$$P_1=\frac{2\sqrt{2}}{\pi} \sin \frac{\pi}{\sqrt{2}}=0.7163756$$
And it's remarkably close to the numerical value. The third one already gives $6$ correct digits:
Is this method useful for evaluating infinite products? It's easy to implement in any algebra system. I know it can be generalized using gamma functions. Can you provide some reference on the topic?