An important difference between floating-point and fixed-point instruction selection, however, involves the inherent rounding of floating-point calculations. In general, the finite floating-point format is an imperfect representation of a real number, hence the need for rounding. The error introduced during rounding can be analyzed, and various approaches have been developed to reproducibly control this error. The IEEE 754 Standard for Binary Floating-Point Arithmetic (IEEE 754) specifies the most common approach. The PowerPC architecture includes the features needed to conform to IEEE 754. The requirements of IEEE 754, however, restrict the optimizations available during code selection.

The PowerPC architecture's floating-point support includes:

• Load and store instructions.
• Elementary arithmetic (addition, subtraction, multiplication and division) instructions.
• Fused multiply-add instructions. The processor maintains the intermediate product at full precision, an attribute that offers considerable advantage to a clever numerical programmer.
• Support for both single-precision and double-precision.
• Rounding and conversion instructions.
• Comparison instructions.
• Instructions to control the Floating-Point Unit and to test floating-point conditions.
• Hardware support for the IEEE 754 standard, including conforming formats, rounding modes, operations, special values and exceptions.

During a load into a Floating-Point Register, the processor converts 32-bit single-precision values to the 64-bit double-precision format. A full complement of single-precision arithmetic operations accept single-precision values as input and round their results to single-precision. Maintaining the single-precision values in double-precision format in the Floating-Point Registers removes the need for these values to undergo conversion to double-precision. In addition, a double-precision operation can trivially perform the multiplication of two single-precision values to produce a double-precision result.

Some implementations perform all floating-point operations (except division) in double-precision format, rounding the final result to the target precision. Some implementations perform single-precision operations with less latency than the comparable double-precision operations. In summary, there may be timing differences between single-precision and double-precision operations, so consult Appendix B and specific implementation documentation for further information.

The 64-bit integer format fills the entire Floating-Point Register, while the 32-bit integer format resides in the low-order half of the register. These integer values result from a direct load into the Floating-Point Register, a conversion from a floating-point value, or the transferal of the contents of the FPSCR by mffs. The PowerPC architecture includes several instructions to support conversion between integer and floating-point formats in the Floating-Point Registers. Detailed algorithms describing their behavior can be found in Appendix B of Book I of The PowerPC Architecture. Code examples that perform these conversions appear in Section 3.3.8 of this book and Appendix E.3 of Book I of The PowerPC Architecture.

An important optimization involves avoiding the conversion of integer induction variables to floating-point values in floating-point computational loops if the integer values take part in floating-point computations. Integer-to-float conversions are time-consuming. One way to avoid this is to create a floating-point value that increments in lock step with the induction variable. Using this floating-point value in the calculation avoids the conversion.

The Floating Round to Single-Precision (frsp) instruction explicitly rounds the value in a Floating-Point Register to single-precision, with all applicable IEEE 754 exceptions. PowerPC single-precision arithmetic operations automatically round their results to single-precision.

The value of the RN field in the Floating-Point Status and Control Register (FPSCR) determines the rounding mode as follows:

• 00—Round to nearest
• 01—Round toward 0
• 10—Round toward +
• 11—Round toward -

3.3.2 Memory Access

3.3.2.1 Single-Precision Loads and Stores

A single-precision store operation does not round the register value to single-precision, but merely performs a format conversion when copying the data. This format conversion may involve denormalizing the value, but does not include rounding. For the storage of a properly rounded single-precision result, it is the compiler's responsibility either to round the result to single-precision with the frsp instruction prior to the store or to ensure that the value in the register falls in the range of the single-precision format. Storing a double-precision value that is not rounded to single-precision can produce unexpected results when the operand is out of range.

3.3.2.2 Double-Precision Loads and Stores

3.3.2.3 Endian Conversion

3.3.2.4 Touch Instructions

3.3.3 Floating-Point Move Instructions

3.3.4 Computation

3.3.4.1 Setting Status Bits

3.3.4.2 Arithmetic

Multiplication followed by a dependent addition is the most common idiom in floating-point code, so the PowerPC architecture includes a set of fused multiply-add instructions. The multiply-add instructions do not round the multiplication result before performing the addition. This property is useful for many algorithms, some of which appear in subsequent sections. On the other hand, this property does not conform to IEEE 754, which requires rounding after every operation, but the results of the fused multiply-add instructions always differ from IEEE 754 results in a way that would be considered more accurate.

