Эпический PosInRange без вещь-чисел и переполнения

mikser · 03.03.2012, 18:40

Даны две точки нужно проверить превышает ли расстояние между ними наше число Radius. Известно что Radius не отрицательно.
Или иными словами попадает ли наша точка в сферу с радиусом Radius или не попадает. То есть на выходе False/True;
Нужно обойтись без вещественных чисел, тоесть всякие float, extendet, real, double и тд запрещено использовать.

Delphi Код:

Type
  TPoint= record 
   x,y,z:integer;
  end;
function isInSphere(Center:TPoint;Radius:integer;P:TPoint):Boolean;
begin
 ...
{result:= Истина если точка внутри сферы иначе Ложь}
end;

Нужно запилить такую функцию что бы не было переполнения, так как координаты x,y,z обещают быть большими. И её потом можно было юзать в боте для ладвы

PS Да кстати простое использование int64 не поможет, даже его может не хватить.

PSS Думаю надо упрощать формулу SQRT ( (x2-x1)^2 + (y2-y1)^2 + (z2-z1)^2 ) < Radius

Добавлено через 1 час 24 минуты
Вот мой черновой набросок

Код:

Type
  TPoint= record 
   x,y,z:integer;
  end;
function isInSphere(Center:TPoint;R:int64;P:TPoint):Boolean;
var a,b,c, AB, Left, Right:int64;
begin
	with Center  do
	begin
		a:= int64(P.x) - int64(x);
		b:= int64(P.y) - int64(y);
		c:= int64(P.z) - int64(z);
	end;//with Center
	{
		SQRT(a*a + b*b + c*c) < R
		Упрощаем:
		//^2 - возведение в квадрат
		слева и справа положительные числа => при возведении в квадрат обоих частей знак не поменяется
		a^2 + b^2 + c^2 < R^2 	
		прибавим к обоим частям +2ab //2 * a * b
		(a^2 +2ab + b^2) + c^2 < R^2 +2ab
		да это же квадрат суммы!!!
		(a+b)^2 + c^2 < R^2 + 2ab
		добавим к обоим частям 2*с(a+b) // Почему именно это, а вот так захотелось
		((a+b)^2 + 2(a+b)c + c^2) < R^2 + 2ab + 2(a+b)c
		омг да это опять квадрат суммы
		((a+b) + c)^2  <  R^2 + 2ab + 2(a+b)c
		Делим обе части на R^2 | знак неравнства не пострадает так как R^2  неотрицателен 
		(a+b+c)^2 / R^2 < 1 + (2ab+2(a+b)c) /R^2
		в левой части выносим квадрат за скобки
		( (a+b+c) / R )^2 < 1 + (2ab+2(a+b)c) /R^2
		вычтем 1^2 из обоих частей //1^2 = 1
		( (a+b+c) / R )^2 - 1^2   <   (2ab+2(a+b)c) /R^2
		Вспоминаем формулу разность квадратов  p^2 - t^2 = (p-t)*(p+t)
		(( (a+b+c) / R ) - 1) * (( (a+b+c) / R ) + 1)   <   (2ab+2(a+b)c) /R^2
		Заталкиваем еденицы в чсилитель умножив их на R - (знаменатель)
		( (a+b+c - R) / R  ) * ( (a+b+c + R ) / R )     <   (2ab+2(a+b)c) /R^2
		что бы умножить две дроби нужно просто перемножить их верх и низ :)
		круто что в знаменателе получится R^2
		 (a+b+c - R) * (a+b+c + R ) / R^2      <   (2ab+2(a+b)c) /R^2
		 Делим обе части на R^2 так как квадрат  не отрицателен всегда то знак не пострадает
		 Молимся что бы R не оказался нулем и делим.
		 (a+b+c - R) * (a+b+c + R )      <   (2ab+2(a+b)c)
		 Выражение упростилось :)
		 Теперь начинается гиморой
		 делим обе части на a*b Но если А* B отрицателен то 
		 знак неравенства поменяется вот мляя :(
		 Придется это отслеживать
			
			Пусть ab>0
				(a+b+c - R) * (a+b+c + R )  /ab     <   (2ab+2(a+b)c) /ab
				правую часть  упростим как можем
				(a+b+c - R) * (a+b+c + R )  /ab     <   2 + (2(a+b)c)/ab
				левую скобку разделим на А а правую на B -
				((a+b+c - R) /a) * ((a+b+c + R )/b )   <   2 + (2(a+b)c)/ab
				сокращаем в левой части что можем
				(1+(b+c-R)/a) * ( 1 +(a+c+R)/b )   <   2 + (2(a+b)c)/ab
				А вот дальше я не могу упростить :(
 	 
			Пусть ab<0  будет тоже самое что и  в предыдущем случа однако знак неравнества сменится
				(1+(b+c-R)/a) * ( 1 +(a+c+R)/b )   >   2 + (2(a+b)c)/ab
		 
			Пусть аb=0
				(a+b+c - R) * (a+b+c + R )      <   (2ab+2(a+b)c)
				Подставляем ab=0 
				(a+b+c - R) * (a+b+c + R )      <   2(a+b)c
				ab=0 значит а или b равно нулю или они оба нули
				пусть а = 0
					(b+c - R) * (b+c + R )      <   2bc
					
				пусть b = 0
					(a+c - R) * (a+c + R )      <   2ac
				пусть a = b = 0
					(c - R) * (c + R )      <   0

	}
	AB:=a*b;
	
	if (AB=0)then //нельзя делить на AB
	begin
		result:= (a+b+c-R) * (a+b+c + R )    <   2*(a+b)*c
	end
	else
	begin
		Left:= (1 + (a+c-R) div a) * (1 + (a+c+R) div b);
		Right:= 2 * (1 + (c*(a+b) div AB));
	
		if AB > 0 then //если делили на + то знак неравнества не меняется
			result:= Left<Right
		else if (AB<0) then// делили на минус знак поменялся
			result:= Left>=Right
	end;

{result:= Истина если точка внутри сферы иначе Ложь}
end;

supernewbie · 03.03.2012, 18:43

юзай длинную арифметику тогда чтоли, формулу можно вот такую сделать:
(x2-x1)^2 + (y2-y1)^2 + (z2-z1)^2) < Radius^2

