Previous Next Contents Index Doc Set Home

The Math Libraries


This chapter describes the math libraries provided with Sun WorkShop Compilers. Some attention is given to IEEE supporting functions and to functions that convert data between IEEE and non-IEEE formats.

Some information can also be obtained from the man page Intro(3).

This chapter has the following organization:

Standard Math Library

page 32

Additional Math Libraries

page 33

Single, Double, and Long Double Precision

page 37

IEEE Support Functions

page 38

Implementation Features of libm and libsunmath

page 47

libc Support Functions

page 50

Standard Math Library

The libm math library contains the functions required by the various standards to which the Solaris operating system conforms. This library is bundled with the Solaris operating system in two forms: libm.a, the static version, and, the shared version.

The default directories for a standard installation of libm are:





The default directories for standard installtion of the header files for libm are:




Table 3-1 lists the functions in libm with the names used for calling them from a C program.

Table  3-1 Contents of libm 

Function Name

Algebraic functions

cbrt, hypot, sqrt

Elementary transcendental


asin, acos, atan, atan2, asinh, acosh, atanh

exp, expm1, pow,

log, log1p, log10,

sin, cos, tan, sinh, cosh, tanh

Higher transcendental functions

bessel(j0, j1, jn, y0, y1, yn),

erf, erfc, gamma, lgamma, gamma_r, lgamma_r

Integral rounding functions

ceil, floor, rint

IEEE standard recommended functions

copysign, fmod, ilogb, nextafter, remainder, scalbn, fabs

IEEE classification functions

isnan, finite1

Old style floating-point functions

logb, scalb, significand

Error handling routine (user-defined)


1 On Solaris systems, this function is provided in the C library libc.

Note that the functions gamma_r and lgamma_r are reentrant versions of gamma and lgamma.

Additional Math Libraries

Sun Math Library

The library libsunmath is part of the libraries supplied with all Sun language products. The library libsunmath contains a set of functions that were incorporated in previous versions of libm from Sun.

The default directories for a standard installation of libsunmath are:



The default directories for standard installation of the header files for libsunmath are:



Table 3-2 lists the functions in libsunmath; the names are those used for calling the double precision version of functions from a C program

Table  3-2 Contents of libsunmath  

Function Name

Functions from Table 3-1

single, extended and quadruple precision available, except for matherr

Elementary transcendental


exp2, exp10,



Trigonometric functions (degree arguments/values)

asind, acosd, atand, atan2d,

sind, cosd, sincosd, tand

Trigonometric functions scaled in

asinpi, acospi, atanpi, atan2pi,

sinpi, cospi, sincospi, tanpi

Trigonometric functions with double precision argument reduction

asinp, acosp, atanp,

sinp, cosp, sincosp, tanp

Financial functions

annuity, compound

Integral rounding functions

aint, anint, irint, nint

IEEE standard recommended functions


IEEE classification functions

fp_class, isinf, isnormal, issubnormal, iszero

Functions that supply useful IEEE values

min_subnormal, max_subnormal,

min_normal, max_normal,

infinity, signaling_nan, quiet_nan

Random number generators

d_addran_, d_addrans_, d_lcran_, d_lcrans_, d_shufrans_, i_addran_, i_addrans_,
i_lcran_, i_lcrans_, i_shufrans_, r_addran_, r_addrans_, r_lcran_, r_lcrans_,
r_shufrans_, u_addrans, u_lcrans_, u_shufrans_

Data conversion


Control rounding mode and floating-point exception flags


Floating-point trap handling


Show status


Toggle hardware between standard and
nonstandard modes (advisory)

standard_arithmetic, nonstandard_arithmetic

Optimized Libraries

Optimized versions of some of the routines in libm are provided in the library libmopt. On SPARC and Intel systems, optimized versions of some of the support routines in libc are provided in the library libcopt. Finally, on SPARC and PowerPC systems, alternate forms of some libc support routines are provided in libcx.

The default directories for a standard installation of libmopt, libcopt, and libcx are:


/opt/SUNWspro/SC4.2/lib/libcopt.a (SPARC and Intel only)

/opt/SUNWspro/SC4.2/lib/libcx.a (SPARC and PowerPC only)

The routines contained in libcopt are not intended to be called by the user directly. Instead, they replace support routines in libc that are used by the compiler.