Many PowerPC implementations have floating-point units designed around the fused multiply-add functionality. Regardless of the arithmetic operation, the actual execution involves both a multiplication and an addition. Addition operations include an effective multiply by 1. Multiplication operations include an effective add to 0. In these implementations, the multiply-add operations are faster than separate multiply and add steps.

For the negative multiply-add or multiply-subtract instructions, the rounding occurs prior to the sign inversion. Therefore, when the rounding mode is round toward ± , the fnmadd instruction, defined as -[A × C + B], is not the same as -A × C - B, and the compiler cannot substitute fnmadd for -A × C - B if IEEE 754 compatibility is required. Similar concerns apply to fnmsub and -A × C + B.

On all existing designs, the division operation requires a large number of cycles to execute. In the special case that the divisor is a power of 2, you may replace the division with multiplication by the reciprocal, which can be represented precisely, unless the divisor is a denormal other than 2^-127 (for single-precision) or 2^-1023 (for double-precision).

3.3.4.3 Floating-Point Comparison

If either operand is a NaN, the Floating Compare Ordered (fcmpo) instruction sets VXVC in FPSCR. If interrupts on Invalid Operations are enabled, the processor generates the corresponding interrupt.

3.3.5 FPSCR Instructions

• All floating-point exceptions caused by previous instructions are recorded in the FPSCR.
• All floating-point interrupts generated by previous instructions have been processed.

FPSCR instructions do not set the FEX and VX bits in FPSCR explicitly, but rather these bits represent the logical OR of the exception bits and the Invalid Operation exception bits, respectively. On some implementations, updating fewer than all eight fields of the FPSCR may result in substantially poorer performance than updating all the fields.

3.3.6 Optional Floating-Point Instructions

3.3.6.1 Square Root

3.3.6.2 Storage Access

3.3.6.3 Reciprocal Estimate

3.3.6.4 Reciprocal Square Root Estimate

3.3.6.5 Selection

3.3.7 IEEE 754 Considerations

The PowerPC architecture's hardware support for the IEEE 754 standard includes:

• Formats—32-bit single-precision and 64-bit double-precision.
• Rounding—All four rounding modes: Round to Nearest, Round toward + , Round toward - , and Round toward 0.
• Operations—Add, subtract, multiply, divide, square root (optionally supported by PowerPC architecture), round to single-precision, convert floating-point value to integer word or doubleword, convert integer doubleword to floating-point value, and compare (result is delivered as a condition code).
• Special values—Infinity, signed zero, QNaN and SNaN.
• Exceptions—Invalid Operation, Divide by Zero, Overflow, Underflow, Inexact. The sub-categories of Invalid Operation include: SNaN, - , ÷ , 0 ÷ 0, × 0, Invalid Compare, Software Request, Invalid Square Root, and Invalid Integer Convert.
• Traps—If the IEEE 754 default response to one of the preceding exceptions is not acceptable to the programmer, an interrupt can be enabled.

• Square root (if not directly supported in hardware).
• IEEE 754 remainder.
• Conversion between floating-point and integer formats.
• Rounding of a floating-point value to a floating-point integer.
• Conversions between binary and decimal.
• IEEE 754 recommended functions.
• Support for extended formats (IEEE 754 recommended).

The operating system must provide:

• A floating-point state that is consistent through context changes.
• Compatible interrupt handling.

3.3.7.1 Relaxations

Multiply-Add

Non-IEEE Mode

3.3.8 Data Format Conversion

The code examples in this section often use the optional fsel instruction. If IEEE 754 compatibility is required or if the instruction is not supported on the target processor, replace this instruction with the equivalent comparison-branch sequence.

3.3.8.1 Floating-Point to Integer

In general, floating-point to integer conversions require rounding. The rounding mode is that indicated in the RN field of the FPSCR unless the "z" form of the convert to integer instruction is used. Then, regardless of the rounding mode, the value is rounded toward zero, effectively truncating the fractional part of the floating-point value.

