Numerical Computation Guide

Appendix E

Standards Compliance

The compilers, header files, and libraries in the Sun WorkShop Compilers products for the Solaris environment support multiple standards: System V Interface Definition (SVID), Edition 3, X/Open and ANSI C. Accordingly, the mathematical library libm and related files have been modified so that C programs comply with the standards. Users' programs are usually not affected, because the differences primarily involve exception handling.

SVID History

To understand the differences between exception handling according to SVID and the point of view represented by the IEEE Standard, it is necessary to review the circumstances under which both developed. Many of the ideas in SVID trace their origins to the early days of UNIX, when it was first implemented on mainframe computers. These early environments have in common that rational floating-point operations +, -, * and / are atomic machine instructions, while sqrt, conversion to integral value in floating-point format, and elementary transcendental functions are subroutines composed of many atomic machine instructions.

Because these environments treat floating-point exceptions in varied ways, uniformity could only be imposed by checking arguments and results in software before and after each atomic floating-point instruction. Because this has too great an impact on performance, SVID does not specify the effect of floating-point exceptions such as division by zero or overflow.

Operations implemented by subroutines are slow compared to single atomic floating-point instructions; extra error checking of arguments and results has little performance impact; so such checking is required by the SVID. When exceptions are detected, default results are specified, errno is set to EDOM for improper operands, or ERANGE for results that overflow or underflow, and the function matherr() is called with a record containing details of the exception. This costs little on the machines for which UNIX was originally developed, but the value is correspondingly small because the far more common exceptions in the basic operations +, -, * and / are completely unspecified.

IEEE 754 History

The IEEE Standard explicitly states that compatibility with previous implementations was not a goal. Instead, an exception handling scheme was developed with efficiency and users' requirements in mind. This scheme is uniform across the simple rational operations (+, -, * and /), and more complicated operations such as remainder, square root, and conversion between formats. Although the Standard does not specify transcendental functions, the framers of the Standard anticipated that the same exception handling scheme would be applied to elementary transcendental functions in conforming systems.

Elements of IEEE exception handling include suitable default results and interruption of computation only when requested in advance.

SVID Future Directions

The current SVID, ("Edition 3" or "SVR4"), identifies certain directions for future development. One of these is compatibility with the IEEE Standard. In particular a future version of the SVID replaces references to HUGE, intended to be a large finite number, with HUGE_VAL, which is infinity on IEEE systems. HUGE_VAL, for instance, is returned as the result of floating-point overflows. The values returned by libm functions for input arguments that raise exceptions are those in the IEEE column in TABLE E-1. In addition, errno no longer needs to be set.

SVID Implementation

The following libm functions provide operand or result checking corresponding to SVID. The sqrt function is the only function that does not conform to SVID when called from a C program that uses the libm in-line expansion templates through -xlibmil, because this causes the hardware instruction for square root, fsqrt[sd], to be used in place of a function call.

TABLE E-1   Exceptional Cases and libm Functions 

