[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
1.3b7 CPX and CPXD bugfixes
|
From: |
Darren J Wilkinson |
|
Subject: |
1.3b7 CPX and CPXD bugfixes |
|
Date: |
Mon, 19 Mar 2001 09:29:38 +0000 |
Dear Keith and Co,
I've downloaded and installed 1.3 beta 7, and it seems to work OK
generally. However, within the Numeric section of the Required
library, the complex number classes, CPX and CPXD are still a bit of a
mess. It's just a hunch, but I'm guessing that complex numbers aren't
Keith's forte! ;-) Anyway, I've fixed up a couple of bugs in the
routines, and I've re-written most of the pre and post conditions, as
most of them were buggy. I've also attached some test classes (which I
actually emailed to Keith a year ago, but maybe they got lost). I'm
not certain that everything is exactly right, but things are a lot
better than they were! Keep up the good work everyone. Cheers,
--
Darren Wilkinson Email: address@hidden
WWW: http://www.btinternet.com/~darrenjwilkinson/-------------------------> GNU Sather - sourcefile <-------------------------
-- Copyright (C) 2000 by K Hopper, University of Waikato, New Zealand --
-- This file is part of the GNU Sather library. It is free software; you may --
-- redistribute and/or modify it under the terms of the GNU Library General --
-- Public License (LGPL) as published by the Free Software Foundation; --
-- either version 2 of the license, or (at your option) any later version. --
-- This library is distributed in the hope that it will be useful, but --
-- WITHOUT ANY WARRANTY without even the implied warranty of MERCHANTABILITY --
-- or FITNESS FOR A PARTICULAR PURPOSE. See Doc/LGPL for more details. --
-- The license text is also available from: Free Software Foundation, Inc., --
-- 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA --
--------------> Please email comments to <address@hidden> <--------------
partial class CPX{STP < $REAL{STP}, ATP} is
-- This class implements the mathematical notion of a complex number
-- within the constraints of the parameter type.
-- Some of the algorithms are taken from:
--
-- Press, Flannery, Teukolsky, and Vettering, "Numerical Recipes in C",
-- 2nd ed, CUP, 1993.
--
-- Some of the choices of branch cut were chosen to be consistent with:
--
-- Guy L. Steele, "Common Lisp, The Language", 2nd ed, Digital 1990
-- Version 1.1 Dec 2000. Copyright K Hopper, U of Waikato
-- Development History
-- -------------------
-- Date Who By Detail
-- ---- ------ ------
-- 8 Aug 97 kh Original from Sather 1.1 Dist
-- 7 Dec 00 kh included is_lt based on magnitude.
-- 18 Mar 01 djw Fixed bug in div, added routine
-- is_similar, and fixed/added numerous
-- pre and post conditions
include COMPARABLE ;
include BINARY ;
include COMPLEX_STR{STP} ;
-- include CPX_FUNCTIONS ; -- include if desired!
const negatable : BOOL := true ;
attr re,
im : STP ; -- Real and imaginary parts.
stub log : SAME ;
-- This routine returns the complex logarithm of self.
create_real(
val : STP
) : SAME is
-- This routine creates a complex number which has a zero imaginary
-- component.
return create(val,STP::create(0.0))
end ;
build(
cursor : BIN_CURSOR
) : SAME
pre ~void(cursor)
and ~cursor.is_done
post true
is
-- This routine returns the complex number contained in the indicated
-- string at the current position.
return create(STP::build(cursor),STP::build(cursor))
end ;
create(
val : CARD
) : SAME is
-- This version of create produces one which has an integral real part
-- but zero imaginary part.
return create_real(STP::create(val))
end ;
create(
val : FIELD
) : SAME is
-- This version of create produces one which has an integral real part
-- but zero imaginary part.
return create_real(STP::create(val))
end ;
create(
val : INT
) : SAME is
-- This version of create produces one which has an integral real part
-- but zero imaginary part.
return create_real(STP::create(val))
end ;
create(
val : INTI
) : SAME
pre true
post (result.re = STP::create(val))
and (result.im = STP::zero)
is
-- This version of create produces one which has an integral real part
-- but zero imaginary part.
return create_real(STP::create(val))
end ;
create(
val : RAT
) : SAME
pre true
post (result.re = STP::create(val))
and (result.im = STP::zero)
is
-- This version of create produces one which has an integral real part
-- but zero imaginary part.
return create_real(STP::create(val))
end ;
create(
val : FLT
) : SAME is
-- This routine creates a complex number with a real part val and
-- zero imaginary part.
return create_real(STP::create(val))
end ;
create(
val : FLTD
) : SAME
pre (STP::maxval.fltd >= val)
and (val >= -STP::maxval.fltd)
post true
is
-- This routine creates a complex number with a real part val and
-- zero imaginary part.
return create_real(STP::create(val))
end ;
zero : SAME is
-- This routine provides a complex zero value.
return create(STP::zero,STP::zero)
end ;
one : SAME is
-- This routine provides a complex number with unit real part and zero
-- imaginary part.
return create(STP::one,STP::zero)
end ;
maxval : SAME is
-- This routine creates a complex number which has the maximum
-- representable real and imaginary parts.
return create(STP::maxval,STP::maxval)
end ;
minval : SAME is
-- This routine creates a complex number which has the minimum
-- representable real and imaginary parts.
return create(STP::minval,STP::minval)
end ;
nil : SAME is
-- This predicate returns a nil complex value.
return create(re.nil,im.nil)
end ;
private absolute : STP
pre true
post create(result.square).is_similar(create(self.magnitude_squared))
is
-- This private routine returns the absolute magnitude of self which is
-- calculated using the algorithm in 'Numerical Recipes in C' p949.
loc_re : STP := re.abs ;
loc_im : STP := im.abs ;
temp : STP ;
if loc_re = STP::zero then
return loc_im
elsif loc_im = STP::zero then
return loc_re
elsif loc_re > loc_im then
temp := loc_im / loc_re ;
return loc_re * (STP::one + temp * temp).sqrt
else
temp := loc_re / loc_im ;
return loc_im * (STP::one + temp * temp).sqrt
end
end ;
abs : SAME
pre true
post (result.re = absolute)
and (result.im = STP::zero)
is
-- This routine returns the absolute value of self. It is here to
-- conform to the interface of $NFE.
return create_real(absolute)
end ;
magnitude : STP
pre true
post result = absolute
is
-- This routine returns the absolute magnitude of self.
return absolute
end ;
magnitude_squared : STP
pre (re / STP::maxval)*re + (im / STP::maxval)*im < STP::one
post true
is
-- This routine returns the square of the absolute magnitude of self.
return re * re + im * im
end ;
conjugate : SAME
pre true
post (self*result).is_similar(create(self.magnitude_squared))
is
-- This routine returns the complex conjugate of self.
return create(re,-im)
end ;
is_eq(
other : SAME
) : BOOL is
-- This predicate returns true if and only if self and other have
-- identical values of real and imaginary parts.
return re = other.re
and im = other.im
end ;
is_lt(
other : SAME
) : BOOL is
-- This predicate returns true if and only if the magnitude of self is
-- less than that of other, otherwise false
return magnitude < other.magnitude
end ;
is_nil : BOOL is
-- This predicate returns true if either component of the number is nil.
return re.is_nil
or im.is_nil
end ;
is_neg : BOOL is
-- This routine returns true if and only if both components are negative.
return (re < STP::zero)
and (im < STP::zero)
end ;
is_zero : BOOL is
-- This routine returns true if and only if both components are zero.
return (re = STP::zero)
and (im = STP::zero)
end ;
is_pos : BOOL is
-- This routine returns true if and only if both components are positive.
return (re > STP::zero)
and (im > STP::zero)
end ;
is_within(
radius : STP,
other : SAME
) : BOOL is
-- This predicate returns true if and only if self is within the given
-- radius of other.
return (self - other).magnitude_squared <= radius*radius
end ;
is_similar(
other : SAME
) : BOOL is
tol::=STP::epsilon.sqrt;
return is_within(tol,other);
end ;
sign : NUM_SIGNS is
-- This routine returns the sign of self which is negative if either
-- component is negative, zero if both are zero - otherwise positive..
if (re < STP::zero)
or (im < STP::zero) then
return NUM_SIGNS::Negative
elsif self = zero then
return NUM_SIGNS::Zero
else
return NUM_SIGNS::Positive
end
end ;
plus(
other : SAME
) : SAME
pre ( ((re / STP::maxval) + (other.re / STP::maxval) < STP::one)
and ((re / STP::maxval) + (other.re / STP::maxval) > -STP::one)
and ((im / STP::maxval) + (other.im / STP::maxval) < STP::one)
and ((im / STP::maxval) + (other.im / STP::maxval) > -STP::one) )
post self.is_similar(result-other)
is
-- This routine returns the sum of self and other.
return create(re + other.re,im + other.im)
end ;
minus(
other : SAME
) : SAME
pre ((self.re.sign = other.re.sign)
or ((STP::maxval - self.re.abs) >= other.re.abs))
and ((self.im.sign = other.im.sign)
or ((STP::maxval - self.im.abs) >= other.im.abs))
post true
is
-- This routine returns the complex difference of subtracting other from
-- self.
return create(re - other.re,im - other.im)
end ;
negate : SAME
pre true
post zero.is_similar(self+result)
is
-- This routine returns the additive inverse of self.
return create(-re,-im)
end ;
times(
other : SAME
) : SAME
pre (
((re / STP::maxval)*other.re - (im / STP::maxval)*other.im
< STP::one) and
((re / STP::maxval)*other.re - (im / STP::maxval)*other.im
> -STP::one) and
((re / STP::maxval)*other.im + (im / STP::maxval)*other.re
< STP::one) and
((re / STP::maxval)*other.im + (im / STP::maxval)*other.re
> -STP::one)
)
post true
is
-- This routine returns the complex product of self and other.
return create(re * other.re - im * other.im,
re * other.im + im * other.re)
end ;
div(
other : SAME
) : SAME
pre true
post self.is_similar(result*other)
is
-- This routine returns the result of complex division of self by other.
denom,
res : STP ;
if other.re.abs >= other.im.abs then
res := other.im/other.re ;
denom := other.re + res * other.im ;
res := res / denom ; -- to make sure no overflow!
return create((re/denom) + (res * im),(im/denom) - (res * re))
else
res := other.re/other.im ;
denom := other.im + res * other.re ;
res := res / denom ; -- to make sure no overflow!
return create((im/denom) + (res * re), (res * im) - (re/denom))
end
end ;
mod(
other : SAME
) : SAME is
-- This routine returns the remainder of the result of dividing self by
-- other. This is zero for a complex number.
return create(0)
end ;
times(
factor : STP
) : SAME
pre (((factor.abs > STP::one)
and ((STP::maxval / factor.abs) <= re.abs))
or ((STP::maxval * factor.abs) >= re.abs))
and (((factor.abs > STP::one)
and ((STP::maxval / factor.abs) <= im.abs))
or ((STP::maxval * factor.abs) >= im.abs))
post result.is_similar(self*create(factor))
is
-- This routine scales both real and imaginary components of self by
-- the given factor.
return create(re * factor,im * factor)
end ;
div(
divisor : STP
) : SAME
pre (divisor.abs >= STP::one)
or (((divisor.abs * STP::maxval) <= re.abs)
and ((divisor.abs * STP::maxval) <= im.abs))
post result.is_similar(self/create(divisor))
is
-- This routine divides both components of self by the given divisor.
return create(re / divisor, im /divisor)
end ;
pow(
other : SAME
) : SAME
pre ~((re = STP::zero)
and (im = STP::zero))
is
-- This routine returns the result of raising self to the power of other.
return (log * other).exp
end ;
reciprocal : SAME
pre true
post one.is_similar(self*result)
is
-- This routine returns the multiplicative inverse of self.
denom,
res : STP ;
if re.abs >= im.abs then
res := im/re ;
denom := re + res * im ;
return create(STP::one/denom,-res/denom)
else
res := re/im ;
denom := im + res * re ;
return create(res/denom,(-STP::one)/denom)
end
end ;
exp : SAME is
-- This routine returns the complex exponential `e^self'.
real_part : STP := re.exp ;
phase : ATP := ATP::radians(im) ;
return create(real_part * phase.cos,real_part * phase.sin)
end ;
sqrt : SAME
pre true
post self.is_similar(result.square)
is
-- This routine returns the complex square root of self. The algorithm
-- is taken from 'Numerical Recipes in C' p949, choosing the branch cut by
--
-- e^((log z)/2)
if re = STP::create(STP::zero) -- zero is special case.
and im = STP::create(STP::zero) then
return create(STP::create(STP::zero),STP::create(STP::zero))
end ;
loc_re : STP := re.abs ;
loc_im : STP := im.abs ;
trial_val : STP ;
loc_half : STP := STP::one / (STP::one + STP::one) ;
if loc_re >= loc_im then
tmp : STP := loc_im / loc_re ;
trial_val := loc_re.sqrt * (loc_half * STP::one) +
(STP::one + tmp * tmp).sqrt.sqrt
else
tmp : STP := loc_re / loc_im ;
trial_val := loc_im.sqrt * (loc_half *
(tmp + (STP::one + tmp * tmp).sqrt)).sqrt
end ;
loc_two : STP := STP::one + STP::one ;
if re >= STP::zero then
return create(trial_val,im / (loc_two * trial_val))
elsif im >= STP::zero then
return create(im / (loc_two * trial_val),trial_val)
else
return create(-im / (loc_two * trial_val),-trial_val)
end
end ;
cube_root : SAME
pre true
post self.is_similar(result.cube)
is
-- This routine returns the complex cube root of self using a preliminary
-- algorithm.
loc_three : STP := STP::one + STP::one + STP::one ;
return self.pow(create_real(STP::one/loc_three))
end ;
square : SAME
-- pre (STP::maxval / re.square < im.square)
post result.is_similar(self.pow(one+one))
is
-- This routine returns the square of self.
return self * self
end ;
cube : SAME
pre (STP::maxval / (re.square * re) < (im.square * im))
post result.is_similar(self.pow(one+one+one))
is
-- This routine returns the cube of self.
return self * self * self
end ;
binstr : BINSTR
pre true
post build(result.cursor) = self
is
-- This routine returns a binary representation of self.
return re.binstr + im.binstr
end ;
end ; -- CPX{T}
-------------------------------------------------------------------------------
immutable class CPX < $COMPLEX{FLT,CPX}, $OPTION, $FLT_FMT is
-- This class implements the class of complex numbers which have real
-- components (of FLT class).
-- Version 1.0 Aug 97. Copyright K Hopper, U of Waikato
-- Development History
-- -------------------
-- Date Who By Detail
-- ---- ------ ------
-- 8 Aug 97 kh Original from Sather 1.1 Dist
-- 18 Mar 01 djw Fixed bug in log, and added
-- post condition to log.
include CPX{FLT,ANGLE} ;
create(
real,
imaginary : FLT
) : SAME is
-- This routine creates a complex number with a real part `re' and
-- imaginary part `im'.
me : SAME ;
return me.re(real).im(imaginary)
end ;
log : CPX
post self.is_similar(result.exp)
is
-- This routine returns the complex logarithm of self. The chosen
-- branch is
--
-- log |self| + i phase(self). See Steele p302.
phase : ANGLE := ANGLE::atan2(im,re) ;
magnitude : FLT := (re * re + im * im).sqrt.log ;
return create(magnitude , phase.radians)
end ;
end ; -- CPX
-------------------------> GNU Sather - sourcefile <-------------------------
-- Copyright (C) 2000 by K Hopper, University of Waikato, New Zealand --
-- This file is part of the GNU Sather library. It is free software; you may --
-- redistribute and/or modify it under the terms of the GNU Library General --
-- Public License (LGPL) as published by the Free Software Foundation; --
-- either version 2 of the license, or (at your option) any later version. --
-- This library is distributed in the hope that it will be useful, but --
-- WITHOUT ANY WARRANTY without even the implied warranty of MERCHANTABILITY --
-- or FITNESS FOR A PARTICULAR PURPOSE. See Doc/LGPL for more details. --
-- The license text is also available from: Free Software Foundation, Inc., --
-- 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA --
--------------> Please email comments to <address@hidden> <--------------
immutable class CPXD < $COMPLEX{FLTD, CPXD}, $FLT_FMT is
-- This class implements the mathematical notion of a 'double length'
-- complex number within the constraints of the parameter type.
-- Some of the algorithms are taken from:
-- Press, Flannery, Teukolsky, and Vettering, "Numerical Recipes in C",
-- second edition, Cambridge University Press, 1993.
--
-- Some of the choices of branch cut were chosen to be consistent with:
-- Guy L. Steele, "Common Lisp, The Language", second edition,
-- 1990, Digital Press.
-- Version 1.1 Oct 98. Copyright K Hopper, U of Waikato
-- Development History
-- -------------------
-- Date Who By Detail
-- ---- ------ ------
-- 25 Aug 97 kh Original from Sather 1.1 Dist
-- 2 Oct 98 kh Factored out formatting - see CPX{?}
-- 18 Mar 01 djw Fixed bugs as for CPX
include CPX{FLTD,ANGLED} ;
create(
real,
imaginary : FLTD
) : SAME is
-- This routine creates a complex number with a real part `re' and
-- imaginary part `im'.
me : SAME ;
return me.re(real).im(imaginary)
end ;
log : CPXD
post self.is_similar(result.exp)
is
-- This routine returns the complex logarithm of self. The chosen
-- branch is
--
-- log |self| + i phase(self). See Steele p302.
phase : ANGLED := ANGLED::atan2(im,re) ;
magnitude : FLTD := (re * re + im * im).sqrt.log ;
return create(magnitude , phase.radians)
end ;
end ; -- CPXD
class TEST_CPX is
-- This class carries out tests on the complex number class CPX. Note
-- that it is not portable for simplicity in building a local library.
-- Version 1.0 Nov 97. Copyright K Hopper, U of Waikato
-- Development History
-- -------------------
-- Date Who By Detail
-- ---- ------ ------
-- 21 Nov 97 kh Original from Sather distribution
-- 3 Apr 00 djw Extended to cover more cases
include TEST ;
main is
class_name("CPX") ;
first : CPX := CPX::create(3.0,4.0) ;
second : CPX := CPX::create(5.0,7.0) ;
third : CPX := CPX::create(8.0,11.0) ;
sum : CPX := first + second ;
test("plus",sum.str,third.str) ;
fourth : CPX := CPX::create(-2.0,-3.0) ;
diff : CPX := first - second ;
test("minus",diff.str,fourth.str) ;
neg : CPX := first + (second.negate) ;
test("negate",neg.str,fourth.str) ;
prod : CPX := first * second ;
fifth : CPX := CPX::create(-13.0,41.0) ;
test("times",prod.str,fifth.str) ;
div : CPX := second / first ;
sixth : CPX := CPX::create(1.72,0.04) ;
test("divide",div.str,sixth.str) ;
recip : CPX := second * (first.reciprocal) ;
test("reciprocal",recip.str,sixth.str) ;
dfirst : CPX := CPX::create(4.0,5.0) ;
dsecond : CPX := dfirst.conjugate ;
test("conjugate",dsecond.str,CPX::create(4.0,-5.0).str) ;
dpow : CPX := dfirst.pow(CPX::create(2.0)) ;
dthird : CPX := CPX::create(-9.0,40.0) ;
test("pow",dpow.re.str(8)+" "+dpow.im.str(6),
dthird.re.str(8)+" "+dthird.im.str(6)) ;
dsqr : CPX := dfirst.square ;
test("square",dsqr.str,dthird.str) ;
dsqrt : CPX := dsqr.sqrt ;
test("sqrt",dsqrt.re.str(8)+" "+dsqrt.im.str(8),
dfirst.re.str(8)+" "+dfirst.im.str(8)) ;
dmag : FLT := dfirst.magnitude_squared ;
test("magnitude_squared",dmag.str,41.0.str) ;
dmagg : FLT := dfirst.magnitude ;
test("magnitude",dmagg.str(8),41.0.sqrt.str(8)) ;
finish
end ;
end ; -- TEST_CPX
class TEST_CPXD is
-- This class carries out tests on the complex number class CPXD. Note
-- that it is not portable for simplicity in building a local library.
-- Version 1.0 Nov 97. Copyright K Hopper, U of Waikato
-- Development History
-- -------------------
-- Date Who By Detail
-- ---- ------ ------
-- 21 Nov 97 kh Original from Sather distribution
-- 3 Apr 00 djw Extended to cover more cases
include TEST ;
main is
class_name("CPXD") ;
first : CPXD := CPXD::create(3.0d,4.0d) ;
second : CPXD := CPXD::create(5.0d,7.0d) ;
third : CPXD := CPXD::create(8.0d,11.0d) ;
sum : CPXD := first + second ;
test("plus",sum.str,third.str) ;
fourth : CPXD := CPXD::create(-2.0d,-3.0d) ;
diff : CPXD := first - second ;
test("minus",diff.str,fourth.str) ;
neg : CPXD := first + (second.negate) ;
test("negate",neg.str,fourth.str) ;
prod : CPXD := first * second ;
fifth : CPXD := CPXD::create(-13.0d,41.0d) ;
test("times",prod.str,fifth.str) ;
div : CPXD := second / first ;
sixth : CPXD := CPXD::create(1.72d,0.04d) ;
test("divide",div.re.str(8)+" "+div.im.str(8),
sixth.re.str(8)+" "+sixth.im.str(8)) ;
recip : CPXD := second * (first.reciprocal) ;
test("reciprocal",recip.re.str(8)+" "+recip.im.str(8),
sixth.re.str(8)+" "+sixth.im.str(8)) ;
dfirst : CPXD := CPXD::create(4.0d,5.0d) ;
dsecond : CPXD := dfirst.conjugate ;
test("conjugate",dsecond.str,CPXD::create(4.0d,-5.0d).str) ;
dpow : CPXD := dfirst.pow(CPXD::create(2.0d)) ;
dthird : CPXD := CPXD::create(-9.0d,40.0d) ;
test("pow",dpow.re.str(8)+" "+dpow.im.str(8),
dthird.re.str(8)+" "+dthird.im.str(8)) ;
dsqr : CPXD := dfirst.square ;
test("square",dsqr.str,dthird.str) ;
dsqrt : CPXD := dsqr.sqrt ;
test("sqrt",dsqrt.re.str(8)+" "+dsqrt.im.str(8),
dfirst.re.str(8)+" "+dfirst.im.str(8)) ;
dmag : FLTD := dfirst.magnitude_squared ;
test("magnitude_squared",dmag.str,41.0d.str) ;
dmagg : FLTD := dfirst.magnitude ;
test("magnitude",dmagg.str(8),41.0d.sqrt.str(8)) ;
finish
end ;
end ; -- TEST_CPXD
- 1.3b7 CPX and CPXD bugfixes,
Darren J Wilkinson <=