Если это ещё актуально, то смотрите что тут нужно делать.
Во-первых, надо понять что вообще требуется. Как вот понять фразу "программа должна работать для различных" - различных чего? Надеюсь, что это различные значения задаваемой точности ε (которая у Вас задаётся как ε=1e-4, но может быть в принципе и другой), а не, например, какие-либо числа помимо пятидесяти - это была бы жесть (а почему - смотрите следующие пункты !).
Во-вторых, √50 = 7•√(1+1/49) - в таком представлении первоначальный корень сводится к вычислению функции √(1+x) с малым параметром x=1/49, для чего есть асимптотические разложения (смотрите следующий пункт !).
В-третьих, чем, извиняюсь за выражение, муждохаться с разложением в ряд функции √(1+x), не лучше ли сразу взять уже готовый для неё ряд Маклорена (в окрестностях нуля, то есть при х≈0 или, по крайней мере когда |х| не превосходит 1) ? Вот этот ряд (смотрите комментарий №1):
В-четвёртых, вычислять значения (ограниченных) сумм этого ряда лучше всего используя рекурренцию, тогда и никаких степеней с факториалами впридачу не понадобится от слова совсем. Здесь это можно проделать так: присмотритесь как устроены слагаемые этого ряда и сразу заметите, что он начинается с единицы, а потом, в самом конце, ещё всю найденную сумму надо будет умножить на 7. Так вот, можно сразу задать переменные n=0, p=7 и s=7, а потом на каждом шаге итерационного цикла прибавлять к n единицу, умножать p на что? Правильно, вот на это:
(-x/4)•2n•(2n-1)/n² или, что то же самое,
на х•(0.5/n-1). Затем служебную переменную eps приравниваем к p/(1-2n) и прибавляем её к значению суммы s, которую и надо вычислить. Цикл прекращается, когда |eps| становится меньше заданной точности ε.
В-пятых, можно в принципе варьировать и задаваемую точность ε, и значение числа а, корень из которого нужно вычислить, только это несколько жестковато, но если варьировать число а в весьма широких пределах, то получается уже полная жесть ! Я же в нижеприведённом коде сделала так, чтобы p и s изначально инициировались семёрками, а это означает, что для того, чтобы модули |х| всегда лежали в пределах единицы, число а должно находиться в диапазоне (0;98), а не где попало. Вот, кстати, и сам код (если C++ Вам, конечно же, подходит):
#include <cmath>
#include <iostream>
using namespace std;
int main()
{
double a, x, p, s, n, ε, eps;
for (;;)
{
cout << "a ε: ";
cin >> a >> ε;
x = a / 49. - 1.;
n = 0.; p = s = 7.;
eps = 1e303;
while (fabs(eps) > ε)
{
n += 1.;
p *= x * (0.5 / n - 1.);
eps = p / (1. - 2. * n);
s += eps;
}
cout.precision(16); cout << s;
cout.precision(6); cout <<
", err=" << s - sqrt(a) <<
", " << n << " terms" << endl;
}
}
Видите какая программа получается небольшая ? А вот что она печатает:
(Смотрите комментарий №2 !)
По-моему, нормально этот код работает. Только мне смешно то, что корни любых (неотрицательных) чисел вычисляются гораздо быстрее и эффективнее вышеизложенного способа при помощи чрезвычайно архаичных алгоритмов, известных с покон веку, а вовсе даже не при помощи разложения в ряды, да и посмотрите лишний раз на скрин с результатами работы программы, на котором видно, что для вычисления √50 с точностью 0,0001 требуется всего-то навсего три слагаемых (terms по-английски), так что, похоже, в данном случае вычисления на компьютере совсем не обязательны - достаточно простого калькулятора или даже ручных вычислений на бумажке. И вообще всё это действительно получается как-то не очень серьёзно...
Приближенные вычисления значений функций с помощью 
