Math, czyli matematyka studenta

Mały przegląd funkcji TP służących do operacji na liczbach rzeczywistych,
zespolonych i macierzach - oraz kilka przydatnych przekształceń.
Wszystkie opisane procedury i funkcje znajdują się w bibliotece math
dostępnej wraz z kodem źródłowym tutaj.

Oto co nam standardowo udostępnia TP 7.0:
sin - zwraca sinus kąta podanego w radianach (2*pi = 360°);
cos - zwraca cosinus podanego kąta;
ln - logarytm naturalny (przy podstawie e) z podanej liczby;
exp - podnosi e do zadanej potęgi;
arctan - zwraca wartość kąta (w radianach!), którego tangens podajemy.

Dodatkowo jeszcze:
sqrt - pierwiastek kwadratowy (x^0.5);
sqr - kwadrat (x^2);
round - zaokrągla;
trunc - ucina część ułamkową;
int - zwraca część całkowitą (zaokrągla w dół);
frac - zwraca część po przecinku;
random - zwraca losową liczbę z przedziału <0,1>.

Niedużo tego, ale zaraz zobaczymy, co z tego można wyczarować.

Na początek zadeklarujemy sobie typ float:

type float = single;

będziemy go używać zamiast typu real, który jest wyjątkowo
nieekonomiczny jeśli chodzi o czas procesora.

Typ single jest mało dokładny (10^-7), ale za to liczenie na nim jest
sensownie szybkie. Jeśli ktoś będzie potrzebował dużej dokładności, to
zamiast single niech podstawi double albo nawet extended. Ale NIE real.

Na pierwszy ogień funkcje sin i cos:

tangens:

function tg(const x : float) : float;
var y : float;
begin
  y := cos(x);
  if y < 0 then tg := sin(x)/y else tg := 3.4e38;
end;

Co do dzielenia przez zero to zakładam, że daje w wyniku nieskończoność.

W naszych warunkach będzie to największa dostępna na danym typie wartość. Dla single jest to 3.4e38.

Teraz cotangens:

function ctg(const x : float) : float;
var y : float;
begin
  y := sin(x);
  if y < 0 then tg := cos(x)/y else ctg := 3.4e38;
end;

(przy okazji - jeśli mamy funkcję tg zdefiniowaną w ten sposób,
że przy dzieleniu przez zero zwraca nieskończoność, to funkcję ctg można
utworzyć tak:

function ctg(const x : float) : float;
var y : float;
begin
  y := tg(x);
  if y < 0 then ctg := 1/y else ctg := 3.4e38;
end;

co jest gorszym, ale ciekawym rozwiązaniem.)

To teoretycznie wszystko w tym dziale funkcji. Kolej na arcusy:

ArcSin:

function ArcSin(const x : float) : float;
begin
  if x < 1 then ArcSin := ArcTan(x/sqrt (1-sqr (x))) else ArcSin := 3.4e38;
end;

ArcCos:

function ArcCos(const x : float) : float;
begin
  if abs(x) < 1 then ArcSin := ArcTan(x/sqrt (1-sqr (x))) else ArcCos := 3.4e38;
end;

A skoro już jesteśmy przy sinusach i im podobnych, to zdefiniujemy sobie
zestaw funkcji hiperbolicznych:

sinh - sinus hiperboliczny:

function Sinh(const x : float) : float;
begin
  sinh := (exp(x) - exp(-x))/2;
end;

cosh - cosinus hiperboliczny (tzw. funkcja łańcuchowa):

function cosh(const x : float) : float;
begin
  cosh := (exp(x) + exp(-x))/2;
end;

th - tangens hiperboliczny:

function th(const x : float) : float;
begin
  th := sinh(x)/cosh(x);
end;

cth - cotangens hiperboliczny:

function cth(const x : float) : float;
begin
  cth := cosh(x)/sinh(x);
end;

Jako ciekawostkę dodam, że chociaż powyższe funkcje nie mają nic
wspólnego z funkcjami sinusoidlanymi, to obowiązują dla nich
podobne prawa (np.: sinh' = cosh; cosh' = sinh; cosh^{2 - sinh}2 = 1).