The code sequences in Figures 3-48 and 3-49 convert a floating-point value to a 32-bit and a 64-bit signed integer value, respectively. The Floating Convert To Integer instructions convert the floating-point value to an integer in the Floating-Point Register. The Store Float-Point Double instructions directly copy the bit pattern of this integer to memory. Then, the integer is loaded into a General-Purpose Register.

Figure 3-48. Convert Floating-Point to 32-Bit Signed Integer Code Sequence

	32-Bit Implementation
	fctiw[z]	FR2,FR1	# convert to integer
	stfd	FR2,disp(R1)	# copy unmodified to memory
	lwz	R3,disp+4(R1)	# load the low-order 32 bits
	64-Bit Implementation
	fctiw[z]	FR2,FR1	# convert to integer
	stfd	FR2,disp(R1)	# copy unmodified to memory
	lwa	R3,disp+4(R1)	# load low-order 32 bits with sign
			# extension

Figure 3-49. Convert Floating-Point to 64-Bit Signed Integer Code Sequence

	64-Bit Implementation Only
	fctid[z]	FR2,FR1	# convert to doubleword integer
	stfd	FR2,disp(R1)	# store float double
	ld	R3,disp(R1)	# load double word

The code sequences in Figures 3-50 and 3-51 convert a floating-point value to a 32-bit and a 64-bit unsigned integer value, respectively. The Floating Convert To Integer instruction converts a floating-point value to a signed integer. Using this instruction to convert to an unsigned integer requires some adjustments. If the floating-point input value is negative, replace it with 0. If the floating-point input value is greater than the maximum for the format (2^M - 1), replace it with 2^M - 1, where M is the number of bits in the resulting integer (32 or 64). If the floating-point input value lies in the range 2^M-1 to 2^M - 1, the range where a signed integer is negative, reduce the input value by 2^M-1. Then, convert the floating-point value to an integer using a Floating Convert To Integer instruction. If the input was greater than 2^M-1, add 2^M-1. Negative values return 0; values exceeding 2^M - 1 return 2^M - 1. The 64-bit implementation of the conversion to a 32-bit unsigned integer uses the fctid[z] instruction to avoid the extra code associated with the 2³¹ to 2³² range.

Figure 3-50. Convert Floating-Point to 32-Bit Unsigned Integer Code Sequence

	# FR0 = 0.0
	# FR1 = value to be converted
	# FR3 = 232 - 1
	# FR4 = 231
	# R3 = returned result
	# disp = displacement from R1
	32-Bit Implementation
	fsel	FR2,FR1,FR1,FR0	# use 0 if < 0
	fsub	FR5,FR3,FR1	# use 232-1 if >= 232
	fsel	FR2,FR5,FR2,FR3
	fsub	FR5,FR2,FR4	# subtract 2**31
	fcmpu	cr2,FR2,FR4
	fsel	FR2,FR5,FR5,FR2	# use diff if >= 231
			# next part same as conversion to
			# signed integer word
	fctiw[z]	FR2,FR2	# convert to integer
	stfd	FR2,disp(R1)	# copy unmodified to memory
	lwz	R3,disp+4(R1)	# load low-order word
	blt	cr2,$+8	# add 231 if input
	xoris	R3,R3,0x8000	# was >= 231
	64-Bit Implementation
	fsel	FR2,FR1,FR1,FR0	# use 0 if < 0
	fsub	FR4,FR3,FR1	# use 232-1 if >= 232
	fsel	FR2,FR4,FR2,FR3
			# next part same as conversion to
			# signed integer word except
			# convert to double
	fctid[z]	FR2,FR2	# convert to doubleword integer
	stfd	FR2,disp(R1)	# copy unmodified to memory
	lwz	R3,disp+4(R1)	# load low-order word, zero extend