Function	errno	error message	SVID	X/Open	IEEE
acos(|x|>1)	EDOM	DOMAIN	0.0	0.0	NaN
acosh(x<1)	EDOM	DOMAIN	NaN	NaN	NaN
asin(|x|>1)	EDOM	DOMAIN	0.0	0.0	NaN
atan2((+-0,+-0)	EDOM	DOMAIN	0.0	0.0	+-0.0,+-pi
atanh(|x|>1)	EDOM	DOMAIN	NaN	NaN	NaN
atanh(+-1)	EDOM/ERANGE	SING	+-HUGE (EDOM)	+-HUGE_VAL (ERANGE)	+-infinity
cosh overflow	ERANGE	-	HUGE	HUGE_VAL	infinity
exp overflow	ERANGE	-	HUGE	HUGE_VAL	infinity
exp underflow	ERANGE	-	0.0	0.0	0.0
fmod(x,0)	EDOM	DOMAIN	x	NaN	NaN
gamma(0 or -integer)	EDOM	SING	HUGE	HUGE_VAL	infinity
gamma overflow	ERANGE	-	HUGE	HUGE_VAL	infinity
hypot overflow	ERANGE	-	HUGE	HUGE_VAL	infinity
j0(|x| > X_TLOSS)	ERANGE	TLOSS	0.0	0.0	correct answer
j1(|x| > X_TLOSS)	ERANGE	TLOSS	0.0	0.0	correct answer
jn(|x| > X_TLOSS)	ERANGE	TLOSS	0.0	0.0	correct answer
lgamma(0 or -integer)	EDOM	SING	HUGE	HUGE_VAL	infinity
lgamma overflow	ERANGE	-	HUGE	HUGE_VAL	infinity
log(0)	EDOM/ERANGE	SING	-HUGE (EDOM)	-HUGE_VAL (ERANGE)	-infinity
log(x<0)	EDOM	DOMAIN	-HUGE	-HUGE_VAL	NaN
log10(0)	EDOM/ERANGE	SING	-HUGE (EDOM)	-HUGE_VAL (ERANGE)	-infinity
log10(x<0)	EDOM	DOMAIN	-HUGE	-HUGE_VAL	NaN
log1p(-1)	EDOM/ERANGE	SING	-HUGE (EDOM)	-HUGE_VAL (ERANGE)	-infinity
log1p(x<-1)	EDOM	DOMAIN	NaN	NaN	NaN
pow(0,0)	EDOM	DOMAIN	0.0	1.0 (no error)	1.0 (no error)
pow(NaN,0)	EDOM	DOMAIN	NaN	NaN	1.0 (no error)
pow(0,neg)	EDOM	DOMAIN	0.0	-HUGE_VAL	+-infinity
pow(neg, non-integer)	EDOM	DOMAIN	0.0	NaN	NaN
pow overflow	ERANGE	-	+-HUGE	+-HUGE_VAL	+-infinity
pow underflow	ERANGE	-	+-0.0	+-0.0	+-0.0
remainder(x,0)	EDOM	DOMAIN	NaN	NaN	NaN
scalb overflow	ERANGE	-	+-HUGE_ VAL	+-HUGE_VAL	+-infinity
scalb underflow	ERANGE	-	+-0.0	+-0.0	+-0.0
sinh overflow	ERANGE	-	+-HUGE	+-HUGE_VAL	+-infinity
sqrt(x<0)	EDOM	DOMAIN	0.0	NaN	NaN
y0(0)	EDOM	DOMAIN	-HUGE	-HUGE_VAL	-infinity
y0(x<0)	EDOM	DOMAIN	-HUGE	-HUGE_VAL	NaN
y0(x > X_TLOSS)	ERANGE	TLOSS	0.0	0.0	correct answer
y1(0)	EDOM	DOMAIN	-HUGE	-HUGE_VAL	-infinity
y1(x<0)	EDOM	DOMAIN	-HUGE	-HUGE_VAL	NaN
y1(x > X_TLOSS)	ERANGE	TLOSS	0.0	0.0	correct answer
yn(n,0)	EDOM	DOMAIN	-HUGE	-HUGE_VAL	-infinity
yn(n,x<0)	EDOM	DOMAIN	-HUGE	-HUGE_VAL	NaN
yn(n, x> X_TLOSS)	ERANGE	TLOSS	0.0	0.0	correct answer

General Notes on Exceptional Cases and libm Functions

TABLE E-1 lists all the libm functions affected by the standards. The value X_TLOSS is defined in <values.h>. SVID requires <math.h> to define HUGE as MAXFLOAT, which is approximately 3.4e+38. HUGE_VAL is defined as infinity in libc. errno is a global variable accessible to C and C++ programs.

<errno.h> defines 120 or so possible values for errno; the two used by the math library are EDOM for domain errors and ERANGE for range errors. See intro(3) and perror(3).

• The ANSI C compiler switches -Xt, -Xa, -Xc, -Xs, among other things, control the level of standards compliance that is enforced by the compiler. Refer to cc(1) for a description of these switches.
• As far as libm is concerned, -Xt and -Xa cause SVID and X/Open behavior, respectively. -Xc corresponds to strict ANSI C behavior.

An additional switch -<x>libmieee, when specified, returns values in the spirit of IEEE 754. The default behavior for libm and libsunmath is to be SVID-compliant on Solaris 2.6, Solaris 7, and Solaris 8.

• For strict ANSI C (-Xc), errno is set always, matherr() is not called, and the X/Open value is returned.
• For SVID (-Xt or -Xs), the function matherr() is called with information about the exception. This includes the value that is the default SVID return value.

A user-supplied matherr() could alter the return value; see matherr(3m). If there is no user-supplied matherr(), libm sets errno, possibly prints a message to standard error, and returns the value listed in the SVID column of TABLE E-1 on page 255.

• For X/Open (-Xa), the behavior is the same as for the SVID, in that matherr() is invoked and errno set accordingly. However, no error message is written to standard error, and the X/Open return values are the same as IEEE return values in many cases.
• For the purposes of libm exception handling, -Xs behaves the same as -Xt. That is, programs compiled with -Xs use the SVID compliant versions of the libm functions listed in TABLE E-1 on page 255.
• For efficiency, programs compiled with inline hardware floating-point do not do the extra checking required to set EDOM or call matherr()if sqrt encounters a negative argument. NaN is returned for the function value in situations where EDOM might otherwise be set.

Thus, C programs that replace sqrt() function calls with fsqrt[sd] instructions conform to the IEEE Floating-Point Standard, but may no longer conform to the error handling requirements of the System V Interface Definition.

Notes on libm

SVID specifies two floating-point exceptions, PLOSS (partial loss of significance) and TLOSS (total loss of significance). Unlike sqrt(-1), these have no inherent mathematical meaning, and unlike exp(+-10000), these do not reflect inherent limitations of a floating-point storage format.

PLOSS and TLOSS reflect instead limitations of particular algorithms for fmod and for trigonometric functions that suffer abrupt declines in accuracy at definite boundaries.

Like most IEEE implementations, the libm algorithms do not suffer such abrupt declines, and so do not signal PLOSS. To satisfy the dictates of SVID compliance, the Bessel functions do signal TLOSS for large input arguments, although accurate results can be safely calculated.

The implementations of sin, cos, and tan treat the essential singularity at infinity like other essential singularities by returning a NaN and setting EDOM for infinite arguments.

Likewise SVID specifies that fmod(x,y) is be zero if x/y overflows, but the libm implementation of fmod, derived from the IEEE remainder function, does not compute x/y explicitly and hence always delivers an exact result.

LIA-1 Conformance

In this section, LIA-1 refers to ISO/IEC 10967-1:1994 Information Technology - Language Independent Arithmetic - Part 1: Integer and floating-point arithmetic.

The C and Fortran-77 compilers (cc and f77) contained in the Sun WorkShop Compilers 6.0 release conform to LIA-1 in the following senses (paragraph letters correspond to those in LIA-1 section 8):

a. TYPES (LIA 5.1): The LIA-1 conformant types are C int and Fortran INTEGER. Other types may conform as well, but they are not specified here. Further specifications for specific languages await language bindings to LIA-1 from the cognizant language standards organizations.

b. PARAMETERS (LIA 5.1):

 #include <values.h> /* defines MAXINT */
 #define TRUE 1
 #define FALSE 0
 #define BOUNDED TRUE
 #define MODULO TRUE
 #define MAXINT 2147483647
 #define MININT -2147483648
 	 logical bounded, modulo
 	 integer maxint, minint
 	 parameter (bounded = .TRUE.)
 	 parameter (modulo = .TRUE.)
 	 parameter (maxint = 2147483647)
 	 parameter (minint = -2147483648)

d. DIV/REM/MOD (LIA 5.1.3):

C / and %, and Fortran / and mod(), provide DIVtI(x,y) and REMtI(x,y). Also, modaI(x,y) is available with this code:

 int modaI(int x, int y) {
 	 int t = x % y;
 	 if (y < 0 && t > 0)
 	 t -= y;
 	 else if (y > 0 && t < 0)
 	 t += y;
 	 	 return t;
 	 }

or this:

 	 integer function modaI(x, y)
 	 integer x, y, t
 	 t = mod(x, y)
 	 if (y .lt. 0 .and. t .gt. 0) t = t - y
 	 if (y .gt. 0 .and. t .lt. 0) t = t + y
 	 modaI = t
 	 return
 	 end

i. NOTATION (LIA 5.1.3): The following table shows the notation by which the LIA integer operations may be realized.

TABLE E-2   LIA-1 Conformance - Notation

LIA	C	Fortran if different
addI(x,y)	x+y
subI(x,y)	x-y
mulI(x,y)	x*y
divtI(x,y)	x/y
remtI(x,y)	x%y	mod(x,y)
modaI(x,y)	see above
negI(x)	-x
absI(x)	#include <stdlib.h> abs(x)	abs(x)
signI(x)	#define signI(x) (x > 0 ? 1 : (x < 0 ? -1 : 0))	see below
eqI(x,y)	x==y	x.eq.y
neqI(x,y)	x!=y	x.ne.y
lssI(x,y)	x<y	x.lt.y
leqI(x,y)	x<=y	x.le.y
gtrI(x,y)	x>y	x.gt.y
geqI(x,y)	x>=y	x.ge.y

The following code shows the Fortran notation for signI(x).

 	 integer function signi(x)
 	 integer x, t
 	 if (x .gt. 0) t=1
 	 if (x .lt. 0) t=-1
 	 if (x .eq. 0) t=0
 	 return
 	 end

j. EXPRESSION EVALUATION: By default, when no optimization is specified, expressions are evaluated in int (C) or INTEGER (Fortran) precision. Parentheses are respected. The order of evaluation of associative unparenthesized expressions such as a + b + c or a * b * c is not specified.

k. METHOD OF OBTAINING PARAMETERS: Include the definitions above in the source code.

n. NOTIFICATION: Integer exceptions are x/0 and x%0 or mod(x,0). By default, these exceptions generate SIGFPE. When no signal handler is specified for SIGFPE, the process terminates and dumps memory.

o. SELECTION MECHANISM: signal(3) or signal(3F) may be used to enable user exception handling for SIGFPE.