File: [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Numbers / multprec_floating_numbers.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:26 2000 UTC (23 years, 7 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.
|
with Standard_Mathematical_Functions; use Standard_Mathematical_Functions;
package body Multprec_Floating_Numbers is
-- NOTES ON THE IMPLEMENTATION :
-- 0) See the notes in the bodies of Multprec_{Natural,Integer}_Numbers.
-- 1) 0.0 must have zero fraction and zero exponent.
-- 2) The truncate and round operations all depend on the decimal radix.
-- 3) The function "Decimal_Places" occupies a central place.
-- 4) The normalized representation has a nonzero least significant digit.
tol : constant double_float := 10.0**(-8);
-- AUXILIARIES :
function Truncate ( f : double_float ) return integer is
i : integer := integer(f);
begin
if i >= 0
then if double_float(i) > f + tol
then i := i-1;
end if;
else if double_float(i) < f - tol
then i := i+1;
end if;
end if;
return i;
end Truncate;
function Max_Size ( i1,i2 : Integer_Number ) return natural is
-- DESCRIPTION :
-- Return the maximal size of i1 and i2.
sz1 : constant natural := Size(i1);
sz2 : constant natural := Size(i2);
begin
if sz1 > sz2
then return sz1;
else return sz2;
end if;
end Max_Size;
procedure Shift_Left ( n : in out Natural_Number; ns : out natural ) is
-- DESCRIPTION :
-- Shifts n until the leading coefficient is larger than Base/10.
-- The number ns on return equals the number of shifts.
begin
if (Empty(n) or else Equal(n,0))
then ns := 0;
else ns := 0;
while (Coefficient(n,Size(n)) < Base/10) loop
Mul(n,10);
ns := ns+1;
end loop;
end if;
end Shift_Left;
procedure Shift_Right ( n : in out Natural_Number; ns : out natural ) is
-- DESCRIPTION :
-- Shifts n until the last digit in the least significant coefficient
-- of n is different from zero.
-- The number ns on return equals the number of shifts.
cff : natural;
begin
if (Empty(n) or else Equal(n,0))
then ns := 0;
else ns := 0;
cff := Coefficient(n,0) mod 10;
while (cff = 0) loop
Div(n,10);
ns := ns+1;
cff := Coefficient(n,0) mod 10;
end loop;
end if;
end Shift_Right;
procedure Shift_Right ( i : in out Integer_Number; ns : out natural ) is
-- DESCRIPTION :
-- Applies the Shift_Right to the number representation of i.
ni : Natural_Number;
begin
if Empty(i)
then ns := 0;
else --ni := Unsigned(i);
Copy(Unsigned(i),ni); -- avoid sharing
Shift_Right(ni,ns);
if Sign(i) < 0
then Clear(i); i := Create(ni); Min(i);
else Clear(i); i := Create(ni);
end if;
Clear(ni);
end if;
end Shift_Right;
procedure Normalize ( f : in out Floating_Number ) is
-- DESCRIPTION :
-- Normalizes the representation of f so that its fraction has as
-- rightmost digit a nonzero number.
ns : natural;
begin
Shift_Right(f.fraction,ns);
Add(f.exponent,ns);
end Normalize;
procedure Adjust_Expo ( f : in out double_float;
expo : in Integer_Number ) is
-- DESCRIPTION :
-- Multiplies f with 10**expo.
cnt : Integer_Number;
begin
Copy(expo,cnt);
while cnt > 0 loop
f := f*10.0;
Sub(cnt,1);
end loop;
while cnt < 0 loop
f := f/10.0;
Add(cnt,1);
end loop;
Clear(cnt);
end Adjust_Expo;
-- CONSTRUCTORS :
function Create ( i : integer ) return Floating_Number is
res : Floating_Number;
begin
res.fraction := Create(i);
res.exponent := Create(0);
Normalize(res);
return res;
end Create;
function Create ( i : Integer_Number ) return Floating_Number is
res : Floating_Number;
begin
res.fraction := +i;
res.exponent := Create(0);
Normalize(res);
return res;
end Create;
function Create ( fraction,exponent : integer ) return Floating_Number is
res : Floating_Number;
begin
res.fraction := Create(fraction);
res.exponent := Create(exponent);
Normalize(res);
return res;
end Create;
function Create ( fraction : Integer_Number;
exponent : integer ) return Floating_Number is
res : Floating_Number;
begin
res.fraction := +fraction;
res.exponent := Create(exponent);
Normalize(res);
return res;
end Create;
function Create ( fraction : integer;
exponent : Integer_Number ) return Floating_Number is
res : Floating_Number;
begin
res.fraction := Create(fraction);
res.exponent := +exponent;
Normalize(res);
return res;
end Create;
function Create ( fraction,exponent : Integer_Number )
return Floating_Number is
res : Floating_Number;
begin
res.fraction := +fraction;
res.exponent := +exponent;
Normalize(res);
return res;
end Create;
function Create ( f : double_float ) return Floating_Number is
-- NOTE :
-- To avoid rounding errors, comparison operations are prefered
-- above relying on intermediate results of calculations, although
-- rounding errors still occur!
res : Floating_Number;
ten : Integer_Number;
wrkflt,fltexp,lowflt,uppflt,incflt : double_float;
intexp : integer;
fraction,exponent,intinc : Integer_Number;
begin
if f = 0.0
then res := Create(0);
return res;
end if;
if f > 0.0
then wrkflt := f;
else wrkflt := -f;
end if;
fltexp := LOG10(wrkflt);
intexp := Truncate(fltexp);
if intexp < 15
then wrkflt := wrkflt*(10.0**(15-intexp));
exponent := Create(intexp-15);
intexp := 15;
else exponent := Create(0);
end if;
lowflt := 10.0**intexp;
uppflt := 10.0**(intexp+1);
while lowflt > wrkflt loop
intexp := intexp-1;
lowflt := 10.0**intexp;
uppflt := 10.0**(intexp+1);
end loop;
ten := Create(10);
fraction := ten**intexp;
for k in 1..16 loop
incflt := 10.0**intexp;
intinc := ten**intexp;
uppflt := lowflt;
for i in 1..10 loop
uppflt := uppflt + incflt;
exit when (uppflt > wrkflt);
lowflt := uppflt;
Add(fraction,intinc);
end loop;
intexp := intexp - 1;
Clear(intinc);
end loop;
Clear(ten);
if f < 0.0
then Min(fraction);
end if;
res := Create(fraction,exponent);
Normalize(res);
return res;
end Create;
-- SELECTORS :
function Fraction ( f : Floating_Number ) return Integer_Number is
begin
return f.fraction;
end Fraction;
function Exponent ( f : Floating_Number ) return Integer_Number is
begin
return f.exponent;
end Exponent;
function Size_Fraction ( f : Floating_Number ) return natural is
begin
return Size(Fraction(f));
end Size_Fraction;
function Size_Exponent ( f : Floating_Number ) return natural is
begin
return Size(Exponent(f));
end Size_Exponent;
function Decimal_Places_Fraction ( f : Floating_Number ) return natural is
begin
return Decimal_Places(Fraction(f));
end Decimal_Places_Fraction;
function Decimal_Places_Exponent ( f : Floating_Number ) return natural is
begin
return Decimal_Places(Exponent(f));
end Decimal_Places_Exponent;
function Decimal_to_Size ( deci : natural ) return natural is
res : natural;
exp : constant natural := Multprec_Natural_Numbers.Exponent;
begin
res := deci/exp;
if res*exp < deci
then res := res+1;
end if;
return res;
end Decimal_to_Size;
function AbsVal ( f : Floating_Number ) return Floating_Number is
res : Floating_Number;
begin
if Negative(f.fraction)
then res.fraction := -f.fraction;
else res.fraction := +f.fraction;
end if;
res.exponent := +f.exponent;
return res;
end AbsVal;
-- COMPARISON AND COPYING :
function Equal ( f1 : Floating_Number; f2 : double_float ) return boolean is
ff2 : Floating_Number := Create(f2);
res : boolean := Equal(f1,ff2);
begin
Clear(ff2);
return res;
end Equal;
function Equal ( f1,f2 : Floating_Number ) return boolean is
res : boolean;
f1deci : constant natural := Decimal_Places(f1.fraction);
f2deci : constant natural := Decimal_Places(f2.fraction);
f1expo : Integer_Number := f1.exponent + f1deci;
f2expo : Integer_Number := f2.exponent + f2deci;
mulfra : Integer_Number;
begin
if Equal(f1expo,f2expo)
then if f1deci = f2deci
then res := Equal(f1.fraction,f2.fraction);
else if f1deci < f2deci
then mulfra := f1.fraction*10;
for i in 1..(f2deci-f1deci-1) loop
Mul(mulfra,10);
end loop;
res := Equal(mulfra,f2.fraction);
else mulfra := f2.fraction*10;
for i in 1..(f1deci-f2deci-1) loop
Mul(mulfra,10);
end loop;
res := Equal(f1.fraction,mulfra);
end if;
Clear(mulfra);
end if;
else res := false;
end if;
Clear(f1expo); Clear(f2expo);
return res;
end Equal;
function "<" ( f1 : double_float; f2 : Floating_Number ) return boolean is
ff1 : Floating_Number := Create(f1);
res : boolean := (ff1 < f2);
begin
Clear(ff1);
return res;
end "<";
function "<" ( f1 : Floating_Number; f2 : double_float ) return boolean is
ff2 : Floating_Number := Create(f2);
res : boolean := (f1 < ff2);
begin
Clear(ff2);
return res;
end "<";
function "<" ( f1,f2 : Floating_Number ) return boolean is
res : boolean;
f1deci : constant natural := Decimal_Places(f1.fraction);
f2deci : constant natural := Decimal_Places(f2.fraction);
f1expo : Integer_Number := f1.exponent + f1deci;
f2expo : Integer_Number := f2.exponent + f2deci;
mulfra : Integer_Number;
begin
if (f1expo < f2expo)
then if Sign(f2.fraction) > 0
then res := true;
elsif Sign(f2.fraction) < 0
then res := false;
else res := (Sign(f1.fraction) < 0);
end if;
elsif (f1expo > f2expo)
then if Sign(f1.fraction) > 0
then res := false;
elsif Sign(f1.fraction) < 0
then res := true;
else res := (Sign(f2.fraction) > 0);
end if;
else if f1deci = f2deci
then res := (f1.fraction < f2.fraction);
else if f1deci < f2deci
then mulfra := f1.fraction*10;
for i in 1..(f2deci-f1deci-1) loop
Mul(mulfra,10);
end loop;
res := (mulfra < f2.fraction);
else mulfra := f2.fraction*10;
for i in 1..(f1deci-f2deci-1) loop
Mul(mulfra,10);
end loop;
res := (f1.fraction < mulfra);
end if;
Clear(mulfra);
end if;
end if;
Clear(f1expo); Clear(f2expo);
return res;
end "<";
function ">" ( f1 : double_float; f2 : Floating_Number ) return boolean is
ff1 : Floating_Number := Create(f1);
res : boolean := (ff1 > f2);
begin
Clear(ff1);
return res;
end ">";
function ">" ( f1 : Floating_Number; f2 : double_float ) return boolean is
ff2 : Floating_Number := Create(f2);
res : boolean := (f1 > ff2);
begin
Clear(ff2);
return res;
end ">";
function ">" ( f1,f2 : Floating_Number ) return boolean is
res : boolean;
f1deci : constant natural := Decimal_Places(f1.fraction);
f2deci : constant natural := Decimal_Places(f2.fraction);
f1expo : Integer_Number := f1.exponent + f1deci;
f2expo : Integer_Number := f2.exponent + f2deci;
mulfra : Integer_Number;
begin
if (f1expo > f2expo)
then if Sign(f1.fraction) > 0
then res := true;
elsif Sign(f1.fraction) < 0
then res := false;
else res := (Sign(f2.fraction) < 0);
end if;
elsif (f1expo < f2expo)
then if Sign(f2.fraction) > 0
then res := false;
elsif Sign(f2.fraction) < 0
then res := true;
else res := (Sign(f1.fraction) > 0);
end if;
else if f1deci = f2deci
then res := (f1.fraction > f2.fraction);
else if f1deci < f2deci
then mulfra := f1.fraction*10;
for i in 1..(f2deci-f1deci-1) loop
Mul(mulfra,10);
end loop;
res := (mulfra > f2.fraction);
else mulfra := f2.fraction*10;
for i in 1..(f1deci-f2deci-1) loop
Mul(mulfra,10);
end loop;
res := (f1.fraction > mulfra);
end if;
Clear(mulfra);
end if;
end if;
Clear(f1expo); Clear(f2expo);
return res;
end ">";
procedure Copy ( f1 : in Floating_Number; f2 : in out Floating_Number ) is
begin
Clear(f2);
Copy(f1.fraction,f2.fraction);
Copy(f1.exponent,f2.exponent);
end Copy;
-- EXPANDING AND SHORTENING THE MANTISSA :
function Expand ( i : Integer_Number; k : natural )
return Array_of_Naturals is
-- DESCRIPTION :
-- Returns the expanded coefficient vector of i.
cff : constant Array_of_Naturals := Coefficients(Unsigned(i));
res : Array_of_Naturals(cff'first..cff'last+k);
begin
res(cff'range) := cff(cff'range);
res(cff'last+1..res'last) := (cff'last+1..res'last => 0);
return res;
end Expand;
procedure Truncate ( i : in Integer_Number; k : in natural;
trc : out Array_of_Naturals; exp : out integer ) is
-- DESCRIPTION :
-- On return is the truncated coefficient vector of i and the number
-- of decimal places that are lost, or gained when exp is positive.
n : Natural_Number;
ns : natural; -- #shifts = #decimal places gained
begin
Copy(Unsigned(i),n);
Shift_Left(n,ns);
exp := k*Multprec_Natural_Numbers.Exponent - ns;
declare
cff : constant Array_of_Naturals := Coefficients(n);
begin
for i in trc'range loop
trc(i) := cff(i+k);
end loop;
Clear(n);
end;
end Truncate;
procedure Truncate ( i : in out Integer_Number; k : in natural;
exp : out integer ) is
n : Natural_Number := Unsigned(i);
neg : constant boolean := Negative(i);
ns : natural;
newint : Integer_Number;
begin
Shift_Left(n,ns);
exp := k*Multprec_Natural_Numbers.Exponent - ns;
declare
cff : constant Array_of_Naturals := Coefficients(n);
trc : Array_of_Naturals(cff'first..cff'last-k);
begin
for i in trc'range loop
trc(i) := cff(i+k);
end loop;
newint := Create(trc);
i := newint; --Copy(newint,i); --Clear(newint);
if neg
then Min(i);
end if;
end;
Clear(n);
end Truncate;
procedure Round ( i : in Integer_Number; k : in natural;
trc : out Array_of_Naturals; exp : out integer ) is
-- DESCRIPTION :
-- On return is the rounded coefficient vector of i and the number
-- of decimal places that are lost, or gained when exp is positive.
n : Natural_Number := Unsigned(i);
ns : natural; -- #shifts = #decimal places gained
begin
Shift_Left(n,ns);
exp := k*Multprec_Natural_Numbers.Exponent - ns;
declare
cff : constant Array_of_Naturals := Coefficients(n);
begin
for i in trc'range loop
trc(i) := cff(i+k);
end loop;
if k > 0
then if cff(k-1) >= Multprec_Natural_Numbers.Base/2
then trc(0) := trc(0)+1;
end if;
end if;
end;
Clear(n);
end Round;
function Expand ( f : Floating_Number; k : natural ) return Floating_Number is
res : Floating_Number;
cff : constant Array_of_Naturals := Expand(f.fraction,k);
begin
res.fraction := Create(Create(cff));
if Negative(f.fraction)
then Min(res.fraction);
end if;
Copy(f.exponent,res.exponent);
Normalize(res);
return res;
end Expand;
procedure Expand ( f : in out Floating_Number; k : in natural ) is
cff : constant Array_of_Naturals := Expand(f.fraction,k);
neg : constant boolean := Negative(f.fraction);
begin
Clear(f.fraction);
f.fraction := Create(Create(cff));
if neg
then Min(f.fraction);
end if;
Normalize(f);
end Expand;
function Round ( f : Floating_Number ) return double_float is
res : double_float;
szf,ns : natural;
n : Natural_Number;
expo : Integer_Number;
begin
if (Empty(f.fraction) or else Equal(f.fraction,0))
then res := 0.0;
else Copy(Unsigned(f.fraction),n);
Shift_Left(n,ns);
expo := f.exponent - ns;
szf := Size(n);
while Coefficient(n,szf) = 0 and szf > 0 loop
szf := szf - 1;
end loop;
res := double_float(Coefficient(n,szf));
if szf >= 1
then res := res*double_float(Multprec_Natural_Numbers.Base)
+ double_float(Coefficient(n,szf-1));
if szf >= 2
then if Coefficient(n,szf-2) > Multprec_Natural_Numbers.Base/2
then res := res+1.0;
end if;
ns := (szf - 1)*Multprec_Natural_Numbers.Exponent;
Add(expo,ns);
end if;
end if;
Adjust_Expo(res,expo);
Clear(n); Clear(expo);
end if;
return res;
end Round;
function Round ( f : Floating_Number; k : natural ) return Floating_Number is
res : Floating_Number;
cff : Array_of_Naturals(0..(Size(f.fraction)-k));
exp : integer;
begin
Round(f.fraction,k,cff,exp);
res.fraction := Create(Create(cff));
if Negative(f.fraction)
then Min(res.fraction);
end if;
Copy(f.exponent,res.exponent);
Add(res.exponent,exp);
Normalize(res);
return res;
end Round;
procedure Round ( f : in out Floating_Number; k : in natural ) is
exp : integer;
n : Natural_Number;
neg : constant boolean := Negative(f.fraction);
ns : natural;
begin
Copy(Unsigned(f.fraction),n);
Shift_Left(n,ns);
exp := k*Multprec_Natural_Numbers.Exponent - ns;
declare
cff : constant Array_of_Naturals := Coefficients(n);
trc : Array_of_Naturals(cff'first..cff'last-k);
begin
for i in trc'range loop
trc(i) := cff(i+k);
end loop;
if k > 0
then if cff(k-1) >= Multprec_Natural_Numbers.Base/2
then trc(0) := trc(0)+1;
end if;
end if;
Clear(f.fraction);
f.fraction := Create(trc);
if neg
then Min(f.fraction);
end if;
end;
Clear(n);
Add(f.exponent,exp);
Normalize(f);
end Round;
function Trunc ( f : Floating_Number ) return double_float is
res : double_float;
szf,ns : natural;
n : Natural_Number;
expo : Integer_Number;
begin
if (Empty(f.fraction) or else Equal(f.fraction,0))
then res := 0.0;
else n := Unsigned(f.fraction);
Shift_Left(n,ns);
expo := f.exponent - ns;
szf := Size(n);
while Coefficient(n,szf) = 0 and szf > 0 loop
szf := szf - 1;
end loop;
res := double_float(Coefficient(f.fraction,szf));
if szf >= 1
then res := res*double_float(Multprec_Natural_Numbers.Base)
+ double_float(Coefficient(f.fraction,szf-1));
if szf >= 2
then ns := (szf - 1)*Multprec_Natural_Numbers.Exponent;
Add(expo,ns);
end if;
end if;
Adjust_Expo(res,expo);
Clear(n); Clear(expo);
end if;
return res;
end Trunc;
function Trunc ( f : Floating_Number; k : natural ) return Floating_Number is
res : Floating_Number;
cff : Array_of_Naturals(0..(Size(f.fraction)-k));
exp : integer;
begin
Truncate(f.fraction,k,cff,exp);
res.fraction := Create(Create(cff));
if Negative(f.fraction)
then Min(res.fraction);
end if;
Copy(f.exponent,res.exponent);
Add(res.exponent,exp);
return res;
end Trunc;
procedure Trunc ( f : in out Floating_Number; k : in natural ) is
exp : integer;
n : Natural_Number;
neg : constant boolean := Negative(f.fraction);
ns : natural;
begin
Copy(Unsigned(f.fraction),n);
Shift_Left(n,ns);
exp := k*Multprec_Natural_Numbers.Exponent - ns;
declare
cff : constant Array_of_Naturals := Coefficients(n);
trc : Array_of_Naturals(cff'first..cff'last-k);
begin
for i in trc'range loop
trc(i) := cff(i+k);
end loop;
Clear(f.fraction);
f.fraction := Create(trc);
if neg
then Min(f.fraction);
end if;
end;
Clear(n);
Add(f.exponent,exp);
end Trunc;
procedure Set_Size ( f : in out Floating_Number; k : in natural ) is
sz : constant natural := Size_Fraction(f);
begin
if sz > k
then Round(f,sz-k);
elsif sz < k
then Expand(f,k-sz);
end if;
end Set_Size;
-- ARITHMETIC OPERATIONS as functions (no data sharing) :
function "+" ( f1 : Floating_Number; f2 : double_float )
return Floating_Number is
ff2 : Floating_Number := Create(f2);
res : Floating_Number := f1+ff2;
begin
Clear(ff2);
return res;
end "+";
function "+" ( f1 : double_float; f2 : Floating_Number )
return Floating_Number is
begin
return f2+f1;
end "+";
function "+" ( f1,f2 : Floating_Number ) return Floating_Number is
res : Floating_Number;
maxsz,diff_size,max_diff : integer;
diff_expo : Integer_Number;
cnt : natural;
begin
if (Empty(f1.fraction) or else Equal(f1.fraction,0))
then Copy(f2,res);
elsif (Empty(f2.fraction) or else Equal(f2.fraction,0))
then Copy(f1,res);
else maxsz := Max_Size(f1.fraction,f2.fraction) + 1;
if Equal(f1.exponent,f2.exponent)
then res.fraction := f1.fraction + f2.fraction;
Copy(f1.exponent,res.exponent);
else diff_expo := f1.exponent - f2.exponent;
max_diff := maxsz*Multprec_Natural_Numbers.Exponent;
if diff_expo > 0
then Copy(f1.fraction,res.fraction);
if diff_expo < 2*max_diff
then cnt := Create(diff_expo);
for i in 1..cnt loop
Mul(res.fraction,10);
end loop;
Add(res.fraction,f2.fraction);
Copy(f2.exponent,res.exponent);
else Copy(f1.exponent,res.exponent);
end if;
else Copy(f2.fraction,res.fraction);
Min(diff_expo);
if diff_expo < 2*max_diff
then cnt := Create(diff_expo);
for i in 1..cnt loop
Mul(res.fraction,10);
end loop;
Add(res.fraction,f1.fraction);
Copy(f1.exponent,res.exponent);
else Copy(f2.exponent,res.exponent);
end if;
end if;
Clear(diff_expo);
end if;
diff_size := Size(res.fraction) + 1 - maxsz;
if diff_size > 0
then Round(res,diff_size);
else Normalize(res);
end if;
end if;
return res;
end "+";
function "+" ( f : Floating_Number ) return Floating_Number is
res : Floating_Number;
begin
Copy(f,res);
return res;
end "+";
function "-" ( f : Floating_Number ) return Floating_Number is
res : Floating_Number;
begin
Copy(f,res);
Min(res.fraction);
return res;
end "-";
function "-" ( f1 : Floating_Number; f2 : double_float )
return Floating_Number is
minf2 : double_float := -f2;
res : Floating_Number := f1+minf2;
begin
return res;
end "-";
function "-" ( f1 : double_float; f2 : Floating_Number )
return Floating_Number is
ff1 : Floating_Number := Create(f1);
res : Floating_Number := ff1 - f2;
begin
Clear(ff1);
return res;
end "-";
function "-" ( f1,f2 : Floating_Number ) return Floating_Number is
res,minf2 : Floating_Number;
begin
minf2.exponent := f2.exponent;
minf2.fraction := -f2.fraction;
res := f1 + minf2;
Clear(minf2.fraction);
return res;
end "-";
function "*" ( f1 : Floating_Number; f2 : double_float )
return Floating_Number is
ff2 : Floating_Number := Create(f2);
res : Floating_Number := f1*ff2;
begin
Clear(ff2);
return res;
end "*";
function "*" ( f1 : double_float; f2 : Floating_Number )
return Floating_Number is
ff1 : Floating_Number := Create(f1);
res : Floating_Number := ff1*f2;
begin
Clear(ff1);
return res;
end "*";
function "*" ( f1,f2 : Floating_Number ) return Floating_Number is
res : Floating_Number;
diffsize : integer;
begin
if not (Empty(f1.fraction) or else Equal(f1.fraction,0))
then if not (Empty(f2.fraction) or else Equal(f2.fraction,0))
then res.fraction := f1.fraction*f2.fraction;
res.exponent := f1.exponent+f2.exponent;
diffsize := Size(res.fraction)
- Max_Size(f1.fraction,f2.fraction);
if diffsize > 0
then Round(res,diffsize);
end if;
Normalize(res);
end if;
end if;
return res;
end "*";
function "**" ( f : double_float; n : Natural_Number )
return Floating_Number is
ff : Floating_Number := Create(f);
res : Floating_Number := ff**n;
begin
Clear(ff);
return res;
end "**";
function "**" ( f : Floating_Number; n : Natural_Number )
return Floating_Number is
res : Floating_Number;
cnt : Natural_Number;
begin
if (Empty(n) or else Equal(n,0))
then res := Create(1);
else Copy(f,res);
cnt := Create(1);
while cnt < n loop
Mul(res,f);
Add(cnt,1);
end loop;
Clear(cnt);
Normalize(res);
end if;
return res;
end "**";
function "**" ( f : Floating_Number; i : integer ) return Floating_Number is
res : Floating_Number;
begin
if i = 0
then res := Create(1.0);
else if i > 0
then Copy(f,res);
for j in 1..(i-1) loop
Mul(res,f);
end loop;
else res := Create(1);
for j in 1..(-i) loop
Div(res,f);
end loop;
end if;
Normalize(res);
end if;
return res;
end "**";
function "**" ( f : double_float; i : Integer_Number )
return Floating_Number is
ff : Floating_Number := Create(f);
res : Floating_Number := ff**i;
begin
Clear(ff);
return res;
end "**";
function "**" ( f : Floating_Number; i : Integer_Number )
return Floating_Number is
res : Floating_Number;
cnt,n : Natural_Number;
begin
if (Empty(i) or else Equal(i,0))
then res := Create(1);
else n := Unsigned(i);
if i > 0
then Copy(f,res);
cnt := Create(1);
while cnt < n loop
Mul(res,f);
Add(cnt,1);
end loop;
else res := Create(1);
cnt := Create(0);
while cnt < n loop
Div(res,f);
Add(cnt,1);
end loop;
end if;
Clear(cnt);
Normalize(res);
end if;
return res;
end "**";
function "/" ( f1 : Floating_Number; f2 : double_float )
return Floating_Number is
ff2 : Floating_Number := Create(f2);
res : Floating_Number := f1/ff2;
begin
Clear(ff2);
return res;
end "/";
function "/" ( f1 : double_float; f2 : Floating_Number )
return Floating_Number is
ff1 : Floating_Number := Create(f1);
res : Floating_Number := ff1/f2;
begin
Clear(ff1);
return res;
end "/";
function Pos_Div ( f1,f2 : Floating_Number ) return Floating_Number is
-- DESCRIPTION :
-- This division only applies when f1 and f2 are both positive.
res : Floating_Number;
rmd_frac,new_rmd_frac : Integer_Number;
diffsize : integer;
maxsz : natural := Max_Size(f1.fraction,f2.fraction);
begin
Div(f1.fraction,f2.fraction,res.fraction,rmd_frac);
res.exponent := f1.exponent - f2.exponent;
diffsize := Size(res.fraction) - maxsz;
if not Equal(rmd_frac,0)
then while diffsize <= 0 loop
while rmd_frac < f2.fraction loop
Mul(rmd_frac,10);
Mul(res.fraction,10);
Sub(res.exponent,1);
end loop;
Div(rmd_frac,f2.fraction,new_rmd_frac);
Add(res.fraction,rmd_frac);
Clear(rmd_frac); rmd_frac := new_rmd_frac;
diffsize := Size(res.fraction) - maxsz;
exit when Equal(rmd_frac,0);
end loop;
end if;
Clear(rmd_frac);
if diffsize > 0
then Round(res,diffsize);
end if;
Normalize(res);
return res;
end Pos_Div;
function "/" ( f1,f2 : Floating_Number ) return Floating_Number is
res,minf1,minf2 : Floating_Number;
begin
if not (Empty(f1.fraction) or else Equal(f1.fraction,0))
then if not (Empty(f2.fraction) or else Equal(f2.fraction,0))
then if Multprec_Integer_Numbers.Positive(f1.fraction)
then if Multprec_Integer_Numbers.Positive(f2.fraction)
then res := Pos_Div(f1,f2);
else minf2.fraction := -f2.fraction;
minf2.exponent := f2.exponent;
res := Pos_Div(f1,minf2);
Clear(minf2.fraction);
Min(res);
end if;
else minf1.fraction := -f1.fraction;
minf1.exponent := f1.exponent;
if Multprec_Integer_Numbers.Positive(f2.fraction)
then res := Pos_Div(minf1,f2);
Min(res);
else minf2.fraction := -f2.fraction;
minf2.exponent := f2.exponent;
res := Pos_Div(minf1,minf2);
Clear(minf2.fraction);
end if;
Clear(minf1.fraction);
end if;
else raise NUMERIC_ERROR;
end if;
end if;
return res;
end "/";
-- ARITHMETIC OPERATIONS as procedures for memory management :
procedure Add ( f1 : in out Floating_Number; f2 : in double_float ) is
ff2 : Floating_Number := Create(f2);
begin
Add(f1,ff2);
Clear(ff2);
end Add;
procedure Add ( f1 : in out Floating_Number; f2 : in Floating_Number ) is
res : Floating_Number;
maxsz,diff_size,max_diff : integer;
diff_expo : Integer_Number;
cnt : natural;
begin
if (Empty(f1.fraction) or else Equal(f1.fraction,0))
then Copy(f2,f1);
elsif (Empty(f2.fraction) or else Equal(f2.fraction,0))
then null;
else maxsz := Max_Size(f1.fraction,f2.fraction) + 1;
if Equal(f1.exponent,f2.exponent)
then Add(f1.fraction,f2.fraction);
else diff_expo := f1.exponent - f2.exponent;
max_diff := maxsz*Multprec_Natural_Numbers.Exponent;
if diff_expo > 0
then if diff_expo < 2*max_diff
then cnt := Create(diff_expo);
for i in 1..cnt loop
Mul(f1.fraction,10);
end loop;
Add(f1.fraction,f2.fraction);
Copy(f2.exponent,f1.exponent);
else null;
end if;
else Copy(f2.fraction,res.fraction);
Min(diff_expo);
if diff_expo < 2*max_diff
then cnt := Create(diff_expo);
for i in 1..cnt loop
Mul(res.fraction,10);
end loop;
Add(f1.fraction,res.fraction);
Clear(res.fraction);
else Copy(f2.exponent,f1.exponent);
Clear(f1.fraction);
f1.fraction := res.fraction;
end if;
end if;
Clear(diff_expo);
end if;
diff_size := Size(f1.fraction) + 1 - maxsz;
if diff_size > 0
then Round(f1,diff_size);
else Normalize(f1);
end if;
end if;
end Add;
procedure Sub ( f1 : in out Floating_Number; f2 : in double_float ) is
minf2 : double_float := -f2;
begin
Add(f1,minf2);
end Sub;
procedure Sub ( f1 : in out Floating_Number; f2 : in Floating_Number ) is
minf2 : Floating_Number;
begin
minf2.exponent := f2.exponent;
minf2.fraction := -f2.fraction;
Add(f1,minf2);
Clear(minf2.fraction);
end Sub;
procedure Min ( f : in out Floating_Number ) is
begin
Min(f.fraction);
end Min;
procedure Mul ( f1 : in out Floating_Number; f2 : in double_float ) is
ff2 : Floating_Number := Create(f2);
begin
Mul(f1,ff2);
Clear(ff2);
end Mul;
procedure Mul ( f1 : in out Floating_Number; f2 : in Floating_Number ) is
diffsize,maxsize : integer;
begin
if Empty(f2.fraction)
then Clear(f1);
elsif Equal(f2.fraction,0)
then Clear(f1);
elsif Empty(f1.fraction)
then null;
elsif Equal(f1.fraction,0)
then null;
else maxsize := Max_Size(f1.fraction,f2.fraction);
Mul(f1.fraction,f2.fraction);
Add(f1.exponent,f2.exponent);
diffsize := Size(f1.fraction) - maxsize;
if diffsize > 0
then Round(f1,diffsize);
end if;
Normalize(f1);
end if;
end Mul;
procedure Div ( f1 : in out Floating_Number; f2 : in double_float ) is
ff2 : Floating_Number := Create(f2);
begin
Div(f1,ff2);
Clear(ff2);
end Div;
procedure Pos_Div ( f1 : in out Floating_Number; f2 : in Floating_Number ) is
-- DESCRIPTION :
-- Does the division in case f1 and f2 are both positive.
rmd_frac,new_rmd_frac : Integer_Number;
diffsize : integer;
maxsz : natural := Max_Size(f1.fraction,f2.fraction);
begin
Div(f1.fraction,f2.fraction,rmd_frac);
Sub(f1.exponent,f2.exponent);
diffsize := Size(f1.fraction) - maxsz;
if not Equal(rmd_frac,0)
then while diffsize <= 0 loop
while rmd_frac < f2.fraction loop
Mul(rmd_frac,10);
Mul(f1.fraction,10);
Sub(f1.exponent,1);
end loop;
Div(rmd_frac,f2.fraction,new_rmd_frac);
Add(f1.fraction,rmd_frac);
Clear(rmd_frac); rmd_frac := new_rmd_frac;
diffsize := Size(f1.fraction) - maxsz;
exit when Equal(rmd_frac,0);
end loop;
end if;
Clear(rmd_frac);
if diffsize > 0
then Round(f1,diffsize);
end if;
Normalize(f1);
end Pos_Div;
procedure Div ( f1 : in out Floating_Number; f2 : in Floating_Number ) is
minf1,minf2 : Floating_Number;
begin
if not (Empty(f1.fraction) or else Equal(f1.fraction,0))
then if not (Empty(f2.fraction) or else Equal(f2.fraction,0))
then if Multprec_Integer_Numbers.Positive(f1.fraction)
then if Multprec_Integer_Numbers.Positive(f2.fraction)
then Pos_Div(f1,f2);
else minf2.fraction := -f2.fraction;
minf2.exponent := f2.exponent;
Pos_Div(f1,minf2);
Clear(minf2.fraction); Min(f1);
end if;
else minf1 := -f1;
if Multprec_Integer_Numbers.Positive(f2.fraction)
then Pos_Div(minf1,f2);
Clear(f1); f1 := minf1; Min(f1);
else minf2.fraction := -f2.fraction;
minf2.exponent := f2.exponent;
Pos_Div(minf1,minf2);
Clear(minf2.fraction);
Clear(f1); f1 := minf1;
end if;
end if;
else raise NUMERIC_ERROR;
end if;
end if;
end Div;
-- DESTRUCTOR :
procedure Clear ( f : in out Floating_Number ) is
begin
Clear(f.fraction);
Clear(f.exponent);
end Clear;
end Multprec_Floating_Numbers;
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to multprec_floating_numbers.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Numbers