Figure 3-51. Convert Floating-Point to 64-Bit Unsigned Integer Code Sequence

	64-Bit Implementation Only
	# FR0 = 0.0
	# FR1 = value to be converted
	# FR3 = 264 - 2048
	# FR4 = 263
	# R4 = 263
	# R3 = returned result
	# disp = displacement from R1
	fsel	FR2,FR1,FR1,FR0	# use 0 if < 0
	fsub	FR5,FR3,FR1	# use max if > max
	fsel	FR2,FR5,FR2,FR3
	fsub	FR5,FR2,FR4	# subtract 263
	fcmpu	cr2,FR2,FR4	# use diff if >= 263
	fsel	FR2,FR5,FR5,FR2
			# next part same as conversion to
			# signed integer doubleword
	fctid[z]	FR2,FR2	# convert to integer
	stfd	FR2,disp(R1)	# copy unmodified to memory
	ld	R3,disp(R1)	# load doubleword integer value
	blt	cr2,$+8	# add 263 if input
	add	R3,R3,R4	# was >= 263

3.3.8.2 Integer to Floating-Point

In a 32-bit implementation, you may convert a 32-bit integer to a floating-point value as follows. Flip the integer sign bit and place the result in the low-order part of a doubleword in memory. Create the high-order part with sign and exponent fields such that the resulting doubleword value interpreted as a hexadecimal floating-point value is 0x1.00000dddddddd×10^D, where 0xdddddddd is the hexadecimal sign-flipped integer value. Then, load the doubleword as a floating-point value. Subtract the hexadecimal floating-point value 0x1.0000080000000×10^D from the previous value to generate the result.

The 64-bit implementations possess the Floating Convert From Integer Doubleword (fcfid) instruction, which converts an integer to a floating-point value. Both 64-bit implementation examples transfer the value to the floating-point unit, where the fcfid instruction is used.

The conversion of a 32-bit integer to a floating-point value is always exact. The conversion of a 64-bit integer to a floating-point value may require rounding, which conforms to the mode indicated in the RN field of the FPSCR.

Figure 3-52. Convert 32-Bit Signed Integer to Floating-Point Code Sequence

	32-Bit Implementation
	# FR2 = 0x4330000080000000
	addis	R0,R0,0x4330	# R0 = 0x43300000
	stw	R0,disp(R1)	# store upper half
	xoris	R3,R3,0x8000	# flip sign bit
	stw	R3,disp+4(R1)	# store lower half
	lfd	FR1,disp(R1)	# float load double of value
	fsub	FR1,FR1,FR2	# subtract 0x4330000080000000
	64-Bit Implementation
	extsw	R3,R3	# extend sign
	std	R3,disp(R1)	# store doubleword integer
	lfd	FR1,disp(R1)	# load integer into FPR
	fcfid	FR1,FR1	# convert to floating-point value

Figure 3-53. Convert 64-Bit Signed Integer to Floating-Point Code Sequence

	64-Bit Implementation Only
	std	R3,disp(R1)	# store doubleword
	lfd	FR1,disp(R1)	# load float double
	fcfid	FR1,FR1	# convert to floating-point integer

The code sequences in Figure 3-54 convert an unsigned 32-bit integer to a floating-point value. These code examples parallel those given for the signed case in Figure 3-52.

In a 32-bit implementation, construct the floating-point value in memory, as before, but do not flip the sign bit. Subtract the hexadecimal floating-point value 0x1.0000080000000×10^D from the loaded value to produce the result.

In a 64-bit implementation, replace the sign extension in the signed version with a zero extension performed by the rldicl instruction.

Figure 3-54. Convert 32-Bit Unsigned Integer to Floating-Point Code Sequence

	32-Bit Implementation
	# FR2 = 0x4330000000000000
	addis	R0,R0,0x4330	# R0 = 0x43300000
	stw	R0,disp(R1)	# store high half
	stw	R3,disp+4(R1)	# store low half
	lfd	FR1,disp(R1)	# float load double of value
	fsub	FR1,FR1,FR2	# subtract 0x4330000000000000
	64-Bit Implementation
	rldicl	R0,R3,0,32	# zero extend value
	std	R0,disp(R1)	# store doubleword value to memory
	lfd	FR1,disp(R1)	# load value to FPU
	fcfid	FR1,FR1	# convert to floating-point integer

The code sequence in Figure 3-55 converts a 64-bit unsigned integer value to a floating-point value in a 64-bit implementation. The first example converts the two 32-bit halves separately and combines the results at the end with a multiply-add instruction.

The second example presents an alternative, shorter sequence that can be used if the rounding mode is Round toward ± , or if the rounding mode does not matter. This example converts the entire integer doubleword in a single step with a subsequent addition to correct for negative values. The addition operation may cause inaccuracy in certain rounding modes.

