!topic LSN-1051 on the Choices for Complex Arithmetic
!reference MS-L.5;4.8
!reference MS-L.5;4.9
!reference 92-1549.a Philippe Kruchten 92-10-5
!from Ken Dritz 92-10-28
!keyword complex, imaginary
!key LSN-1051
!discussion
This LSN discusses the current state of proposals for complex arithmetic, both
as independent ISO standards and in Ada 9X. It focuses on the relative merits
of competing alternatives and makes recommendations for which I am seeking
concurrence.
1. HISTORY AND POLITICS
The long (five-year) history of the SIGAda NumWG's work on generic packages for
complex arithmetic is one of recurring compromises. It has been impossible to
satisfy everyone's conflicting goals of functionality and efficiency without
some impact on style and readability.
Although the NumWG's proposals haven't yet been accepted (or even submitted to)
WG9, a month ago they were considered "ready enough for standardization"
(in view of their anticipated submittal to the NRG in October, as the final
step before WG9) to serve as the basis for features I proposed in version 4.8
of the Numerics Annex. Perhaps naively, I assumed that the response would be
about as favorable as it was earlier for the elementary and primitive
functions, which have a similar provenance. Instead, lots of "better ideas"
came forth. It is interesting to speculate what the response would have been,
had the ancestors of the Ada 9X proposals for complex arithmetic reached the
same level of acceptance in the standardization community as their cousins. On
the other hand, it is pointless to make any assumptions in that regard, since
there is no real assurance that the NumWG's proposals would survive scrutiny by
WG9 intact.
The NRG received the NumWG's proposals on October 13. Normally, the NRG might
be expected to make very few changes in a proposal it receives from the NumWG,
especially when five years have already been spent on the effort. The NRG
wants to see all these numeric standards reach fruition without further delay.
(In addition to the need to get the standards in place, so that practitioners
can start reaping the promised benefits of portability, there is the very real
danger that the committees that have been doing the work will run out of steam
before finishing. With job transfers and even retirements, the committees have
dwindled to a bare minimum in size and only questionably reflect a consensus of
the numeric user community now.) However, at this particular NRG meeting I
played a somewhat different role from my usual one. I conveyed to the NRG the
comments of the DRs. I suggested that certain changes were expected in my
proposals for Ada 9X, and that to retain compatibility with the separate ISO
standards, the NRG should seriously consider making some changes in the
proposals now before it. Unfortunately, even this appeal failed to rally the
NRG to make any bold changes. Thus, the choice now seems to be between
adopting in the Numerics Annex a set of facilities for complex arithmetic that
is unremarkable but compatible with (what might become) an ISO standard, and a
solution that seems better suited to Ada and follows through on contemporary
recommendations but that will have the effect of delaying the work of the
NumWG and NRG even further. Neither prospect is welcome.
While it is appropriate for me to bring to the Ada 9X Numerics Annex the
collective suggestions of the numerics community, as reflected for example in
the work of the NumWG and NRG, I also have an obligation as an annex author to
do what I feel is right for Ada 9X. In the present instance, I must choose to
reject some of the advice of my colleagues on the NumWG and NRG, which has, in
some cases perhaps, been overly influenced by their Fortran background.
I now set politics behind (except for occasional comments later) and turn to
technical matters.
2. THE NUMERICS ANNEX VERSION 4.8 PROPOSAL
I repeat here the specification of Generic_Complex_Types from version 4.8 of
the Numerics Annex, to serve as the basis of further discussion. (Actually,
this is from version 4.9, where I used mixed case, as was recommended. There
were no other changes in version 4.9 except the addition of a note labeling the
specification as preliminary and subject to change.) This is referred to
subsequently as "the existing proposal," "my proposal," etc., to distinguish it
from alternatives.
with Elementary_Functions_Exceptions;
generic
type Float_Type is digits <>;
package Generic_Complex_Types is
subtype Float_Base is Float_Type'Base;
type Complex is
record
Re, Im : Float_Base;
end record;
function Real_Part (X : Complex) return Float_Base;
function Imag_Part (X : Complex) return Float_Base;
procedure Set_Real_Part (X : in out Complex; Real_Part : in Float_Base);
procedure Set_Imag_Part (X : in out Complex; Imag_Part : in Float_Base);
function Compose_From_Cartesian (Real_Part : Float_Base) return Complex;
function Compose_From_Cartesian (Real_Part : Float_Base;
Imag_Part : Float_Base) return Complex;
function Modulus (X : Complex) return Float_Base;
function "abs" (X : Complex) return Float_Base renames Modulus;
function Argument (X : Complex) return Float_Base;
function Argument (X : Complex; Cycle : Float_Base) return Float_Base;
function Compose_From_Polar (Modulus : Float_Base;
Argument : Float_Base) return Complex;
function Compose_From_Polar (Modulus : Float_Base;
Argument : Float_Base;
Cycle : Float_Base) return Complex;
function "+" (X : Complex) return Complex;
function "-" (X : Complex) return Complex;
function Conjugate (X : Complex) return Complex;
function "+" (X, Y : Complex) return Complex;
function "+" (X : Float_Base; Y : Complex) return Complex;
function "+" (X : Complex; Y : Float_Base) return Complex;
function "-" (X, Y : Complex) return Complex;
function "-" (X : Float_Base; Y : Complex) return Complex;
function "-" (X : Complex; Y : Float_Base) return Complex;
function "*" (X, Y : Complex) return Complex;
function "*" (X : Float_Base; Y : Complex) return Complex;
function "*" (X : Complex; Y : Float_Base) return Complex;
function "/" (X, Y : Complex) return Complex;
function "/" (X : Float_Base; Y : Complex) return Complex;
function "/" (X : Complex; Y : Float_Base) return Complex;
function "**" (X : Complex; N : Integer) return Complex;
Argument_Error : exception
renames Elementary_Functions_Exceptions.Argument_Error;
end Generic_Complex_Types;
There was also a specification for Generic_Complex_Elementary_Functions, which
I will discuss (only briefly) in Section 6.
The functions, procedures, and exceptions declared above are also present in
the proposed ISO standard for complex arithmetic packages for Ada 83. Their
names, and those of the formal parameters, are also the same, except that
Float_Base doesn't exist in the separate ISO proposal, and Float_Type is used
there wherever Float_Base appears above as a component, parameter, or result
subtype. Of course, this was motivated by the desire to remain as compatible
as possible with the proposed ISO standard.
(Political digression: The specification for Generic_Complex_Types in the
proposed ISO standard is about three times as long; it contains declarations
for vectors and matrices of complex components, along with relevant operations
on those vector and matrix types. There is also a parallel generic package
called Generic_Real_Types in that proposal, which provides vectors and matrices
of real components and operations thereon. (The scalar operations are not
defined because they are already available in Standard.) Yet another generic
package takes care of mixed real and complex vector and matrix operations.
Early external comments from three years ago questioned why vector and matrix
types weren't defined as a separate composable abstraction (in a separate
generic package that imports the component type as a private type). The NumWG
had an answer that did not satisfy everyone, but it is not particularly
relevant here. This is one of the main issues that is expected to be raised
again by WG9. I elected to decouple the features proposed in the Numerics
Annex from that controversy by not including vector and matrix types in any
predefined unit in the Numerics Annex, and the NumWG and NRG support me in that
decision; the less highly developed nature of the core facilities for array
manipulation in Ada, as compared to Fortran 90, was another factor in that
decision. In any case, the important point is that the Numerics Annex features
for complex arithmetic are a small subset of the features proposed in the
separate ISO standard, which therefore has a strong reason to continue with a
separate existence. Some might question the wisdom of targeting the separate
ISO standards for Ada 83 at this point in history, but that, too, is not
relevant to the discussion at hand. The version in the Numerics Annex makes
the appropriate adaptations to Ada 9X that a future Ada 9X version of the
separate ISO standards would be expected to make.)
The proposal above elicited the following reactions at the last DRs' meeting:
o The names of some operations, particularly Compose_From_Cartesian and
Compose_From_Polar, are too long. Short names are desired in expressions in
mathematical software, because that is the style of the underlying
mathematics (and also because the expressions are sometimes lengthy).
o The component names (Re, Im) of the complex type are not too attractive,
despite being short.
o The complex type should be private. (Eventually, the visibility of this
type was accepted. It has been argued that programmers should have the
aggregate notation available for constructing complex constants, as in other
languages; that higher-level abstractions that import the complex
abstraction should know the representation of Complex for efficiency
reasons; and that making the type visible was one way of specifying and
enforcing a Cartesian representation, which is necessary for well-defined
and portable accuracy requirements.)
o There are too many different ways of expressing the same thing (too much
redundancy). One can compose or decompose complex values in several
different ways, for instance.
o Some complex constants, like i and Zero, should be defined.
o The operands of the infix operators should be named Left and Right, for
uniformity with the rest of Ada.
o The complex elementary functions should be bundled with the other operations
in the complex types generic package, so as to form a complete abstraction.
o The subtype of the components of the Complex type should be called Real,
which is a better contrast to Complex than Float_Base.
The presentation on complex arithmetic motivated Kruchten, Bardin, and Mathis
to think about organizational and stylistic issues. They shared with me their
previous approaches to defining Complex, tailored slightly to the present
context.
3. ALTERNATE APPROACHES
All three of the alternate approaches had in common the goal of allowing the
programmer to write, for example, 3.0 + 5.0*i (with surrounding parentheses, if
needed), where the existing proposal requires either (3.0, 5.0) or
Compose_From_Cartesian(3.0, 5.0). All of them define the constant i. They
differ primarily in the type they give to i, and in other minor respects.
All three of the alternate approaches define Complex as a record containing two
real components, as did I and the NumWG. They all used the name Real for this
type, following through with a suggestion from the DRs' meeting. In all three
cases, Real was the name given to the generic formal parameter (instead of
Float_Type). Until recently I planned to retain Float_Type as the name of the
generic formal parameter and instead use Real for the name of the unconstrained
subtype thereof. (I have argued that it is inappropriate for the components of
the results of complex operations to be limited by arbitrary range constraints
imported with the generic formal parameter.) My present thoughts on this issue
can be found in Section 5.
Kruchten's proposal, which was posted to ada9x-mrt and which is reproduced in
Appendix A, was the least "radical," differing only modestly from the existing
proposal. Kruchten defined I as a Complex constant having value (0.0, 1.0).
Thus, 5.0*i is of type Complex and is accommodated by the mixed-mode
(real/complex) "*" operator already provided in my proposal (and, of course,
the NumWG's proposal). Then, 3.0 + 5.0*i is, of course, of type Complex and is
accommodated by the mixed-mode "+" operator, which is also already provided.
An appropriate body for the latter operator (assuming the names that Kruchten
used) is, of course,
function "+" (Left : Real; Right : Complex) return Complex is
begin
return (Left + Right.Real_Part, Right.Imag_Part);
end;
It is important, for reasons of efficiency and for other reasons that I will
cover later, not to implement this as
function "+" (Left : Real; Right : Complex) return Complex is
begin
return Make(Left) + Right;
end;
("Make" was the name Kruchten used for Compose_From_Cartesian.)
An appropriate body for the mixed-mode "*" operator used in this example is
function "*" (Left : Real; Right : Complex) return Complex is
begin
return (Left*Right.Real_Part, Left*Right.Imag_Part);
end;
which is preferable to
function "*" (Left : Real; Right : Complex) return Complex is
begin
return Make(Left)*Right;
end;
on several counts.
Unfortunately, even with the preferred implementation shown here, in the
context of 3.0 + 5.0*i a redundant (degenerate) multiplication, of 5.0 by 0.0
to obtain the real part of the intermediate result, and a redundant
(degenerate) addition, of 3.0 to 0.0 to obtain the real part of the final
result, will be performed, unless of course the compiler recognizes
opportunities for algebraic optimizations.
The second proposal (reproduced in Appendix B), by Bardin, seeks to do slightly
better by overcoming the need for such optimizations. Bardin coincidentally
defined the Complex type as a record type with component names X and Y, so that
one could extend the customary view of a complex value Z as being decomposable
as X + iY to the decomposition of a Complex object Z into the components Z.X
and Z.Y; I'll use his names for the following discussion.
Bardin defined an auxiliary pure-imaginary type, Imaginary. For reasons that
he didn't discuss (but which I'll supply later), he defined it as a private
type, rather than pursuing the more obvious alternative of visibly deriving it
from Real. He gave the full type declaration as
type Imaginary is
record
Value : Real:
end record;
rather than as
type Imaginary is new Real;
for subtle reasons that merely affect the style of what the implementer of the
complex abstraction will write in the bodies of the mixed-mode operations.
(Compare
function "+" (Left : Imaginary: Right : Complex) return Complex is
begin
return (Right.X, Left.Value + Right.Y);
end;
with
function "+" (Left : Imaginary: Right : Complex) return Complex is
begin
return (Right.X, Real(Left) + Right.Y);
end;
for example.)
Bardin then went on to define the constant i as a deferred constant,
i : constant Imaginary;
with full type declaration
i : constant Imaginary := (Value => 1.0);
Unlike Kruchten, he also defined an identical constant j, which is preferred by
engineers. And, unlike Kruchten, he made a point of departing from the
initial-capital convention for these particular names (of course, this affects
only the appearance of the specification and doesn't affect what users can do,
but it might plant a useful suggestion in the mind of a user who follows
stylistic guidelines too rigidly).
Bardin's proposal provides operators not only for mixed real/complex
operations, but also for appropriate mixed imaginary/complex and imaginary/real
operations; since Imaginary is a private type, it also must explicitly provide
the expected operations on pure Imaginary types, like addition and subtraction.
Thus, 5.0*i is of type Imaginary (in contrast to Complex in Kruchten's
proposal), the multiplication being accommodated by a mixed-mode operation for
which an appropriate body is
function "*" (Left : Real; Right : Imaginary) return Imaginary is
begin
return (Value => Left*Right.Value);
end;
Then, 3.0 + 5.0*i is, of course, of type Complex, the addition being
accommodated by a mixed-mode operation for which an appropriate body is
function "+" (Left : Real; Right : Imaginary) return Complex is
begin
return (Left, Right.Value);
end;
which, as expected, is just as efficient as Compose_From_Cartesian. (They are
not renamings of each other, unfortunately, because their right operands have
different types.)
What has Bardin achieved in contrast to the more modest proposal from Kruchten?
No degenerate arithmetic of the kind that begs for compiler optimizations is
performed in his case, either in the multiplication or the addition, for the
example given.
I did not classify the multiplication of 5.0 by 1.0, to obtain the imaginary
part of the intermediate result, as degenerate in either proposal, though in
reality it is (and, of course, it can be optimized away). The third proposal
(reproduced in Appendix C), from Mathis, seeks to eliminate even that
multiplication! Mathis defines a pure-imaginary type, Imaginary, as well; he
visibly derives it from Real. He defines i (and j) as literals of an auxiliary
(visible) enumeration type, Imaginary_Indicator. Thus, he provides an
overloaded "*" for multiplying an operand of type Real by an operand of type
Imaginary_Indicator to obtain a result of type Imaginary. The body of this
operation just returns the first operand. With appropriate inlining, but
without the need for any algebraic optimizations, Z := 3.0 + 5.0*i becomes as
efficient, in his proposal, as Z := (3.0, 5.0).
This is certainly attractive, but I think the drive for efficiency in this case
makes it impossible to do other "obvious" things. For example, although one
can declare variables of type Imaginary, one cannot assign the value i to such
a variable with an assignment as simple as Y := i. This is a reasonable thing
to expect to do, and is possible with Bardin's proposal. With Mathis's
proposal, one would have to write Y := 1.0*i. Mathis did not fully flesh out
his proposal. Perhaps he intended to include overloadings of unary "+" as
appropriate coercions from Real or Imaginary to Complex, as does Bardin; if so,
he could also include coercions from Imaginary_Indicator to Imaginary and to
Complex. (Y := +i is better than Y := 1.0*i, but not as good as Y := i. In
either case, for Complex Z, Z := +i is required.) There is no problem with
something like 3.0 + i, in Mathis's proposal, provided that an overloading of
binary "+" for Real + Imaginary_Indicator -> Complex is included. All in all,
I slightly favor Bardin's proposal because it needs to define fewer overloaded
operators, clearly allows i to be used *wherever* any pure imaginary value is
allowed, and gives the user no unnecessary or irrelevant "stuff" (e.g., useless
attributes and variables of the enumeration type Imaginary_Indicator). I,
personally, am willing to accept multiplications by 1.0 as the cost of the more
straightforward design--or, alternatively, the need for a compiler optimization
to remove them. (How many compilers already perform this optimization?)
Actually, this potential cost would be borne only by those who use the notation
x + y*i. The notation (x, y) is still available; I do not propose to take it
away. Bardin made the point that multiple styles and both long and short names
should be allowed, because the cost of allowing the user a choice is very
small. His implicit answer to the complaint that there is already too much
redundancy in the existing proposal is that redundancy is not necessarily bad,
e.g. when it supports a variety of styles all of which have their adherents. I
agree. His proposal provides even more redundancy than mine.
Recall that Bardin made the Imaginary type private, but Mathis didn't. Why
does this type have to be private? If it is visibly derived from Real, then
one could write Y := 5.0, where Y is of type Imaginary; furthermore,
expressions like X + 5.0 become ambiguous (where X is of type Real). (Is the
"+" operator the one for Real + Real -> Real, or is it the one for Real +
Imaginary -> Complex?) We really don't *want* to allow a simple real literal
to appear anywhere an imaginary value is expected or allowed. We want to force
the notation 5.0*i (for example) where an imaginary literal is desired as a
.
Of course, one can deny the implicit conversion of a real literal to an
Imaginary value without making the Imaginary type private. For example, one
could use the record type that Bardin used for Imaginary, but simply make it
visible. But this allows Y := (Value => 5.0) and X + (Value => 5.0), which is
not the kind of notation we wish to encourage. There is insufficient advantage
in making Imaginary visible in this way. It is best to leave it private.
The idea of defining the Complex type to contain one component of type Real and
one of type Imaginary, rather than two of type Real, comes to mind. If this
were done, the Imag_Part function should also be redefined to yield an
Imaginary result, and Compose_From_Cartesian should take an Imaginary second
argument. Then, if Y were a variable of type Imaginary, one could write
Y := Imag_Part(Z), in contrast to Y := i*Imag_Part(Z), which all three
alternatives presently require. However, defining Complex to have one Real and
one Imaginary component would no longer allow complex "literals" to be written
as (3.0, 5.0); they would have to be written as (3.0, 5.0*i) or as 3.0 + 5.0*i.
Furthermore, defining it that way departs too much from the usual convention in
mathematics, and in other programming languages, whereby the real and imaginary
parts of a complex quantity are real. Lastly, about the only choice available
for the I/O format of a complex value on an external medium (see Section 7) is
(3.0, 5.0), if compatibility with existing data files is to be preserved. All
of these reasons point to representing Complex as a pair of Reals.
Armed with the Lexington reaction to the proposal in the Numerics Annex, and
the alternate proposals, I went to the NRG meeting hoping to secure some
changes in the proposed ISO standards, so that I could modify my proposal for
the Numerics Annex and still remain compatible with the secondary standards.
4. THE NRG MEETING
The following decisions were reached at the NRG meeting in London, October
13-14, and apply to its proposals for complex arithmetic. I report them in the
order in which I reported the "problems" with the Numerics Annex proposal in
Section 2, above. My response to the NRG's decisions is presented in Section
5.
o A few other names were considered for Compose_From_Cartesian and
Compose_From_Polar, but they were rejected. The NRG found disfavor with
Cartesian and Polar, since those names would, by usual conventions, be taken
to imply the form of the result, not the source. Of course, the form of the
result is always Cartesian, and what we need is a way to distinguish the
form of the source. Complex is out; it is already used (for the type), and
anyway is only one of two needed names (by itself, it doesn't distinguish
the form of the source). No further candidates were considered, and the NRG
wished to stick with its two long names.
o The component names Re and Im were retained for the Complex type. An
element of the committee's reasoning was that short names are required so
that the bodies of operations like "*" and "/", which are non-trivial, can
be made readable; component names like Real_Part and Imaginary_Part (e.g.,
from Mathis's proposal) would wreak havoc with the source formatting of
lengthy expressions. X and Y (from Bardin's proposal) were not mnemonic
enough to suit the committee. It was pointed out that the names Re and Im
are often used in mathematics as operators on complex values. To capitalize
on that, the committee decided not only to retain Re and Im as component
names, but also to change the names of the selector functions Real_Part and
Imag_Part to Re and Im; the names of the procedures Set_Real_Part and
Set_Imag_Part to Set_Re and Set_Im; and the names of the Real_Part and
Imag_Part parameters of Compose_From_Cartesian, Set_Real_Part, and
Set_Imag_Part to Re and Im. Thus, Re and Im are now used uniformly in all
contexts to denote the real and imaginary parts of a complex value or
object.
o The NRG was not asked to make the Complex type private; they would surely
not have done so. This does not appear to be a current issue.
o The issue of redundancy (too many ways of accomplishing the same thing) was
considered hand in hand with using the style of the three alternate
proposals (i.e. 3.0 + 5.0*i) as a means of denoting the composition of a
complex value from a pair of reals. In this specific case, it was
recognized that Compose_From_Cartesian cannot be dispensed with entirely,
since the one-parameter form (for promotion from real to complex) needs to
be available as a potential generic actual parameter in an instantiation of
a higher-level abstraction, for example a linear equation solver that works
for either real or complex types.
o Constants like i and Zero will not be included in Generic_Complex_Types.
A NumWG member on the NRG reminded me that the NumWG's proposal originally
included the constant i, but that it was dropped when it was realized that
such an item poses problems if multiple instantiations of the generic
package are used simultaneously (the name i would not be directly visible).
This argument fails to acknowledge that the same problem exists for the type
name Complex itself, and that it is a familiar problem that users of
generics have ways to accommodate. (This problem is further discussed in
Section 9.) There were also concerns that the name i might conflict with
common uses of I as a loop index. Implicit in this decision is the decision
not to define a pure-imaginary type, which requires the constant i in order
to be useful.
o The parameter names of the infix operators will not be changed to Left and
Right. (It was assumed that the request to use Left and Right applied only
to predefined operators that are being overloaded in this generic package,
not to any other function or procedure being declared therein.) Cited as
rationale for keeping X and Y were the following reasons: (1) X and Y are
"familiar mathematical notation" (I personally feel this is bogus, since
mathematicians do not normally think in terms of names for the operands of
the familiar mathematical operators); (2) it seems somewhat arbitrary to use
Left and Right for operators (including "abs") but not functions like
Modulus or Conjugate (the renaming of Modulus as "abs" would also rename the
parameter); (3) the committee didn't like the precedent set in Ada 83, and
would change it if it could (John Barnes explicitly voiced this opinion);
(4) even the RM doesn't insist on using Left and Right for the predefined
operators everywhere (cf. an example in RM 12.1.3); (5) bodies for some of
the more complicated operators, like complex multiplication and division,
would have awkward and hard-to-format expressions with parameter names as
long, even, as Left and Right (clearly this affects only implementers, not
users, and should not have been given so much weight); (6) parameter names
of X and Y were accepted for "**" in Generic_Elementary_Functions (they were
approved by WG9, SC22, and JTC1, and are about to be cast in concrete in the
final edit of DIS 11430; reviewers have not previously complained about the
parameters of "**" in the version of Generic_Elementary_Functions in the
Numerics Annex); (7) several vendors have already implemented
Generic_Elementary_Functions as proposed, using the parameter names X and Y
for "**".
o The complex elementary functions should remain separate from the generic
package defining the complex type and arithmetic operations. Not all users
of complex arithmetic will require the complex elementary functions.
Furthermore, the ISO version of the generic complex types package, which is
already huge because of the vector and matrix types and operations that it
contains, would become just too unwieldy. We do not need to go out of our
way to ensure that all operations of Complex are primitive, since there is
little value to deriving from Complex.
o The proposals before the NRG, based as they are on Ada 83, do not declare
the unconstrained subtype, Float_Base, of the generic formal parameter
Float_Type. The NRG did, however, consider changing the name of the generic
formal parameter itself to Real, but rejected it. All of the NRG's proposed
secondary ISO standards are consistent in using the name Float_Type for a
generic formal parameter, and the NRG felt that this consistency should be
maintained. However, it is not possible to change Float_Type to Real
everywhere, because in the proposed ISO standard for complex types, there is
the aforementioned generic package parallel to Generic_Complex_Types called
Generic_Real_Types, which defines and exports a subtype of the generic
formal parameter called Real, intended to be used by the user. Another
reason to avoid globally changing Float_Type to Real is that it is used, of
course, as the name of the generic formal parameter of
Generic_Elementary_Functions too, and that standard has been approved up
through JTC1, is in final editing, and has been implemented by several
vendors, as has already been pointed out.
5. SYNTHESIZING AN AGREEABLE PROPOSAL FOR COMPLEX TYPES IN ADA 9X
On my return from the NRG meeting, a copy of a paper by W. Kahan and J. Thomas,
"Augmenting a Programming Language with Complex Arithmetic," was waiting for
me. This paper makes a strong case for a pure-imaginary type.
The paper first points out that prematurely promoting a real quantity to
complex (with an imaginary part of zero) is not only wasteful, given that
multiplications by zero or additions to zero may ensue, but in IEEE
environments is actually harmful. For example, if in R*(X, Y) a real R is
first coerced to (R, 0.0), and the multiplication is then performed as a
complex multiplication, yielding (R*X - 0.0*Y, R*Y + 0.0*X), then
3.0*(Infinity, 5.0) will not yield (Infinity, 15.0) but rather (Infinity, NaN).
This is because, on IEEE systems, zero times infinity is an invalid operation
and yields a NaN ("Not-a-Number"), if not trapped. (If invalid operations are
trapped, the premature coercion is even worse, because then we may get no
result.) It is for this reason that I emphasized in Section 3 that the
mixed-mode multiplication R*(X, Y) should be performed as (R*X, R*Y), without
coercing R to complex first. Similarly, if R is first coerced to complex in
R + (X, Y), and the addition is then performed as a complex addition, yielding
(R + X, 0.0 + Y), then 3.0 + (5.0, -0.0) will not yield (8.0, -0.0) but rather
(8.0, +0.0). The reason for this is that, on IEEE systems, addition of two
quantities differing only in sign yields +0.0 (unless the rounding mode is
"round towards minus infinity"). This was the reason for emphasizing, in
Section 3, that the mixed-mode addition R + (X, Y) should be performed as
(R + X, Y), without coercing R to complex first.
Similar problems involving mixed imaginary/complex operations can obviously be
displayed. However, if a pure-imaginary type is not offered, the programmer
has no choice but to use, in place of an imaginary value, a complex value with
a zero real part, which represents an explicit, premature coercion of the
imaginary quantity to complex. Problems like those demonstrated above are
unavoidable in that case. For example, the multiplication of two
pure-imaginary quantities, one of which is infinite and the other nonzero,
should yield an infinite real result, but if the multiplication is performed
as a complex multiplication then the imaginary part of the result will be a NaN
instead of zero. By offering a pure-imaginary type, we make it possible for
the programmer to exercise equal control over the premature coercion of either
pure-real or pure-imaginary values to the complex type.
I present below my response to the NRG actions and the suggestions I have
received. It features the inclusion of a pure-imaginary type, as proposed by
Bardin and Mathis (but following more closely the details of Bardin's proposal)
and discussed by Kahan and Thomas. Again, I discuss the issues in the same
order as was used in Sections 2 and 4.
o Regarding the long names Compose_From_Cartesian and Compose_From_Polar, I
propose to retain those names and to provide renamings with shorter names,
giving the user a choice of styles. I would rename the former with Cmplx,
which is short, sufficiently mnemonic, and the name used by Fortran for the
same function. Since the representation of a complex value is Cartesian, I
see no need to reiterate that in the short form. Many applications will
have no need for polar components, and for them to use names containing or
implying the word Cartesian is plainly redundant. For those that do have a
need for constructing (Cartesian) complex values from polar components,
Polar_To_Cmplx seems like a suitably short alternative renaming of
Compose_From_Polar. It uses an accepted style for coercion functions (e.g.,
a_To_b) and shares the fragment Cmplx with the other coercion function.
Admittedly it is only four letters shorter, but it is hard to do any better
while saying something about both the source and target types. I will
entertain suggestions to omit Polar_To_Complx entirely and have a short
renaming (as Cmplx) only for Compose_From_Cartesian. I do not like the idea
of retaining only Cmplx and Polar_To_Cmplx, discarding
Compose_From_Cartesian and Compose_From_Polar, since the explicit
renaming of Compose_From_Cartesian as Cmplx helps to document, in the
specification, semantics that the name Cmplx doesn't convey all by itself.
o Since users cannot shorten component names of record types by renaming, and
since component selection will be done by some users using expanded names,
rather than function calls, it is important to have short names for the
components of the Complex type, and I think Re and Im are appropriate. I
will retain the relatively long names Real_Part and Imag_Part (but not
Imaginary_Part) for the selector functions, but I will rename them as Re and
Im for those who prefer the functional style and need short names.
(Imag_Part allows for better text layout than Imaginary_Part when it is in
apposition to Real_Part, for those who are fanatic about lining up their
source. Furthermore, Imag is familiar from Fortran--where, of course, it
has to be called Aimag!) Although Set_Real_Part and Set_Imag_Part are not
expected to be used often enough to warrant shorter names, I will
nevertheless provide renamings of them as Set_Re and Set_Im, for uniformity;
the renamed versions will also have Re and Im, instead of Real_Part and
Imag_Part, as parameter names. The renaming of Compose_From_Cartesian as
Cmplx, discussed above, will also involve a renaming of the parameters to Re
and Im.
o The Complex type continues to be a visible type.
o I have increased the redundancy, rather than reducing it. The minutes of
the last DRs' meeting reflect the fact that there was discussion about the
redundancy (and I think it was mostly Dewar's), but there was no motion or
vote on it. Since the alternative proposals have added yet another way to
create a complex value from its Cartesian components, I judge that some
reviewers think that redundancy is OK when it supports a variety of
different styles, all of which may be preferred by some user or other.
Thus, I have included a declaration of a pure-imaginary type, along with
appropriate operations on the type and on mixtures of that type and either
the real or the complex type. (My primary motivation for doing so was not
simply to offer more choices, however; I was strongly influenced by the
Kahan and Thomas paper. In contrast to other languages in which they might
be interested, in Ada the inclusion of such a type can be accomplished quite
elegantly, and it would be nice to see them admire that.) I have also
included a broader range of renamings for those who prefer short names,
while retaining longer names for those who prefer them, as noted above.
o The constants i and j are now defined, as in Bardin's proposal. I have not
defined Zero and One, since those can be easily enough defined by the user.
(The user cannot define i, on the other hand, because it is a constant of a
private type.) I have resisted using longer names than i and j, and thus
have not addressed the problem of conflict with loop indices, since nothing
else has the desired notational elegance of 3.0 + 5.0*i.
o Motivated mostly by the status achieved by X and Y as parameter names for
"**" in DIS 11430, and by the desire not to change that DIS or to impose a
change on the vendors who currently support that DIS (or, alternatively, not
to introduce a gratuitous difference between the Ada 9X version of
Generic_Elementary_Functions and that DIS), I will keep the names currently
proposed for the parameters of the overloaded predefined arithmetic
operators. These are X and Y in most cases (or just X), and X and N in the
case of complex exponentiation by an integer. I note also that parameter
names other than Left and Right have been used for overloadings of "*" and
"+" in the string-handling facilities of the Information Systems Annex. I
think it is appropriate for application-oriented extensions in the form of
overloadings of the predefined operators to use parameter names relevant to
the application area. (If I am overruled by an excessively rigid and
paramount stylistic requirement to use Left and Right, then the adverse
consequences for DIS 11430, or the vendors who have already implemented it,
must be accepted.)
o Generic_Complex_Types and Generic_Complex_Elementary_Functions remain
separate generic packages. Being able to write Rho * Exp(i*Theta) as an
alternate form of Compose_From_Polar(Rho, Theta) does not hinge on merging
the two generic packages. It does, however, hinge on providing in
Generic_Complex_Elementary_Functions a version of Exp for pure-imaginary
arguments; otherwise, one would have to write Rho * Exp((0.0, Theta) or
Rho * Exp(0.0 + i*Theta)). I return to this discussion later.
o Since the version of Generic_Complex_Types in the Numerics Annex defines the
unconstrained subtype Float_Base of the generic formal parameter Float_Type,
it seems easy to arrange for Real to be used for the subtype of the
components of Complex and the formal parameters and results of the included
operations without changing the name Float_Type; one merely substitutes Real
for Float_Base. Actually, this corrects a slight problem that went
undetected in versions 4.8 and 4.9 of the Numerics Annex: instantiations of
both Generic_Elementary_Functions and Generic_Complex_Types export the name
Float_Base, which is consequently not directly visible when both are used
simultaneously. While Float_Base is mostly for convenience within the
generic package specifications, and is not particularly intended for the
user to use, it nevertheless seems best to avoid the cancellation effect.
But, if this change were made, then instantiations of Generic_Complex_Types
and the NumWG's/NRG's Generic_Real_Types would both export the subtype name
Real, which we must avoid. In thinking about various other alternatives, I
reevaluated my own supposition that using the unconstrained subtype of
Float_Type for the components of the complex type and certain subprogram
parameters and results was appropriate. I now conclude that it is not. It
was intended to allow arithmetic operations and elementary functions
yielding Complex results to escape the limiting effect of range constraints
on (the components of) the result, inherited from Float_Type, in analogy to
the predefined arithmetic operations and to the real elementary functions
(which are limited only by the capabilities of the hardware and not by the
range constraints of any subtype). But there is a major difference in the
case of the composite type, Complex: with the existing proposal, there is no
way for a user to impose range constraints on the components of an object of
type Complex. It seems reasonable to want to do that. For example, it
should be possible to define a type (and objects of the type) representing
a particular rectangular subset of the complex plane, and to be notified by
the raising of Constraint_Error whenever the result of a calculation lies
outside its boundary. The important thing is to give the user the ability
to define such a type, if needed, or to opt instead for no limitations on
the range of the components except for those imposed by the hardware. The
user exercises this choice by deciding whether to instantiate
Generic_Complex_Types with a constrained or an unconstrained floating-point
subtype. The necessary control cannot be exercised in any other way. Thus,
I have abandoned the idea of defining the unconstrained subtype of the
generic formal parameter in Generic_Complex_Types. It is only one small
further step to change the name of that parameter to Real; in so doing, I am
explicitly departing from the naming convention preferred by the NumWG and
NRG, rationalizing that other considerations are more important. It should
be noted that this approach is identical to what Kruchten and Bardin
instinctively did, and it may be what Mathis had in mind for the proposal
that he only sketched.
Taking all these considerations together, but omitting a few of Bardin's
mixed-mode operators, I come up with the following:
with Elementary_Functions_Exceptions;
generic
type Real is digits <>;
package Generic_Complex_Types is
-- ***********************************************************************
-- * *
-- * Types and subtypes *
-- * *
-- ***********************************************************************
type Imaginary is private; -- to preclude implicit conversion
-- of real literals to Imaginary
i : constant Imaginary;
j : constant Imaginary;
type Complex is
record
Re, Im : Real;
end record;
-- ***********************************************************************
-- * *
-- * Selector functions for extracting real or imaginary parts *
-- * of complex values *
-- * *
-- ***********************************************************************
function Real_Part (X : Complex) return Real;
function Imag_Part (X : Complex) return Real;
function Re (X : Complex) return Real renames Real_Part;
function Im (X : Complex) return Real renames Imag_Part;
-- ***********************************************************************
-- * *
-- * Procedures for replacing real or imaginary parts *
-- * of complex values, leaving the other part untouched *
-- * *
-- ***********************************************************************
procedure Set_Real_Part (X : in out Complex;
Real_Part : in Real);
procedure Set_Imag_Part (X : in out Complex;
Imag_Part : in Real);
procedure Set_Re (X : in out Complex;
Re : in Real) renames Set_Real_Part;
procedure Set_Im (X : in out Complex;
Im : in Real) renames Set_Imag_Part;
-- ***********************************************************************
-- * *
-- * Constructor functions for creating a complex value *
-- * from its real and imaginary parts *
-- * *
-- ***********************************************************************
function Compose_From_Cartesian (Real_Part : Real) return Complex;
function Compose_From_Cartesian (Real_Part : Real;
Imag_Part : Real) return Complex;
function Cmplx (Re : Real) return Complex renames Compose_From_Cartesian;
function Cmplx (Re : Real;
Im : Real) return Complex renames Compose_From_Cartesian;
-- ***********************************************************************
-- * *
-- * Functions to obtain components of a complex value in polar form *
-- * *
-- ***********************************************************************
function Modulus (X : Complex) return Real;
function Argument (X : Complex) return Real;
function Argument (X : Complex; Cycle : Real) return Real;
function "abs" (X : Complex) return Real renames Modulus;
-- ***********************************************************************
-- * *
-- * Constructor functions for creating a complex value *
-- * from polar components *
-- * *
-- ***********************************************************************
function Compose_From_Polar (Modulus : Real;
Argument : Real) return Complex;
function Compose_From_Polar (Modulus : Real;
Argument : Real;
Cycle : Real) return Complex;
function Polar_To_Cmplx (Modulus : Real;
Argument : Real) return Complex
renames Compose_From_Polar;
function Polar_To_Cmplx (Modulus : Real;
Argument : Real;
Cycle : Real) return Complex
renames Compose_From_Polar;
-- ***********************************************************************
-- * *
-- * Complex conjugation *
-- * *
-- ***********************************************************************
function Conjugate (X : Complex) return Complex;
function Conj (X : Complex) return Complex renames Conjugate;
-- ***********************************************************************
-- * *
-- * Pure-imaginary arithmetic operations *
-- * *
-- ***********************************************************************
function "+" (X : Imaginary) return Imaginary;
function "-" (X : Imaginary) return Imaginary;
function "+" (X, Y : Imaginary) return Imaginary;
function "-" (X, Y : Imaginary) return Imaginary;
-- ***********************************************************************
-- * *
-- * Pure-complex arithmetic operations *
-- * *
-- ***********************************************************************
function "+" (X : Complex) return Complex;
function "-" (X : Complex) return Complex;
function "+" (X, Y : Complex) return Complex;
function "-" (X, Y : Complex) return Complex;
function "*" (X, Y : Complex) return Complex;
function "/" (X, Y : Complex) return Complex;
-- ***********************************************************************
-- * *
-- * Mixed real/imaginary arithmetic operations *
-- * *
-- ***********************************************************************
function "*" (X, Y : Imaginary) return Real;
function "*" (X : Imaginary; Y : Real) return Imaginary;
function "*" (X : Real; Y : Imaginary) return Imaginary;
function "/" (X, Y : Imaginary) return Real;
function "/" (X : Imaginary; Y : Real) return Imaginary;
function "/" (X : Real; Y : Imaginary) return Imaginary;
-- ***********************************************************************
-- * *
-- * Mixed real/complex arithmetic operations *
-- * *
-- ***********************************************************************
function "+" (X : Real; Y : Complex) return Complex;
function "+" (X : Complex; Y : Real) return Complex;
function "-" (X : Real; Y : Complex) return Complex;
function "-" (X : Complex; Y : Real) return Complex;
function "*" (X : Real; Y : Complex) return Complex;
function "*" (X : Complex; Y : Real) return Complex;
function "/" (X : Real; Y : Complex) return Complex;
function "/" (X : Complex; Y : Real) return Complex;
-- ***********************************************************************
-- * *
-- * Mixed imaginary/complex arithmetic operations *
-- * *
-- ***********************************************************************
function "+" (X : Imaginary; Y : Complex) return Complex;
function "+" (X : Complex; Y : Imaginary) return Complex;
function "-" (X : Imaginary; Y : Complex) return Complex;
function "-" (X : Complex; Y : Imaginary) return Complex;
function "*" (X : Imaginary; Y : Complex) return Complex;
function "*" (X : Complex; Y : Imaginary) return Complex;
function "/" (X : Imaginary; Y : Complex) return Complex;
function "/" (X : Complex; Y : Imaginary) return Complex;
-- ***********************************************************************
-- * *
-- * Mixed real/imaginary/complex arithmetic operations *
-- * *
-- ***********************************************************************
function "+" (X : Real; Y : Imaginary) return Complex;
function "+" (X : Imaginary; Y : Real) return Complex;
function "-" (X : Real; Y : Imaginary) return Complex;
function "-" (X : Imaginary; Y : Real) return Complex;
-- ***********************************************************************
-- * *
-- * Exponentiation by an integer (other exponenetiation operations *
-- * are in Generic_Complex_Elementary_Functions) *
-- * *
-- ***********************************************************************
function "**" (X : Imaginary; N : Integer) return Complex;
function "**" (X : Complex; N : Integer) return Complex;
-- ***********************************************************************
-- * *
-- * Exceptions *
-- * *
-- ***********************************************************************
Argument_Error : exception
renames Elementary_Functions_Exceptions.Argument_Error;
private
-- ***********************************************************************
-- * *
-- * Full type declarations of private types and deferred constants *
-- * *
-- ***********************************************************************
type Imaginary is
record
Value : Real;
end record;
i : constant Imaginary := (Value => 1.0);
j : constant Imaginary := (Value => 1.0);
end Generic_Complex_Types;
I have omitted the following operators that were included by Bardin (see
Appendix B): +Real -> Complex, -Real -> Complex, +Imaginary -> Complex,
-Imaginary -> Complex, Real*Imaginary -> Complex, Imaginary*Real -> Complex.
Bardin wanted the operators that promote pure-reals to complex so that one
could write assignments that initialize complex variables to pure-real values,
as in Z := +1.0. Unfortunately, inclusion of these operators leads to an
ambiguity: in -X + Y*i, the "+" could be either Real + Imaginary -> Complex or
Complex + Imaginary -> Complex. Bardin wanted the operators that promote
pure-imaginaries to complex so that one could write assignments that initialize
complex variables to pure-imaginary values, as in Z := -i. Unfortunately,
inclusion of these operators leads to the following ambiguity: in -i + X, the
"+" could be either Imaginary + Real -> Complex or Complex + Real -> Complex.
Finally, Bardin wanted the operators that promote products of reals and
imaginaries directly to complex so that one could write other assignments that
initialize complex variables to pure-imaginary values, as in Z := 5.0*i.
Inclusion of these operators leads, unfortunately, to yet another ambiguity,
about which Bardin seemed to be aware: in X + Y*i, the "+" could be either Real
+ Imaginary -> Complex or Real + Complex -> Complex. By precluding the
ambiguities, we are requiring that the examples of the simple assignments cited
above be written in the alternate forms Z := 1.0 + 0.0*i, Z := 0.0 - 1.0*i, and
Z := 0.0 + 5.0*i, which actually "look" more type-safe and are no less
efficient (if the compiler inlines and optimizes multiplication by one). Of
course, the programmer can use aggregate notation, if desired.
Also, I have not followed Bardin's suggestion of using "not" to denote complex
conjugation, since it is not customary mathematical notation (i.e., on
encountering a use of this operator to denote conjugation, the reader of a
program, other than the author, will likely not understand it).
6. COMPLEX ELEMENTARY FUNCTIONS REVISITED
The paper by Kahan and Thomas discusses the value of having versions of some of
the complex elementary functions for pure-imaginary arguments. Mostly this is
a matter of efficiency; i.e., when Z = (0.0, Y) is pure-imaginary, instead of
computing Sin(Z) as (Sin(Z.Re)*Cosh(Z.Im), Cos(Z.Re)*Sinh(Z.Im)), it is less
costly to just compute it as (0.0, Sinh(Z.Im)). I believe that following this
suggestion to the maximum extent possible has only marginal gain, since most of
the unnecessary computation can be avoided merely by testing for degenerate
cases. However, I propose that the function
function Exp (X : Imaginary) return Complex;
be added to Complex_Generic_Elementary_Functions. The efficiency gain is not
the real motivation here, since it is negligible. Rather, inclusion of this
overloading of Exp allows one to construct (Cartesian) complex values from
their polar components Rho and Theta by writing Rho * Exp(i*Theta). The
efficiency of this should be comparable to the efficiency of
Polar_To_Cmplx(Rho, Theta), particularly if the compiler optimizes
multiplication by one.
7. COMPLEX I/O
Versions 4.8 and 4.9 of the Numerics Annex neglected to include facilities for
complex I/O. It is curious that the NumWG has not thought at all about complex
I/O, although most vendors who provide complex arithmetic packages find it
necessary to offer some capability for input and output of complex values.
Kruchten (Appendix A) provided in his version of Generic_Complex_Types for Get
and Put conversions between String and Complex in the form of generic Get and
Put procedures parameterized by Get and Put for reals. He did not address the
problem, associated with his Put, of allocating space within the given string
between the real and imaginary parts (the Put procedure for reals fills the
string it is given). Instead, I define a generic package Generic_Complex_IO
that exports, like the generic packages nested within Text_IO, a full set of
Get and Put procedures. Like Generic_Complex_Elementary_Functions, it
is a library unit that imports an instantiation of Generic_Complex_Types; it
also imports an instantiation of Text_IO.Float_IO (for the same floating-point
type), from which it obtains the necessary building blocks. Its specification
is as follows:
with Text_IO; use Text_IO;
with Generic_Complex_Types;
generic
with package Complex_Types is new Generic_Complex_Types (<>);
with package Real_IO is new Float_IO (<>);
package Generic_Complex_IO is
use Complex_Types;
Default_Fore : Field := Real_IO.Default_Fore;
Default_Aft : Field := Real_IO.Default_Aft;
Default_Exp : Field := Real_IO.Default_Exp;
procedure Get (File : in File_Type;
Item : out Complex;
Width : in Field := 0);
procedure Get (Item : out Complex;
Width : in Field := 0);
procedure Put (File : in File_Type;
Item : in Complex;
Fore : in Field := Default_Fore;
Aft : in Field := Default_Aft;
Exp : in Field := Default_Exp);
procedure Put (Item : in Complex;
Fore : in Field := Default_Fore;
Aft : in Field := Default_Aft;
Exp : in Field := Default_Exp);
procedure Get (From : in String;
Item : out Complex:
Last : out Positive);
procedure Put (To : out String;
Item : in Complex;
Aft : in Field := Default_Aft;
Exp : in Field := Default_Exp);
end Generic_Complex_IO;
The format used for output of a complex value is that of an Ada aggregate. (I
did not think it necessary to offer all the choices enumerated by Kruchten.)
The values given or defaulted for Fore, Aft, and Exp are used for both the real
part and the imaginary part. For output to a file, the exact format is a left
parenthesis, followed by the output of the real part (using the Put procedure
for reals), followed by a comma and a space, followed by the output of the
imaginary part (using the Put procedure for reals), followed by a right
parenthesis. For output to a string, the left parenthesis, real part (with a
Fore of 0), and comma are left-adjusted in the given string without intervening
blanks, and the imaginary part (with a Fore of 0) and right parenthesis are
right-adjusted in the given string without intervening blanks; any blanks
necessary to fill the given string are placed between the comma and the
imaginary part. (Other formats can be imagined, but it is not clear that they
are worth the additional complexity in description and implementation.) On
input, the left parenthesis, real part, comma, imaginary part, and right
parenthesis are expected, with arbitrary leading blanks before each part;
exactly Width characters are consumed if Width is specified and nonzero.
Exceptional conditions are reported using the exceptions defined in Text_IO; in
particular, Data_Error is raised if input is ill-formed, and Layout_Error is
raised by an attempt to Put too many characters to a string.
8. (FRAGMENTS OF) AN EXAMPLE
I examined a Fortran application of complex arithmetic that solves the
three-dimensional Ginzburg-Landau equations for a superconducting cube in a
fixed external magnetic field. The revised proposal for Ada 9X meets the needs
of this application quite well, especially insofar as the application makes
moderately heavy use of exponentials of pure-imaginary values. Thus, in the
update of the Fourier transform of the vector potential throughout the cube for
a new time step, three intermediate quantities are first computed (inside the
main loop):
x4a = 2.d0*zi*dimag(ps(ix,iy,iz)*dconjg(ux(ix,iy,iz)*ps(ix+1,iy,iz)))
x4b = 2.d0*zi*dimag(ps(ix,iy,iz)*dconjg(ux(ix,iy,iz)*ps(ix,iy+1,iz)))
x4c = 2.d0*zi*dimag(ps(ix,iy,iz)*dconjg(ux(ix,iy,iz)*ps(ix,iy,iz+1)))
The variables x4a, x4b, x4c, zi and the arrays ps and ux are (double precision)
complex; zi is a "constant" (in a common block) initialized to (0.d0, 1.d0),
that is, the constant i. The version of Fortran used went beyond the standard
in offering a double-precision complex type (which explains the nonstandard
intrinsic function names dimag and dconjg). It is clear that the three
assignments give pure-imaginary values to x4a, x4b, and x4c. (Although time is
wasted on computing the real part of the temporary value passed do dimag,
presumably the form of the formula matches the physics or mathematics well, so
I won't criticize it.) The obvious translation to Ada 9X would be
x4a := 2.0*i*Im(ps(ix,iy,iz)*Conj(ux(ix,iy,iz)*ps(ix+1,iy,iz)))
x4b := 2.0*i*Im(ps(ix,iy,iz)*Conj(ux(ix,iy,iz)*ps(ix,iy+1,iz)))
x4c := 2.0*i*Im(ps(ix,iy,iz)*Conj(ux(ix,iy,iz)*ps(ix,iy,iz+1)))
where x4a, x4b, and x4c are of type Imaginary and ps and ux are arrays of
Complex. Other intermediate quantities are computed as well, for example x1a
and x1b (these are also pure-imaginary). Eventually, the new value for the
Fourier transform of the x-component of the vector potential at the point
(ix,iy,iz), unx(ix,iy,iz), is computed:
unx(ix,iy,iz) = ux(ix,iy,iz)*exp(coeff*(x4a - x1a/hlx2 - x1b/hlz2))
(The other two components, uny and unz, are computed in like fashion, and then
the arrays ux, uy, and uz are updated from unx, uny, and unz.) This line is
converted to Ada just by changing the assignment symbol, assuming Exp for
pure-imaginary arguments is added to Generic_Complex_Elementary_Functions, as I
have proposed.
Very few complex constants are used in this example. In fact, only zero and i
are used, and they are assigned to variables in common (zzero and zi,
respectively). Nor is the Fortran intrinsic function Cmplx used at all. For
these reasons, plus the fact that Fortran defines neither a pure-imaginary type
nor adding operators for combining a pure-real and a pure-imaginary value into
a complex value, this example affords no real opportunity to evaluate the
utility of permitting X + Y*i as a means for constructing a complex value.
9. FEEDBACK SOUGHT
I strongly urge reviewers to consider my revised proposal and post comments to
ada9x-mrt. If at all possible, I would like to determine the consensus of
opinion before the WG9 meeting. Please do keep in mind, however, that the
history of work on the standards for complex arithmetic has revealed that it
seems impossible to please everyone, and that some compromise is necessary.
I am aware that the names i and j of the constants exported by instantiations
of Generic_Complex_Types pose problems when multiple instantiations (e.g., for
different precisions) are used simultaneously: the names are not directly
visible. Whereas the strategies for dealing with this problem (i.e., renaming
or the use of expanded names) in the case of the type names Complex and
Imaginary are also available for the constants i and j, of course, their
application affects expressions and not just declarations; those strategies
seem less satisfactory in the case of i and j. I see no practical solution to
this problem, but am open to suggestions if it is thought serious enough to
solve. I don't think that the obvious approach of making i and j functions is
acceptable on efficiency grounds. The only other solution I can foresee is to
define i and j as literals of an enumeration type (e.g., Mathis's
Imaginary_Indicator), but in a separate, non-generic package that the user
would also have to use (only once, regardless of the number of instantiations
of Generic_Complex_Types). Operations like Real*Imaginary_Indicator ->
Imaginary would be defined in Generic_Complex_Types. But having both to use
the non-generic package and to instantiate Generic_Complex_Types seems too
cumbersome for the user, especially in the common case of a single
instantiation; besides, it is subject to the criticisms that I leveled on the
enumeration type in Section 3.
APPENDIX A (Alternative proposal of Kruchten)
topic Small style issues on package Generic_Complex_Types
reference MS-L.5; 4.8
from Philippe Kruchten 92-10-05
keywords complex, style
<>
discussion
Based on our discussion last week in Boston, I "merged" my own package complex
with the current proposal, making a certain number of changes. At the end, you
will see the resulting package specification. Most of it, I am afraid, is
a matter of personal taste. :-)
1. Capitalization is done according to the new Ada 9X rule
2. Use Real instead of Float_Base
Float seems to me related to the representation in machine; REAL is
the effective counterpart of COMPLEX.
The suffix _Type is not necessary and inconsistent: you do not
have Complex_Type, for instance.
3. Parameters for operators are named according to LRM rules (as defined
in Standard):
Instead of X, Y, use Left and Right.
4. Other arbitrary complex parameters are named C rather than X, to remind
of their type
5. Introduction of a constant I to allow an alternate style for
expressing complex numbers:
C1 := 1.0 + 2.0 * I;
C1 := R1 + R2 * I;
6. Introduction of a constant Zero, for better legibility:
C1 := Zero;
instead of
C1 := (0.0, 0.0);
or
C1 := Compose_From_Cartesian (0.0); -- :-(
7. Replace Compose_From_Cartesian by Make and Compose_From_Polar by make_polar
because they are too long to be practical.
Alternative could be simply Cartesian and Polar, as in
C1 := Cartesian (R1, R2);
C2 := Polar (Modulus, Argument);
If this does not read well one can always qualify:
C1 := Complex'(Cartesian (R1, R2));
C2 := Complex'(Polar (Modulus, Argument));
People like me who usually use prefix notation (and minimize use
clauses) would prefix by the package name (assume Cplx is the instance name):
C1 := Cplx.Cartesian (R1, R2);
C2 := Cplx.Polar (Modulus, Argument);
This works also for Cplx.I and Cplx.Zero ...
8. Introduce the generic procedure "Put to a string", which is intended to
work like the corresponding Text_io.Float_Io.Put, with similar rules and same
exceptions.
If necessary only one display style need be implemented (I guess it
would be the Fortran style. I remember that a Fortran complex
literal is like an ada aggragate, but I do not remember how it is
output; using 2 Gw.d formats?).
9. Same thing for "Get from a string"
Those 2 procedures is what I found was the minimal support needed to
do some form of I/O with complex numbers.
10. Personally I would be in favor of initializing the Complex (with outrageous
values like Real'Large) but this would certainly start a religious war,
so I left this one out.
Applying all this gives the following alternate proposal (Ada 83 version):
with Elementary_Functions_Exceptions;
generic
type Real is digits <>;
package Generic_Complex_Support is
type Complex is
record
Re, Im : Real;
end record;
I : constant Complex := (0.0, 1.0);
Zero : constant Complex := (0.0, 0.0);
function Make (Real_Part : Real) return Complex;
function Make (Real_Part, Imag_Part : Real) return Complex;
function Modulus (C : Complex) return Real;
function "abs" (C : Complex) return Real renames Modulus;
function Argument (C : Complex) return Real;
function Argument (C : Complex;
Cycle : Real) return Real;
function Make_Polar (Modulus : Real;
Argument : Real)
return Complex;
function Make_Polar (Modulus : Real;
Argument : Real;
Cycle : Real)
return Complex;
function "+" (Right : Complex) return Complex;
function "-" (Right : Complex) return Complex;
function Conjugate (C : Complex) return Complex;
function "+" (Left,
Right : Complex) return Complex;
function "+" (Left : Real;
Right : Complex) return Complex;
function "+" (Left : Complex;
Right : Real) return Complex;
function "-" (Left,
Right : Complex) return Complex;
function "-" (Left : Real;
Right : Complex) return Complex;
function "-" (Left : Complex;
Right : Real) return Complex;
function "*" (Left,
Right : Complex) return Complex;
function "*" (Left : Real;
Right : Complex) return Complex;
function "*" (Left : Complex;
Right : Real) return Complex;
function "/" (Left,
Right : Complex) return Complex;
function "/" (Left : Real;
Right : Complex) return Complex;
function "/" (Left : Complex;
Right : Real) return Complex;
function "**" (Left : Complex;
Right : Integer) return Complex;
type Display_Style is (Ada_Aggregate, -- (-34.2,23.6)
Re_Plus_Im_I, -- -34.2+23.6i
Bracketed -- [-34.2;23.6]
);
generic
with procedure Put (To : out String; Item : Real);
Style : Display_Style := Re_Plus_Im_I;
procedure Put (To : out String; Item : in Complex);
generic
with procedure Get (From : in String;
Item : out Real;
Last : out Positive);
procedure Get (From : in String; Item : out Complex; Last : out Positive);
Argument_Error : exception
renames Elementary_Functions_Exceptions.Argument_Error;
end Generic_Complex_Support;
APPENDIX B (Alternative proposal of Bardin)
Ken,
The following is not finished. It is my intent to incorporate the complex
elementary functions package into this package, because it then forms a
"complete" abstraction which supports expressions like: Rho * Exp (i * Theta),
for instance. It is not necessary to have a separate package for elementary
functions as is the case for type Real. The package does not yet address the
Text_Io question.
A bit about philosophy. My approach is non-traditional in Ada terms. My
overriding concern is with the end user's view, particularly that of a
mathematician, scientist, or engineer coming to Ada from FORTRAN. He is used
to relatively terse identifiers, and those are accommodated. For the
long-identifier Ada purist, those are also retained. (The cost for this
flexibility is trivial: a few renaming declarations.) I provide a sample of
the use of each construct at the point of definition for clarity; this may be
stripped out later if desired. Further examples are given in the sample
application program. It is my goal to make the use of Ada for this class of
application seem to be more natural and more expressive in Ada than in Fortran.
I agree with Philippe that Real is the natural type name for the components of
type Complex.
This package compiles. I have also looked at some sample bodies. The use of
a distinct type for imaginary values allows some potential optimazations for
some classes of operations, such as Imaginary x Imaginary -> Real and for
functions like:
function "**" (Left : Imaginary; Right : Integer) return Complex;
The main remaining problem with this approach is that:
- +
must be explicitly converted (or qualified) as:
Real(- ) +
I have seen this phenomenon before, in defining a rational arithmetic package
for my design course, and I believe it is a(n unecessary) defect in the
Ada overload resolution rules. It should be possible to resolve using
multiple implicit conversions. Perhaps this is something that can be fixed in
Ada 9X.
This needs more work, but I think it is close to meeting my personal goals.
Feel free to show it around and accumulate suggestions for rework.
Regards,
Bryce
package Elementary_Functions_Exceptions is
Argument_Error : exception;
-- ...
end;
with Elementary_Functions_Exceptions;
generic
type Real is digits <>;
package Generic_Complex_Types is
-- *** This is not done yet. ***
type Complex is
record
X, Y : Real;
-- These component names allow direct references to an object Obj of
-- type Complex of the form: Obj.X and Obj.Y
end record;
type Imaginary is private;
-- The use of lower case is intended.
i : constant Imaginary; -- for mathematicians and scientists
j : constant Imaginary; -- for engineers
One : constant Complex := (1.0, 0.0);
Zero : constant Complex := (0.0, 0.0);
-- Pure Complex operations:
function "+" (Right : Complex) return Complex;
-- Cplx := + Cplx;
function "-" (Right : Complex) return Complex;
-- Cplx := - Cplx;
function "+" (Left,
Right : Complex) return Complex;
-- Cplx := Cplx + Cplx;
function "-" (Left,
Right : Complex) return Complex;
-- Cplx := Cplx - Cplx;
function Conjugate (Z : Complex) return Complex;
-- Cplx := Conjugate (Cplx);
function "not" (Right : Complex) return Complex;
-- Cplx := not Cplx;
-- If "not" is allowed, then this is unecessary.
function Conj (Z : Complex) return Complex renames Conjugate;
-- Cplx := Conj (Cplx);
function "*" (Left,
Right : Complex) return Complex;
-- Cplx := Cplx * Cplx;
function "/" (Left,
Right : Complex) return Complex;
-- Cplx := Cplx / Cplx;
function "**" (Left : Complex;
Right : Integer) return Complex;
-- Cplx := Cplx ** Int;
-- Mixed Real and Complex operations:
function Complex_From_Cartesian (X : Real) return Complex;
-- Cplx := Complex_From_Cartesian (Re);
function Cart (X : Real) return Complex renames Complex_From_Cartesian;
-- Cplx := Cart (Re);
function Complex_From_Cartesian (X,
Y : Real) return Complex;
-- Cplx := Complex_From_Cartesian (Re, Re);
function Cart (X,
Y : Real) return Complex renames Complex_From_Cartesian;
-- Cplx := Cart (Re, Re);
function Complex_From_Polar (Modulus : Real;
Argument : Real) return Complex;
-- Cplx := Complex_From_Polar (Re, Re);
function Polar (Magnitude : Real;
Angle : Real) return Complex renames Complex_From_Polar;
-- Cplx := Polar (Re, Re);
function Complex_From_Polar (Modulus : Real;
Angle : Real;
Cycle : Real) return Complex;
-- Cplx := Complex_From_Polar (Re, Re, Re);
function Polar (Magnitude : Real;
Angle : Real;
Cycle : Real) return Complex renames Complex_From_Polar;
-- Cplx := Polar (Re, Re, Re);
function Real_Part (Z : Complex) return Real;
-- Re := Real_Part(Cplx);
function X (Z : Complex) return Real renames Real_Part;
-- Re := X(Cplx);
function Imaginary_Part (Z : Complex) return Real;
-- Re := Imaginary_Part(Cplx);
function Y (Z : Complex) return Real renames Real_Part;
-- Re := Y(Cplx);
procedure Set_Real_Part (Z : in out Complex;
X : Real);
-- Set_Real_Part (Cplx, Re);
procedure Set_X (Z : in out Complex;
X : Real) renames Set_Real_Part;
-- Set_X (Cplx, Re);
procedure Set_Imaginary_Part (Z : in out Complex;
X : Real);
-- Set_Imaginary_Part (Cplx, Re);
procedure Set_Y (Z : in out Complex;
Y : Real) renames Set_Imaginary_Part;
-- Set_Y (Cplx, Re);
function Modulus (Z : Complex) return Real;
-- Re := Modulus (Cplx);
function "abs" (Right : Complex) return Real renames Modulus;
-- Re := abs Cplx;
function Argument (Z : Complex) return Real;
-- Re := Argument (Cplx);
function Angle (Z : Complex) return Real renames Argument;
-- Re := Angle (Cplx);
function Argument (Z : Complex;
Cycle : Real) return Real;
-- Re := Argument (Cplx, Re);
function Angle (Z : Complex;
Cycle : Real) return Real renames Argument;
-- Re := Angle (Cplx, Re);
function "+" (Right : Real) return Complex;
-- Cplx := + Re;
function "-" (Right : Real) return Complex;
-- Cplx := - Re;
function "+" (Left : Real;
Right : Complex) return Complex;
-- Cplx := Re + Cplx;
function "+" (Left : Complex;
Right : Real) return Complex;
-- Cplx := Cplx + Re;
function "-" (Left : Real;
Right : Complex) return Complex;
-- Cplx := Re - Cplx;
function "-" (Left : Complex;
Right : Real) return Complex;
-- Cplx := Cplx - Re;
function "*" (Left : Real;
Right : Complex) return Complex;
-- Cplx := Re * Cplx;
function "*" (Left : Complex;
Right : Real) return Complex;
-- Cplx := Cplx * Re;
function "/" (Left : Real;
Right : Complex) return Complex;
-- Cplx := Re / Cplx;
function "/" (Left : Complex;
Right : Real) return Complex;
-- Cplx := Cplx / Re;
-- Mixed Real and Imaginary operations:
function "*" (Left : Real;
Right : Imaginary) return Imaginary;
-- Im := Re * Im;
function "*" (Left : Imaginary;
Right : Real) return Imaginary;
-- Im := Im * Re;
function "*" (Left,
Right : Imaginary) return Real;
-- Re := Im * Im;
function "/" (Left,
Right : Imaginary) return Real;
-- Re := Im / Im;
function "/" (Left : Real;
Right : Imaginary) return Imaginary;
-- Im := Re / Im;
function "/" (Left : Imaginary;
Right : Real) return Imaginary;
-- Im := Im / Re;
-- Pure Imaginary operations:
function "+" (Right : Imaginary) return Imaginary;
-- Im := + Im;
function "-" (Right : Imaginary) return Imaginary;
-- Im := - Im;
function "+" (Left,
Right : Imaginary) return Imaginary;
-- Im := Im + Im;
function "-" (Left,
Right : Imaginary) return Imaginary;
-- Im := Im - Im;
-- Mixed Imaginary and Complex operations:
function "+" (Right : Imaginary) return Complex;
-- Cplx := + Im;
function "-" (Right : Imaginary) return Complex;
-- Cplx := - Im;
function "+" (Left : Complex;
Right : Imaginary) return Complex;
-- Cplx := Cplx + Im;
function "+" (Left : Imaginary;
Right : Complex) return Complex;
-- Cplx := Im + Cplx;
function "-" (Left : Complex;
Right : Imaginary) return Complex;
-- Cplx := Cplx - Im;
function "-" (Left : Imaginary;
Right : Complex) return Complex;
-- Cplx := Im - Cplx;
function "*" (Left : Imaginary;
Right : Complex) return Complex;
-- Cplx := Im * Cplx;
function "*" (Left : Complex;
Right : Imaginary) return Complex;
-- Cplx := Cplx * Im;
function "/" (Left : Imaginary;
Right : Complex) return Complex;
-- Cplx := Im / Cplx;
function "/" (Left : Complex;
Right : Imaginary) return Complex;
-- Cplx := Cplx / Im;
function "**" (Left : Imaginary;
Right : Integer) return Complex;
-- Cplx := Im ** Int;
-- Mixed Real, Imaginary, and Complex operations:
function "+" (Left : Real;
Right : Imaginary) return Complex;
-- Cplx := Re + Im;
function "+" (Left : Imaginary;
Right : Real) return Complex;
-- Cplx := Im + Re;
function "-" (Left : Real;
Right : Imaginary) return Complex;
-- Cplx := Re - Im;
function "-" (Left : Imaginary;
Right : Real) return Complex;
-- Cplx := Im - Re;
-- function "*" (Left : Real;
-- Right : Imaginary) return Complex;
-- Cplx := Re * Im;
-- function "*" (Left : Imaginary;
-- Right : Real) return Complex;
-- Cplx := Im * Re;
Argument_Error : exception
renames Elementary_Functions_Exceptions.Argument_Error;
private
type Imaginary is
record
Value : Real;
end record;
i : constant Imaginary := (Value => 1.0);
j : constant Imaginary := (Value => 1.0);
end Generic_Complex_Types;
with Generic_Complex_Types;
procedure Use_Complex_Types is
subtype Real is Float;
package C is new Generic_Complex_Types (Float);
use C;
Re, Rho, Theta, X, Y : Real;
Im, Im1, Im2 : Imaginary;
Cplx, Cplx1, Cplx2 : Complex := Zero;
begin
-- *** This is not done yet. ***
-- Pure Complex operations:
Cplx := Zero;
Cplx := One;
Cplx := + One;
Cplx := - One;
Cplx := One + One;
Cplx := Cplx1 - Cplx2;
Cplx := Cplx1 * Cplx2;
Cplx := Cplx1 / Cplx2;
Cplx := Cplx**3;
Cplx := (-1.0, 2.0);
Cplx := Complex'(X => 2.0, Y => -2.0);
Cplx := not Cplx;
-- Pure Real operations:
Re := 1.0;
Re := + Re;
Re := + 3.0;
Re := - Re;
Re := - 1.0;
Re := Re + Re;
Re := 1.0 + 1.0;
Re := - 1.0 + 1.0;
Re := Re - Re;
Re := 1.0 - 2.0;
Re := - 3.0 - 5.0;
Re := Re * Re;
Re := - 2.0 * 3.0;
X := 1.0;
Y := 2.0;
Re := X / Y;
Re := 2.0 / 3.0;
Re := Re ** 2;
Re := 2.0 ** 3;
-- Pure Imaginary operations:
Im := + i;
Im := - i;
Im := Im + Im;
Im := Im - j;
-- Mixed Imaginary and Complex operations:
Cplx := j ** 2;
Cplx := + j;
Cplx := - i;
-- Mixed Real and Imaginary operations:
Im := 1.0 * i;
Im := 2.0 * j;
Re := j * j;
Re := Im1 / Im2;
Im := Re / Im;
Im := 0.1 / j;
Im := Im / Re;
Im := j / 2.0;
-- Mixed Real and Complex operations:
X := Cplx.X;
X := C.X(Cplx);
Y := Cplx.Y;
Y := C.Y(Cplx);
Set_X(Cplx, 1.0);
Set_Y(Cplx, -1.0);
Cplx.X := 3.0;
Cplx.Y := -3.0;
Re := abs Cplx;
Cplx := + Re;
Cplx := + 3.0;
Cplx := - Re;
Cplx := - 4.0;
Cplx := Re + Cplx;
Cplx := 1.0 + One;
Cplx := Real(- 2.0) + (2.0, 3.0);
Cplx := Cplx + Re;
Cplx := Zero + 1.0;
-- Mixed Real, Imaginary, and Complex operations:
Cplx := Re + Im;
Cplx := X + j * Y;
Cplx := Real(- 1.0) + 2.0 * i;
Cplx := j + Re;
Cplx := 2.0 * i + 1.0;
Cplx := X - i * Y;
Cplx := Re - Im;
Cplx := 3.0 - i;
Cplx := + (Re * Im);
Cplx := + (2.0 * i);
Cplx := + (j * Y);
Cplx := + (i * 0.5);
end Use_Complex_Types;
APPENDIX C (Sketch of alternative proposal by Mathis)
Here is the file I showed you at the DR meeting in Lexington.
Note: based on discussion with Ken Dritz and Bryce Barden at
dinner October 1, 1992, while at Ada9X DR meeting in Lexington.
type imaginary_indicator is (i, j);
type imaginary is new real;
type complex is
record
real_part, imaginary_part: real;
end record;
function "*" (mark: imaginary_indicator; y: real)
return imaginary is
begin
return imaginary(y);
end "*";
function "*" (y: real; mark: imaginary_indicator)
return imaginary is
begin
return imaginary(y);
end "*";
function "+" (x: real; y: imaginary)
return complex is
begin
return (x, real(y));
end "+";
function "+" (y: imaginary; x: real)
return complex is
begin
return (x, real(y));
end "+";
function "-" (x: real; y: imaginary)
return complex is
begin
return (x, -real(y));
end "-";
function "-" (y: imaginary; x: real)
return complex is
begin
return (-x, real(y));
end "+";
function "-" (y: imaginary)
return imaginary is
begin
return imaginary( -real(y));
end "-";
function e_i (theta: real) return complex is
begin
return complex(cos(theta), sin(theta));
end e_i;
function "*" (x: real; z: complex) return complex is
begin
return complex(x*z.real_part, x*z.imaginary_part);
end "*";
-- above gives us possibility of saying
-- z:= modulus * e_i (angle);
-- we were also trying to get a way to invoke
-- the polar form constructor
type complex_angle_type is new real;
function e_2_i (theta: real)
return complex_angle_type is
begin
return complex_angle_type(theta);
end e_2_i;
function "*" (r: real;
strange_thing: complex_angle_type)
return polar_complex is
begin
return (r, real(strange_thing));
end "*";
-- then we could also say
-- z_polar := modulus * e_2_i (angle)
-- this does it in polar form
-- trig functions would only be invoked in
-- z := convert_polar_to_cartesian (z);
We should work out most of the constructors and functions and
then do a couple of example computations to illustrate use.
I am also surprised that this hasn't been done before, but I
certainly don't know of any.
-- Bob