Sqrt от "отрицательного нуля", как лучше обработать?

Песочница. Раздел для новичков
1

Всем добра!

Если кратко, иногда при вычислениях с float получается -0.0 в выражениях, которые математически не могут принимать отрицательные значения при заданных граничных условиях. Далее по алгоритму, функция sqrt() возвращает NaN от аргумента -0.0. Какие существуют изящные способы решение?

Если вдруг интересны подробности:

Спойлер

Кусок кода:

// это в объявлениях
static const float eps = 0.000001;
// радиус, за пределы которого не должны выходить конечности
float max_radius = 35;
// текущий вектор скорости
static float current_velocity_dx = 0;
static float current_velocity_dy = 0;
// координаты находящейся на поверхности конечности
static float landed_limb_x = 0;
static float landed_limb_y = 0;
// координаты поднятой конечности
static float floating_limb_x = 0;
static float floating_limb_y = 0;

// ЭТО В ЦИКЛЕ (часть кода)

float floating_velocity_dx, floating_velocity_dy;
// текущая скорость
float current_velocity_value = hypot(current_velocity_dx, current_velocity_dy);

if (current_velocity_value > eps) {
  // определение целевой позиции поднятой конечности - на краю окружности в направлении current_velocity
  float floating_target_x = (current_velocity_dx / current_velocity_value) * max_radius;
  float floating_target_y = (current_velocity_dy / current_velocity_value) * max_radius;
  // текущее расстояние от центра до позиции приземленной конечности
  float center_to_landed_limb_distance = hypot(landed_limb_x, landed_limb_y);
  float landed_to_out_distance;
  if (center_to_landed_limb_distance > eps) {
    // конечность движется в направлении, обратном движению - вычисляем косинус угла между векторами landed_limb и -current_velocity
    float cosinus = (landed_limb_x * current_velocity_dx + landed_limb_y * current_velocity_dy) / (center_to_landed_limb_distance * current_velocity_value);
    // расстояние - сторона треугольника напротив угла с найденным косинусом, к которому прилегают стороны center_to_landed_limb_distance и max_radius
    
    float a_trouble = center_to_landed_limb_distance * center_to_landed_limb_distance + max_radius * max_radius + 2 * cosinus * center_to_landed_limb_distance * max_radius;
    landed_to_out_distance = sqrt(a_trouble);

  } else landed_to_out_distance = max_radius;
  // определение вектора из текущего положения поднятой конечности в ее целевую позицию
  float to_floating_target_dx = floating_target_x - floating_limb_x;
  float to_floating_target_dy = floating_target_y - floating_limb_y;
  float to_floating_target_dl = hypot(to_floating_target_dx, to_floating_target_dy);
  if (to_floating_target_dl > eps) {
    // определение значения скорости поднятой конечности, необходимой чтобы успеть переместиться в целевую позицию
    float floating_velocity_value = to_floating_target_dl * current_velocity_value / landed_to_out_distance;
    // определение вектора скорости поднятой конечности
    floating_velocity_dx = (to_floating_target_dx / to_floating_target_dl) * floating_velocity_value;
    floating_velocity_dy = (to_floating_target_dy / to_floating_target_dl) * floating_velocity_value;
  } else {
    // если to_floating_target_dl слишком мало - просто перемещаемся в to_floating_target
    floating_velocity_dx = to_floating_target_dx;
    floating_velocity_dy = to_floating_target_dy;
  }

} else {
  floating_velocity_dx = 0;
  floating_velocity_dy = 0;
}

floating_limb_x += floating_velocity_dx;
floating_limb_y += floating_velocity_dy;

Есть точка A (landed_limb_x, landed_limb_y), движется по вектору V (current_velocity_dx, current_velocity_dy) внутри окружности с центром в начале координат и радиусом R max_radius. Необходимо найти расстояние от точки A, до пересечения промой, проходящей через нее с окружностью, в направлении V.

image
image677×604 28.4 KB

(извиняюсь за уродскую зарисовку)
Решаю через треугольник со сторонами R, OA и углом между векторами OA и V. В общем, иногда в строке 33 a_trouble приобретает значение -0.0 и геометрия “падает”.

Я понимаю, что можно как минимум сравнивать a_trouble > 0, или вовсе брать sqrt(a_trouble + eps), но как красиво сделать?))

2

А можно код сделать коротким (в несколько строк), но полным? Чтобы я мог у себя его запустить и он заработал и показал Вашу проблему и, главное, чтобы мне было что поменять для показа решения?

3

я не знаю, как сделать имитацию вычисления, которое даст в результате -0.0

void setup() {

Serial.begin(9600);

}

void loop() {

  float negative_zero = ???
  float value = sqrt(negative_zero);

  Serial.println(value);

}

если подскажете, что вписать вместо ???
а то в моем коде до вышеупомянутого места идет прием данных с пульта, вычисление setpoint вектора скорости, перемещение текущего вектора скорости в setpoint с заданным ускорением, и только потом расчеты под спойлером. А ошибка с NaN может 10 минут не появляться.

4

хотя, наверное, вот так?

void setup() {

Serial.begin(9600);

}

void loop() {

  float negative_zero = (float)-1/0;
  float value = sqrt(negative_zero);

  Serial.println(value);

}
5

Serial.println(value, 6);

6

