Использование формулы Герона для вычисления квадратного корня в C

Question

Я реализовал эту функцию:

double heron(double a)
{
    double x = (a + 1) / 2;
    while (x * x - a > 0.000001) {
        x = 0.5 * (x + a / x);
    }
    return x;
}

Эта функция работает как задумано, однако я хотел бы улучшить ее. Он должен использовать и бесконечный цикл while, чтобы проверить, что-то похожее на x * x является a, a число, которое пользователь должен ввести.

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

double heron(double a)
{
    double x = (a + 1) / 2;
    while (x * x != a) {
        x = 0.5 * (x + a / x);
    }
    return x;
}

Это мой первый пост, так что если есть что-то неясное или что-то, что я должен добавить, пожалуйста, дайте мне знать.

Неудачная попытка № 2:

double heron(double a)
{
    double x = (a + 1) / 2;
    while (1) {
        if (x * x == a){
            break;
        } else {
            x = 0.5 * (x + a / x);
        }
    }
    return x;
}

Формула цапли

Answer 1

Предполагается использовать и бесконечно while цикл, чтобы проверить, если что-то похожее на x * x является a

Проблемы:

Медленная сходимость

Когда начальный x совершенно неправильно, улучшенный |x - sqrt(a)| ошибка все еще может быть только вдвое меньше. Учитывая широкий спектр double, может потребоваться сотни итераций, чтобы приблизиться.

Ссылка: формула Герона.

Для нового 1-го метода оценки: быстрый обратный квадратный корень.

перелив

x * x в x * x != a склонен к переполнению. x != a/x позволяет подобный тест без этой проблемы диапазона. Если произойдет переполнение, x может заразиться "бесконечностью" или "не числом" и не достичь конвергенции.

колебания

однажды x "близко" к sqrt(a) (с коэффициентом 2) сходимость ошибок квадратична - число битов "вправо" удваивает каждую итерацию. Это продолжается до x == a/x или из-за особенностей double математика, x будет бесконечно колебаться между двумя значениями, как и частное.

Попадание в это колебание приводит к тому, что цикл OP не прерывается

---

Соединение этого вместе с помощью испытательного жгута демонстрирует адекватную сходимость.

#include <assert.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>

double rand_finite_double(void) {
  union {
    double d;
    unsigned char uc[sizeof(double)];
  } u;
  do {
    for (unsigned i = 0; i < sizeof u.uc; i++) {
      u.uc[i] = (unsigned char) rand();
    }
  } while (!isfinite(u.d));
  return u.d;
}

double sqrt_heron(double a) {
  double x = (a + 1) / 2;
  double x_previous = -1.0;
  for (int i = 0; i < 1000; i++) {
    double quotient = a / x;
    if (x == quotient || x == x_previous) {
      if (x == quotient) {
        return x;
      }
      return ((x + x_previous) / 2);
    }
    x_previous = x;
    x = 0.5 * (x + quotient);
  }
  // As this code is (should) never be reached, the `for(i)`
  // loop "safety" net code is not needed.
  assert(0);
}

double test_heron(double xx) {
  double x0 = sqrt(xx);
  double x1 = sqrt_heron(xx);
  if (x0 != x1) {
    double delta = fabs(x1 - x0);
    double err = delta / x0;
    static double emax = 0.0;
    if (err > emax) {
      emax = err;
      printf("    %-24.17e %-24.17e %-24.17e %-24.17e\n", xx, x0, x1, err);
      fflush(stdout);
    }
  }
  return 0;
}

int main(void) {
  for (int i = 0; i < 100000000; i++) {
    test_heron(fabs(rand_finite_double()));
  }
  return 0;
}

---

улучшения

• sqrt_heron(0.0) работает.

• Измените код для лучшего начального предположения.

---

double sqrt_heron(double a) {
  if (a > 0.0 && a <= DBL_MAX) {
    // Better initial guess - halve the exponent of `a`
    // Could possible use bit inspection if `double` format known.  
    int expo;
    double significand = frexp(a, &expo);
    double x = ldexp(significand, expo / 2);

    double x_previous = -1.0;
    for (int i = 0; i < 8; i++) {  // Notice limit moved from 1000 down to < 10
      double quotient = a / x;
      if (x == quotient) {
        return x;
      }
      if (x == x_previous) {
        return (0.5 * (x + x_previous));
      }
      x_previous = x;
      x = 0.5 * (x + quotient);
    }
    assert(0);
  }
  if (a >= 0.0) return a;
  assert(0);  // invalid argument.
}