freepascal.ru

Всем доброго времени суток.
Предлагаю корректную реализацию функций arcsin(X) и arccos(X), а также функцию вычисления угла по синусу и косинусу либо по длине сторон.
Для того, чтобы функции работали корректно есть ряд констант:

Код: Выделить всё

var
    Matem : packed record
            PI,             { Число Пи                                        }
            PI12,           { Число Пи * 1/2                                  }
            PI2,            { Число Пи * 2                                    }
            Rad,            { Число 180 / Пи                                  }
            Q2_2,           { Число Корень из 2 / на 2                        }
    end;

procedure Init;
begin
  Matem.PI:=pi;             { Число Пи                                        }
  Matem.PI12:=pi/2;         { Расчет числа Пи * 1/2                           }
  Matem.PI2:=pi*2;          { Расчет числа Пи * 2                             }
  Matem.Rad:=180/pi;        { Расчет числа 180 / Пи                           }
  Matem.Q2_2:=sqrt(2)/2;    { Расчет корня из 2 / на 2                        }
end;


Ну и, собственно, сами функции:

Код: Выделить всё

function ArcSin(const X:extended)  :extended;
{     вычисление угла в пределах -pi/2<ArcSin(X)<pi/2                      }
begin
  if abs(X)>Matem.Q2_2
    then
      begin
        if (1-sqr(X))<0
          then
            ArcSin:=Matem.Pi12*Sign11(X)
          else
            begin
              if X<0
                then
                  ArcSin:=-(Matem.Pi12+ArcTan(sqrt(1-sqr(X))/X))
                else
                  ArcSin:=Matem.Pi12-ArcTan(sqrt(1-sqr(X))/X);
            end;
      end
    else
      ArcSin:=ArcTan(X/sqrt(1-sqr(X)));
end;{ArcSin}
{--------------------------------------------------------------------------}
   
function ArcCos(const X:extended)  :extended;
{ вычисление угла в пределах 0<ArcCos(X)<pi                                }
begin
  ArcCos:=Matem.PI12-ArcSin(X);
end;{ArcCos}


{--------------------------------------------------------------------------}
function Ugl(const S,C:extended):extended;
{ вычисление угла в пределах -pi<Ugl(S,C)<pi                               }
{ S=Sin, C=Cos - входные параметры                                         }
var
F : extended;

begin
  if abs(S) > abs(C)
    then
      begin
        F:=C/S;
        if F>=0
          then
            begin
              if C>=0
                then
                  begin
                    if (C=0) and (S<0)
                      then
                        Ugl:=-(Matem.PI12+ArcTan(F))
                      else
                        Ugl:=Matem.PI12-ArcTan(F);
                  end
                else
                 begin
                   if S>=0
                     then
                       Ugl:=(Matem.PI12+Matem.PI)-ArcTan(F)
                     else
                       Ugl:=-ArcTan(F)-Matem.PI12;
                 end;
            end
          else
            begin
              if C>=0
                then
                  Ugl:=-(Matem.PI12+ArcTan(F))
                else
                  begin
                    if S>=0
                      then
                        Ugl:=Matem.PI12-ArcTan(F)
                      else
                        Ugl:=-((Matem.PI12+Matem.PI)+ArcTan(F));
                  end;
            end;
      end
    else
      begin
        if C=0
          then
            Ugl:=0
          else
            begin
              if C>0
                then
                  Ugl:=ArcTan(S/C)
                else
                  begin
                    if S>=0
                      then
                        Ugl:=ArcTan(S/C)+Matem.PI
                      else
                        Ugl:=ArcTan(S/C)-Matem.PI
                  end;
            end;
      end;
end;{Ugl}

{--------------------------------------------------------------------------}
function Sign11(const X:extended)  :extended;
{ возвращает значение '1' со знаком X. Sign11(0)=1                         }
{--------------------------------------------------------------------------}
begin
  if (X<0) then Sign11:=-1 else Sign11:=1;
end;{Sign11}



Если кому-нибудь пригодятся, буду рад

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

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