The routines contained in libmopt replace corresponding routines in libm. The libmopt versions are generally noticeably faster. Note that unlike the libm versions, which can be configured to provide any of ANSI/POSIX, SVID, X/Open, or IEEE-style treatment of exceptional cases, the libmopt routines only support IEEE-style handling of these cases. (See Appendix F, "Standards Compliance.") To see which libm routines have counterparts in libmopt, use the nm command to search for all global symbols defined in libmopt.a.

Both libmopt and libcopt may be linked with a program by specifying the
-xlibmopt flag at link time. Note that the SunOS 5.x C compiler, however, does not accept the -xlibmopt flag. To use both libmopt and libcopt with the Solaris 2.x C compiler, link with both -lmopt and -lcopt.

SPARC and Power PC: On Solaris 2.x, libcx contains faster versions of the 128-bit quadruple precision floating point arithmetic support routines. These routines are not intended to be called directly by the user; instead, they are called by the compiler. On Solaris 2.x, the Fortran 77 compiler links with libcx automatically unless the -nocx flag is used, but the C compiler does not automatically link with libcx. To use libcx with C programs, link with
-lcx. (Note: A shared version of libcx, called, is also provided. This version can be preloaded at run time by setting the environment variable LD_PRELOAD to the full path name of the file.)

Vector Math Library (SPARC only)

On SPARC systems running the Solaris 2.x operating environment, the library libmvec provides routines that evaluate common mathematical functions for an entire vector of arguments.

The default directory for a standard installation of libmvec is:



Table 3-3 lists the functions in libmvec. Note that only double precision versions of these functions are available.

Table  3-3 Contents of libmvec

Function Name

Algebraic functions

vhypot_, vz_abs_

Exponential and related functions

vexp_, vlog_, vpow_, vz_exp_, vz_log_, vz_pow_

Trigonometric functions

vatan_, vatan2_, vcos_, vsin_, vsincos_

Note that libmvec_mt.a provides parallel versions of the vector functions that rely on multiprocessor parallelization. To use libmvec_mt.a, you must have an Sun WorkShop license and link with -xparallel.

See the libmvec(3m) manual page for more information.

Single, Double, and Long Double Precision

Most numerical functions are available in single, double, and long-double precision. Examples of calling different precision versions from different languages are shown in Table 3-4.

Table  3-4 Calling Single, Double, and Quadruple Functions


C, C++

#include <sunmath.h>
float x,y,z;
x = sinf(y);
x = fmodf(y,z);
x = max_normalf();
x = r_addran_();

#include <math.h>
double x,y,z;
x = sin(y);
x = fmod(y,z);

#include <sunmath.h>
double x,y,z;
x = max_normal();
x = d_addran_();

#include <sunmath.h>
long double x,y,z;
x = sinl(y);
x = fmodl(y,z);
x = max_normall();


REAL x,y,z
x = sin(y)
x = r_fmod(y,z)
x = r_max_normal()
x = r_addran()

REAL*8 x,y,z
x = sin(y)
x = d_fmod(y,z)
x = d_max_normal()
x = d_addran()

REAL*16 x,y,z
x = sin(y)
x = q_fmod(y,z)
x = q_max_normal()

In C, names of single-precision functions are formed by appending f to the double-precision name, and names of quadruple-precision functions are formed by adding l. Because FORTRAN calling conventions differ, libsunmath provides r_..., d_..., and q_... versions for single, double, and quadruple precision functions, respectively. FORTRAN intrinsic functions can be called by the generic name for all three precisions.

Not all functions have q_... versions. Refer to math.h and sunmath.h for names and definitions of libm and libsunmath functions.

In FORTRAN programs, remember to declare r_... functions as real, d_... functions as double precision, and q_... functions as real*16. Otherwise, type mismatches may result.

Note - The Intel version of FORTRAN supports real and double precision only; it does not support REAL*16. The Intel version of C, however, supports long double.

IEEE Support Functions

This section describes the IEEE recommended functions, the functions that supply useful values, ieee_flags, ieee_retrospective, and standard_arithmetic and nonstandard_arithmetic. Refer to Chapter 4, "Exceptions and Exception Handling," for more information on the functions ieee_flags and ieee_handler.

