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

Question

Я знаю, это звучит как домашнее задание, но это не так. В последнее время меня интересуют алгоритмы, используемые для выполнения определенных математических операций, таких как синус, квадратный корень и т. Д. В настоящее время я пытаюсь написать вавилонский метод вычисления квадратных корней в C #.

Пока у меня есть это:

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

Он работает просто отлично и каждый раз выдает тот же самый ответ, что и метод Math.Sqrt () .NET Framework. Как вы, вероятно, можете догадаться, однако, он медленнее, чем собственный метод (примерно на 800 тиков). Я знаю, что этот конкретный метод никогда не будет быстрее, чем собственный метод, но мне просто интересно, есть ли какие-либо оптимизации, которые я могу сделать.

Единственной оптимизацией, которую я сразу увидел, был тот факт, что вычисление будет выполнено 100 раз, даже после того, как ответ уже будет определен (в этот момент значение r всегда будет одним и тем же значением). Итак, я добавил быструю проверку, чтобы увидеть, совпадает ли вновь вычисленное значение с ранее вычисленным значением и выходит из цикла. К сожалению, это не имело большого значения для скорости, но казалось, что это правильно.

И прежде чем вы скажете: «Почему бы просто не использовать вместо этого Math.Sqrt ()?» ... Я делаю это в качестве учебного упражнения и не собираюсь фактически использовать этот метод в любом производственном коде.

Answer 1

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

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

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

Answer 2

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/

Answer 3

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

Answer 4

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

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

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

mantissa *= 2^K
exponent -= K

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

EDIT. Ниже я приведу код 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
}

Answer 5

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

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

Answer 6

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

Answer 7

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

Answer 8

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

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

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

Answer 9

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

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

Answer 10

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