Как я могу улучшить этот метод квадратного корня?

public static double SquareRoot(double x) {
    if (x == 0) return 0;

    double r = x / 2; // this is inefficient, but I can't find a better way
                      // to get a close estimate for the starting value of r
    double last = 0;
    int maxIters = 100;

    for (int i = 0; i < maxIters; i++) {
        r = (r + x / r) / 2;
        if (r == last)
            break;
        last = r;
    }

    return r;
}

eps = 1e-10  // pick any small number
if (Math.Abs(r-last) < eps) break;

Ответ 2 · рейтинг 5 · Chris H

float InvSqrt (float x){
  float xhalf = 0.5f*x;
  int i = *(int*)&x;
  i = 0x5f3759df - (i>>1);
  x = *(float*)&i;
  x = x*(1.5f - xhalf*x*x);
  return x;}

Это мой любимый быстрый квадратный корень. На самом деле это обратный квадратный корень, но вы можете инвертировать его после, если хотите .... Я не могу сказать, будет ли это быстрее, если вам нужен квадратный корень, а не обратный квадратный корень, но это чертовски круто все равно .
http://www.beyond3d.com/content/articles/8/

Ответ 3 · рейтинг 4 · Jasiu

Здесь вы выполняете метод Ньютона для поиска корня. Таким образом, вы можете просто использовать более эффективный алгоритм поиска корня. Вы можете начать поиск здесь.

Ответ 4 · рейтинг 2 · Richard

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

У вас гораздо больше шансов получить улучшение, выйдя из цикла раньше, поэтому либо сохраните новый r в переменной и сравните со старым r, либо сохраните x/r в переменной и сравните его с r, прежде чем выполнять сложение и деление. .

Ответ 5 · рейтинг 2 · AlbertoPL

Вместо разрыва цикла и последующего возврата r вы можете просто вернуть r. Может не обеспечить какого-либо заметного увеличения производительности.

Ответ 6 · рейтинг 2 · Eric Bainville

С вашим методом каждая итерация удваивает количество правильных битов.

Используя таблицу для получения начальных 4 битов (например), у вас будет 8 бит после 1-й итерации, затем 16 бит после второй и все биты, которые вам нужны, после четвертой итерации (поскольку double хранит 52+1 бит мантиссы).