mikser · 03.03.2012, 18:53

Добавлено через 1 минуту

Цитата:

Сообщение от supernewbie

юзай длинную арифметику тогда чтоли, формулу можно вот такую сделать:
(x2-x1)^2 + (y2-y1)^2 + (z2-z1)^2) < Radius^2

Длинную арифметику пускай юзают астрономы и криптографы.
У нас не такие астрономические вычисления, это как из базуки по воробьям стрелять...

mira · 05.03.2012, 10:44

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

mikser · 05.03.2012, 13:09

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

Yegor · 12.03.2012, 10:57

mikser, самый простой и быстрый способ - это заменить сферу на куб

и проверять попадание в него банальными сравнениями. Или бот должен быть прецензионным?

mikser · 12.03.2012, 14:17

Цитата:

Сообщение от Yegor

mikser, самый простой и быстрый способ - это заменить сферу на куб

и проверять попадание в него банальными сравнениями. Или бот должен быть прецензионным?

Да Да я уже пару дней назад сам до этого дошёл.
На самом деле это решает всё. Оказалось всё настолько просто что
мне даже стало стыдно как я раньше до этого не додумался.
Алгоритм:

Опишем вокруг сферы куб (Сфера вписана в куб, 2D аналог окружность вписана в квадрат)
Если точка лежит за пределами куба -то лежит за пределами сферы //тут отсеятся 90% всех "плохих" точек
Если точка лежит внутри куба {тут начинается самое интересное} Стоп! А что такое Куб? Если взглянуть на "левый нижний угол куба" - это же начало декартовой системы координат. Точнее его зона + + + (где x y z положительны или 0) Как мы только это осознали сразу переводим координаты сферы в новую систему координат {По сути мы переносим точку отсчёта в левый нижний угол куба}
Переводим координаты точки в новую систему координат (так как она лежит внутри куба)
Приводим тип всех координат в Longint и считаем обычной формулой.
Переполнения не произойдёт потому что куб мал

Добавлено через 16 минут
Правда есть ещё шанс что куб будет занимать всю карту если сфера размером во всю карту или её половину.
Надо подумать как лучше поступить в этом случае.
Масштабирование может помочь но потеряется точность
и надо правильно подобрать коэффицент на который делить //желательно степень двойки
Как его правильно подобрать?
Кто силён в вычислительной геометрии?

Aries · 12.03.2012, 14:33

Чего-то ты перемудрил)))
Пункты 1,2 из алгоритма - хорошее решение.
Пункты 3,4 - извращение)
У тебя переполнение возникает в случае большого расстояния между точками и перемещение в др. систему координат никак не улучшит эту ситуацию)
В итоге алгоритм: 1, 2, 5.

В остальном хз... Хотя я немного сомневаюсь, что будут такие расстояния, если речь конечно идет о ла2.

mikser · 12.03.2012, 19:43

Цитата:

У тебя переполнение возникает в случае большого расстояния между точками и перемещение в др. систему координат никак не улучшит эту ситуацию)

Точно тут ты прав

Надо думать дальше
случай если радиус на пол карты тоже надо предусмотреть.

Yegor · 14.03.2012, 05:41

mikser, ты как математик простое превратил в сложное. Берешь 32 битный тип int и не будет никакого переполнения. В L2 координаты вписываются даже в +-1.000.000 с большим запасом.