Pora na funkcje wykładnicze. Mamy funkcję exp, czyli e^x, gdzie x jest
rzeczywiste (e to liczba Eulera równa około 2.71). Czyli jak widać można
podnieść sobie e do potęgi niecałkowitej, np.: do 1/2 (czyli spierwiastkować).

A czy można podnieść dowolną liczbę rzeczywistą do takiej niecałkowitej potęgi?

No jasne! Trzeba tylko uzależnić ją od e, tzn zapisać w postaci x = e^ln(x).

Wtedy x^{y = (e}ln(x))^{y = e}(ln(x)*y):

power - potęga rzeczywista (z liczby dodatniej!!!):

function power(const x,y : float) : float;
begin
  power := exp(ln(abs(x))*y);
end;

Jak widać x znajduje się pod logarytmem, czyli musi być dodatnie.
Natomiast y może być dowolny. Teraz można sobie sprawdzić czy to
działa: policzmy sqr(2)=4=power(2,2) i sqrt(2)=1.41=power(2,0.5).

Teraz logarytm. Turbo Pascal udostępnia nam funkcję ln, czyli
logarytm naturalny (logarytm przy podstawie e). Jeśli chcemy
policzyć algorytm przy innej podstawie, to sięgamy do wzoru z początku
liceum: log_a(b) = ln(a)/ln(b).

log - logarytm przy podstawie a (a > 0 i a <> 1) z liczby b (b > 0):

function log(const a,b : float) : float;
begin
  if (a > 0) and (b > 0) then
    log := ln(a)/ln(b)
  else
    log := 3.4e38;
end;

Jeśli piszesz program, który bardzo intensywnie korzysta z którejś z
wymienionych funkcji (najczęściej jest to sin, rzadziej cos), to warto
zbudować tablicę, do niej wrzucić wartości funkcji, a potem korzystać
z tablicy pobierając przybliżone wartości funkcji. Przyspiesza to
bardzo działanie programu.

Teraz mocniejsze czary: implementacja liczb zespolonych. Nie będę dokładnie tłumaczyć, co to są liczby
zespolone, bo jeśli tego nie wiesz, to znaczy że nie i potrzebujesz tej wiedzy.

Ogólnie liczba zespolona to suma dwóch liczb rzeczywistych, z tym, że druga jest pomnożona przez tzw. jednostkę urojoną - i. Wygląda to tak: x + i*y.

Jednostka urojona ma takie miłe właściwości, że i^2 = -1 (stąd wniosek, że istnieją pierwiastki liczb ujemnych itp).

Implementacja samego typu zespolonego wygląda tak:

type
  complex = record
  case byte of
    0: (re, im : float);
    1: (z,  fi : float);
  end;

gdzie jednocześnie są użyte dwa typu zapisu liczb zespolonych: klasyczna
x + i*y oraz eulerowska |z|exp(ifi) (zawsze nieujemny moduł z, oraz
faza fi-przedział od 0 do 2pi). Tu przydadzą się funkcje do przeliczania
jednej postaci na drugą oraz do pobierania części rzeczywistej
(Realis) i urojonej (Imaginarius).

function Re(const z : complex) : float; {część rzeczywista}
begin Re := z.Re end;

function Im(const z : complex) : float; {część urojona}
begin Im := z.Im end;

function zabs(const z : complex) : float; {moduł liczby zespolonej}
begin
  zabs := sqrt(sqr(z.re)+sqr(z.im));
end;

procedure efi(const z : complex;var y : complex);
{przeliczanie z postaci klasycznej na eulerowską}
begin
  y.z := zabs(z);
  if (z.re = 0) and (z.im < 0) then y.fi := pi/2 else y.fi := abs(arctan(z.im/z.re));

  if (z.re < 0) and (z.im > 0) then y.fi := y.fi + pi/2 else

  if (z.re < 0) and (z.im <= 0) then y.fi := y.fi + pi else

  if (z.re > 0) and (z.im < 0) then y.fi := y.fi + 3*pi/2;

end;

