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: loga(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(i*fi) (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;