ieee_functions(3m) and ieee_sun(3m)

The functions described by ieee_functions(3m) and ieee_sun(3m) provide capabilities either required by the IEEE standard or recommended in its appendix. These are implemented as efficient bit mask operations.

Table  3-5 ieee_functions(3m)  



Header file


x with y's sign bit


Absolute value of x

fmod(x, y)

Remainder of x with respect to y


Base 2 unbiased exponent of x in integer format


Next representable number after x, in the direction y


Remainder of x with respect to y


Table  3-6 ieee_sun(3m)  



Header file


Classification function


Classification function


Classification function


Classification function


Classification function


Classification function


Toggle hardware


Toggle hardware



remainder(x,y) is the operation specified in IEEE Standard 754-1985. The difference between remainder(x,y) and fmod(x,y) is that the sign of the result returned by remainder(x,y) might not agree with the sign of either x or y, whereas fmod(x,y) always returns a result whose sign agrees with x. Both functions return exact results and do not generate inexact exceptions.

Table  3-7 Calling ieee_functions from FORTRAN  

IEEE Function
Single Precision
Double Precision
Quadruple Precision





















Table  3-8 Calling ieee_sun from FORTRAN

IEEE Function
Single Precision
Double Precision
Quadruple Precision





Note - You must declare d_function as double precision and q_function as REAL*16 in the FORTRAN program that uses them.


IEEE values like infinity, NaN, maximum and minimum positive floating-point numbers are provided by the functions described by the ieee_values(3m) man page. Table 3-9, Table 3-10, Table 3-11, and Table 3-12 show the decimal values and hexadecimal IEEE representations of the values provided by ieee_values(3m) functions.

Table  3-9 IEEE Values: Single Precision 

IEEE value
Decimal value
hexadecimal representation

C, C++

max normal


r = max_normalf();
r = r_max_normal()

min normal


r = min_normalf();
r = r_min_normal()

max subnormal


r = max_subnormalf();
r = r_max_subnormal()

min subnormal


r = min_subnormalf();
r = r_min_subnormal()


r = infinityf();
r = r_infinity()

quiet NaN


r = quiet_nanf(0);
r = r_quiet_nan(0)

signaling NaN


r = signaling_nanf(0);
r = r_signaling_nan(0)

Table  3-10 IEEE Values: Double Precision 

IEEE value
Decimal Value
hexadecimal representation

C, C++

max normal


7fefffff ffffffff

d = max_normal();
d = d_max_normal()

min normal


00100000 00000000

d = min_normal();
d = d_min_normal()

max subnormal


000fffff ffffffff

d = max_subnormal();
d = d_max_subnormal()

min subnormal


00000000 00000001

d = min_subnormal();
d = d_min_subnormal()


7ff00000 00000000

d = infinity();
d = d_infinity()

quiet NaN


7fffffff ffffffff

d = quiet_nan(0);
d = d_quiet_nan(0)

signaling NaN


7ff00000 00000001

d = signaling_nan(0);
d = d_signaling_nan(0)

Table  3-11 IEEE Values: Quadruple Precision (SPARC and PowerPC)

IEEE value
Decimal value
hexadecimal representation

C, C++

max normal


7ffeffff ffffffff ffffffff ffffffff

q = max_normall();
q = q_max_normal()

min normal


00010000 00000000 00000000 00000000

q = min_normall();
q = q_min_normal()

max subnormal


0000ffff ffffffff ffffffff ffffffff

q = max_subnormall();
q = q_max_subnormal()

min subnormal


00000000 00000000 00000000 00000001

q = min_subnormall();
q = q_min_subnormal()


7fff0000 00000000 00000000 00000000

q = infinityl();
q = q_infinity()

quiet NaN


7fff8000 00000000 00000000 00000000

q = quiet_nanl(0);
q = q_quiet_nan(0)

signaling NaN


7fff0000 00000000 00000000 00000001

q = signaling_nanl(0);
q = q_signaling_nan(0)

Table  3-12

IEEE value
Decimal value
hexadecimal representation (80 bits)

C, C++

max normal


7ffe ffffffff ffffffff

x = max_normall();

min positive normal


0001 80000000 00000000

x = min_normall();

max subnormal


0000 7fffffff ffffffff