procedure compl(const z : complex;var y : complex);
{przeliczanie z postaci Eulera na klasyczną}
begin
  y.re := z.z*cos(z.fi);
  y.im := z.z*sin(z.fi);
end;

Przydałoby się też kilka procedurek obsługujących najprostsze działania na typie zespolonym:

procedure zsqr(const z : complex;var y : complex); {kwadrat}
begin
  y.re := sqr(z.re) - sqr(z.im);
  y.im := 2*z.im*z.re;
end;

procedure zmul(const a,b : complex;var y : complex); {dzielenie a/b}
begin
  y.re := a.re*b.re - a.im*b.im;
  y.im := a.re*b.im + a.im*b.re;
end;

procedure zadd(const a,b : complex;var y : complex); {dodawanie}
begin
  y.re := a.re + b.re;
  y.im := a.im + b.im;
end;

procedure zsub(const a,b : complex;var y : complex); {odejmowanie}
begin
  y.re := a.re - b.re;
  y.im := a.im - b.im;
end;

procedure zdiv(const a,b : complex;var y : complex); {dzielenie}
var
  c : complex;
begin
  c.re := b.re;
  c.im := -b.im;
  zmul(a,c,y);
  c.z := sqr(b.im) + sqr(b.re);
  y.re := y.re/c.z;
  y.im := y.im/c.z;
end;

Potrzebne też będzie pierwiastkowanie; trzeba pamiętać, że pierwiastek
stopnia n w dziedzinie zespolonej ma n rozwiązań, symetrycznie
rozłożonych wokół punktu (0,0).

procedure zsqrt(const z : complex;var x1,x2 : complex);
{pierwiastek kwadratowy - 2 rozwiązania}
var
  a : complex;
begin
  efi(z,a);

  a.re := sqrt(a.re);

  x1.re := a.re*sin(a.fi);
  x1.im := a.re*cos(a.fi);

  x2.re := -x1.re;
  x2.im := -x1.im;
end;

procedure zsqrn(const z : complex;n : byte;var X : array of complex);
{pierwiastek dowolnego naturalnego stopnia - n rozwiązań!!!}
var
  a : complex;
  i : byte;
begin
  efi(z,a);
  a.z := power(a.z,1/n);

  for i := 0 to n-1 do
  begin
    x[i].re := a.z*cos(a.fi+ (i/n)*2*pi);
    x[i].im := a.z*sin(a.fi+ (i/n)*2*pi);
  end;
end;

Na koniec zostało jeszcze potęgowanie rzeczywiste - można je wykorzystać do
policzenia pierwszego pierwiastka zespolonego dowolnego stopnia:

procedure zpower(const z : complex;const n : float;var y : complex);
var
  a : complex;
begin
  efi(z,y);
  a.z  := power(y.z,n);
  a.fi := y.fi*n;
  y.re := a.z*cos(a.fi);
  y.im := a.z*sin(a.fi);
end;

To chyba wszystko, co najważniejsze i najniezbędniejsze w temacie liczb zespolonych.

Kolej na wektory i macierze. Najpierw definicja typu macierzowego:

const
  matrixsize = 255 div SizeOf(float);

type
  PVector = ^TVector;
  TVector = array[1..matrixsize] of float;
  PMatrix = ^TMatrix;
  TMatrix = array[1..matrixsize] of TVector; {macierz kwadratowa}

To dość prostackie rozwiązanie, bo nie pozwala na użycie tylko niezbędnej
ilości pamięci, więc budowa dużych macierzy jest niemożliwa (zajmują one
bardzo dużo pamięci, zbyt dużo). Częściowo można to obejść stosując wskaźniki

PVector i PMatrix, ale to połowiczne załatwienie problemu. Jednak w tym artykule
nie chodzi o oszczędność i prędkość, ale w ogóle o umożliwienie wykonywania danej operacji.

To taka mała dygresja ~:-) Teraz konkrety - podstawowe operacje na macierzach:

procedure MAdd(const A,B : TMatrix;const x,y : byte;var N : TMatrix);
{dodawanie macierzy: N = A+B. x i y to wymiary macierzy A (takie same jak B)}
var
  i,j : byte;