Figure 3-55. Convert 64-Bit Unsigned Integer to Floating-Point Code Sequence

	64-Bit Implementation (All Rounding Modes)
	# FR4 = 232
	rldicl	R2,R3,32,32	# isolate high half
	rldicl	R0,R3,0,32	# isolate low half
	std	R2,disp(R1)	# store high half
	std	R0,disp+8(R1)	# store low half
	lfd	FR2,disp(R1)	# load high half
	lfd	FR1,disp+8(R1)	# load low half
	fcfid	FR2,FR2	# convert each half to floating-
	fcfid	FR1,FR1	# point integer (no round)
	fmadd	FR1,FR4,FR2,FR1	# (232)*high + low
			# (only add can round)
	Alternate Version (Only for Round toward ± )
	# FR2 = 264
	std	R3,disp(R1)	# store doubleword
	lfd	FR1,disp(R1)	# load float double
	fcfid	FR1,FR1	# convert to floating-point integer
	fadd	FR4,FR1,FR2	# add 264
	fsel	FR1,FR1,FR1,FR4	# if R3 < 0

3.3.8.3 Rounding to Floating-Point Integer

You may round a floating-point value to a 32-bit integer as follows. For non-negative values, add 2⁵² and then subtract it, allowing the floating-point hardware to perform the rounding using the rounding mode given by the RN field. For negative values, add -2⁵² and then subtract it. If you require a rounding mode different from that specified in the RN field, this technique would require modification of RN.

In a 64-bit implementation, you may round a floating-point value to a 64-bit floating-point integer by converting to a 64-bit integer and then converting back to a floating-point value. If VXCVI is set, indicating an invalid integer conversion, the value does not fit into the 64-bit integer format (i.e., the value is greater than or equal to 2⁶⁴). These large values have no fractional part.

Figure 3-56. Round to Floating-Point Integer Code Sequence

	32-Bit Implementation
	# FR0 = 0.0
	# FR2 = 0x43300000 = 252
	# FR3 = 0xC3300000 = -252
	fcmpu	cr6,FR1,FR0
	bt	cr6[lt],lab	# branch if value < 0.0
	fcmpu	cr7,FR1,FR2
	bt	cr7[gt],exit	# input was floating-point integer
	fadd	FR4,FR1,FR2	# add 252
	fsub	FR1,FR4,FR2	# subtract 252
	b	exit
lab:
	fcmpu	cr7,FR1,FR3
	bt	cr7[lt],exit	# input was floating-point integer
	fadd	FR4,FR1,FR3	# add -252
	fsub	FR1,FR4,FR3	# subtract -252
exit:
	64-Bit Implementation
	mtfsb0	23	# clear VXCVI
	fctid[z]	FR3,FR1	# convert to fixed-point integer
	fcfid	FR3,FR3	# convert back again
	mcrfs	7,5	# transfer VXCVI to CR
	bf	31,$+8	# skip if VXCVI was 0
	fmr	FR3,FR1	# input was floating-point integer

3.3.9 Floating-Point Branch Elimination

The code examples in Figure 3-57 perform a greater than or equal to 0.0 comparison for both comparison-branch and fsel object code. The fsel code does not cause an exception for a NaN, unlike the comparison-branch code.

The use of the Condition Register Logical instruction reduces the number of branches from two to one. Unlike the integer case, greater than or equal to is not equivalent to not less than because the result of the comparison could be unordered.

Figure 3-57. Greater Than or Equal to 0.0 Code Example

	Fortran Source Code
	if (a .ge. 0.0) then x = y
	else x = z
	Branching Assembly Code
	# (FR0) = 0.0
	fcmpo	cr7,FRa,FR0	# causes exception if a is NaN
	cror	cr5[fe],cr7[fe],cr7[fg]	# temp = (a .gt. 0.0)
			# .or. (a .eq. 0.0)
	bf	cr5[fe],lab1	# if (.not.temp) then branch
	fmr	FRx,FRy	# x = y
	b	lab2
lab1:
	fmr	FRx,FRz	# x = z
lab2:
	fsel Assembly Code
	fsel	FRx,FRa,FRy,FRz	# if (a .ge. 0.0) then x = y
			# else x = z
			# no exception if a is NaN

