File: [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Numbers / generic_complex_numbers.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:25 2000 UTC (23 years, 8 months ago) by maekawa
Branch: PHC, MAIN
CVS Tags: v2, maekawa-ipv6, RELEASE_1_2_3, RELEASE_1_2_2_KNOPPIX_b, RELEASE_1_2_2_KNOPPIX, RELEASE_1_2_2, RELEASE_1_2_1, HEAD Changes since 1.1: +0 -0
lines
Import the second public release of PHCpack.
OKed by Jan Verschelde.
|
package body Generic_Complex_Numbers is
TWO : constant number := one+one;
-- CREATORS :
function Create ( i : integer ) return Complex_Number is
res : Complex_Number;
begin
if i = 0
then res.re := +zero;
elsif i = 1
then res.re := +one;
else res.re := Create(i);
end if;
res.im := +zero;
return res;
end Create;
function Create ( f : number ) return Complex_Number is
res : Complex_Number;
begin
res.re := +f;
res.im := +zero;
return res;
end Create;
function Create ( re,im : number ) return Complex_Number is
res : Complex_Number;
begin
res.RE := +re;
res.IM := +im;
return res;
end Create;
function Conjugate ( c : Complex_Number ) return Complex_Number is
res : Complex_Number;
begin
res.RE := +c.RE;
res.IM := -c.IM;
return res;
end Conjugate;
-- COMPARISON and COPYING :
function Equal ( x,y : Complex_Number ) return boolean is
begin
return (Equal(x.RE,y.RE) and Equal(x.IM,y.IM));
end Equal;
procedure Copy ( x : in Complex_Number; y : in out Complex_Number ) is
begin
Copy(x.RE,y.RE);
Copy(x.IM,y.IM);
end Copy;
function "<" ( x,y : Complex_Number ) return boolean is
avx : number := AbsVal(x);
avy : number := AbsVal(y);
res : boolean := (avx < avy);
begin
Clear(avx);
Clear(avy);
return res;
end "<";
function ">" ( x,y : Complex_Number ) return boolean is
avx : number := AbsVal(x);
avy : number := AbsVal(y);
res : boolean := (avx > avy);
begin
Clear(avx);
Clear(avy);
return res;
end ">";
-- SELECTORS :
function REAL_PART ( x : Complex_Number ) return number is
begin
return +x.RE;
end REAL_PART;
function IMAG_PART ( x : Complex_Number ) return number is
begin
return +x.IM;
end IMAG_PART;
function AbsVal ( x : Complex_Number ) return number is
res : number := AbsVal(x.RE);
acc : number := AbsVal(x.IM);
begin
Add(res,acc);
Clear(acc);
return res;
end AbsVal;
function AbsVal ( x : Complex_Number ) return Complex_Number is
abx : number := AbsVal(x);
res : Complex_Number := Create(abx);
begin
Clear(abx);
return res;
end AbsVal;
-- ARITHMETIC OPERATIONS AS FUNCTIONS :
function "+" ( x : Complex_Number; y : number ) return Complex_Number is
begin
return (x.RE+y,+x.IM);
end "+";
function "-" ( x : Complex_Number; y : number ) return Complex_Number is
begin
return (x.RE-y,+x.IM);
end "-";
function "*" ( x : Complex_Number; y : number ) return Complex_Number is
begin
return (x.RE*y,x.IM*y);
end "*";
function "/" ( x : Complex_Number; y : number ) return Complex_Number is
begin
return (x.RE/y,x.IM/y);
end "/";
function "+" ( x : number; y : Complex_Number ) return Complex_Number is
begin
return (x+y.RE,+y.IM);
end "+";
function "-" ( x : number; y : Complex_Number ) return Complex_Number is
begin
return (x-y.RE,-y.IM);
end "-";
function "*" ( x : number; y : Complex_Number ) return Complex_Number is
begin
return (x*y.RE,x*y.IM);
end "*";
function "/" ( x : number; y : Complex_Number ) return Complex_Number is
res : Complex_Number;
acc,avyim,avyre : Number;
begin
if Equal(y.IM,zero)
then res.RE := x/y.RE;
res.IM := +zero;
elsif Equal(y.RE,zero)
then res.RE := +zero;
res.IM := x/y.IM; Min(res.IM);
else avyim := AbsVal(y.IM); avyre := AbsVal(y.RE);
if avyim < avyre
then acc := y.IM/y.RE;
res.RE := x*acc; Mul(acc,y.IM); Add(acc,y.RE); Div(res.RE,acc);
res.IM := x/acc; Min(res.IM);
Clear(acc);
elsif avyim > avyre
then acc := y.RE/y.IM;
res.IM := x*acc; Min(res.IM); Mul(acc,y.RE); Add(acc,y.IM);
res.RE := x/acc;
Clear(acc);
elsif Equal(y.IM,y.RE)
then acc := TWO*y.IM;
res.RE := x/acc;
res.IM := -res.RE;
Clear(acc);
else -- y.IM = -y.RE then
acc := TWO*y.IM;
res.RE := x/acc; Min(res.RE);
res.IM := x/acc; Min(res.IM);
Clear(acc);
end if;
Clear(avyim); Clear(avyre);
end if;
return res;
end "/";
function "+" ( x : Complex_Number ) return Complex_Number is
begin
return (+x.RE,+x.IM);
end "+";
function "-" ( x : Complex_Number ) return Complex_Number is
begin
return (-x.RE,-x.IM);
end "-";
function "+" ( x,y : Complex_Number ) return Complex_Number is
begin
return (x.RE+y.RE,x.IM+y.IM);
end "+";
function "-" ( x,y : Complex_Number ) return Complex_Number is
begin
return (x.RE-y.RE,x.IM-y.IM);
end "-";
function "*" ( x,y : Complex_Number ) return Complex_Number is
res : Complex_Number;
acc,avyim,avyre : number;
begin
if Equal(y.IM,zero)
then res.RE := x.RE*y.RE;
res.IM := x.IM*y.RE;
elsif Equal(y.RE,zero)
then res.RE := x.IM*y.IM; Min(res.RE);
res.IM := x.RE*y.IM;
else avyre := AbsVal(y.RE); avyim := AbsVal(y.IM);
if avyre < avyim
then acc := y.RE/y.IM;
res.RE := x.RE*acc; Sub(res.RE,x.IM); Mul(res.RE,y.IM);
res.IM := x.IM*acc; Add(res.IM,x.RE); Mul(res.IM,y.IM);
Clear(acc);
elsif avyre > avyim
then acc := y.IM/y.RE;
res.RE := x.IM*acc; Sub(res.RE,x.RE); Min(res.RE);
Mul(res.RE,y.RE);
res.IM := x.RE*acc; Add(res.IM,x.IM); Mul(res.IM,y.RE);
Clear(acc);
elsif Equal(y.RE,y.IM)
then res.RE := x.RE - x.IM; Mul(res.RE,y.RE);
res.IM := x.IM + x.RE; Mul(res.IM,y.RE);
else -- y.RE = -y.IM then
res.RE := x.RE + x.IM; Mul(res.RE,y.RE);
res.IM := x.IM - x.RE; Mul(res.IM,y.RE);
end if;
Clear(avyre); Clear(avyim);
end if;
return res;
end "*";
function "/" ( x,y : Complex_Number ) return Complex_Number is
res : Complex_Number;
acc,avyre,avyim : number;
begin
if Equal(y.IM,zero)
then res.RE := x.RE/y.RE;
res.IM := x.IM/y.RE;
elsif Equal(y.RE,zero)
then res.RE := x.IM/y.IM;
res.IM := x.RE/y.IM; Min(res.IM);
else avyre := AbsVal(y.RE); avyim := AbsVal(y.IM);
if avyre < avyim
then acc := y.RE/y.IM;
res.RE := x.RE*acc; Add(res.RE,x.IM);
res.IM := x.IM*acc; Sub(res.IM,x.RE);
Mul(acc,y.RE); Add(acc,y.IM);
Div(res.RE,acc);
Div(res.IM,acc);
Clear(acc);
elsif avyre > avyim
then acc := y.IM/y.RE;
res.RE := x.IM*acc; Add(res.RE,x.RE);
res.IM := x.RE*acc; Sub(res.IM,x.IM); Min(res.IM);
Mul(acc,y.IM); Add(acc,y.RE);
Div(res.RE,acc);
Div(res.IM,acc);
Clear(acc);
elsif Equal(y.RE,y.IM)
then acc := TWO*y.RE;
res.RE := x.RE + x.IM; Div(res.RE,acc);
res.IM := x.IM - x.RE; Div(res.IM,acc);
Clear(acc);
else -- y.RE = -y.IM then
acc := TWO*y.RE;
res.RE := x.RE - x.IM; Div(res.RE,acc);
res.IM := x.IM + x.RE; Div(res.IM,acc);
Clear(acc);
end if;
Clear(avyre); Clear(avyim);
end if;
return res;
end "/";
function "**" ( x : Complex_Number; m : integer ) return Complex_Number is
res : Complex_Number;
begin
if m = 0
then res := Create(1);
elsif m > 0
then res := +x; --Copy(x,res); <- Copy CAUSED BUS ERROR!!
for j in 2..m loop
Mul(res,x);
end loop;
else res := Create(1);
for j in 1..(-m) loop
Div(res,x);
end loop;
end if;
return res;
end "**";
-- ARITHMETIC OPERATIONS AS PROCEDURES :
procedure Add ( x : in out Complex_Number; y : in number ) is
begin
Add(x.RE,y);
end Add;
procedure Sub ( x : in out Complex_Number; y : in number ) is
begin
Sub(x.RE,y);
end Sub;
procedure Mul ( x : in out Complex_Number; y : in number ) is
begin
Mul(x.RE,y);
Mul(x.IM,y);
end Mul;
procedure Div ( x : in out Complex_Number; y : in number ) is
begin
Div(x.RE,y);
Div(x.IM,y);
end Div;
procedure Add ( x : in out Complex_Number; y : in Complex_Number ) is
begin
Add(x.RE,y.RE);
Add(x.IM,y.IM);
end Add;
procedure Sub ( x : in out Complex_Number; y : in Complex_Number ) is
begin
Sub(x.RE,y.RE);
Sub(x.IM,y.IM);
end Sub;
procedure Min ( x : in out Complex_Number ) is
begin
Min(x.RE);
Min(x.IM);
end Min;
procedure Mul ( x : in out Complex_Number; y : in Complex_Number ) is
res : Complex_Number;
acc,avyim,avyre : number;
begin
if Equal(y.IM,zero)
then Mul(x.RE,y.RE);
Mul(x.IM,y.RE);
elsif Equal(y.RE,zero)
then res.RE := x.IM*y.IM; Min(res.RE);
res.IM := x.RE*y.IM;
Clear(x); x := res;
else avyre := AbsVal(y.RE); avyim := AbsVal(y.IM);
if avyre < avyim
then acc := y.RE/y.IM;
res.RE := x.RE*acc; Sub(res.RE,x.IM); Mul(res.RE,y.IM);
res.IM := x.IM*acc; Add(res.IM,x.RE); Mul(res.IM,y.IM);
Clear(acc);
elsif avyre > avyim
then acc := y.IM/y.RE;
res.RE := x.IM*acc; Sub(res.RE,x.RE); Min(res.RE);
Mul(res.RE,y.RE);
res.IM := x.RE*acc; Add(res.IM,x.IM); Mul(res.IM,y.RE);
Clear(acc);
elsif Equal(y.RE,y.IM)
then res.RE := x.RE - x.IM; Mul(res.RE,y.RE);
res.IM := x.IM + x.RE; Mul(res.IM,y.RE);
else -- y.RE = -y.IM then
res.RE := x.RE + x.IM; Mul(res.RE,y.RE);
res.IM := x.IM - x.RE; Mul(res.IM,y.RE);
end if;
Clear(avyre); Clear(avyim);
Clear(x); x := res;
end if;
end Mul;
procedure Div ( x : in out Complex_Number; y : in Complex_Number ) is
res : Complex_Number;
acc,avyre,avyim : number;
begin
if Equal(y.IM,zero)
then Div(x.RE,y.RE);
Div(x.IM,y.RE);
elsif Equal(y.IM,zero)
then res.RE := x.IM/y.IM;
res.IM := x.RE/y.IM; Min(res.IM);
Clear(x); x := res;
else avyre := AbsVal(y.RE); avyim := AbsVal(y.IM);
if avyre < avyim
then acc := y.RE/y.IM;
res.RE := x.RE*acc; Add(res.RE,x.IM);
res.IM := x.IM*acc; Sub(res.IM,x.RE);
Mul(acc,y.RE); Add(acc,y.IM);
Div(res.RE,acc);
Div(res.IM,acc);
Clear(acc);
elsif avyre > avyim
then acc := y.IM/y.RE;
res.RE := x.IM*acc; Add(res.RE,x.RE);
res.IM := x.RE*acc; Sub(res.IM,x.IM); Min(res.IM);
Mul(acc,y.IM); Add(acc,y.RE);
Div(res.RE,acc);
Div(res.IM,acc);
Clear(acc);
elsif Equal(y.RE,y.IM)
then acc := TWO*y.RE;
res.RE := x.RE + x.IM; Div(res.RE,acc);
res.IM := x.IM - x.RE; Div(res.IM,acc);
Clear(acc);
else -- y.RE = -y.IM then
acc := TWO*y.RE;
res.RE := x.RE - x.IM; Div(res.RE,acc);
res.IM := x.IM + x.RE; Div(res.IM,acc);
Clear(acc);
end if;
Clear(avyre); Clear(avyim);
Clear(x); x := res;
end if;
end Div;
-- DESTRUCTOR :
procedure Clear ( x : in out Complex_Number ) is
begin
Clear(x.RE);
Clear(x.IM);
end Clear;
end Generic_Complex_Numbers;