да тоже не совсем то, там бесконечность получается, а не -0.0

7

Как вы фиксируете, что результат вычисления равен -0.0?

8

все числа, по модулю меньше эпсилон, считай нулём:
if (abs(x) < Epsilon) x = 0;
Времени это сравнение не займет особо. Вызывай только перед корнем. Или корень себе перепиши с этим учетом… или вообще его перепиши, он не самый быстрый. Кто ж знает, что именно тебе надо?

1 лайк
9

вывод на экран пульта, и да, там два знака всего, но какая разница? я положительные числа перемножаю и делю, единственное место, где может получиться минус:

float cosinus = (landed_limb_x * current_velocity_dx + landed_limb_y * current_velocity_dy) / (center_to_landed_limb_distance * current_velocity_value);
if (isnan(cosinus) && err == "") err = "err 7.1";
/ расстояние - сторона треугольника напротив угла с найденным косинусом, к которому прилегают стороны center_to_landed_limb_distance и max_radius
float a_trouble = center_to_landed_limb_distance * center_to_landed_limb_distance + max_radius * max_radius + 2 * cosinus * center_to_landed_limb_distance * max_radius;
landed_to_out_distance = sqrt(a_trouble + eps);
if (isnan(landed_to_out_distance) && err == "") err = "e72 " + String(cosinus) + " " + String(center_to_landed_limb_distance) + " " + String(a_trouble);

это если cosinus получился отрицательным, и 2 * cosinus * center_to_landed_limb_distance * max_radius получилось по абсолютному значению больше, чем center_to_landed_limb_distance * center_to_landed_limb_distance + max_radius * max_radius. То есть сумма двух квадратов меньше, чем 2 * на эти же числа не в квадратах и умноженное на косинус, который не больше единицы (или чуть больше).

Там ноль-то можно получить, только если center_to_landed_limb_distance == max_radius (оба положительные причем). Т.е. вероятно, дело именно в погрешностях при вычислении float

10

спасиб, так и понял, что подобное будет решением… кстати перед корнем можно без abs(), наверное

11

цитата просто, чтоб сослаться.
Приведенный код настолько далек от оптимизации, что просто добавь там модуль и не парься вообще. Если тебе хватает скорости в таком случае, то даже и голову ерундой не забивай.

12

Оцените исполнительный механизм для которого делается такой точный расчет. Шаг сервопривода, длина конечности - все это отражается на минимальном расстоянии перемещения. Если это не общая библиотека пригодная и для детской игрушки и для хирургического робота.
Возможно, на некоторых шагах вычислений достаточно добавить небольшую погрешность (да, тот eps), которая напрочь убирает проблемы связаные с малыми величинами, но вообще не влияет на воплощенное в реальность перемещение.
Осторожно с вычислениями, в которых погрешность может набегать (и в прямых-обратных преобразованиях). А для прямых расчетов это может сократить несколько проверок внутри, просто добавив величину перед началом расчета.

1 лайк
13

Как же я люблю “оверскилл” советы! :wink:
ТС считает третью сторону треугольника во float, не в double, просто по теореме косинусов.
А ты, точно как в старом анекдоте: “…а если бы он вез патроны!?”

14

по всей видимости - мой вариант

240мгц, 2 ядра)) но код чутка причешется, разумеется, сейчас просто пилю и обкатываю, всплывающие проблемы пытаюсь решать

15

если бы у меня акселерометр горизонт на полетном контроллере удерживал, то парился бы. А у меня конечности просто должны от края к краю допустимой зоны перемещаться и переставляться, с учетом изменения вектора скорости. Там при каждом приближении к краю области координаты вычисляются по сути заново, основываясь на показаниях с джойстика, погрешности некогда накапливаться

16

По факту расчет перемещения ноги известной длины и с приводом от сервы с известным шагом предварительно загонятеся в табличку. Расчет ее идет при задании геометрии механизма на этапе калибровки. Это избыточно на больших мощностях, но помогает при оптимизации для слабых вычислителей.

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

17

ни х… себе

18

вторая часть “камасутры” :grinning_face_with_smiling_eyes:

19

Да, конечно, надо вписать -FLT_MIN (только не забыть включить "float.h")

Самое эффективное решение проблемы – использовать функцию copysign. Если Вы заглянете в math.h, то увидите что она состоит всего

из двух машинных команд 
/**
    The copysign() function returns \a __x but with the sign of \a __y.
    They work even if \a __x or \a __y are NaN or zero.
*/
__ATTR_CONST__ static inline double copysign (double __x, double __y)
{
    __asm__ (
	"bst	%D2, 7	\n\t"
	"bld	%D0, 7	"
	: "=r" (__x)
	: "0" (__x), "r" (__y) );
    return __x;
}

Вот так выглядит решение задачи:

#include <float.h>

void setup(void) {
	Serial.begin(9600);
	//
	float negative_zero = - FLT_MIN;
	Serial.print("negative_zero: ");
	Serial.println(negative_zero);
	// 
	float value = sqrt(copysign(negative_zero, 1));
	Serial.print("value: ");
	Serial.println(value);
}

void loop(void) {}

Печатает

negative_zero: -0.00
value: 0.00
1 лайк
20

сомневаюсь, чтобы тс имел в виду вот то что на хабре. Заголовок пугает, не стал читать, завтра)