The code examples in Figure 3-58 perform a greater than 0.0 comparison. The fsel code uses the comparison (-a 0), while reversing the assignments. If a is a NaN, the result of the assignment is reversed from the branch case. Again, an SNaN or a QNaN does not cause an exception in the fsel case.

Figure 3-58. Greater Than 0.0 Code Example

	Fortran Source Code
	if (a .gt. 0.0) then x = y
	else x = z
	Branching Assembly Code
	# (FR0) = 0.0
	fcmpo	cr7,FRa,FR0	# causes exception if a is NaN
	bf	cr7[fg],lab1	# if .not.(a .gt. 0.0) then branch
	fmr	FRx,FRy	# x = y
	b	lab2
lab1:
	fmr	FRx,FRz	# x = z
			# if (a is NaN) then x = z
lab2:
	fsel Assembly Code
	fneg	FRs,FRa	# s = -a
	fsel	FRx,FRs,FRz,FRy	# if (s .ge. 0.0) then x = z
			# else x = y
			# if (a is NaN) then x = y
			# no exception if a is NaN

The code examples in Figure 3-59 perform an equal to 0.0 comparison. Two fsel instructions are used in series so that x is set to y only if a 0.0 and a 0.0. Otherwise, x is set to z. Again, a NaN does not cause an exception.

Figure 3-59. Equal to 0.0 Code Example

	Fortran Source Code
	if (a .eq. 0.0) then x = y
	else x = z
	Branching Assembly Code
	# (FR0) = 0.0
	fcmpo	cr7,FRa,FR0	# causes exception if a is NaN
	bf	cr7[fe],lab1	# if (a .ne. 0) then branch
	fmr	FRx,FRy	# x = y
	b	lab2
lab1:
	fmr	FRx,FRz	# x = z
lab2:
	fsel Assembly Code
	fsel	FRx,FRa,FRy,FRz	# if (a .ge. 0.0) then x = y
			# else x = z
	fneg	FRs,FRa	# s = -a
	fsel	FRx,FRs,FRx,FRz	# if (s .ge. 0.0) then x = x
			# else x = z
			# no exception if a is NaN

The code examples in Figure 3-60 perform the min(a,b) function. You may compute this function by comparing 0.0 to the difference between a and b, and using the result to select a or b conditionally. If either a or b is a NaN, a different result from the comparison-branch version may occur. Moreover, the subtraction operation may cause exceptions that do not occur for the comparison-branch case. The code for the maximum function is analogous.

Figure 3-60. Minimum Code Example

	Fortran Source Code
	x = min(a,b)
	Branching Assembly Code
	fcmpo	cr7,FRa,FRb	# causes exception if a or b is NaN
	bf	cr7[fl],lab1	# if (.not.(a .lt. b)) then branch
	fmr	FRx,FRa	# min = a
	b	lab2
lab1:
	fmr	FRx,FRb	# min = b
			# if (a or b is NaN) then min = b
lab2:
	fsel Assembly Code
	fsub	FRs,FRa,FRb	# s = a - b
			# causes exception if
			# a or b is SNaN
	fsel	FRx,FRs,FRb,FRa	# if (s .ge. 0.0) then min = b
			# else min = a
			# if (a or b is NaN) then min = a
			# no exception if a or b is QNaN

The code examples in Figure 3-61 compare a and b for equality. Two fsel instructions are used in series so that x is set to y only if (a - b) 0.0 and (a - b) 0.0. Otherwise, x is set to z. This comparison between two values is similar to the minimum function, having similar incompatibilities with the comparison-branch approach.

Figure 3-61. a Equal To b Code Example

	Fortran Source Code
	if (a .eq. b) then x = y
	else x = z
	Branching Assembly Code
	fcmpo	cr7,FRa,FRb	# causes exception if a or b is NaN
	bf	cr6[fe],lab1	# if (a .eq. b) then branch
	fmr	FRx,FRy	# x = y
	b	lab2
lab1:
	fmr	FRx,FRz	# x = z