begin
  for i := 1 to x do
  for j := 1 to y do
  N[i,j] := A[i,j] + B[i,j];
end;

procedure MSub(const A,B : TMatrix;const x,y : byte;var N : TMatrix);
{odejmowanie macierzy N = A-B}
var
  i,j : byte;
begin
  for i := 1 to x do
  for j := 1 to y do
  N[i,j] := A[i,j] - B[i,j];
end;

procedure MMul(const A,B : TMatrix;const Ax,By : byte;var N : TMatrix);
{mnożenie macierzy N = A*B. Ax - szerokość A, By - wysokość B; z założenia Ax = By}
var
  i,x,y : byte;

begin
  if Ax < By then exit;

  for x:=1 to Ax do
  for y:=1 to By do
  begin
    N[x,y] := 0;
    for i := 1 to Ax do
    N[x,y] := N[x,y] + A[i,y]*B[x,i];
  end;
end;

procedure MMulS(const A : TMatrix;const skalar : float;const x,y : byte;var N : TMatrix);
{mnożenie macierzy przez skalar: N = A*k}
var
  i,j : byte;

begin
  for i:=1 to x do
  for j:=1 to y do
  N[j,i] := A[i,j]*skalar;
end;

procedure MTrans(const A : TMatrix;const x,y : byte;var N : TMatrix);

{transponowanie macierzy: N = AT}
var
  i,j : byte;

begin
  for i:=1 to x do
  for j:=1 to y do
  N[j,i] := A[i,j];
end;

procedure MZeros(var A : TMatrix;const x,y : byte);
{zerowanie macierzy}
var
  i : byte;

begin
  for i:=1 to y do
  Fill16(A[1,i],x*sizeof(float),0);
end;

To tyle jeśli chodzi o najważniejsze rzeczy. Potęgowanie naturalne można zrealizować
poprzez pętlę mnożenia, gorzej z dzieleniem. Za jakiś czas pojawi się procedura
do odwracania macierzy (A^-1), na razie jednak trzeba się jakoś obejść
bez tego. Oto przykład: rozwiązywanie układu równań A*x = b metodą Gaussa
(analogiczna do x = Ab).

function MResolve(var A : TMatrix;d : word;var x,b : TVector) : boolean;
{rozwiązuje metodĄ Gaussa układ równań A*x=b, gdzie:
  A - macierz kwadratowa współczynników;
  d - wymiar macierzy A
  x - wektor wyników [x1,x2,..,xn]
  b - wektor współczynników (A[n,1]*x[1]+A[n,2]*x[2]+..=b[n]).
UWAGA!!! Procedura modyfikuje wszystkie podane macierze - A, x i b}
var
  c        : float;
  n,m,i,j  : word;
  ok       : boolean;

procedure reduce(var a,n : TVector;var ab,bb : float;i : word); {i-ta kolumna}
var
  y : float;
  w : word;

begin
  if a[i] = 0 then begin ok := false; exit; end;
  y := n[i]/a[i];
  for w := i to d do n[w] := n[w] - y*a[w];
  bb := bb - y*ab;
end;

begin
  ok := true;
  for i := 1 to d-1 do
  for j := i+1 to d do
  begin
    reduce(a[i],a[j],b[i],b[j],i);
    if not ok then break;
  end;

  if ok then
  for i := d downto 1 do
  begin
    if a[i,i] = 0 then begin ok := false; break end;
    x[i] := b[i]/a[i,i];

    for j := 1 to i-1 do
    begin
      c := a[j,i]/a[i,i];
      b[j] := b[j] - c*b[i];
      a[j,i] := 0;
    end;
  end;

  MResolve := ok;
end;

Funkcja zwróci wartość true, jeśli układ ma rozwiązanie bądź wiele rozwiązań,
jeśli układ nie ma rozwiązania, to zwrócona będzie wartość false. Po wywołaniu
funkcji w zmiennej x znajduje się wektor rozwiązań, macierz A oraz wektor b są niszczone.

Kod źródłowy biblioteki math znajduje się tutaj.