John Barnes

# 6.3 Numerics

Although Ada 95 introduced unsigned integer types in the form of modular types, nevertheless, the strong typing rules of Ada have not made it easy to get unsigned and signed integers to work together. The following discussion using Ada 95 is based on that in AI-340.
Suppose we wish to implement a simulation of a typical computer which has addresses and offsets. We make it a generic
generic
type Offset_Type is range <>;
...
package Simulator is
...
end Simulator;
Addresses are represented as unsigned integers (a modular type), whereas offsets are signed integers. The function Calc_Address aims to add an offset to a base address and return an address. The offset could be negative.
Naïvely we might hope to write
begin
return Base_Add + Offset;    -- illegal
but this is plainly illegal because Base_Add and Offset are of different types.
We can try a type conversion thus
or perhaps, since Address_Type might have a constraint,
but in any case the conversion is doomed to raise Constraint_Error if Offset is negative.
We then try to be clever and write
but this raises Constraint_Error if Address_Type'Modulus > Offset_Type'Base'Last which it often will be. To see this consider for example a 32-bit machine with
type Offset_Type is range –(2**31) .. 2**31–1;
in which case Address_Type'Modulus is 2**32 which is greater than Offset_Type'Base'Last which is 2**31–1.
So we try an explicit test for a negative offset
if Offset >= 0 then
else
end if;
But if Address_Type'Base'Last < Offset_Type'Last then this will raise Constraint_Error for some values of Offset. Unlikely perhaps but this is a generic and so ought to work for all possible pairs of types.
If we attempt to overcome this then we run into problems in trying to compare these two values since they are of different types and converting one to the other can raise the Constraint_Error problem once more. One solution is to use a bigger type to do the test but this may not exist in some implementations. We could of course handle the Constraint_Error and then patch up the answer. The ruthless programmer might even think of Unchecked_Conversion but this has its own problems. And so on – 'tis a wearisome tale.
The problem is neatly overcome in Ada 2005 by the introduction of a new functional attribute
function S'Mod(Arg: universal_integerreturn S'Base;
S'Mod applies to any modular subtype S and returns
Arg mod S'Modulus
In other words it converts a universal_integer value to the modular type using the corresponding mathematical mod operation. We can then happily write
begin
and this always works.
The next topic in the numerics area concerns rounding. One of the problems in the design of any programming language is getting the correct balance between performance and portability. This is particularly evident with numeric types where the computer has to implement only a crude approximation to the mathematician's integers and reals. The best performance is achieved by using types and operations that correspond exactly to the hardware. On the other hand, perfect portability requires using types with precisely identical characteristics on all implementations.
An interesting example of this problem arises with conversions from a floating point type to an integer type when the floating type value is midway between two integer values.
In Ada 83 the rounding in the midway case was not specified. This upset some people and so Ada 95 went the other way and decreed that such rounding was always away from zero. As well as this rule for conversion to integer types, Ada 95 also introduced a functional attribute to round a floating value. Thus for a subtype S of a floating point type T we have
function S'Rounding(X: T) return T;
This returns the nearest integral value and for midway values rounds away from zero.
Ada 95 also gives a bit more control for the benefit of the statistically minded by introducing
function S'Unbiased_Rounding(X: T) return T;
This returns the nearest integral value and for midway values rounds to the even value.
However, there are many applications where we don't care which value we get but would prefer the code to be fast. Implementers have reported problems with the elementary functions where table look-up is used to select a particular polynomial expansion. Either polynomial will do just as well when at the midpoint of some range. However on some popular hardware such as the Pentium, doing the exact rounding required by Ada 95 just wastes time and the resulting function is perhaps 20% slower. This is serious in any comparison with C.
This problem is overcome in Ada 2005 by the introduction of a further attribute
function S'Machine_Rounding(X: T) return T;
This does not specify which of the adjacent integral values is returned if X lies midway. Note that it is not implementation defined but deliberately unspecified. This should discourage users from depending upon the behaviour on a particular implementation and thus writing non-portable code.
Zerophiles will be pleased to note that if S'Signed_Zeros is true and the answer is zero then it has the same sign as X.
It should be noted that Machine_Rounding, like the other rounding functions, returns a value of the floating point type and not perhaps universal_integer as might be expected. So it will typically be used in a context such as
X: Some_Float;
Index: Integer;
...
Index := Integer(Some_Float'Machine_Rounding(X));
...    -- now use Index for table look-up
Implementations are urged to detect this case in order to generate fast code.
The third improvement to the core language in the numerics area concerns fixed point arithmetic. This is a topic that concerns few people but those who do use it probably feel passionately about it.
The trouble with floating point is that it is rather machine dependent and of course integers are just integers. Many application areas have used some form of scaled integers for many decades and the Ada fixed point facility is important in certain applications where rigorous error analysis is desirable.
The model of fixed point was changed somewhat from Ada 83 to Ada 95. One change was that the concepts of model and safe numbers were replaced by a much simpler model just based on the multiples of the number small. Thus consider the type
Del: constant := 2.0**(–15);
type Frac is delta Del range –1.0 .. 1.0;
In Ada 83 small was defined to be the largest power of 2 not greater than Del, and in this case is indeed 2.0**(–15). But in Ada 95, small can be chosen by the implementation to be any power of 2 not greater than Del provided of course that the full range of values is covered. In both languages an aspect clause can be used to specify small and it need not be a power of 2. (Remember that representation clauses are now known as aspect clauses.)
A more far reaching change introduced in Ada 95 concerns the introduction of operations on the type universal_fixed and type conversion.
A minor problem in Ada 83 was that explicit type conversion was required in places where it might have been considered quite unnecessary. Thus supposing we have variables F, G, H of the above type Frac, then in Ada 83 we could not write
H := F * G;    -- illegal in Ada 83
but had to use an explicit conversion
H := Frac(F * G);    -- legal in Ada 83
In Ada 83, multiplication was defined between any two fixed point types and produced a result of the type universal_fixed and an explicit conversion was then required to convert this to the type Frac.
This explicit conversion was considered to be a nuisance so the rule was changed in Ada 95 to say that multiplication was only defined between universal_fixed operands and delivered a universal_fixed result. Implicit conversions were then allowed for both operands and result provided the type resolution rules identified no ambiguity. So since the expected type was Frac and no other interpretation was possible, the implicit conversion was allowed and so in Ada 95 we can simply write
H := F * G;    -- legal in Ada 95
Similar rules apply to division in both Ada 83 and Ada 95.
Note however that
F := F * G * H;    -- illegal
is illegal in Ada 95 because of the existence of the pervasive type Duration defined in Standard. The intermediate result could be either Frac or Duration. So we have to add an explicit conversion somewhere.
One of the great things about Ada is the ability to define your own operations. And in Ada 83 many programmers wrote their own arithmetic operations for fixed point. These might be saturation operations in which the result is not allowed to overflow but just takes the extreme implemented value. Such operations often match the behaviour of some external device. So we might declare
function "*"(Left, Right: Frac) return Frac is
begin
return Standard."*"(Left, Right);
exception
when Constraint_Error =>
if (Left>0.0 and Right>0.0) or (Left<0.0 and Right<0.0) then
return Frac'Last;
else
return Frac'First;
end if;
end "*";
and similar functions for addition, subtraction, and division (taking due care over division by zero and so on). This works fine in Ada 83 and all calculations can now use the new operations rather than the predefined ones in a natural manner.
Note however that
H := Frac(F * G);
is now ambiguous in Ada 83 since both our own new "*" and the predefined "*" are possible interpretations. However, if we simply write the more natural
H := F * G;
then there is no ambiguity. So we can program in Ada 83 without the explicit conversion.
However, in Ada 95 we run into a problem when we introduce our own operations since
H := F * G;
is ambiguous because both the predefined operation and our own operation are possible interpretations of "*" in this context. There is no cure for this in Ada 95 except for changing our own multiplying operations to be functions with identifiers such as mul and div. This is a very tedious chore and prone to errors.
It has been reported that because of this difficulty many projects using fixed point have not moved from Ada 83 to Ada 95.
This problem is solved in Ada 2005 by changing the name resolution rules to forbid the use of the predefined multiplication (division) operation if there is a user-defined primitive multiplication (division) operation for either operand type unless there is an explicit conversion on the result or we write Standard."*" (or Standard."/").
This means that when there is no conversion as in
H := F * G;
then the predefined operation cannot apply if there is a primitive user-defined "*" for one of the operand types. So the ambiguity is resolved. Note that if there is a conversion then it is still ambiguous as in Ada 83.
If we absolutely need to have a conversion then we can always use a qualification as well or just instead. Thus we can write
F := Frac'(F * G) * H;
and this will unambiguously use our own operation.
On the other hand if we truly want to use the predefined operation then we can always write
H := Standard."*"(F, G);
Another example might be instructive. Suppose we declare three types TL, TA, TV representing lengths, areas, and volumes. We use centimetres as the basic unit with an accuracy of 0.1 cm together with corresponding consistent units and accuracies for areas and volumes. We might declare
type TL is delta 0.1 range –100.0 .. 100.0;
type TA is delta 0.01 range –10_000.0 .. 10_000.0;
type TV is delta 0.001 range –1000_000.0 .. 1000_000.0;
for TL'Small use TL'Delta;
for TA'Small use TA'Delta;
for TV'Small use TV'Delta;
function "*"(Left: TL; Right: TL) return TA;
function "*"(Left: TL; Right: TA) return TV;
function "*"(Left: TA Right: TL) return TV;
function "/"(Left: TV; Right: TL) return TA;
function "/"(Left: TV; Right: TA) return TL;
function "/"(Left: TA; Right: TL) return TL;
XL, YL: TL;
XA, YA: TA;
XV, YV: TV;
These types have an explicit small equal to their delta and are such that no scaling is required to implement the appropriate multiplication and division operations. This absence of scaling is not really relevant to the discussion below but simply illustrates why we might have several fixed point types and operations between them.
Note that all three types have primitive user-defined multiplication and division operations even though in the case of multiplication, TV only appears as a result type. Thus the predefined multiplication or division with any of these types as operands can only be considered if the result has a type conversion.
As a consequence the following are legal
XV := XL * XA;    -- OK, volume = length × area
XL := XV / XA;    -- OK, length = volume ÷ area
but the following are not because they do not match the user-defined operations
XV := XL * XL;    -- no, volume ≠ length × length
XV := XL / XA;    -- no, volume ≠ length ÷ area
XL := XL * XL;    -- no, length ≠ length × length
But if we insist on multiplying two lengths together then we can use an explicit conversion thus
XL := TL(XL * XL);    -- legal, predefined operation
and this uses the predefined operation.
If we need to multiply three lengths to get a volume without storing an intermediate area then we can write
XV := XL * XL * XL;
and this is unambiguous since there are no explicit conversions and so the only relevant operations are those we have declared.
It is interesting to compare this with the corresponding solution using floating point where we would need to make the unwanted predefined operations abstract as discussed in an earlier chapter (see 2.7).
It is hoped that the reader has not found this discussion to be too protracted. Although fixed point is a somewhat specialized area, it is important to those who find it useful and it is good to know that the problems with Ada 95 have been resolved.
There are a number of other improvements in the numerics area but these concern the Numerics annex and are discussed in Section 7.6.

© 2005, 2006, 2007 John Barnes Informatics.