This section describes implementation features of libm and libsunmath:
Argument reduction using infinitely precise ??, and trigonometric functions scaled in ??
Data conversion routines for converting floating-point data between IEEE and non-IEEE formats
Random number generators
The elementary functions in libm and libsunmath on SPARC-based systems are implemented with table-driven and polynomial/rational approximation algorithms. These algorithms are subject to change between releases for better performance or accuracy. Some elementary functions in libm and libsunmath on x86-based systems are implemented using the elementary function kernel instructions provided in the x86 instruction set; other functions are implemented using the same table-driven or polynomial/rational approximation algorithms used on SPARC-based systems.
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-based systems, 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. (The common functions include the exponential, logarithm, and power functions and circular trigonometric functions of radian arguments. Other functions, such as the hyperbolic trig functions and higher transcendental functions, are less accurate.) These error bounds have been obtained by direct analysis of the algorithms. Users can also test the accuracy of these routines using BeEF, the Berkeley Elementary Function test programs, available from netlib in the ucbtest package http://www.netlib.org/fp/ucbtest.tgzhttp://www.netlib.org/fp/ucbtest.tgz.
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 might be dominated by the argument reduction error, while even for fairly large arguments, the error due to argument reduction might 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 be 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, and tanpi (see Table 16) scales the input argument by ?? to avoid inaccuracies introduced by range reduction.
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-based systems.
There are three facilities for generating uniform pseudo-random numbers in 32-bit integer, single precision floating-point, and double precision floating-point formats:
The functions described in the addrans(3m) manual page are based on a family of table-driven additive random number generators.
The functions described in the lcrans(3m) manual page are based on a linear congruential random number generator.
The functions described in the mwcrans(3m) manual page are based on multiply-with-carry random number generators. These functions also include generators that supply uniform pseudo-random numbers in 64-bit integer formats.
In addition, the functions described on the shufrans(3m) manual page can be used in conjunction with any of these generators to shuffle an array of pseudo-random numbers, thereby providing even more randomness for applications that need it. Note that there is no facility for shuffling arrays of 64-bit integers.
Each of the random number facilities includes routines that generate one random number at a time, i.e., one per function call, as well as routines that generate an array of random numbers in a single call. The functions that generate one random number at a time deliver numbers that lie in the ranges shown in Table 30.
|
The functions that generate an entire array of random numbers in a single call allow the user to specify the interval in which the generated numbers will lie. Examples gives several examples that show how to generate arrays of random numbers uniformly distributed over different intervals.
Note that the addrans and mwcrans generators are generally more efficient than the lcrans generators, but their theory is not as refined. “Random Number Generators: Good Ones Are Hard To Find”, by S. Park and K. Miller, Communications of the ACM, October 1988, discusses the theoretical properties of linear congruential algorithms. Additive random number generators are discussed in Volume 2 of Knuth's The Art of Computer Programming.