Хабр Курсы для всех
РЕКЛАМА
Практикум, Хекслет, SkyPro, авторские курсы — собрали всех и попросили скидки. Осталось выбрать!
Хабр доступен 24/7 благодаря поддержке друзей
Комментарии 73
А зачем максимум и минимум берутся от модулей чисел? В треугольнике все стороны как правило больше нуля.
Хм… хороший вопрос, о котором я не подумал пока переводил статью.
Видимо автор решил обобщить формулу, которая служит для нахождения гипотенузы на более широкий класс входных данных.
Здесь же более важно не то, как правильно находить гипотенузу, а какими принцам нужно следовать при реализации математических формул.
Видимо автор решил обобщить формулу, которая служит для нахождения гипотенузы на более широкий класс входных данных.
Здесь же более важно не то, как правильно находить гипотенузу, а какими принцам нужно следовать при реализации математических формул.
Не уверен в этом, в английской wiki по ссылке en.wikipedia.org/wiki/Hypot#Implementation это трактуется именно как знак.
Кстати там же интересно посмотреть в каких версиях каких языков именно этот способ вычисления реализован.
Кстати там же интересно посмотреть в каких версиях каких языков именно этот способ вычисления реализован.
Какой неожиданный поворот событий. Осталось только понять как это работает www.netlib.org/fdlibm/e_hypot.c
Не могу понять, почему используется
а не просто
Неужели при положительных a и b выражение a-(-b) может дать не такой результат, как (a+b)?
w = sqrt(t1*t1-(b*(-b)-t2*(a+t1)));
а не просто
w = sqrt(t1*t1+(b*b+t2*(a+t1)));
Неужели при положительных a и b выражение a-(-b) может дать не такой результат, как (a+b)?
Знак double хранится в одном выделенном бите (msb). Диапазоны положительных и отрицательных полностью сопавдают.
Не годится — для вещественных чисел используется обратный код (знаковый бит и абсолютное значение), так что представление x и -x отличается ровно одим битом, диапазоны чисел с одинаковым порядком, точностью и т.п. полностью симметричны. Возможно, что как-то по-разному работает округление, и это важно — ведь здесь они борются за ошибку не больше sqrt(1/2) младшего бита!
nan вообще нечто странное и nan+(-nan) это не ноль вроде как.
простой пример на коленке:
вывод:
Собственно теже грабли как и с целыми. MIN_INT != — MAX_INT. а если к ним применить минус и сложить, то ноль как раз получится.
простой пример на коленке:
#include <iostream>
#include <stdint.h>
int main()
{
uint64_t a_pos= 0x7FFFFFFFFFFFFFFFUL;
uint64_t a_neg= 0xFFFFFFFFFFFFFFFFUL;
double *a, *b;
a = (double*)&a_pos;
b = (double*)&a_neg;
std::cout << *a << " " << *b << " " << *a+*b << std::endl;
a_pos = - a_pos;
a_neg = - a_neg;
std::cout << *a << " " << *b << " " << *a+*b << std::endl;
a_pos = - a_pos;
a_neg = - a_neg;
std::cout << *a << " " << *b << " " << *a+*b << std::endl;
return 0;
}
вывод:
$ ./a.out
nan -nan nan
-4.94066e-324 4.94066e-324 0
nan -nan nan
$
Собственно теже грабли как и с целыми. MIN_INT != — MAX_INT. а если к ним применить минус и сложить, то ноль как раз получится.
Если мы рисуем диагональ на экране (стрелочку из точки X,Y в точку X',Y'), то очень удобно считать разницу расстояний между координатами как X'-X и Y'-Y. Если стрелочка идёт «вправо, вниз», стороны треугольника будут положительными, а если влево-вверх, то обе стороны будут отрицательными. Можно писать функцию, которая это «прожует», а можно передавать в функцию модули. Выбор за вами.
«есть большая вероятность, что сделано это не просто так и есть какие-то глубокие причины для такого решения» :-)
Вычислительная геометрия зачастую оперирует векторами, а не отрезками, потому длина может быть условно отрицательной. Это позволяет написать что-то вроде hypot(x1-x2, y1-y2) и не парится с направлением вектора.
Из справки std::hypot:
Т.е. ф-ия чуть более общая, чем указал автор.
This is the length of the hypotenuse of a right-angled triangle with sides of length x and y, or the distance of the point (x,y) from the origin (0,0), or the magnitude of a complex number x+iy
Т.е. ф-ия чуть более общая, чем указал автор.
Не знаю, ответили ли уже в комментариях, но как x, так и y могут быть отрицательными.
Вдруг берется система координат и треугольник отражен вертикально и/или горизонтально?
Тогда х = 0 + (-5) = -5, а у = 0 + (-6) = -6.
Вдруг берется система координат и треугольник отражен вертикально и/или горизонтально?
Тогда х = 0 + (-5) = -5, а у = 0 + (-6) = -6.
Если они имеют знак, то какой знак должен быть у гипотенузы?
Но мы же вычисляем длину гипотенузы :) она положительна.
Если брать формулу sqrt(x^2+y^2) — знаки не важны.
Но если идти по статье, то в случае отрицательных чисел
Даст иной результат, нежели
Если брать формулу sqrt(x^2+y^2) — знаки не важны.
Но если идти по статье, то в случае отрицательных чисел
max = maximum(|x|, |y|);
min = minimum(|x|, |y|);
r = min / max;
return max*sqrt(1 + r*r);
Даст иной результат, нежели
max = maximum(x, y);
min = minimum(x, y);
r = min / max;
return max*sqrt(1 + r*r);
Как-то сравнивал расстояния между точками, потом понял, что собственно корень извлекать мне не нужно, сойдут и квадраты расстояний для сравнения. Это я к тому, что иногда не нужно вычислять гипотенузу.
P.S. Случай автора натолкнул на мысль. Ничего против его метода не имею.
P.S. Случай автора натолкнул на мысль. Ничего против его метода не имею.
Если одна из сторон
0, то алгоритм работает правильно, но если обе стороны по нулям, то алгоритм должен вернуть 0, разве нет? Или предполагается что это очевидная вещь и проверки должны написать сами?Если мы используем функцию hypot(x,y) для определения расстояния между точками (а это случается в сотни раз чаще, чем вычисление «гипотенузы треугольника»), то 0 (и даже два нуля) на входе — это вполне корректные данные. И мы вправе ожидать от алгоритма правильного результата.
Вот что должен выдавать atan2(0,0), не знаю. Но что не exception, это точно.
Вот что должен выдавать atan2(0,0), не знаю. Но что не exception, это точно.
Замечательно. И как же называется функция, вычисляющая длину вектора в двумерной евклидовой плоскости?
Про atan2 в MSDN написано: «If both parameters of atan2 are 0, the function returns 0. » Так что Ваше утверждение, что «Для atan2 же заявлены одновременно ненулевые аргументы, поэтому вариант atan2(0,0) по умолчанию неправомерен.» — не вполне достоверно :(
Про hypot они пишут: «The _hypot function calculates the length of the hypotenuse of a right triangle, given the length of the two sides x and y (in other words, the square root of x^2 + y^2). » Вы хотите сказать, что это описание является внутренне противоречивым? Про x>0, y>0 нет ни слова.
А «простой if», если писать его перед каждым вызовом функции — это очень дорого, учитыая, что почти всегда вызов будет происходить. Лучше уж сама функция проверит аргументы наиболее эффективным способом.
Про atan2 в MSDN написано: «If both parameters of atan2 are 0, the function returns 0. » Так что Ваше утверждение, что «Для atan2 же заявлены одновременно ненулевые аргументы, поэтому вариант atan2(0,0) по умолчанию неправомерен.» — не вполне достоверно :(
Про hypot они пишут: «The _hypot function calculates the length of the hypotenuse of a right triangle, given the length of the two sides x and y (in other words, the square root of x^2 + y^2). » Вы хотите сказать, что это описание является внутренне противоречивым? Про x>0, y>0 нет ни слова.
А «простой if», если писать его перед каждым вызовом функции — это очень дорого, учитыая, что почти всегда вызов будет происходить. Лучше уж сама функция проверит аргументы наиболее эффективным способом.
Да, в системах символьных вычислений вполне можно вернуть atan2(0,0) как Undefined. Впрочем, для них atan2(1,1) может оказаться выражением Pi/4, а atan2(1,-2) — величиной Pi/2-atan(2), так что им можно. Хорошо, что в численных библиотеках это не так.
А «обернуть if в новую функцию»… мне этих длин векторов и арктангенсов приходится вычислять по несколько миллионов в секунду в реальном времени, так там каждая проверка с переходом на счету, не то, что вызовы функций. Функция, конечно, заинлайнится, но лишних проверок лучше избежать.
А «обернуть if в новую функцию»… мне этих длин векторов и арктангенсов приходится вычислять по несколько миллионов в секунду в реальном времени, так там каждая проверка с переходом на счету, не то, что вызовы функций. Функция, конечно, заинлайнится, но лишних проверок лучше избежать.
— The _hypot function calculates the length of the hypotenuse of a right triangle, given the length of the two sides x and y (in other words, the square root of x^2 + y^2).У меня в
— Треугольник это термин. Треугольник это геометрическая фигура, образованная тремя отрезками, которые соединяют три не лежащие на одной прямой точки. Из этого определение следует, что нет ситуации, когда у треугольника хотя бы одна сторона имеет нулевую длину, а когда такая ситуация наступает, то треугольник перестает быть таковым, даже сам смысл термина гипотенуза в таком случае теряется.
man hypot написано прямо в заголовке: «hypot, hypotf, hypotl - Euclidean distance function». И ниже:The hypot() function returns sqrt(x*x+y*y). This is the length of the hypotenuse of a right-angled triangle with sides of length x and y, or the distance of the point (x,y) from the origin.Только в man это продублировали в описании, а в стандарте написано
hypot, hypotf, hypotl — Euclidean distance functionи
These functions shall compute the value of the square root of x²+y² without undue overflow or underflow., но ниже только
Upon successful completion, these functions shall return the length of the hypotenuse of a right-angled triangle with sides of length x and y.
В RSA шифровании тоже есть такая проблема
Дело в том, что в RSA значения N и P обычно очень большие (~200 цифр), так что N^P будет состоять уже из 40 000 цифр, в то время как итоговый результат будет состоять из стольки цифр, из скольки состоит M (за счет
N ^ P mod M. В криптографические библиотеки для этого добавляют специально оптимизированную функцию modexp.Дело в том, что в RSA значения N и P обычно очень большие (~200 цифр), так что N^P будет состоять уже из 40 000 цифр, в то время как итоговый результат будет состоять из стольки цифр, из скольки состоит M (за счет
mod M). Оптимизация позволяет этого избежать.
С этим справляется Быстрое возведение в степень, если результат каждого умножения брать mod M.
Код
__int64 FastPow (__int64 a, __int64 x, __int64 p)
{
__int64 res = 1, s = a;
while (x) {
if (x % 2 == 1)
res = res * s % p;
s = s * s % p;
x /= 2;
}
return res;
}
прошу прощения, я дэбил, поторопился :)
Кто-нибудь может объяснить, почему стандартный
Кстати, вариант из статьи замедляет вычисления раза примерно в 2.
hypot(x, y) в 10-20 раз (тестировал на Visual Studio C++ 2010 и MinGW 4.7.2) медленнее наивного рукописного варианта?inline double naive_hypot(double x, double y)
{
return sqrt(x * x + y * y);
}
Кстати, вариант из статьи замедляет вычисления раза примерно в 2.
Выше в комментариях mrTuborg приведена ссылка на исходник www.netlib.org/fdlibm/e_hypot.c
Вариант из статьи как-то плохо себя поведет около нуля. Вероятно, в окрестности нуля надо переходить к вычислению простым способом:
sqrt(x * x + y * y).Меня тоже смущает это место:
Если min — очень маленькое число, а max — очень большое число, а потом ещё это и в квадрат и прибавить к единице…
Не будет ли здесь потеря точности?
С другой стороны, здесь
тоже получается в этом случае, что мы очень большое складываем с очень маленьким и где будет больше потеря точности с ходу непонятно.
Нужно подумать над конкретными числовыми примерами.
r = min / max;Если min — очень маленькое число, а max — очень большое число, а потом ещё это и в квадрат и прибавить к единице…
Не будет ли здесь потеря точности?
С другой стороны, здесь
sqrt(x * x + y * y); тоже получается в этом случае, что мы очень большое складываем с очень маленьким и где будет больше потеря точности с ходу непонятно.
Нужно подумать над конкретными числовыми примерами.
тоже получается в этом случае, что мы очень большое складываем с очень маленьким и где будет больше потеря точности с ходу непонятно.
x=1, y=123123123123sqrt(x * x + y * y);
здесь потери точности нету
r = min / max;
здесь может быть
У вас 53 бита мантиссы при делении.
В любом варианте сложение в формуле есть, различие порядков слагаемых одинаково. Либо вы складываете a2 и b2, либо 1 и a2/b2:
(комментарий к последнему равенству: любой логарифм 1 равен нулю)
Так что как вы ни вычисляйте, потери точности произойдут одинаковые.
В любом варианте сложение в формуле есть, различие порядков слагаемых одинаково. Либо вы складываете a2 и b2, либо 1 и a2/b2:
(комментарий к последнему равенству: любой логарифм 1 равен нулю)
Так что как вы ни вычисляйте, потери точности произойдут одинаковые.
Я вот на JS проверил:
function hypo1 (x,y) {
return Math.sqrt(x*x+y*y);
}
function hypo2 (x,y) {
var i, a, r;
x = Math.abs(x);
y = Math.abs(y);
a = Math.max(x,y);
i = Math.min(x,y);
r = i / a;
return a*Math.sqrt(1 + r*r);
}
var
x = 100000013,
y = 13;
console.log( hypo1(x,y) ); // 100000013.00000083
console.log( hypo2(x,y) ); // 100000013.00000085
var
x = 221412142513,
y = 1242353;
console.log( hypo1(x,y) ); // 221412142516.48547
console.log( hypo2(x,y) ); // 221412142516.48544
Задача на сообразительность: в каких случаях вычисление гипотенузы сводится к формуле
Не видевшие реализации, ни в жизнь не догадаются.
Версия для float: code.metager.de/source/xref/glibc/sysdeps/ieee754/flt-32/e_hypotf.c
(для справки: hb &= 0x7fffffff — взятие мантиссы и экспоненты без учёта знака, ha == 0x7f800000 — проверка на бесконечность)
Ссылку на версию для double дали выше, но это уже совсем другая история.
|x| * |y|?Не видевшие реализации, ни в жизнь не догадаются.
Версия для float: code.metager.de/source/xref/glibc/sysdeps/ieee754/flt-32/e_hypotf.c
(для справки: hb &= 0x7fffffff — взятие мантиссы и экспоненты без учёта знака, ha == 0x7f800000 — проверка на бесконечность)
Ссылку на версию для double дали выше, но это уже совсем другая история.
Когда одно из чисел x,y равно бесконечности или NAN?
А точнее, именно NaN, но не бесконечность (ha == 0x7f800000 — бесконечность, ha > 0x7f800000 — NaN). А делается это, вероятно, для того, чтобы корректно обработать сигнальность NaN. Если я ошибаюсь, прошу поправить.
Вот еще хороший вопрос — а как найти гипотенузу если нет функции sqrt. Например в если нет возможности подключить math lib для sqlite а гипотенуза нужна? в наличии только +,-,*/ :)?
И log и exp нет?
Можно извлечь корень за 3 деления (предполагается sizeof(int)=4):
Гипотенузу вычислять как обычно.
double FSqrt(double x){
double y=x;
((int*)&y)[1]=(((int*)&x)[1]>>1)+0x1FF76CF6;
y=(y+x/y)*0.5;
y=(y+x/y)*0.5;
y=(y+x/y)*0.5;
return y;
}
Гипотенузу вычислять как обычно.
0x1FF76CF6 — это специальная магическая константа для получения хорошего первого приближения?
Да, подобрана экспериментально.
Вот тут есть подробное описание, как получается аналогичная константа для вычисления обратной функции.
Плюс, ей посвящен известный комикс.
Плюс, ей посвящен известный комикс.
Если считать sqrt(x) как x*rsqrt(x), то можно обойтись без делений вообще (но за 16 умножений):
Случаи INF и NAN, а также result>MAX_DOUBLE не обрабатываются.
double Hypo(double a,double b){
int ca=((int*)&a)[1]&0x7ff00000,cb=((int*)&b)[1]&0x7ff00000;
const int maxdif=0x1e00000; // 2^-30
if(ca-cb>maxdif) return fabs(a);
if(cb-ca>maxdif) return fabs(b);
if(ca<cb) ca=cb;
ca-=0x40000000;
((int*)&a)[1]-=ca; ((int*)&b)[1]-=ca;
double x=a*a+b*b;
double x2=x*0.5,y=0;
((int*)&y)[1]=0x5fe6eb50-(((int*)&x)[1]>>1);
y*=1.5-x2*y*y;
y*=1.5-x2*y*y;
y*=1.5-x2*y*y;
y*=1.5-x2*y*y;
y*=x;
((int*)&y)[1]+=ca;
return y;
}
Случаи INF и NAN, а также result>MAX_DOUBLE не обрабатываются.
Интересно, как быстро вычислить хотя бы 16 бит гипотенузы, если деления тоже нет? x и y целые (не превосходят по модулю 2^30.5).
А то же самое, но для катета от гипотенузы и второго катета?)
Зарегистрируйтесь на Хабре, чтобы оставить комментарий
Что сложного может быть в вычислении гипотенузы?