wct@po.CWRU.Edu (William C. Thompson) (04/17/91)
Here's a MATH unit I've spent quite a good deal of time
working on, but I "suppose" I can share it....
unit math;
interface
const e=2.7182818284905;
function factorial(n:word):extended;
{ computes factorial }
function stirling(n:real):extended;
{ computes factorial according to Stirling's approximation }
function combination(n,r:word):word;
{ computes nCr }
function permutation(n,r:word):word;
{ computes nPr }
function adjust(a:real):real;
{ adjusts an angle to fit in (-pi,pi] }
procedure quad(a,b,c: real; var x1,x2: real; var im:real);
{ solves a quadratic equation a*x^2 + b*x +c
returns real parts of x1 and x2
stat returns 0 if solutions are real, imaginary part if complex }
function relerror(observed,actual: real):real;
{ computes relative error between the observed and the actual
If the actual is 0, it returns the absolute error }
function pow(a,b:real):real;
{ computes a^b - allows for 0^0=1 }
function sign(r:real):integer;
{ returns sign of number }
function rmax(r,s:real):real;
{ returns larger of two reals }
function rmin(r,s:real):real;
{ returns smaller of two reals }
function imax(m,n:integer):integer;
{ returns larger of two integers }
function imin(m,n: integer):integer;
{ returns smaller of two integers }
function log(x,base: real): real;
{ computes log of x base "base" }
function tan(a:real):real;
{ computes tan of angle - infinity is represented by 1e+38 }
function arcsin(x:real):real;
{ computes arcsin of x - returns [-pi/2,pi/2] }
function arccos(x:real):real;
{ computes arccos of x - returns [0,pi] }
function arctan2(y,x:real):real;
{ computes arctangent of y/x }
function gcd(m,n:integer):integer;
{ greatest common divisor }
function lcm(m,n:integer):integer;
{ least common multiple }
function prime(n:longint):boolean;
{ returns true if n is a prime number
note: 0 and 1 are NOT primes }
function rnd(a,b:real):real;
{ generates a uniformly distributed variable between [a,b]
note: if you want (a,b] or something similar, you must
use a+eps and b as parameters, where eps is some small value }
function gaussian(mu,sigma:real):real;
{ generates a gaussian random variable with parameters mu (mean)
and sigma (standard deviation) }
function distance(x1,y1,x2,y2:real):real;
{ computes distance from point 1 to point 2 }
function findline(x1,y1,x2,y2: real; var a,b:real):boolean;
{ finds equation of line going through (x1,y1) & (x2,y2)
a = slope b = constant - return TRUE if successful, FALSE o/w }
function hero(a,b,c:real):real;
{ Finds area of triangle with sides of length s }
implementation
function adjust(a:real):real;
begin
repeat
if a<=-pi then a:=a+2*pi;
if a>pi then a:=a-2*pi
until (a>-pi) and (a<=pi);
adjust:=a
end;
procedure quad(a,b,c:real; var x1,x2: real; var im:real);
var d,q: real;
begin
im:=0.0;
if a=0 then begin
x1:=-c/b; x2:=-c/b end
else begin
b:=b/a; c:=c/a; a:=1;
d:=b*b-4*c;
q:=-b/2;
if d<0 then begin
x1:=q;
x2:=q;
im:=sqrt(-d)/2
end
else begin
if b<0 then q:=q+sqrt(d)/2
else q:=q-sqrt(d)/2;
x1:=q;
if q=0 then x2:=-b
else x2:=c/q;
if x2<x1 then begin
d:=x1;
x1:=x2;
x2:=d
end
end { d>=0 }
end { a<>0 }
end;
function relerror(observed, actual: real):real;
begin
if actual=0.0 then relerror:=abs(observed)
else relerror:=abs(observed/actual-1)
end;
function pow(a,b: real):real;
begin
if a=0 then
if b=0 then pow:=1
else if b<0 then pow:=1/a
else pow:=0
else if a<0 then
if abs(b)<13-10 then pow:=1
else if relerror(b,round(b))<1e-8 then
pow:=(1-2*ord(odd(round(b))))*exp(b*ln(abs(a)))
else if (relerror(1/b,round(1/b))<1e-8) and odd(round(1/b)) then
pow:=-exp(b*ln(abs(a)))
else pow:=exp(b*ln(a))
end;
function sign(r:real):integer;
begin
if r<0 then sign:=-1
else if r>0 then sign:=1
else sign:=0
end;
function rmax(r,s:real):real;
begin if s>=r then rmax:=s else rmax:=r end;
function rmin(r,s:real):real;
begin if s<=r then rmin:=s else rmin:=r end;
function imax(m,n:integer):integer;
begin if m>=n then imax:=m else imax:=n end;
function imin(m,n:integer):integer;
begin if m<=n then imin:=m else imin:=n end;
function log(x,base:real):real;
begin
if (base=0) and (x=0) then log:=1
else if (base=0) and (x=1) then log:=0
else log:=ln(x)/ln(base)
end;
function tan(a:real):real;
begin
a:=adjust(a);
if abs(a-pi/2)<1e-10 then tan:=9.99e+37
else if abs(a-3*pi/2)<1e-10 then tan:=-9.99e+37
else tan:=sin(a)/cos(a)
end;
function arcsin(x:real):real;
begin
if x<0 then arcsin:=-arcsin(-x)
else
if x=1 then arcsin:=pi/2
else arcsin:=arctan(x/sqrt(1-x*x))
end;
function arccos(x:real):real;
begin
if x<0 then arccos:=pi-arccos(-x)
else
if x=0 then arccos:=pi/2
else arccos:=arctan(sqrt(1-x*x)/x)
end;
function arctan2(y,x:real):real;
var at2:real;
begin
if x=0 then
if y>0 then arctan2:=pi/2
else arctan2:=-pi/2
else begin
at2:=arctan(y/x);
if x<0 then at2:=adjust(at2+pi);
arctan2:=at2
end
end;
function gcd(m,n:integer):integer;
var k,l,r:integer;
begin
m:=abs(m); n:=abs(n);
if m=0 then gcd:=n else begin
k:=m;
l:=n;
while n<>0 do begin
r:=m mod n;
m:=n;
n:=r
end;
gcd:=m
end
end;
function lcm(m,n:integer):integer;
begin
if (m=0) and (n=0) then lcm:=0
else lcm:=abs((m*n) div gcd(m,n))
end;
function prime(n:longint):boolean;
var
p: boolean;
i, limit:longint;
begin
n:=abs(n);
if n in [0..1] then prime:=false
else if n in [2..3] then prime:=true
else begin
p:=false;
i:=2;
limit:=round(sqrt(n));
while (i<=limit) and not p do begin
if (n mod i=0) then p:=true;
i:=i+1
end;
prime:=not p
end
end;
function rnd(a,b: real):real;
begin
rnd:=a+(65535.0*random(65535)+random(65535))/(sqr(65535.0)-1)*(b-a)
end;
function gaussian(mu,sigma:real):real;
var
r,v1,v2: real;
begin
repeat
v1:=rnd(-1,1);
v2:=rnd(-1,1);
r:=sqr(v1)+sqr(v2)
until (r<1) and (r>0);
gaussian:=v1*sqrt(-2*ln(r)/r)*sigma+mu
end;
function distance(x1,y1,x2,y2:real):real;
begin
distance:=sqrt(sqr(x1-x2)+sqr(y1-y2))
end;
function findline(x1,y1,x2,y2: real; var a,b:real):boolean;
begin
if x2=x1 then begin
findline:=false;
exit
end;
findline:=true;
a:=(y2-y1)/(x2-x1);
b:=(x2*y1-x1*y2)/(x2-x1);
end;
function hero(a,b,c:real):real;
{ Finds area of triangle with sides of length s }
var
s: real;
begin
s:=a/2+b/2+c/2;
hero:=sqrt(s*(s-a)*(s-b)*(s-c));
end;
function factorial(n:word):extended;
{ computes factorial }
var
i: word;
f: extended;
begin
f:=1;
for i:=1 to n do f:=f*i;
factorial:=f
end;
function stirling(n:real):extended;
{ computes factorial according to Stirling's approximation }
var
s: extended;
begin
s:=1/n;
s:=s*exp(n*ln(n));
stirling:=s*exp(-n)/sqrt(2*pi);
end;
function combination(n,r:word):word;
{ computes nCr }
begin
combination:=round((factorial(n)/factorial(n-r))/factorial(r))
end;
function permutation(n,r:word):word;
{ computes nPr }
begin
permutation:=round(factorial(n)/factorial(n-r))
end;
END.
--
Ticking away, the moments that make up a dull day. | William C. Thompson
You fritter and waste the hours in an offhand way. | Michelson 620D (wct)
Kicking around on a piece of ground in your hometown,| a.k.a. Master of Time,
waiting for someone or something to show you the way.| Space, and Dimension