x = max_subnormall();

min positive subnormal


0000 00000000 00000001

x = min_subnormall();


7fff 80000000 00000000

x = infinityl();

quiet NaN


7fff c0000000 00000000

x = quiet_nanl(0);

signaling NaN


7fff 80000000 00000001

x = signaling_nanl(0);

IEEE Values: Double Extended Precision (Intel)


ieee_flags (3m) is the recommended interface to:

The syntax for a call to ieee_flags(3m) is:

i = ieee_flags(action, mode, in, out);

The ASCII strings that are the possible values for the parameters are shown in Table 3-13:

Table  3-13 Parameter Values for ieee_flags

C or C++ Type
All Possible Values


char *

get, set, clear, clearall


char *

direction, precision, exception


char *

nearest, tozero, negative, positive, extended, double, single, inexact, division, underflow, overflow, invalid, all, common


char **

nearest, tozero, negative, positive, extended, double, single, inexact, division, underflow, overflow, invalid, all, common

The ieee_flags(3m) man page describes the parameters in complete detail.

Some of the arithmetic features that can be modified by using ieee_flags are covered in the following paragraphs. Chapter 4 contains more information on ieee_flags and IEEE exception flags.

When mode is ``direction'', the specified action applies to the current rounding direction. The possible rounding directions are: round towards nearest, round towards zero, round towards + or round towards -. The IEEE default rounding direction is round towards nearest. This means that when the mathematical result of an operation lies halfway between two representable numbers, the one nearest to the mathematical result is delivered as the result. If the mathematical result lies halfway between the two representable numbers, then round to nearest mode specifies round to the nearest even number (final bit is zero.)

Rounding towards zero is the way many pre-IEEE computers work, and corresponds mathematically to truncating the result. For example, if 2/3 is rounded to 6 decimal digits, the result is .666667 when the rounding mode is round towards nearest, but .666666 when the rounding mode is round towards zero.

When using ieee_flags to examine, clear, or set the rounding direction, possible values for the four input parameters

Table  3-14 ieee_flags Input Values for the Rounding Direction

Possible value (mode is direction)


get, set, clear, clearall


nearest, tozero, negative, positive


nearest, tozero, negative, positive

are shown in Table 3-14.

When mode is ``precision'', the specified action applies to the current rounding precision. On Intel, the possible rounding precisions are: single, double, and extended. The default rounding precision is extended; in this mode, arithmetic operations that deliver a result to a floating point register round their result to the full 64-bit precision of the extended double register format. When the rounding precision is single or double, arithmetic operations that deliver a result to a floating point register round their result to 24 or 53 significant bits, respectively. Although most programs produce results that are at least as accurate, if not more so, when extended rounding precision is used, some programs that require strict adherence to the semantics of IEEE arithmetic will not work correctly in extended rounding precision mode and must be run with the rounding precision set to single or double as appropriate.

Rounding precision cannot be set on systems using SPARC or PowerPC processors. On these systems, calling ieee_flags with mode = ``precision'' has no effect on computation.

Finally, when mode is ``exception'', the specified action applies to the current IEEE exception flags. See Chapter 4, "Exceptions and Exception Handling", for more information about using ieee_flags to examine and control the IEEE exception flags.


The libsunmath function ieee_retrospective prints information about unrequited exceptions and nonstandard IEEE modes. It reports:

The necessary information is obtained from the hardware floating-point status register.

