Multics Technical Bulletin MTB-687
Fortran HFP mode.
To: Distribution
From: Hal Hoover
Date: 19 October 1984
Subject: Fortran Hexadecimal Floating Point Support.
1. Abstract
This MTB discusses the implementation of Fortran's hexadecimal
floating point mode. The changes that have been made to the
Fortran compiler, the Fortran runtime library, and the Operating
system in general, are discussed.
Comments on this MTB should be sent to the author -
via Multics mail to:
Hoover.Multics on System M
via posted mail to:
Hal Hoover
Advanced Computing Technology Centre
Foothills Professional Building
Room #301, 1620 - 29th Street N.W.
Calgary, Alberta T2N-4L7
CANADA
via telephone to:
(403)-270-5400
(403)-270-5422
________________________________________
Multics project internal documentation; not to be reproduced or
distributed outside the Multics project.
MTB-687 Multics Technical Bulletin
Fortran HFP mode.
TABLE OF CONTENTS
Section Page Subject
======= ==== =======
1 i Abstract
2 1 Introduction
3 3 Compiler Changes
3.1 4 . . The Hexadecimal Floating Point option
3.2 5 . . Parsing changes
3.3 6 . . Compiler Output Changes
3.3.1 6 . . . . Symbol Table Changes
3.3.2 6 . . . . Object Code Changes
3.3.2.1 6 . . . . . . Managing HFP Mode
3.3.2.2 8 . . . . . . Integer to Floating Point
3.3.2.3 9 . . . . . . Floating Point to Integer
3.3.2.4 11 . . . . . . Calls to 'pl1_operators_'
3.3.2.5 12 . . . . . . DU address modification
3.4 13 . . Other Compiler Changes
4 14 Runtime Library Changes
4.1 15 . . I/O Routine Changes
4.2 16 . . External Builtin Routine Changes
5 19 System Support Changes
5.1 20 . . 'any_to_any_' Changes
5.2 21 . . 'pl1_operators_' Changes
5.3 26 . . Math Routines
5.3.1 27 . . . . ALM-Coded Math Routines
5.3.1.1 27 . . . . . . 'arc_sine_'
5.3.1.2 28 . . . . . . '(double_)arc_tangent_'
5.3.1.3 28 . . . . . . 'call_math_error_'
5.3.1.4 29 . . . . . . '(double_)exponential_'
5.3.1.5 30 . . . . . . 'fort_math_ops_'
5.3.1.6 30 . . . . . . 'integer_power_integer_'
5.3.1.7 30 . . . . . . '(double_)logarithm_'
5.3.1.8 31 . . . . . . 'math_constants_'
5.3.1.9 31 . . . . . . 'power_'
5.3.1.10 31 . . . . . . 'power_integer_'
5.3.1.11 32 . . . . . . '(double_)sine_'
5.3.1.12 32 . . . . . . '(double_)square_root_'
5.3.1.13 33 . . . . . . '(double_)tangent_'
5.3.2 34 . . . . Non-ALM Coded Math Routines
5.3.2.1 34 . . . . . . 'cabs_'
5.3.2.2 35 . . . . . . 'ccos_'
5.3.2.3 35 . . . . . . 'cexp_'
5.3.2.4 35 . . . . . . 'clog_'
5.3.2.5 35 . . . . . . '(d)cosh_'
5.3.2.6 36 . . . . . . 'csin_'
5.3.2.7 36 . . . . . . 'csqrt_'
5.3.2.8 36 . . . . . . 'cxp2_'
5.3.2.9 36 . . . . . . '(d)sinh_'
Multics Technical Bulletin MTB-687
Fortran HFP mode.
5.3.2.10 37 . . . . . . '(double_)tanh_'
Multics Technical Bulletin MTB-687
Fortran HFP mode.
2. Introduction
Prior to the introduction of the 870M series of processors,
Multics hardware supported two types of floating point number:
Binary Floating Point (BFP) and Decimal Floating Point (DFP). A
BFP number consists of a binary fraction (i.e. a binary number
between -1 and +1) times a power of two. A DFP number consists
of a decimal integer times a power of ten. BFP numbers can be
stored in two formats: single precision (which has 27 bits or
about 8.1 decimal digits of precision) and double precision
(which has 63 bits or about 18.9 decimal digits of precision).
DFP numbers can be stored with a precision from 1 through 59
decimal digits. The largest BFP number is about 1.7e38; while
the largest DFP number is about 1e187. Although BFP numbers have
both less range and less precision than DFP numbers, they are
used much more because they occupy less space for the same
precision and because binary arithmetic operations are very much
faster.
The 870M series of processors added hardware support for a third
type of floating point number: Hexadecimal Floating Point (HFP).
An HFP number consists of a hexadecimal fraction times a power of
sixteen. The precision of an HFP number depends on the value of
the number. In single precision, an HFP number has from 7.2 to
8.1 decimal digits of precision; while in double precision, it
has from 18 to 18.9 decimal digits of precision. The largest HFP
number is about 8.4e152. HFP numbers occupy the same amount of
storage (for equivalent precision) as BFP numbers and arithmetic
operations on HFP numbers are at the same speed as BFP numbers.
Thus at the cost of a small loss in precision HFP numbers give
the same performance as BFP numbers, but with a substantially
greater range.
There has long been grumbling among scientists and engineers that
the range of numbers available on Multics FORTRAN (which does not
have the ability to use DFP numbers) is too small. Adding
support for DFP numbers would not help much because the types of
programs that need the larger range are usually very
computation-intensive and so could not afford the performance
penalties associated with the much slower decimal arithmetic.
Adding support to FORTRAN for HFP numbers would be very helpful
for these types of programs because the range would be
substantially increased at only the cost of a small loss in
precision.
MTB-687 Multics Technical Bulletin
Fortran HFP mode.
At the request of Marketing, a project was struck to add HFP
support to Multics FORTRAN. The aim was to provide HFP as a
natural extension of the existing compiler and runtime, with the
restriction that the changes should not adversely affect
performance of nonHFP programs. This MTB describes the changes
made to accomplish that goal.
Multics Technical Bulletin MTB-687
Fortran HFP mode.
3. Compiler Changes
There are obviously changes required in the user interface areas
of the compiler, since there must be a way to specify HFP and the
object generated for HFP will be different than for BFP. The bit
pattern of the HFP representation of a number is usually
different than that for the BFP representation of the same
number. Therefore, changes are also needed wherever the compiler
must deal with floating point constants -- such as constant
folding in the optimizer.
MTB-687 Multics Technical Bulletin
Fortran HFP mode.
3.1. The Hexadecimal Floating Point option
The most obvious change required of the compiler is the ability
to specify when to use BFP and when to use HFP. This choice
could be made on any of three levels: by variable, by routine,
or by compilation unit. Selection of BFP or HFP on a by-variable
basis is by far the most powerful approach. It is also the
hardest to implement as it requires support for mixed base
arithmetic. Selection by routine is not much more powerful than
is selection by compilation unit since any Fortran subroutine or
function can be compiled as a separate unit.
After weighing the relative difficulties of implementing the
three selection methods against their utility -- and the fact
that Marketing only required that HFP be supported for whole
programs, it was decided to allow BFP/HFP selection by
compilation unit only.
Restricting BFP/HFP selection to whole compilation units allows a
very easy way for the user to specify his selection: control
arguments on the 'fortran' command. Thus we have introduced two
new 'fortran' control arguments, '-binary_floating_point' ('-bfp'
for short) and '-hexadecimal_floating_point' ('-hfp' for short).
As most users will want to always compile a given program in the
same mode, 'binary_floating_point' ('bfp' for short) and
'hexadecimal_floating_point' ('hfp' for short) options are
permitted on the nonstandard '%global' statement. If a mode is
specified on the command line and the conflicting mode is
specified in a '%global' statement in the source, a warning
diagnostic is issued and the mode specified on the command line
is honoured. If the mode is not specified on the command line
nor in the source in a '%global' statement, BFP mode is assumed
in order to maintain compatibility with previous versions.
Multics Technical Bulletin MTB-687
Fortran HFP mode.
3.2. Parsing changes
There is a compiler change required in the lexical analyzer. The
lexical analyzer splits the source program into tokens for the
benefit of the parser. In particular, the lexical analyzer is
responsible for converting numeric constants into their internal
representation. Thus it must be changed to allow the larger
range of HFP constants and to generate the correct bit pattern
for these constants (since that will usually be different than
the bit pattern for BFP constants).
MTB-687 Multics Technical Bulletin
Fortran HFP mode.
3.3. Compiler Output Changes
There are basically three kinds of output from the compiler: the
listing, the object code and the symbol table. No changes are
needed to the listing. The only change needed to the symbol
table is where data types are specified: new codes for HFP
variables must be defined and used. Only minor changes are
needed in the object code.
3.3.1. Symbol Table Changes
The symbol table is basically a description of the layout of the
object (i.e. where the code for each source statement begins)
and the type and location of the variables. It supplies the
information needed to 'probe' a program. The only effect on the
symbol table of adding HFP support is the addition of new data
type codes to describe HFP variables. The codes for the data
types used in the various Multics languages are defined in the
system include file 'std_descriptor_types.incl.pl1'. The codes
for BFP data are: 3 (real single precision), 4 (real double
precision), 7 (complex single precision) and 8 (complex double
precision). Through arrangement with the custodian of this
include file, we have added the following codes for HFP data: 47
(real single precision), 48 (real double precision), 49 (complex
single precision) and 50 (complex double precision). Code 50 is
not used by FORTRAN, since we do not support complex double
precision, but is provided for future expansion.
'Probe' has been updated by its maintainer to handle the new
codes used to describe HFP data.
3.3.2. Object Code Changes
3.3.2.1. Managing HFP Mode
When a floating point instruction is executed, a processor mode
determines whether it is executed according to the rules of BFP
arithmetic or those of HFP arithmetic. There are two bits, one
in the Indicator Register (IR) and one in the Mode Register (MR)
that determine whether the processor is in BFP or HFP mode. The
bit in the MR says whether the user is allowed to enter HFP mode
and the bit in the IR says whether the user wants BFP or HFP
mode. The processor is in HFP mode if and only if both bits are
on.
Multics Technical Bulletin MTB-687
Fortran HFP mode.
The HFP bit in the IR can be set or cleared at the user's
discretion by executing an 'ldi' instruction. It takes a
privileged instruction to alter the MR. Fortunately, the
hardcore maintainers have provided a system routine,
'hcs_$set_hexfp_control', to manage the HFP bit in the MR.
However, this routine may not necessarily grant a user's request
to set the HFP bit in the MR, since the hardware may not support
HFP arithmetic (only 870M processors do) or the user may not be
allowed (by the System Administrator) to use HFP mode. Thus we
must be prepared to signal an error indicating the inability to
enter HFP mode.
This scheme implies that we must be very careful about how we
handle calls to external procedures. We must guard that a call
from an HFP routine to a BFP routine does not leave the BFP
routine in HFP mode. We must also guard that a return from a BFP
routine called from an HFP routine restores HFP mode. These are
not new problems; another bit in the IR which controls whether
arithmetic overflows and underflows signal a hardware fault
suffers from similar problems. Thus the standard calling and
return sequences already address this problem.
The standard calling sequence (either that used for high-level
languages such as FORTRAN or PL/I, or the 'call' and 'short_call'
operators used in ALM) store the caller's IR along with the
return address. The standard return sequence (either that used
for high-level languages or the 'return' and 'short_return' ALM
operators) restore the caller's IR upon return. Thus we need not
worry that, subsequent to a call, we will be in the wrong mode.
Things are not quite so good for the called procedure, though.
The standard entry sequence used for a high-level language clears
the IR, so a high-level language procedure cannot inherit the
wrong mode. However, an entry sequence is not required of an ALM
procedure (it can be entered through a 'segdef' entry point) and
the standard ALM 'entry' operator does not clear the IR. Thus an
ALM procedure inherits the mode of its caller. This is really
not much of a problem because there is so little use of ALM in
application programs and most use does not involve floating point
operations. Moreover, an ALM procedure that does floating point
operations would probably have to be modified to work in HFP mode
anyway, since the bit representation of the floating point
constants is different in HFP mode.
MTB-687 Multics Technical Bulletin
Fortran HFP mode.
From the above description, it may seem like the only change
needed to the object (other than to get the proper bit
representation for HFP constants) is to add code to manage the
HFP bits in the IR and MR appropriately. Unfortunately, this is
not quite true because the instructions needed to convert between
fixed point and floating point are different in HFP mode than in
BFP mode.
3.3.2.2. Integer to Floating Point
The mantissa of an HFP number in the EAQ register does not have
an integral number of hexadecimal digits. It has a sign bit
followed by 15.75 hex digits. In order to understand how
conversions are done between HFP numbers and integers (the only
kind of fixed point numbers in Fortran), we must analyze the code
used in the BFP case.
First let's consider the problem of floating an integer. Suppose
that we have an integer in the AQ and we would like to get the
floating point equivalent in the EAQ. In BFP mode, this is done
by loading an exponent of 71, and then normalizing:
lde =71b25,du
fad =0.0,du
The exponent is loaded with 71 to indicate that the implied
binary point is 71 bits (i.e. the number of data bits in the AQ)
after the sign bit. Another way to look at this is through an
example. Suppose that the AQ contains the integer 1, which we
want to float. We know that the normalized floating point
representation of 1.0 is 0.5*2**1, so what do we have to do to
the E and the AQ to end up with that value? First, we have to
shift the AQ left 70 bits to change it from 71 zeroes and a one
to a zero, a one and 70 zeroes, which is the mantissa of 0.5. We
must decrement the E one for each left shift of the mantissa if
we are to preserve the same floating point value. We want to end
up with an exponent of 1, so we must start with an exponent of
1+70 = 71.
Now let's consider the same example in HFP mode. In HFP mode,
the normalized floating point representation of 1.0 is 1/16 *
16**1. Thus we must shift the AQ left 67 bits (since a mantissa
of 1/16 has a zero sign bit, 3 leading zeroes, a one and 67
trailing zeroes). In HFP mode, the exponent is to the base 16,
so when we normalize, a change of 1 in the exponent requires a
shift of 4 in the mantissa to balance it. Obviously, loading the
exponent with a suitable value and then normalizing will not work
in the HFP case, because we need the normalization to result in a
left shift of 67, which is impossible since that is not a
Multics Technical Bulletin MTB-687
Fortran HFP mode.
multiple of 4. The closest we can come with normalization is a
left shift of 68, which would give us a mantissa of 1/8 instead
of 1/16. A normalization that yielded a left shift of 68 would
decrement the exponent by 68/4 = 17. Thus if we start with the
integer 1 in the AQ and 18 in the E, a normalization will give us
a mantissa of 1/8 and an exponent of 1. The exponent is right,
but the mantissa is too big by a factor of 2. We can fix the
mantissa by multiplying by 0.5. Therefore, the HFP mode code to
convert an integer in the AQ to the equivalent floating point
number in the EAQ is:
lde =18b25,du
fad =0.0,du
fmp =0.5,du
3.3.2.3. Floating Point to Integer
Now let's consider the inverse operation: truncating a floating
point number in the EAQ to an integer in the AQ. For
concreteness, suppose that the EAQ contains 1.5 and we wish to
truncate it, thus ending up with the integer 1 in the AQ. In BFP
mode, this can be accomplished with one instruction, an
unnormalized floating add of a floating point number with an
exponent of 71 and a mantissa of zero: 'ufa =71b25,du'.
In order to perform the addition, the hardware must first get
both addends to have the same exponent. It does this by
incrementing the exponent and right shifting the mantissa of the
addend with smaller exponent until the exponents are the same.
In our example, the EAQ starts off containing a normalized 1.5,
which means that E contains 1 and AQ contains 0.75 (i.e. a zero,
two ones and 69 zeroes). We are adding a number with an exponent
of 71 and a mantissa of zero, so we must add 70 to E to get the
exponents the same. This is balanced by a right shift of 70,
which leaves the AQ with 71 zeroes and a one. Then we add the
mantissae, which has no effect on the AQ since the other addend
has a mantissa of zero. Thus we are left with an integer one in
the AQ, as desired.
It looks like all that is needed to truncate the EAQ is a 'ufa
=71b25,du'. Unfortunately, this is not quite the case; numbers
in the range [-1, 0) have exponents of zero or less (-1.0 is
represented as -1.0*2**0 rather than -0.5*2**1, as might be
expected), which causes the method to breakdown because such
numbers are insignificant with respect to a number with an
exponent of 71 (i.e. the result will be 0 rather than the
desired -1).
MTB-687 Multics Technical Bulletin
Fortran HFP mode.
This problem can be compensated for by preceding the 'ufa' with
the instructions:
dufa almost_one
dufs almost_one
where the location 'almost_one' contains the largest double
precision floating point number less than 1.0 (i.e. the number
with an exponent of 0 and a mantissa consisting of a sign of zero
followed by 63 ones). These instructions have no effect if E is
greater than zero. If E is zero or less, these instructions
change E to 1 and right shift the AQ appropriately.
A more obvious (but less elegant) solution to the problem is to
split the conversion into two cases, according to whether the
number is positive or negative. If the number is positive, we
don't have a problem. If the number is negative, we truncate its
absolute value and then negate. The code for this looks like:
fad =0.0,du " set indicators for EAQ
tmi 3,ic " if EAQ is positive:
ufa =71b25,du
tra 4,ic " else:
fneg
ufa =71b25,du
negl
Now, let's consider the same example in HFP mode. In HFP mode
the normalized representation of 1.5 is 3/32 * 16**1. Thus E
contains 1 and AQ contains a zero sign bit followed by 3 zeroes,
2 ones and 66 more zeroes. We want to end up with the integer 1
in the AQ, so we must arrange for the AQ to be shifted right 67
places. We can't do that with a 'ufa' of an appropriate
unnormalized zero, though, because 67 is not a multiple of 4.
However, if we first multiply the EAQ by 2, we would then have to
shift right 68 bits, and that we can do with a 'ufa'. Since the
exponent is 1, we will need an exponent 17 larger to force a
right shift of 68. Thus the desired code is:
fmp =o002100,du
ufa =18b25,du
Multics Technical Bulletin MTB-687
Fortran HFP mode.
Note that we cannot use '=2.0,du' in the multiply because 2.0 has
a different bit pattern in HFP mode. Since ALM does not support
HFP constants, we had to specify the desired value in octal.
Also note that the 'fmp' will result in a right shift of 3 of the
AQ and an increment of E if the magnitude of the mantissa was
greater than or equal to a half. This doesn't matter though,
since either the number was too big to represent as a 71-bit
integer or the discarded bits will be in the fraction part (in
which case they don't matter).
We have a similar kind of problem with the 'ufa =18b25,du' as in
the BFP case: Numbers in the range [-0.5, 0) will give us a
result of 0 instead of -1. We fix this in exactly the same way
as in the BFP case, so the final code to truncate the EAQ is:
fmp =o002100,du
dufa almost_one
dufs almost_one
ufa =18b25,du
Or, we can solve the problem, as with the BFP case, by using the
negative of the truncation of the absolute value of a negative:
fmp =o002100,du
tmi 3,ic " if positive:
ufa =18b25,du
tra 4,ic " else:
fneg
ufa =18b25,du
negl
3.3.2.4. Calls to 'pl1_operators_'
Another area which may have to differ in HFP mode is where the
compiler generates calls to operators in 'pl1_operators_' rather
than generating inline code. Some of the operators will work
equally well in either mode (e.g. the operator to multiply two
complex numbers), while others require different code in HFP mode
(e.g. the operator that implements the FORTRAN DMOD builtin
function). There are two obvious ways to address this problem:
define new operators to handle the cases where existing operators
won't work in HFP mode, or overload the existing operators in a
similar way to the way the hardware overloads the floating point
opcodes. The advantage of adding operators is that that ensures
that we will not degrade performance of BFP programs. The
advantage of overloading the existing operators is that we then
don't have to make the compiler aware of any new operators, thus
reducing the changes needed to the compiler.
MTB-687 Multics Technical Bulletin
Fortran HFP mode.
It turns out that it is possible to overload the existing
operators in such a way that there is no performance penalty on
BFP programs, so we have opted for that choice. How this is done
is described in section 6.2, which details the changes needed to
'pl1_operators_' to support HFP mode. Thus, we end up adding
only two operators to 'pl1_operators_' ('enter_BFP_mode' and
'enter_HFP_mode') which have the dual function of putting the
processor in the right mode and telling 'pl1_operators_' which
(BFP or HFP) operator is desired when an overloaded operator is
subsequently called.
3.3.2.5. DU address modification
There is one more area of the object code that needs changing in
HFP mode: the use of 'du' modifications to supply floating point
constants in some of the boilerplate code sequences generated by
the compiler. This turns out to be not much of a problem because
zero and numbers in the ranges [-1.0, -0.5) and [0.5, 1.0) have
exactly the same bit patterns in both modes. In fact, there is
only one 'du' constant outside these ranges that is used, and
that is the constant 1.0. Moreover, it only occurs in floating
point add or subtract operations, so we need only use the
constant -1.0 with the inverse operation (e.g. we replace 'fad
=1.0,du' by 'fsb =-1.0,du').
Multics Technical Bulletin MTB-687
Fortran HFP mode.
3.4. Other Compiler Changes
Most of the compiler deals only with integer and character data.
However, there are a few parts of the compiler that deal with
floating point data derived from the source being compiled. This
manipulation of floating point data falls into two classes,
constant coercion and constant expression folding. For example,
consider the statement:
parameter(pi_by_2 = 3.14159/2)
This statement requires the compiler to coerce the integer
constant 2 to the floating point constant 2.0 and then fold the
constant expression 3.14159/2.0 into the constant 1.570795 which
is then stored as the value of the parameter 'pi_by_2'.
Constant coercion and constant expression folding occur in the
parser for nonexecutable statements and in the optimizer for
executable statements. Constant coercion also occurs in both
code generators when the code for mixed mode expressions is
generated.
The compiler does not handle constant coercion and constant
expression folding in a uniform way. In the parser, most of such
operations are handled by calls to routines in a segment called
'fort_parm_math'. 'form_parm_math' is a collection of routines
that are coded in FORTRAN and which perform constant coercion and
the basic mathematical operations of addition, subtraction,
multiplication, division, exponentiation and negation. In the
rest of the compiler, such operations are performed by inline
PL/I code.
It is easy to alter the parts of the compiler that manipulate
floating point constants through calls to the 'fort_parm_math'
routines so that they handle both BFP and HFP data correctly.
All that is needed is to provide two copies of these FORTRAN
coded routines, one compiled in BFP mode and the other compiled
in HFP mode. It turns out to be especially simple to select the
right set of routines because these routines are not referenced
directly but rather through arrays of entry variables that have
been initialized to these routines. Thus all that is needed to
accomplish BFP or HFP routine selection is to change the entry
variable initialization code.
Since PL/I does not support HFP arithmetic, the parts of the
compiler that manipulate floating point data via inline PL/I code
must be changed to call the same FORTRAN coded routines as used
in the parser.
MTB-687 Multics Technical Bulletin
Fortran HFP mode.
4. Runtime Library Changes
The only routines in the FORTRAN runtime library that have
anything to do with floating point (and hence may need changes)
are the I/O routines and some external builtin functions (such as
'amod_').
Multics Technical Bulletin MTB-687
Fortran HFP mode.
4.1. I/O Routine Changes
The I/O routines are all in the module 'fortran_io_.pl1'. There
are both obvious and subtle changes required to these routines.
The obvious changes are that formatted, list-directed and
namelist I/O must be altered to handle the different bit
representations of BFP and HFP numbers. The subtle changes occur
because of the way routines in 'fortran_io_' are invoked.
Although 'fortran_io_' is coded as a normal PL/I procedure, it is
not accessed (except for the very first call) in the normal way.
This is because it was discovered that when the I/O routines were
called in the standard way, 30% of the CPU time was spent in the
call and return code. Unfortunately, the method used to call and
return from I/O procedures does not save the IR before the call
and restore the IR upon return, nor does it set BFP mode for the
I/O routines. In order to fix these problems, we must make
changes in 'pl1_operators_', 'fortran_io_' and 'return_to_user'
(a small ALM module in the runtime library which is used to
effect the nonstandard return from the I/O routines).
'pl1_operators' must be changed to save and then clear the IR
before calling the I/O routines. Code in 'fortran_io_' that
changes the return address from the place in 'pl1_operators_'
where the I/O routines are called to the place in the user's
program where 'pl1_operators_' was called must be moved into
'pl1_operators_'. This is because the way 'fortran_io_' does it
clears the location where 'pl1_operators_' stores the IR.
'return_to_user' must be changed to restore the IR upon return.
(For more information on the nonstandard calling/return sequence
used for FORTRAN I/O routines, consult the comments following the
journalization information in 'fortran_io_.pl1'.)
MTB-687 Multics Technical Bulletin
Fortran HFP mode.
4.2. External Builtin Routine Changes
The compiler tries to deal with calls to builtin functions in the
most efficient manner possible. If a builtin function is fairly
simple (e.g. ABS), the compiler evaluates the function with
inline code. The more complicated functions (e.g. SIN) are
evaluated by calling the appropriate operator in
'pl1_operators_'. This handling of the builtin functions
presents the compiler with a problem when it finds the name of a
builtin is being passed as an argument to a subprogram. The
compiler resolves this problem by passing the entry address of a
FORTRAN or PL/I procedure that implements that function. Such a
procedure is referred to as an External Builtin Routine.
The External Builtin Routines needed by FORTRAN are kept in one
of two places: 'bound_math_' and 'bound_fort_runtime_'.
Builtins that are common to both FORTRAN and PL/I (e.g. SIN) are
generally found in 'bound_math_'. The rest are found in
'bound_fort_runtime_'. The routines in 'bound_math_' are all
coded in PL/I, while those in 'bound_fort_runtime_' are coded in
FORTRAN.
There are three types of external builtin routines in
'bound_math_'. The first type are those for which the PL/I
compiler generates inline code (e.g. 'abs'). The second type
are those which PL/I evaluates by a call to an operator in
'pl1_operators_' (e.g. 'sin'). The third type are those for
which PL/I generates a standard procedure call (e.g. 'sinh').
The routines in 'bound_math_' which are of the first two types
are merely shells that reference the appropriate PL/I builtin to
evaluate the function. For example, the routine that corresponds
to the single precision sine builtin is just:
sin_: proc (arg) returns (float bin);
dcl arg float bin;
return (sin (arg));
end sin_;
The routines in 'bound_math_' of the third kind contain actual
PL/I code to evaluate the function. They are, in fact, the sole
definition of those functions on Multics. This may seem to
contradict the previous assertion that FORTRAN replaces builtin
calls either by inline code or a call to 'pl1_operators_', since
FORTRAN uses some of these routines that are coded in PL/I (e.g.
SINH).
Multics Technical Bulletin MTB-687
Fortran HFP mode.
The conflict is resolved as follows. There are operators in
'pl1_operators_' for all the routines which are defined as PL/I
procedures. FORTRAN uses these operators, but PL/I makes direct
calls to the procedures. These operators merely call (via the
routine 'fort_math_ops_') the appropriate PL/I procedure to
evaluate the function.
This raises the question: If some of the operators in
'pl1_operators_' that evaluate FORTRAN builtins call external
PL/I procedures to do the work, what do the other operators do?
The answer is that they transfer to ALM routines which evaluate
the functions very efficiently. Thus, another way to look at the
three types of builtin is to say that type 1 is evaluated by
compiler-generated inline code, type 2 is evaluated by calling an
ALM coded routine, and type 3 is evaluated by calling a PL/I
coded routine.
From the above discussion it should be clear that we must do
quite a lot of work to add HFP versions of the FORTRAN builtins.
Let's examine what needs to be done for the HFP version according
to which of the three types the BFP version fits into.
Builtins for which the compiler generates inline code are easy.
All we have to do is code a FORTRAN function that just calls the
builtin. Since these external functions will only be used by
FORTRAN, we should put them into 'bound_fort_runtime_'. For
consistency, it makes sense that we ought to provide FORTRAN
functions in 'bound_fort_runtime_' for the corresponding BFP
builtins, rather than have some of the BFP versions in
'bound_math_' and the rest of the BFP versions and all the HFP
versions in 'bound_fort_runtime_'.
Since the FORTRAN compiler requires that all routines in a
compilation unit be compiled in the same floating point base, we
cannot put both the BFP and the HFP external builtin routines
into 'fort_math_builtins_' (the object in 'bound_fort_runtime_'
that currently contains all the FORTRAN coded external builtin
routines.) Thus we will replace 'fort_math_builtins_' by three
objects: 'fort_bfp_builtins_' (which contains all the floating
point routines to be compiled in BFP mode), 'fort_hfp_builtins_'
(which contains the HFP version of the routines in
'fort_bfp_builtins_') and 'fort_int_builtins_' (which contain all
the integer builtins, which happen to also be the only
nonfloating point builtins).
MTB-687 Multics Technical Bulletin
Fortran HFP mode.
Builtins for which the compiler generates operator calls are less
easy because we must implement the routines that the operators
call to work in HFP mode. These routines could be implemented
all in ALM, or all in FORTRAN (but not PL/I, since it doesn't
support HFP) or some in ALM and some in FORTRAN. It makes sense
to try to keep the HFP implementation as close as possible to the
BFP one, so we will implement the HFP version of each BFP routine
that is currently coded in ALM as another entry point in the same
ALM module. Similarly, we will code the HFP version of each BFP
routine that is currently coded in PL/I in FORTRAN and place
these in 'fort_hfp_builtins_'.
There is one more change that we must make to support the
movement of all FORTRAN external builtin functions into
'bound_fort_runtime_': We must change references in
'pl1_operators_' to these routines from the form 'entry' to the
form 'segment$entry' so that we can distinguish the BFP versions
from the HFP versions, and also so we can distinguish the
'bound_fort_runtime_' BFP versions from those in 'bound_math_'.
This means that we need to update the bind file for
'bound_fort_runtime_' to include add-names for
'fort_bfp_builtins_', 'fort_hfp_builtins_' and
'fort_int_builtins_'. We can also discard the add-names for the
routines currently in 'fort_math_builtins_', thus cleaning up the
name space associated with FORTRAN.
Multics Technical Bulletin MTB-687
Fortran HFP mode.
5. System Support Changes
There are a number of routines in 'bound_library_wired_' on which
FORTRAN programs depend, and some of these must be changed to
support HFP programs. Some of these changes have been alluded to
in previous sections. In this section, we will formalize what
these changes are.
MTB-687 Multics Technical Bulletin
Fortran HFP mode.
5.1. 'any_to_any_' Changes
We have already mentioned that the compiler must be updated to
support input of HFP numbers and the runtime I/O routines must be
updated to support both input and output of such numbers. We did
not specify how this should be done. We shall now do so.
The compiler currently generates the binary representation of the
floating point numbers it finds in the source via the PL/I
'convert' builtin function. The runtime routines use PL/I
picture variable assignments and calls to the system 'assign_'
routine to accomplish its conversion between binary and decimal
representation of floating point number. All of these methods
resolve to calls to the system routine 'any_to_any_' to perform
the actual conversions.
The task of 'any_to_any_' is to convert data from one format to
another (e.g. from 'fixed bin (35)' to 'float bin (63)'). The
input to 'any_to_any_' is a pointer to the data to be converted;
the type, scale and precision of the data; a pointer to the
result buffer; and the type, scale and precision of the result.
The type is specified by the codes defined in
'std_descriptor_types.incl.pl1' (i.e. the same codes that are
used in symbol tables).
Since all type conversions are currently handled in
'any_to_any_', it makes sense that that is where we ought to
place the type conversions we need to support HFP numbers. If we
do that, we can access these conversions through the 'assign_'
system routine, but not through the PL/I 'convert' builtin nor
through PL/I picture assignments, since PL/I does not support
HFP. Thus we need to change the compiler and runtime to perform
conversions through explicit calls to 'assign_'. The upgrade to
any_to_any_ to support HFP numbers has been done and is discussed
in MTB121.
Multics Technical Bulletin MTB-687
Fortran HFP mode.
5.2. 'pl1_operators_' Changes
As previously seen a substantial number of new operators are
required to support HFP mode. Most of these were added by
overloading existing operators (i.e. by making existing
operators work differently according to whether we are in BFP or
HFP mode). This is very convenient because it maintains the
natural parallelism of the two modes and it minimizes the number
of new operators needed -- we added only two: 'enter_BFP_mode'
and 'enter_HFP_mode'. The best part, though, is that we were
able to do this without degrading the performance of nonHFP
programs.
Actually, the concept of having overloaded operators is not new
to 'pl1_operators_'; it is already used to support the 'trace'
command. We have merely extended the concept. In order to see
how it works, we must first review how the operators are invoked.
A look at the code generated by FORTRAN or PL/I will show that
the operators are called by a 'tsx0 pr0|N', where the offset 'N'
determines which operator is called. For example, 'tsx0 pr0|68'
invokes 'fl2_to_fx2', the operator which truncates a floating
point number in the EAQ to an integer in the AQ. The contents of
PR0 is the address of a vector of transfer instructions inside of
'pl1_operators_'. The Nth element of the vector is a transfer to
the code that implements operator N. For example, element 68 of
the vector is 'tra fl2_to_fx2'. The operators are referenced
indirectly through the transfer vector so that changes in
'pl1_operators_' that alter the location of code for an operator
will not be visible outside of 'pl1_operators_'.
It seems very reasonable to keep the address of the
'pl1_operators_' transfer vector in a register, considering how
often the operators are called, but how does that address get
into PR0 in the first place? A look at the code generated by
FORTRAN or PL/I will show that PR0 is referenced in many places,
but it will not appear to be set anywhere. Actually, it is set
as a result of the very first executable code in the routine,
which looks like this:
eax7 M
epp2 pr7|28,*
tsp2 pr2|N
This code is a call to one of the entry operators in
'pl1_operators_'. The first instruction loads X7 with the number
of words of stack space required by the routine. The second
instruction loads PR2 with the address of the 'pl1_operators_'
transfer vector from word 28 of the stack header. (Part of the
Multics stack discipline is that PR7 points to the base of the
MTB-687 Multics Technical Bulletin
Fortran HFP mode.
stack after a call instruction and that word 28 of the stack
points to the 'pl1_operators_' transfer vector.) The third line
invokes an entry operator which then creates the stack frame and
fills in the stack frame header. While doing this, PR0 and word
28 of the stack frame header are set to the address of the
'pl1_operators_' transfer vector.
So far we have been talking as though there was a single transfer
vector inside of 'pl1_operators_'. In fact, there are two. The
second one is almost the same as the first. They differ only in
the address of the transfer instructions for a few of the
operators which have to do with external procedure calls. The
purpose of these differences is to support the 'trace' command.
When the 'trace' command is invoked, it changes the pointer at
word 28 of the stack header to point to the second transfer
vector. Then when a procedure is invoked, it calls a different
entry operator than usual. This entry operator does all the
things that the ordinary one does (except that it sets PR0 and
the pointer at word 28 of the stack frame header to point to the
second transfer vector). However, it also calls a routine in
'trace' so that the trace information can be gathered.
Thus the operators associated with external procedure calls are
overloaded. Whether the normal or trace version of these
operators is executed depends on which transfer vector PR0 points
to. We have generalized this scheme to now include four
different transfer vectors corresponding to the four possible
modes of operation: untraced BFP, untraced HFP, traced BFP and
traced HFP.
There are a couple of important details about the transfer
vectors that we neglected above for simplicity. The first is
that there are tables of constants and some simple code that are
referenced from FORTRAN and PL/I code by negative offsets from
PR0. These tables and code are duplicated in the same relative
order before each of the transfer vectors.
The other important detail is the placement of the various tables
and transfer vectors. The system keeps part of 'pl1_operators_'
in wired memory and part in paged memory. This is because
'pl1_operators_' is called by the paging system and so we must
insure that anything that could be used while processing a page
fault is in wired memory. The tables and transfer vector
associated with untraced BFP operation together with the code and
data they reference are in the wired portion. The tables and
transfer vectors associated with traced or HFP operation,
together with the code and data used only in those modes is
placed in the paged portion in order to minimize the size of the
Multics Technical Bulletin MTB-687
Fortran HFP mode.
wired portion. This means than when an HFP program first starts,
a few page faults will occur to get the needed parts of the paged
portion of 'pl1_operators_' into memory. After that, the
operators should be accessed sufficiently often by the code that
the HFP portion will stay in memory and performance will be the
same as though it were wired.
Now that we have described how the appropriate version of an
overloaded operator is selected, we will describe which operators
need to be overloaded. Since it is intended to restrict support
of HFP to only FORTRAN, we could restrict our attention to only
those operators which are used by FORTRAN. However, we decided
that it would not take much more work to do all the operators
(which would be of future benefit if FORTRAN started using more
operators, or if more languages were made to support HFP), so we
added HFP support for all operators.
For each operator, we had to ask two questions: Would HFP mode
affect the operation of the operator, and would the operator
affect HFP mode (i.e. clear the HFP bit in the caller's IR)? If
the answer to either question was yes, then the operator had to
be overloaded, or protected or both.
There are basically two kinds of operator: Those for which all
the work is done inside of 'pl1_operators_' (e.g. 'fl2_to_fx2')
and those for which some or all of the work is done by calling
routines outside of 'pl1_operators_'. By finding each use of an
'ldi' and each use of a floating point instruction and then
examining the context of that use, we were able to determine that
the only operators of the first kind that were affected by or
affected HFP mode were: 'fx2_to_fl2', 'divide2', 'fl2_to_fx2',
'fl2_to_fxscaled', 'fort_mdfl1', 'fort_dmod', 'rfb2_to_cflb1',
'mdfl1', 'trunc_fl', 'floor_fl', 'ceil_fl',
'nearest_whole_number' and 'nearest_integer'. The problem with
all of these operators, except for 'divide2', was that they
converted between fixed and floating point, which must be done
differently in HFP mode (i.e. these operators need to be
overloaded).
The 'divide2' operator divides one double precision integer by
another. The problem in this operator was that it used the 'fno'
instruction to count the leading zeroes in the AQ, which doesn't
work in HFP mode. Rather than overload this rather large
operator, we decided to make the operator work in both modes by
replacing the single 'fno' instruction with:
sti sp|temp_indicators
ldi =0,dl
fno
ldi sp|temp_indicators
MTB-687 Multics Technical Bulletin
Fortran HFP mode.
These instructions save the current mode, force BFP mode, do the
normalize and then restore the current mode. The three extra
instructions are very fast and this operator is very long so the
extra time is negligible.
There are two kinds of operator that call outside of
'pl1_operators_'. The first kind are those that implement math
functions such as 'sin'. These operators consist of only a
transfer to an approriate entry point in one of the other
routines in 'bound_library_wired_'. Thus the only change for HFP
mode for these operators was to use a different entry point,
namely the one for the HFP version of the math routine. (The
next section describes the changes needed to the math routines to
support these new entry points.)
The second kind of operator that calls external routines is just
an interface between the calling program and the routine. The
operator merely takes care of setting up the call. Operators
that implement I/O operations are examples of this type. The
usual problem with this type of operator was that the discipline
of saving the caller's IR along with the return address was not
followed. Thus when the external routine tried to restore the
caller's IR upon return, zero was loaded and HFP mode would get
turned off. The FORTRAN I/O operators not only had that trouble,
but they also required that the IR be cleared (BFP mode set)
after saving the IR, since they were entered in a nonstandard way
that bypassed the entry sequence code that would ordinarily clear
the IR for the called routine.
Besides protecting and overloading existing operators, we needed
to add two new operators to get us into and out of HFP mode:
'enter_HFP_mode' and 'enter_BFP_mode'. Actually, FORTRAN only
requires 'enter_HFP_mode' since we can rely on the restoration of
the caller's indicators upon return from a procedure to revert
HFP mode when appropriate. We would only need 'enter_BFP_mode'
if we wanted to swap between HFP and BFP mode inside a routine.
However, we decide to add 'enter_BFP_mode' for the sake of
completeness.
The 'enter_BFP_mode' operator is very simple. It clears the
indicators, then resets PR0 and the pointer at word 24 of the
stack frame header to the address of the appropriate (normal BFP
or traced BFP) transfer vector.
The 'enter_HFP_mode' operator is a little more complicated
because, unlike BFP mode, it is not always possible to enter HFP
mode. First the indicators are cleared, except for the HFP bit,
which is set. This will put us into HFP mode if the HFP bit is
also set in the MR. Next we check if we are in HFP mode. This
is done by loading an HFP 1.0 into the EAQ and normalizing. If
Multics Technical Bulletin MTB-687
Fortran HFP mode.
we are in HFP mode, the normalize will do nothing, but if we are
in BFP mode, it will shift the AQ left 3 bits and decrement the
exponent by 3 because in BFP mode an HFP 1.0 looks like an
unnormalized 0.125. If the normalize shows that we are still in
BFP mode, we call 'hcs_$set_hexfp_control' to set the HFP bit in
the MR. This may not be possible since the user may not be
allowed by the System Administrator to use HFP mode or the
hardware may not support it. Thus, when control is returned, we
check if we are now in HFP mode. If not, we signal the
restartable condition 'unable_to_enter_HFP_mode'. If we are
restarted, we assume that the user has done something (e.g.
pleaded with the System Administrator) that will now allow us to
enter HFP mode, so we repeat the call to 'hcs_$set_hexfp_control'
and again check to see if we are in HFP mode. If not, we again
signal 'unable_to_enter_HFP_mode'. This process continues ad
nauseum until we are able to enter HFP mode or the user gives up
restarting us.
If we ever establish that we are in HFP mode, we then proceed as
with the 'enter_BFP_mode' operator and set PR0 and the pointer at
word 24 of the stack frame header to the address of the
appropriate (normal HFP or traced HFP) transfer vector.
Note that switching between BFP and HFP mode does not affect
trace mode: If you were tracing before switching, you will still
be tracing afterward, and if you weren't tracing, you won't be
either after the switch.
MTB-687 Multics Technical Bulletin
Fortran HFP mode.
5.3. Math Routines
New versions of the standard math routines are required to
support the larger range and different bit representation of HFP
numbers. There are presently two kinds of math routines: those
coded in ALM and those coded in PL/I.
Multics Technical Bulletin MTB-687
Fortran HFP mode.
5.3.1. ALM-Coded Math Routines
All the math routines that are currently coded in ALM are in
'bound_library_wired_'. Related routines are entry points into
the same module. For example, the module 'sine_' contains the
entry points 'cosine_degrees_', 'cosine_radians_',
'sine_degrees_' and 'sine_radians_'. We shall extend this
principle by placing the HFP versions of routines in the same
module as the corresponding BFP version. The entry point names
will be derived from the corresponding BFP version by adding the
prefix 'hfp_' (e.g. 'hfp_cosine_degrees_').
Most of the existing math routines are very old and are either
not as fast or not as accurate as they should be. The GCOS math
routines (from which the Multics versions were originally
translated) have been improved considerably. In the near future,
we will be replacing the BFP math routines with translations of
the current GCOS versions, which also support HFP mode.
Therefore, the current addition of HFP math routines will be only
temporary and we won't worry about making them fast. Instead, we
will take the simplest approach possible for their
implementation.
In most cases, the simplest approach is to use mathematical
properties of a function to derive a formula for evaluating the
function in terms of the existing BFP version of the function.
This implies the need to convert between the BFP and HFP
representations of a number. We will answer this need by adding
two new functions, 'bfp_to_hfp_' and 'hfp_to_bfp_' to
'bound_library_wired_'. These new functions will convert a BFP
number in the EAQ to HFP format, and vice versa, over the range
of BFP numbers. If you attempt to convert an HFP number to BFP
format and the number is outside the range of BFP numbers,
'hfp_to_bfp_' will return the largest BFP number of the same
sign.
We will now look at each of BFP math modules and specify how the
HFP entry points will be calculated.
5.3.1.1. 'arc_sine_'
The 'arc_sine_' module has entry points to calculate the single
or double precision value, in degrees or radians, of the inverse
sine and inverse cosine functions. The calculations are actually
performed by the inverse tangent function of two arguments
according to the identities:
acos(x) = atan2(sqrt(1 - x*x), x)
asin(x) = atan2(x, sqrt(1 - x*x))
MTB-687 Multics Technical Bulletin
Fortran HFP mode.
The names 'acos', 'asin', 'atan' and 'sqrt' above should be taken
to be generic. For example, the entry point that calculates the
double precision inverse cosine function in degrees uses the
double precision inverse tangent in degrees and the double
precision square root routines.
We will use the same strategy to calculate the HFP equivalents of
the BFP entry points, except we will use the HFP versions of the
required 'atan2' and 'sqrt' functions.
5.3.1.2. 'double_arc_tangent_'
The 'arc_tangent_' routine has entry points to calculate the
single precision inverse tangent function of 1 or 2 arguments in
degrees or radians. The 'double_arc_tangent_' routines has entry
points to calculate the double precision equivalents of the
'arc_tangent_' entry points.
The inverse tangent function has the very useful properties that
atan(x) = x for very small values of x, and
atan(x) = pi or -pi for very large values of x.
It turns out that these properties hold for values of 'x' that
are well within the range of BFP numbers. Thus we can calculate
the HFP versions of the various entry points by:
if abs(x) < 2**(-35)
then hfp_atan (x) = x;
else hfp_atan (x) = hfp (atan (bfp (x));
where the 'hfp' function returns the HFP representation of the
BFP number which is its argument, and the 'bfp' function returns
the BFP representation of the nearest BFP number to the HFP
number which is its argument.
5.3.1.3. 'call_math_error_'
This is not really a math routine, but rather a subroutine to the
various math routines. It is called to signal the appropriate
error condition when a math routine is called with an invalid
argument. We include it in this list because it needs to be
changed to work correctly in HFP mode. The problem with it is
Multics Technical Bulletin MTB-687
Fortran HFP mode.
that it overwrites the location in the user's stack frame header
where the IR is stored, so that in the user is no longer in HFP
mode after a math error has occurred. It turns out that the
solution is trivial. All that is needed is to remove the
offensive store because the information being stored is never
used!
5.3.1.4. 'double_exponential_'
These modules calculate the exponential function to single and
double precision, respectively. We evaluate the HFP versions in
terms of the corresponding BFP versions by careful use of the
mathematical identity 'exp(a + b) = exp(a)*exp(b)', as follows.
Suppose that we want to calculate 'exp(x)' for some HFP number
'x'. The above identity suggests that we look for a value 'y'
such that 'exp(x-y)' and 'exp(y)' are both easily calculated,
since 'exp(x)' is the same as 'exp(x-y)*exp(y)'. Now, it is very
easy in HFP mode to multiply by a power of 16; we need only add
the power to the exponent of the multiplicand. Thus, it makes
sense to choose 'y' so that 'exp(y)' will be a power of 16. In
other words, we want 'y' to be the logarithm of a power of 16: y
= log(16**N), where 'N' is an integer. From the properties of
logarithms, we know that 'log(16**N)' is the same as 'N*log(16)',
so we have the identity:
exp(x) = exp(x - N*log(16)) * 16**N
The question is: How do we choose 'N' so that 'exp(x -
N*log(16))' is easy to evaluate? The obvious answer is: Choose
it so that 'x - N*log(16)' falls in the domain of the BFP version
of the exponential function (i.e. from -88 to 88,
approximately). This gives us great lattitude in our choice.
However, we also want to get the best possible accuracy, which
means that we want to choose 'N' so that 'x - N*log(16)' is as
close as possible to zero, since that is where the BFP
exponential function is most accurate. If we rewrite 'x -
N*log(16)' as '(x/log(16) - N)*log(16)', it becomes apparent that
we want to choose 'N' to be the nearest integer to 'x/log(16)'.
We want to avoid floating point divisions, since they are both
slower and less accurate than floating point multiplications, so
we are pleased to note that '1/log(16)' is the same as
'log16(e)'.
In summary, if 'x' is an HFP number then 'exp(x)' can be
calculated in terms of the BFP version of the exponential
function by:
hfp(exp(bfp(x - N*log(16)))) * 16**N
MTB-687 Multics Technical Bulletin
Fortran HFP mode.
where 'N = int(x*log16(e))'.
5.3.1.5. 'fort_math_ops_'
This module is an interface used by FORTRAN to call the math
routines that are not coded in ALM. It must be modified to save
the caller's indicators in the standard location in the stack
frame. It also must be changed to include entry points for all
the new non-ALM coded math routines (which will be described in
6.3.2). Finally, it must be changed to reflect the movement from
'bound_math_' to 'bound_fort_runtime_' (described in 5.2) of all
the non-ALM coded math routines used by FORTRAN.
5.3.1.6. 'integer_power_integer_'
This routine raises an integer to an integer power, and so is
unaffected by HFP mode. However, its use has a disastrous
side-effect: It turns off HFP mode. The problem is that the
code does an 'ldi 0,dl' to turn off the overflow and exponent
overflow indicators prior to execution of a repeat instruction
that is set to terminate on overflow. We could fix this problem
by saving the IR and later restoring it, but a simpler solution
is to replace the 'ldi 0,dl' with:
teo 1,ic clear exponent overflow indicator
tov 1,ic clear overflow indicator
5.3.1.7. 'double_logarithm_'
These modules respectively calculate the single and double
precision versions of the logarithm function. There are three
entry points in each, corresponding to logarithms to the bases 2,
e and 10. The HFP versions are based on the BFP base 2 logarithm
routine and the following mathematical identities:
log(x) = log2(x)*log(2)
log10(x) = log2(x)*log10(2)
log2(a*b) = log2(a) + log2(b)
The first two identities above show us that we need only figure
out how to calculate HFP logarithms to the base 2. The third
identity shows us how to reduce the calculation of an HFP base 2
logarithm to a BFP base 2 logarithm, as follows.
Suppose that we want to calculate the base 2 logarithm of the HFP
number 'x'. Let 'M' be the mantissa of 'x', and let 'E' be the
exponent of 'x' in HFP notation. Then:
Multics Technical Bulletin MTB-687
Fortran HFP mode.
log2(x) = log2(M*16**E)
= log2(M) + log2(16**E)
= log2(M) + E*log2(16)
= log2(M) + 4*E
Since the magnitude of 'M' is less than or equal to one, we can
calculate 'log2(M)' by converting 'M' to BFP, using the BFP base
2 logarithm routine, and converting the result back to HFP.
5.3.1.8. 'math_constants_'
This module centralizes the definition of various commonly used
math constants (e.g. 'pi'). It needs to be augmented by the
definitions of the HFP version of each constant. The name of the
HFP version is derived by adding the prefix 'hfp_' to the name of
the corresponding BFP constant (e.g. 'hfp_pi').
5.3.1.9. 'power_'
This module calculates the result of raising a floating point
number to a floating point power. If the floating point power is
integral, it calls 'power_integer_' to calculate the power;
otherwise it uses the mathematical identity 'a**b =
exp(b*log(a))' to calculate the result. We implement the HFP
version of the module the same way, except using HFP versions of
the called routines.
5.3.1.10. 'power_integer_'
This module raises a floating point number to an integer power.
Except for two minor problems, this code will work equally well
in either BFP or HFP mode. The first problem is that the
instruction 'fld =1.0,du' occurs in three places. This will load
8.0 in the HFP case, so alternate code must be used. Since this
instruction is executed only once or twice per invocation of the
routine, we solve the problem by replacing the instruction with:
fld =0.5,du
fad =0.5,du
MTB-687 Multics Technical Bulletin
Fortran HFP mode.
which solves the problem since the 0.5 has the same
representation in HFP mode as in BFP mode. The other problem is
that the instruction 'ldi =0,dl' is used to clear the overflow
and exponent overflow indicators prior to execution of a reapeat
instruction which terminates upon overflow. This is disastrous
because it turns off HFP mode. This problem could be solved by
first saving the IR and later restoring it, but a simpler
approach is to replace the 'ldi =0,dl' with:
teo 1,ic clear exponent overflow indicator
tov 1,ic clear overflow indicator
5.3.1.11. 'double_sine_'
The module 'sine_' has entry points to compute the single
precision sine or cosine of an angle given in degrees or radians.
The module 'double_sine_' has entry points for the double
precision versions of these routines. Since the sine and cosine
functions are periodic with a period of 2*pi, we can evaluate the
HFP versions in terms of the BFP versions by:
hfp_cos (x) = hfp (cos (bfp (mod (x, 2*pi))))
hfp_sin (x) = hfp (cos (bfp (mod (x, 2*pi))))
However, we must be careful when taking the sine of angles close
to zero, since the above will give us 'hfp_sin (x) = 0' for 'x'
that is too small to represent in BFP mode. Thus for small 'x',
we use the approximation:
hfp_sin (x) = x
5.3.1.12. 'double_square_root_'
These modules calculate the square root of a floating point
number. If we think of an HFP number in terms of its hardware
representation, we see that it has the form M*16**E, where M is
either zero or has a magnitude in the range [1/16, 1), and E is
an integer in the range [-128, 127]. Since the square root
operation distributes over multiplication, we have the following
two identities which we can use to reduce the problem of taking
the square root of an HFP number to that of taking the square
root of a related BFP number:
If E is even: sqrt(M*16**E) = sqrt(M)*16**(E/2)
If E is odd: sqrt(M*16*E) = 4*sqrt(M)*16**((E-1)/2)
Multics Technical Bulletin MTB-687
Fortran HFP mode.
Thus we can take HFP square roots via:
if E is even
then factor = 1.0;
else factor = 4.0;
hfp_sqrt = factor*hfp (sqrt (bfp (M)))*16**floor(E/2);
Of course, we don't really multiply by 16**floor(E/2). Instead,
we add floor(E/2) to the exponent of 'factor*hfp (sqrt (bfp
(M)))'.
5.3.1.13. 'double_tangent_'
These modules calculate the single and double precision tangents
of an angles given in degrees or radians. We implement the HFP
equivalents in terms of the HFP sine and cosine routines, using
the mathematical identity:
tan (x) = sin (x)/cos (x)
MTB-687 Multics Technical Bulletin
Fortran HFP mode.
5.3.2. Non-ALM Coded Math Routines
Several of the current BFP math routines are coded in PL/I.
These are kept in 'bound_math_'. Since PL/I does not support HFP
mode, we must code the HFP versions in either FORTRAN or ALM. It
is a lot easier to code in FORTRAN, so we have opted to do so for
this routines. The resulting objects will be kept in
'bound_fort_runtime_'.
The following sections describe the algorithms used for these
routines.
5.3.2.1. 'cabs_'
This module computes the absolute value of a complex number. If
the real part of the complex number is 'x' and the imaginary part
is 'y', then the absolute value is defined as 'sqrt (x*x + y*y)'.
This formula is not a good one to use because it will overflow
too easily. For example, if 'x' is 1e100 and 'y' is zero, the
absolute value is 1e100, but we cannot calculate it with this
formula because 'x*x' overflows.
There is a better way to calculate the absolute value. Suppose
that 'x' is bigger than 'y' in absolute value. Then 'y/x' will
be between -1 and 1. If we factor 'x*x' out of our original
formula, we get:
cabs(x+iy) = abs(x)*sqrt(1 + (y/x)**2)
Now, we are taking the square root of a number between 1 and 2,
so there can be no overflow there. There still can be an
overflow when we multiply by 'abs(x)', but if there is, it is a
legitimate overflow due to the absolute value being outside the
range of HFP numbers.
If the magnitude of 'y' is bigger than that of 'x', we factor
'y*y' out of the original formula, giving:
cabs(x+iy) = abs(y)*sqrt((x/y)**2 + 1)
If the magnitudes of 'x' and 'y' are the same, we can factor out
either 'x*x' or 'y*y', except if they are zero, in which case we
can just return zero for the absolute value.
5.3.2.2. 'ccos_'
This module calculates the cosine of a complex number. We do
this by using the following identity:
Multics Technical Bulletin MTB-687
Fortran HFP mode.
cos (x+iy) = cos(x)*cosh(y) - i*sin(x)*sinh(y)
5.3.2.3. 'cexp_'
This module calculates the value of the number 'e' raised to a
complex power. We do this by using the following identity:
exp(x+iy) = exp(x)*(cos(y) + i*sin(y))
5.3.2.4. 'clog_'
This module calculates the natural logarithm of a complex number.
We do this by using the following identity:
log(x+iy) = log(abs(x+iy)) + i*atan2(y, x)
5.3.2.5. 'dcosh_'
These modules calculate the single and double precision
hyperbolic cosine, which is defined by:
cosh(x) = 0.5*(exp(abs(x)) + 1/exp(abs(x)))
From this definition, we see that if the magnitude of 'x' is
large enough (specifically, if abs(x) > 352.8119), the result is
outside the range of HFP numbers. In that case, we cause an
exponent overflow fault and (if restarted) return the largest HFP
number.
If the magnitude of 'x' is moderately large (i.e. > 9.704 in
single precision, or 22.18 in double precision), 1/exp(abs(x))
will be insignificant with respect to exp(abs(x)), so we can use
the approximation:
cosh(x) = 0.5*exp(abs(x))
For other values of 'x', we can just use the definition:
cosh(x) = 0.5*((exp(abs(x)) + 1/(exp(abs(x)))
5.3.2.6. 'csin_'
This module calculates the sine of a complex number. We do this
by the identity:
sin(x+iy) = sin(x)*cosh(y) + i*sinh(x)*cos(y)
MTB-687 Multics Technical Bulletin
Fortran HFP mode.
5.3.2.7. 'csqrt_'
This module calculates the square root of a complex number. We
do this by the identity:
sqrt(x+iy) = sqrt(abs(x+iy) + x) + i*sign(sqrt(abs(x+iy) -
x), y)
5.3.2.8. 'cxp2_'
This module calculates the value of a complex number raised to a
complex power. We do this using the identity:
(u+iv)**(x+iy) = exp((x+iy)*log(u+iv))
5.3.2.9. 'dsinh_'
These modules calculate the single and double precision
hyperbolic sine functions, which are defined by:
sinh(x) = sign(0.5*(exp(abs(x)) - 1/exp(abs(x))), x)
From this definition, we see that for 'x' of sufficient magnitude
(specifically, if abs(x) > 352.8119), sinh(x) will be outside the
range of HFP numbers. When that is the case, we cause an
exponent overflow fault and (if restarted) returned the largest
HFP number of the same sign as 'x'.
If the magnitude of 'x' is moderately large (i.e. > 9.704 in
single precision or 22.18 in double precision), 1/exp(abs(x))
will be insignificant with respect to exp(abs(x)), so we can use
the approximation:
sinh(x) = sign(0.5*exp(abs(x)), x)
If the magnitude of 'x' is sufficiently small (i.e. < 2.11e-4 in
the single precision case, or 8.06e-10 in the double precision
case), exp(abs(x)) will be so close to 1 that exp(abs(x) -
1/exp(abs(x)) will give us zero. Thus for 'x' of this small a
magnitude, we must find some other way to approximate sinh(x).
We do this by taking the first few terms of the Taylor series
expansion of sinh(x):
sinh(x) = x + x**3/3! + x**5/5! + x**7/7! + x**9/9! + ...
in the single precision case, the first 4 terms are enough, while
in the double precision case, we need the first 7 terms.
Multics Technical Bulletin MTB-687
Fortran HFP mode.
For x whose magnitude is between "sufficiently small" and
"moderately large" (as defined above), we can just use the
definition:
sinh(x) = sign(0.5*(exp(abs(x)) - 1/exp(abs(x))), x)
5.3.2.10. 'double_tanh_'
These modules calculate the single and double precision
hyperbolic tangent function, which is defined by:
tanh(x) = sinh(x)/cosh(x)
We calculate this function exactly as specified by the above
definition.