Для поиска в таблице вы можете извлечь мантисса в [0,5,1[ и показатель степени из входных данных (используя функцию, например frexp), а затем нормализовать мантисса в [64,256[, используя умножение на подходящую степень 2.

mantissa *= 2^K
exponent -= K

После этого ваш входной номер по-прежнему будет mantissa*2^exponent. K должен быть равен 7 или 8, чтобы получить четный показатель степени. Вы можете получить начальное значение для итераций из таблицы, содержащей все квадратные корни целой части мантиссы. Выполните 4 итерации, чтобы получить квадратный корень r из мантиссы. Результатом является r*2^(exponent/2), построенный с использованием функции, подобной ldexp.

РЕДАКТИРОВАТЬ. Ниже я привожу код C++, чтобы проиллюстрировать это. Функция OP sr1 с улучшенным тестом занимает 2,78 с для вычисления 2 ^ 24 квадратных корней; моя функция sr2 занимает 1,42 с, а аппаратный sqrt - 0,12 с.

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

double sr1(double x)
{
  double last = 0;
  double r = x * 0.5;
  int maxIters = 100;
  for (int i = 0; i < maxIters; i++) {
    r = (r + x / r) / 2;
    if ( fabs(r - last) < 1.0e-10 )
      break;
    last = r;
  }
  return r;
}

double sr2(double x)
{
  // Square roots of values in 0..256 (rounded to nearest integer)
  static const int ROOTS256[] = {
    0,1,1,2,2,2,2,3,3,3,3,3,3,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,6,6,6,6,6,6,6,6,6,6,6,6,
    7,7,7,7,7,7,7,7,7,7,7,7,7,7,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,9,9,9,9,9,9,9,9,9,9,9,9,9,
    9,9,9,9,9,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,11,11,11,11,11,
    11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,12,12,12,12,12,12,12,12,12,12,12,12,
    12,12,12,12,12,12,12,12,12,12,12,12,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,
    13,13,13,13,13,13,13,13,13,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,
    14,14,14,14,14,14,14,14,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,
    15,15,15,15,15,15,15,15,15,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16 };

  // Normalize input
  int exponent;
  double mantissa = frexp(x,&exponent); // MANTISSA in [0.5,1[ unless X is 0
  if (mantissa == 0) return 0; // X is 0
  if (exponent & 1) { mantissa *= 128; exponent -= 7; } // odd exponent
  else { mantissa *= 256; exponent -= 8; } // even exponent
  // Here MANTISSA is in [64,256[

  // Initial value on 4 bits
  double root = ROOTS256[(int)floor(mantissa)];

  // Iterate
  for (int it=0;it<4;it++)
    {
      root = 0.5 * (root + mantissa / root);
    }

  // Restore exponent in result
  return ldexp(root,exponent>>1);
}

int main()
{
  // Used to generate the table
  // for (int i=0;i<=256;i++) printf(",%.0f",sqrt(i));

  double s = 0;
  int mx = 1<<24;
  // for (int i=0;i<mx;i++) s += sqrt(i); // 0.120s
  // for (int i=0;i<mx;i++) s += sr1(i);  // 2.780s
  for (int i=0;i<mx;i++) s += sr2(i);  // 1.420s
}

Ответ 7 · рейтинг 1 · Sinan Ünür

Определите допуск и вернитесь раньше, когда последующие итерации попадают в этот допуск.

Ответ 8 · рейтинг 1 · Unknown

Поскольку вы сказали, что приведенный ниже код недостаточно быстр, попробуйте следующее:

    static double guess(double n)
    {
        return Math.Pow(10, Math.Log10(n) / 2);
    }

Он должен быть очень точным и, надеюсь, быстрым.

Вот код для начальной оценки, описанной здесь. Кажется, это неплохо. Используйте этот код, а затем вы также должны выполнять итерацию, пока значения не сойдутся в пределах эпсилон разницы.

    public static double digits(double x)
    {
        double n = Math.Floor(x);
        double d;

        if (d >= 1.0)
        {
            for (d = 1; n >= 1.0; ++d)
            {
                n = n / 10;
            }
        }
        else
        {
            for (d = 1; n < 1.0; ++d)
            {
                n = n * 10;
            }
        }


        return d;
    }

    public static double guess(double x)
    {
        double output;
        double d = Program.digits(x);

        if (d % 2 == 0)
        {
            output = 6*Math.Pow(10, (d - 2) / 2);
        }
        else
        {
            output = 2*Math.Pow(10, (d - 1) / 2);
        }

        return output;
    }

Ответ 9 · рейтинг 1 · ccook

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

Первый заключался в использовании аппроксимации ряда Тейлора первого порядка по x0:

    Func<double, double> fNewton = (b) =>
    {
        // Use first order taylor expansion for initial guess
        // http://www27.wolframalpha.com/input/?i=series+expansion+x^.5
        double x0 = 1 + (b - 1) / 2;
        double xn = x0;
        do
        {
            x0 = xn;
            xn = (x0 + b / x0) / 2;
        } while (Math.Abs(xn - x0) > Double.Epsilon);
        return xn;
    };

Во-вторых, попробовать третий заказ (более дорогой), повторить

    Func<double, double> fNewtonThird = (b) =>
    {
        double x0 = b/2;
        double xn = x0;
        do
        {
            x0 = xn;
            xn = (x0*(x0*x0+3*b))/(3*x0*x0+b);
        } while (Math.Abs(xn - x0) > Double.Epsilon);
        return xn;
    };

Я создал вспомогательный метод для определения времени функций

public static class Helper
{
    public static long Time(
        this Func<double, double> f,
        double testValue)
    {
        int imax = 120000;
        double avg = 0.0;
        Stopwatch st = new Stopwatch();
        for (int i = 0; i < imax; i++)
        {
            // note the timing is strictly on the function
            st.Start();
            var t = f(testValue);
            st.Stop();
            avg = (avg * i + t) / (i + 1);
        }
        Console.WriteLine("Average Val: {0}",avg);
        return st.ElapsedTicks/imax;
    }
}

Оригинальный метод был быстрее, но опять же, может быть интересным :)

Ответ 10 · рейтинг 1 · Kristof Verbiest

Замена "/ 2" на "* 0,5" делает это примерно в 1,5 раза быстрее на моей машине, но, конечно, не так быстро, как нативная реализация.

Ответ 11 · Chad Birch

Что ж, нативная функция Sqrt(), вероятно, не реализована на C#, она, скорее всего, будет реализована на низкоуровневом языке и наверняка будет использовать более эффективный алгоритм. Поэтому пытаться соответствовать его скорости, вероятно, бесполезно.

Однако, что касается просто попытки оптимизировать вашу функцию для heckuvit, на странице Википедии, на которую вы ссылаетесь, рекомендуется, чтобы "начальное предположение" было 2 ^ этаж (D/2), где D представляет количество двоичных цифр в числе. Вы могли бы попробовать, я не вижу ничего другого, что можно было бы значительно оптимизировать в вашем коде.

Ответ 12 · рейтинг -2 · Noam Gal

Вы можете попробовать r = x >> 1;

вместо / 2 (также в другом месте у вас устройство на 2). Это может дать вам небольшое преимущество. Я бы также переместил 100 в цикл. Наверное, ничего, но мы говорим здесь о клещах.

просто сейчас проверю.

РЕДАКТИРОВАТЬ: исправлено > в >>, но это не работает для двойников, так что неважно. встраивание 100 не дало мне увеличения скорости.