алгоритм квадратного корня C++

Question

Я не могу понять, почему этот алгоритм входит в бесконечный цикл, если введенное число более 12 цифр. Кто-нибудь может понять, почему это никогда не закончится? Благодарю. Я просто обновил алгоритм, чтобы использовать функцию fabs(), и все равно получаю бесконечный цикл.

double squareroot(double x) 

{ /* computes the square root of x */ 

/* make sure x is not negative .. no math crimes allowed! */ 
assert(x >= 0); 
if (x==0) return 0; 

/* the sqrt must be between xhi and xlo */ 
double xhi = x; 
double xlo = 0; 
double guess = x/2; 

/* We stop when guess*guess-x is very small */ 

while (abs(guess*guess-x) > 0.00001) 
{ 
    if (guess*guess > x){ 
     xhi = guess; 
    } 

    else { 
     xlo = guess; 
    } 

    guess = (xhi + xlo)/2; 
} 
return guess; 
} 

Answer 1

Я считаю, что вы должны использовать относительную погрешность для прекращения, а не абсолютная погрешность.

while (abs((guess*guess-x)/guess) > 0.00001) 

В противном случае это займет очень много времени (это не бесконечный цикл) для вычисления квадратного корня из очень длинных значений.

http://en.wikipedia.org/wiki/Approximation_error

Ура!

EDIT: кроме того, как указано ниже в комментариях, он достоин того, чтобы проверить, если guess уже догадались, с тем чтобы избежать бесконечного цикла с некоторыми конкретными случаями угловых.

Answer 2

http://floating-point-gui.de/errors/comparison/

Хех, похоже, вы не должны использовать абс(), как это. Есть некоторые случаи, когда он должен остановиться, но я не буду, поскольку это ограничена точность.

Вместо этого используйте FABS()

http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm

Answer 3

Я бы сказал, что, когда числа достаточно велики, вы не можете использовать абсолютное значение эпсилона, потому что оно не вписывается в точность.

Попробуйте вместо этого использовать относительное сравнение. Рассмотрим следующую функцию, чтобы проверить, если 2 двойники равны:

bool are_eql(double n1, double n2, double abs_e, double rel_e) 
{ 
    // also add check that n1 and n2 are not NaN 
    double diff = abs(n1 - n2); 
    if (diff < abs_e) 
    return true; 
    double largest = max(abs(n1), abs(n2)); 
    if (diff < largest * rel_e) 
    return true; 
    return false; 
} 

Answer 4

Это не является прямым ответом на ваш вопрос, но альтернативное решение.

Вы можете использовать Newton's method for finding roots:

assert(x >= 0); 
if (x == 0) 
    return 0; 

double guess = x; 
for (int i=0; i<NUM_OF_ITERATIONS; i++) 
    guess -= (guess*guess-x)/(2*guess); 
return guess; 

24 итераций вы должны получить достаточно хорошее приближение, но вы также можете проверить абсолютную разницу.

Answer 5

Я предлагаю подождать, пока вы не получили стабильный ответ, а не возиться со значениями эпсилон:

double squareroot(double x) 
{ 
    double xhi = x; 
    double xlo = 0; 
    double guess = x/2; 

    while (guess * guess != x) 
    { 
     if (guess * guess > x) 
      xhi = guess; 
     else 
      xlo = guess; 

     double new_guess = (xhi + xlo)/2; 
     if (new_guess == guess) 
      break; // not getting closer 
     guess = new_guess; 
    } 
    return guess; 
}