ieee_retrospective prints information about exception flags that are raised, and exceptions for which a trap is enabled. These two distinct, if related, pieces of information should not be confused. If an exception flag is raised, then that exception occurred at some point during program execution. If a trap is enabled for an exception, then the exception may not have actually occurred (but if it had, a SIGFPE signal would have been delivered). The ieee_retrospective message is meant to alert you about exceptions that may need to be investigated (if the exception flag is raised), or to remind you that exceptions may have been handled by a signal handler (if the exception's trap is enabled.) Chapter 4, "Exceptions and Exception Handling," discusses exceptions, signals, and traps, and shows how to investigate the cause of a raised exception.

ieee_retrospective can be called anytime, but it is usually called before exit points. Programs compiled with f77 call ieee_retrospective on exit by default.

The syntax for calling this function is:

C, C++

call ieee_retrospective()
For the C function, the argument fp specifies the file to which the output will be written. The FORTRAN function always prints output on stderr.

The following example shows four of the six ieee_retrospective warning messages:

 Note: IEEE floating-point exception flags raised:  
    Inexact; Underflow; 
 Rounding direction toward zero 
 IEEE floating-point exception traps enabled: 
 See the Numerical Computation Guide, ieee_flags(3M),
    ieee_handler(3M), ieee_sun(3m) 

A warning message appears only if trapping is enabled or an exception was raised.

You can suppress ieee_retrospective messages from FORTRAN programs by one of three methods. One approach is to clear all outstanding exceptions, disable traps, and restore round-to-nearest, extended precision, and standard modes before the program exits. To do this, call ieee_flags, ieee_handler, and standard_arithmetic as follows:

character*8 out 
i = ieee_flags('clearall', '', '', out) 
call ieee_handler('clear', 'all', 0)
call standard_arithmetic() 

Note - Clearing outstanding exceptions without investigating their cause is not recommended.
Another way to avoid seeing ieee_retrospective messages is simply to redirect stderr to a file. Of course, this method should not be used if the program sends output other than ieee_retrospective messages to stderr.

The third approach is to include a dummy ieee_retrospective function in the program, for example:

subroutine ieee_retrospective


As discussed in Chapter 2, "IEEE Arithmetic," IEEE arithmetic handles underflowed results using gradual underflow. On some SPARC systems, gradual underflow is often implemented partly with software emulation of the arithmetic. If many calculations underflow, this may cause performance degradation.

To obtain some information about whether this is a case in a specific program, you can use ieee_retrospective or ieee_flags to determine if underflow exceptions occur, and check the amount of system time used by the program. If a program spends an unusually large amount of time in the operating system, and raises underflow exceptions, gradual underflow may be the cause. In this case, using non-IEEE arithmetic may speed up program execution.

The function nonstandard_arithmetic causes underflowed results to be flushed to zero on those SPARC implementations that have a mode in hardware in which flushing to zero is faster. The trade-off for speed is accuracy, because the benefits of gradual underflow are lost.

The function standard_arithmetic resets the hardware to use the default IEEE arithmetic. Both functions have no effect on implementations that provide only the default IEEE 754 style of arithmetic--SuperSPARC is such an implementation.

Implementation Features of libm and libsunmath

This section describes implementation features of libm and libsunmath:

About the Algorithms

The elementary functions in libm and libsunmath on SPARC and PowerPC systems are implemented with an ever-changing combination of table-driven and polynomial/rational approximation algorithms. Some elementary functions in libm and libsunmath on Intel systems are implemented using the elementary function kernel instructions provided in the Intel instruction set; other functions are implemented using the same table-driven or polynomial/rational approximation algorithms used on SPARC and PowerPC.

Both the table-driven and polynomial/rational approximation algorithms for the common elementary functions in libm and the common single precision elementary functions in libsunmath deliver results that are accurate to within one unit in the last place (ulp). On SPARC and PowerPC, the common quadruple precision elementary functions in libsunmath deliver results that are accurate to within one ulp, except for the expm1l and log1pl functions, which deliver results accurate to within two ulps. These error bounds have been obtained by direct analysis of the algorithms. Users may also test the accuracy of these routines using BeEF, the Berkeley Elementary Function test programs, developed by Z. Alex Liu and available from netlib in the ucbtest package.

Argument Reduction for Trigonometric Functions

Trigonometric functions for radian arguments outside the range [-/4,/4] are usually computed by reducing the argument to the indicated range by subtracting integral multiples of /2.

Because is not a machine-representable number, it must be approximated. The error in the final computed trigonometric function depends on the rounding errors in argument reduction (with an approximate as well as the rounding), and approximation errors in computing the trigonometric function of the reduced argument. Even for fairly small arguments, the relative error in the final result may be dominated by the argument reduction error, while even for fairly large arguments, the error due to argument reduction may be no worse than the other errors.

There is widespread misapprehension that trigonometric functions of all large arguments are inherently inaccurate, and all small arguments relatively accurate. This is based on the simple observation that large enough machine-representable numbers are separated by a distance greater than .

There is no inherent boundary at which computed trigonometric function values suddenly become bad, nor are the inaccurate function values useless. Provided that the argument reduction is done consistently, the fact that the argument reduction is performed with an approximation to is practically undetectable, because all essential identities and relationships are as well preserved for large arguments as for small.

libm and libsunmath trigonometric functions use an "infinitely" precise for argument reduction. The value 2/ is computed to 916 hexadecimal digits and stored in a lookup table to use during argument reduction.

The group of functions sinpi, cospi, tanpi (see Table 3-2 on page 34), scales the input argument by to avoid inaccuracies introduced by range reduction.

Data Conversion Routines

In libm and libsunmath, there is a flexible data conversion routine, convert_external, used to convert binary floating-point data between IEEE and non-IEEE formats.

Formats supported include those used by SPARC (IEEE), IBM PC, VAX®, IBM S/370, and Cray®.

Refer to the man page on convert_external(3m) for an example of taking data generated on a Cray, and using the function convert_external to convert the data into the IEEE format expected on SPARC systems.

Random Number Facilities

There are two facilities for generating uniform pseudo-random numbers, addrans(3m) and lcrans(3m). addrans is a table-driven additive random number generator, and lcrans is a linear congruential random number generator.

In addition, shufrans(3m) shuffles a set of pseudo-random numbers to provide even more randomness for applications that need it.

The functions in addrans (3m) are generally more efficient but their theory is not as refined as those in lcrans (3m). Refer to the article "Random Number Generators: Good Ones Are Hard To Find", by S. Park and K. Miller, Communications of the ACM, October 1988, for a discussion of the theoretical properties of these algorithms. Additive random number generators are discussed in Volume 2 of Knuth's The Art of Computer Programming.

The random number generators are available in four versions, for the four data types: signed integer, unsigned integer, single precision floating-point, and double precision floating-point. These functions can be used in two ways: to generate random numbers one at a time (one per function call), or to generate arrays of random numbers (one array per function call.)

When using addrans and lcrans to provide variates one at a time, they are constrained to the intervals shown in Table 3-15.

Table  3-15 Interval for Single-Value Random Number Generators


Lower Bound
Upper Bound

Signed integer



Unsigned integer



Single precision
floating point



Double precision
floating point



That is, these are the upper and lower bounds over which the data is uniformly distributed. These bounds cannot be altered.

However, if addrans or lcrans are used to generate arrays of data, the upper and lower bounds of the interval over which the data is uniformly distributed can be controlled.

Appendix A, "Examples," on page 91, shows an example that uses addrans to generate an array of 1000 random numbers, uniformly distributed over the interval, [0, 3.7e-9].

libc Support Functions

This section describes the base conversion feature for libc support (applicable only to the Solaris operating environment).

Base Conversion

Base conversion is used by I/O routines, like printf and scanf in C, and read, write, and print in FORTRAN. For these functions you need conversions between numbers representations in bases 2 and 10:

The algorithms for base conversion are now table driven as well. These algorithms yield correctly-rounded base conversion routines. In addition to improved accuracy, table-driven algorithms reduce the worst-case times for correctly-rounded base conversion.

Base conversion is not too difficult for integer formats but gets more complicated when floating point is involved; as with other floating-point operations, rounding errors occur.

For instance, on conversion from internal floating point to an external decimal representation, the mathematical problem to be solved in default rounding mode is as follows:

In the case of printf %f or FORTRAN F format with n digits after the decimal point:

Given integers i and e, find an integer J such that

i × 2**e × 10**n - J 0.5

In the case of printf %e or FORTRAN E format with n significant digits:

Given integers i and e, find integers J and E such that

i × 2**e - J × 10**E 0.5 × 10**E


10**(n-1) J 10**n - 1

In both cases, distinguish exact conversions, and round halfway cases to nearest even.

Correct rounding is required by the IEEE Standard 754. The standard requires correct rounding for typical numbers whose magnitudes range from 10 -44 to 10+44, but permits slightly incorrect rounding for larger exponents--the Sun table-driven algorithm for rounding is accurate through the entire range. (See section 5.6 of IEEE Standard 754.)

See Appendix G, "References," on page 237, for references on base conversion. Particularly good references are Coonen's thesis and Sterbenz's book.

Previous Next Contents Index Doc Set Home