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

Во-первых, вместо проверки на равенство (r == last) вы должны проверять сходимость, где r близко к last, где close определяется произвольным эпсилоном:

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

Как упоминается в статье в Википедии, на которую вы ссылаетесь, вы не эффективно вычисляете квадратные корни с помощью метода Ньютона - вместо этого вы используете логарифмы.

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/

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

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

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

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

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

Используя таблицу для получения начальных 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
}

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

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

    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;
    }

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

Первый заключался в использовании аппроксимации ряда Тейлора первого порядка по 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;
    }
}

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

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

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

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

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

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

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

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