Код: Выделить всё

{--------------------------------------------------------------------------}
function Sign11(const X:extended)  :extended;
{ возвращает значение '1' со знаком X. Sign11(0)=1                         }
{--------------------------------------------------------------------------}
begin
  if (X<0) then Sign11:=-1 else Sign11:=1;
end;{Sign11}


а if (X<-eps) не корректней будет?

Откройте для себя функцию ArcTan2 (модуль Math)

Мне для определения угла полярных координат хватает вот такой поделки (тип аргументов можно изменить на extended по вкусу):

Код: Выделить всё

function Arctan3r(x,y: Integer): Extended;
begin
  if x = 0 then
  begin
    if y <= 0 then
      Result := PI
    else
      Result := 0.0;
  end
  else
    Result := Arctan2(x,y);
end;


alexrayne
Нет у нас в проекте целых... То есть они, конечно есть Но не в углах и математике.
Код согласен, выглядит громоздко, но это тоже, к сожалению, требования...
Функции эти считают чуть точнее, чем стандартные функции мат. сопроцессора. Правда проверялись на точность они крайний раз во времена P-1.
А для нашего проекта точность, блин, важна... Иначе за 20 часов полёта можно не туда прилететь

Sergei I. Gorelkin
Сергей, спасибо, поизучаю.
Ибо получается как обычно. Написанный "низ" обычно не переписывается, а только дописывается. Крайний раз математику перекапывали в 2006 году, когда из Delphi+Win32 переползали в Lazarus+Linux (МСВС), а до этого... А до этого, наверное, когда в Win32 из DOS переползали...

P.S. Есть ещё модуль векторной алгебры... Если нужно или интересно - готов поделиться...

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

сопроцессор семейства х86 работает в родном режиме с числами extended, более точного числа по стандартувроде пока нету (как минимум фрюха не ведает), (амд позволяет работать с огромными числами мантиссой 128бит), и при соответсвующей настройке (о которой тоже надо позаботится) можно загрубить вычислитель до double или real.

alexrayne
Некорректно выразился. Точнее и корректнее, чем стандартные функции модуля math. Однако, после сообщения Сергея буду заново смотреть модуль fpc'шной математики.
Разница в функции sign в модуле math и у нас умышленна.

Добавлено спустя 18 часов 58 минут 40 секунд:
Маленький Disclaimer, то есть "отмазка".

Проект над которым я работаю имеет давнюю историю. Начало было ажно в 1973 году, когда и персоналок то ещё не было... Вычислениями занималась БЭСМ-6. Потом были, естественно, переработки кода в связи с появлениями персоналок, потом появилась Винда, потом Винды стало мало и началось "переползание" под Linux.

Но всегда основной идеей была гарантированная работа функций независимо от платформы. Пусть даже и в небольшой ущерб времени выполнения той или иной функции. Поэтому выше описанные функции будут работать везде. И для intel и для amd и для PowerPC и для всего, где работает (уже) fpc. Естественно, это накладывает некоторые ограничения на код. Половину математики можно реализовать с помощью ассемблера. Но нет никаких гарантий, что когда надо будет запустить что-нибудь из проекта на авиационном крейте, всё будет хорошо...

Вот.

Добавлено спустя 6 минут 3 секунды:
Чуть-чуть ещё...
Есть функция вычисления высоты над поверхностью при переходе от геоцентрической к географической СК. Изначально она была реализована через разложение в ряд, потому как не хватало производительности машин. Когда Пром. PC пошли с относительно быстрыми CPU её переписали для вычисления по точным формулам. И так со всем...

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

Добавлено спустя 23 часа 52 минуты 11 секунд:

Sergei I. Gorelkin писал(а):Откройте для себя функцию ArcTan2 (модуль Math)

Сергей, ещё раз спасибо за наводку по поводу Arctan2.
Погонял тесты. Работает в 3..8 раз быстрее по сравнению с нашей функцией Ugl при неизменности результата.

А вот Синус и Косинус из модуля Math работают чуть медленнее и результат различается... Несильно, но всё же.