Power-Programmierung/CD1.mdf/pascal/library/dos/wct/math.pas

Pascal/Delphi Source File | 1992-01-01 | 7.5 KB | 333 lines

unit math;

{ Written by William C. Thompson - [e-mail] }

{ This unit was written to perform several basic mathematical calculations }
{ This unit automatically generates a random seed }

interface

const e=2.7182818284905;

function adjust(a:real):real;
procedure quad(a,b,c: real; var x1,x2: real; var im:real);
function relerror(observed,actual: real):real;
function pow(a,b:real):real;
function sign(r:real):integer;
function rmax(r,s:real):real;
function rmin(r,s:real):real;
function imax(m,n:integer):integer;
function imin(m,n: integer):integer;
function log(x,base: real): real;
function tan(a:real):real;
function arcsin(x:real):real;
function arccos(x:real):real;
function arctan2(y,x:real):real;
function degtorad(d:real):real;
function radtodeg(r:real):real;
function gcd(m,n:longint):word;
function lcm(m,n:integer):integer;
function prime(n:longint):boolean;
function rnd(a,b:real):real;
function rnd2(a,b:real):real;
function gaussian(mu,sigma:real):real;
function distance(x1,y1,x2,y2:real):real;
function findline(x1,y1,x2,y2: real; var a,b:real):boolean;
function hero(a,b,c:real):real;
function factorial(n:word):extended;
function stirling(x:real):extended;
function combination(n,r:word):word;
function permutation(n,r:word):word;

implementation

function adjust(a:real):real;
{ Adjusts angle to fit into (-pi,pi] }
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);
{ Solves any quadratic equation.  If im=0, x1 and x2 are
  two real solutions.  If im<>0, x1 and x2 are real parts
  of two solutions and im is imaginary part.  The two
  solutions would be x1 + im I and x2 - im I }
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;
{ Relative error }
begin
  if actual=0.0 then relerror:=abs(observed)
  else relerror:=abs(observed/actual-1)
end;

function pow(a,b: real):real;
{ Computes a^b }
begin
  if a=0 then
    if b=0 then pow:=1                           { 0^0 = 1 }
    else if b<0 then pow:=exp(b*ln(a))           { force error }
    else pow:=0                                  { 0^x = 0 }
  else if a<0 then
    if abs(b)<1e-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))                       { force error }
  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;
{ Computes log of x base base }
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;
{ Computes angle of point (x,y) }
var at2:real;
begin
  if x=0 then
    if y>=0 then arctan2:=pi/2
    else arctan2:=-pi/2
  else begin
    at2:=arctan(abs(y/x));
    if x>0 then
      if y>0 then arctan2:=at2
      else arctan2:=adjust(-at2)
    else
      if y>0 then arctan2:=adjust(pi-at2)
      else arctan2:=adjust(-pi+at2);
    end
end;

function degtorad(d:real):real;
{ Convert degrees to radians }
begin
  degtorad:=d*pi/180
end;

function radtodeg(r:real):real;
{ Convert radians to degrees }
begin
  radtodeg:=r*180/pi
end;

function gcd(m,n:longint):word;
{ Greatest Common Divisor }
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;
{ Least common multiple }
begin
  if (m=0) and (n=0) then lcm:=0
  else lcm:=abs(round((m*1.0)*n) div gcd(m,n))
end;

function prime(n:longint):boolean;
{ Tests a number for primeness }
var
  p: boolean;
  i, limit: word;
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;
{ divides interval from a to b, inclusive, into 65535 values and
  chooses one }
begin
  rnd:=a+random(65535)/65534*(b-a)
end;

function rnd2(a,b: real):real;
{ Similar to rnd, but MANY more divisions (over 4 million) }
begin
  rnd2:=a+(65535.0*random(65535)+random(65535))/(sqr(65535.0)-1)*(b-a)
end;

function gaussian(mu,sigma:real):real;
{ Produces a Gaussian distributed random variable with mean mu
  and standard deviation sigma }
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;
{ Finds equation for line y=ax+b between (x1,y1) and (x2,y2) }
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(x:real):extended;
{ computes factorial according to Stirling's approximation -
  this may not be correct...  it should be, but ask your math professor }
var
  s: extended;
  i: word;
begin
  s:=1;
  for i:=1 to round(x) do s:=s/e*round(x);
  stirling:=s*sqrt(round(x))*sqrt(2.0*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;

begin
  randomize
end.