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.