С++ вычисление с точностью

Какая проблема вычислить факториал? Это почти так же просто как вычисление целой степени.

#include <iostream>
#include <iomanip>
using namespace std;
pair<double, int> sum(double x, double eps = 1e-6) {
auto k = 0;
auto s = 0.0;
auto m = 0.0;
do {
++k;
auto fk = 1ULL;
for (auto i = 1; i <= k + 1; ++i) fk *= i;
auto x2k = x;
for (auto i = 1; i < k + k; ++i) x2k *= x;
x2k *= k & 1 ? -1 : 1;
m = x2k / (fk * k);
s += m;
} while (fabs(m) > eps);
return { s, k };
}
int main() {
double x;
do {
cout << "x: ";
cin >> x;
} while (x <= 0 || x > 1);
auto[value, quantity] = sum(x);
cout << fixed << setprecision(15)
<< "value: " << value << '\n'
<< "quantity: " << quantity << '\n';
system("pause");
}

Я тоже не вижу тут проблемы. Вот правильный цикл:

y=-x*x; p=1; s=0; k=0; do { k++; p*=y/(k+1); e=p/k; s+=e; } while (abs(e) > 1e-6); cout.precision(16); cout << s << ' ' << k << '\n';

Проблема есть только в представлении чисел в ЭВМ с ограниченной точностью, из за чего даже теоретически абсолютно сходящиеся ряды вдруг при некоторых не таких уж и больших аргументах х вдруг начинают считаться с дикими ошибками. Это происходит и с синусом, и с косинусом, и иже с ними. И тут, по-моему, та же фишка! В написанном выше цикле условие выхода из него - последний считаемый член ряда по модулю не превосходит установленной ошибки в одну миллионную. Но это всё работает только в теории! А здесь чисто фактически при |x|=8..10 происходит чёрт знает что! А наоборот, при, например, -π≤х≤π ряд ведёт очень хорошо и результаты получаются правильными.