lab2:
	fesl Assembly Code
	fsub	FRs,FRa,FRb	# s = a - b
			# causes exception if
			# a or b is SNaN
	fsel	FRx,FRs,FRy,FRz	# if (s .ge. 0.0) then x = y
			# else x = z
	fneg	FRs,FRs	# s = -s
	fsel	FRx,FRs,FRx,FRz	# if (s .ge. 0.0) then x = z
			# else x = z
			# no exception if a or b is QNaN

3.3.10 DSP Filters

Figure 3-62. Matrix Product: C Source Code

	for(i=0;i<10;i++){
	for(j=0;j<10;j++){
	c[i][j]=0;
	for(k=0;k<10;k++){
	c[i][j] = c[i][j] + a[i][k] * b[k][j];
	}}}

The central fragment/inner loop code is:

c[i][j] = c[i][j] + a[i][k] × b[k][j];

assume:

Ra -> points to array a - 8

Rb-> points to array b - 80

Rc-> points to array c - 8

Figure 3-63 shows the object code for the double-precision floating-point case. The multiply-add instructions and update forms of the load and store instructions combine to form a tight loop. This example represents an extremely naive compilation that neglects loop transformations.

Figure 3-63. Double-Precision Matrix Product: Assembly Code

dloop:
	lfdu	FR0,8(Ra)	# load a[i][k], bump to next element
	lfdu	FR1,80(Rb)	# load b[k][j], bump to next element
	fmadd	FR2,FR0,FR1,FR2	# c[i][j] + a[i][k] * b[k][j]
	bdnz	dloop
	stfdu	FR2,8(Rc)	# store c[i][j] element

3.3.11 Replace Division with Multiplication by Reciprocal

Figure 3-64. Convert Division to Multiplication by Reciprocal Code Example

	Fortran Source Code
	do 10 j = 1,n
10	x(j) = x(j)/y
	Assembly Code
	# FR1 = 1.0
	# FR5 = y
	# R3 = address of x(j) - 8
	# CTR = n
	fdiv	FR2,FR1,FR5	# calculate the reciprocal
			# FR2 = rec = 1/y
loop:
	lfd	FR3,8(R3)	# load x(j)
	fmul	FR4,FR3,FR2	# multiply by the reciprocal
			# FR4 = q = x(j)*rec
	fneg	FR6,FR4
	fmadd	FR7,FR6,FR5,FR3	# calculate residual
			# FR5 = res = x(j) - q*y
	fmadd	FR8,FR7,FR2,FR4	# correct x(j)*rec
			# x(j)/y = x(j)*rec + (res * rec)
	stfdu	FR8,8(R3)	# store x(j)/y with update
	bdnz	loop

3.3.12 Floating-Point Exceptions

These software alternatives involve placing the processor into Ignore Exceptions mode, but enabling the desired exceptions through the Floating-Point Status and Control Register (FPSCR). Interrogation of the FEX bit in the FPSCR reveals whether an enabled exception has occurred. This interrogation may be carried out using run-time library functions, such as those described in the AIX API. Another approach uses a recording floating-point instruction form followed by a conditional branch to test the copy of FEX placed in CR1. If FEX is set indicating that an enabled exception has occurred, control is transferred to the appropriate trap handling routine. This control transfer can be managed with an unconditional trap instruction that generates an interrupt, which transfers control to the usual floating-point interrupt handler or with a simple call to a service routine. Figure 3-65 shows an example that uses a recording floating-point instruction, a conditional branch to interrogate FEX, and an unconditional trap to transfer control to the interrupt handler. The interrupt handler must fully diagnose the cause of the interrupt, which involves interrogating the FPSCR register.

Figure 3-65. Precise Interrupt in Software Code Example

	fma.	FR1,FR1,FR2,FR3	# recording floating-point operation
	bt	cr1[fex],lab2	# branch if there is an exception
lab1:
	...
lab2:
	trap		# unconditional trap to handler
	b	lab1	# return

Depending on the frequency of floating-point operations and the cost of running in Precise mode for the implementation under consideration, the overhead associated with the application of software-controlled trapping to all floating-point instructions may degrade performance more than running in Precise mode. The potential benefit of these approaches may occur only when a single test for exceptions at the end of a series of floating-point operations suffices. In this case, if an exception occurs during the execution of the series, you may not be able identify the operation that caused it.