Войти
ПрограммированиеФорумОбщее

быстрый квадратный корень

Страницы: 1 2 Следующая »
#0
23:13, 14 янв. 2006

Есть функция которая с низкой точностью вычисляет квадратный корень:

double fsqrt(double n)
{
	double f;
	unsigned long *i;

	f = n;
	i = ((unsigned long*)&f) + 1;
	*i = (0xbfcd4600 - *i) >> 1;
	n *= f;

	return (n * (1.5 - 0.5 * f * n));
}

Малость не понятно выражение i = ((unsigned long*)&f) + 1;
Получаем целую чась n, но почему прибавляем 1?

Как преобразовать функцию для float вместо double?


#1
23:23, 14 янв. 2006

Проблема решена :)

float fsqrtf(float n)
{
	float f;
	unsigned long *i;

	f = n;
	i = ((unsigned long*)&f);
	*i = (0xbe6f0000 - *i) >> 1;
	n *= f;

	return (n * (1.5 - 0.5 * f * n));
}

Может у кого есть идеи как еще ускорить функцию?

#2
23:24, 14 янв. 2006

прикольно добавлять 1 к указателю )
а это вообще работает?

#3
23:25, 14 янв. 2006

Вообще Ламот на мой взгляд написал самую быструю функцию извлечения квадратного корня

float fast_sqrt(float f)
{
    float result;
    _asm
    {
        mov eax, f
        sub eax, 0x3f800000
        sar eax, 1
        add eax, 0x3f800000
        mov result, eax
    }

    return result;
}

#4
23:29, 14 янв. 2006

Это формула разложения в степенной ряд, причем взято только 2 первых членов. Еденица там возникла из-за того, что в качествен начального приближения x0=1 принята еденица (это сделано потому, что из нее удобно вычислять квадратный корень, ну и любую степень тоже).
С тремя членами эта формула будет иметь вид

x=1+(x-1)/2-(x-1)^2/4

#5
23:39, 14 янв. 2006

2 DissDoc
точность по сравнению с моим вариантом не очень
случайно нет такой же только для нахождения 1.0f/sqrt(x)?

#6
23:42, 14 янв. 2006

bmax
Бери формулу Тейлора и вычисляй для любой желательной функции. Нет проблем.

#7
23:46, 14 янв. 2006

2 Федор
странно, но в варианте с float, эта единица дает ошибочный результат

#8
23:53, 14 янв. 2006

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

#9
0:00, 15 янв. 2006

bmax
Чего странного, double - 8 байт, float - 4 байта, unsigned long тут тоже 4 байта. +1 просто перемещает указатель на unsigned long во вторую половину double - туда, где  экспонента.

#10
0:19, 15 янв. 2006

Но я не понял твою функцию. Там i не попадает на выход и на самом деле результат это

return (n*n * (1.5 - 0.5 * n *n* n));

Поскольку у тебя
f = n;
n *= f;

Я не исключаю что (n*n * (1.5 - 0.5 * n *n* n)); может давать для некоторого дипазона правильный результат, но это надо еще выяснить.

Сейчас выяснил. Наиболее точен результат при x в районе 1, но очень не точно.
Видемо, здесь что-то не так.

#11
0:21, 15 янв. 2006

Короче бери мою формулу, а еще лучше функцию DissDoc'а.

#12
1:27, 15 янв. 2006

Кому интересна тема, рекомендую почитать тут.

#13
2:44, 15 янв. 2006

bmax
>случайно нет такой же только для нахождения 1.0f/sqrt(x)?

#define ONE_AS_INTEGER ((DWORD)(0x3F800000))
float __fastcall InvQSqrt(const float & x)
{
  DWORD   tmp = ((ONE_AS_INTEGER << 1) + ONE_AS_INTEGER - *(DWORD*)&x) >> 1;
  float y = *(float*)&tmp;
  return y * (1.47f - 0.47f * x * y * y);
}
#14
11:51, 15 янв. 2006

После некоторых экспериментов получились такие функции:

typedef float real;
typedef unsigned long uint;

real math_sqrt(const real x)
{
	real r = x;
	uint *i = ((uint*)&r);

	*i = (0xbe6f0000 - *i) >> 1;

	return (x * r * (1.5f - 0.5f * x * r * r));
}

real math_isqrt(const real x)
{
	real r = x;
	uint *i = (uint*)&r;

	*i = 0x5f375a86 - (*i >> 1);

	return (r * (1.5f - 0.5f * x * r * r));
}

Вычисляют sqrt, 1/sqrt с точностью 2 значений после запятой

Страницы: 1 2 Следующая »
ПрограммированиеФорумОбщее

Тема в архиве.