MATH(3M) RISC/os Reference Manual MATH(3M)
NAME
math - introduction to mathematical library functions
DESCRIPTION
These functions constitute the C math library libm. There
are two versions of the math library libm.a and libm43.a.
The first, libm.a, contains routines written in MIPS assem-
bly language and tuned for best performance and includes
many routines for the float data type. The routines in
there are based on the algorithms of Cody and Waite or those
in the 4.3 BSD release, whichever provides the best perfor-
mance with acceptable error bounds. Those routines with
Cody and Waite implementations are marked with a `*' in the
list of functions below.
The second version of the math library, libm43.a, contains
routines all based on the original codes in the 4.3 BSD
release. The difference between the two version's error
bounds is typically around 1 unit in the last place, whereas
the performance difference may be a factor of two or more.
The link editor searches this library under the "-lm" (or
"-lm43") option. Declarations for these functions may be
obtained from the include file <math.h>. The Fortran math
library is described in intro(3F).
LIST OF FUNCTIONS
The cycle counts of all functions are approximate; cycle
counts often depend on the value of argument. The error
bound sometimes applies only to the primary range.
Error Bound (ULPs) Cycles
Name Appears on Page Description libm.a libm43.a libm.a libm43.a
acos sin(3M) inverse trigonometric function 3 3? ?
acosh asinh(3M) inverse hyperbolic function 3 3? ?
asin sin(3M) inverse trigonometric function 3 3? ?
asinh asinh(3M) inverse hyperbolic function 3 3? ?
atan sin(3M) inverse trigonometric function 1 152 260
atanh asinh(3M) inverse hyperbolic function 3 3? ?
atan2 sin(3M) inverse trigonometric function 2 2? ?
cabs hypot(3M) complex absolute value 1 1? ?
cbrt sqrt(3M) cube root 1 1? ?
ceil floor(3M) integer no less than 0 0? ?
copysign ieee(3M) copy sign bit 0 0? ?
cos sin(3M) trigonometric function 2 128 243
cosh sinh(3M) hyperbolic function ? 1
342 294
drem ieee(3M) remainder 0 0? ?
erf erf(3M) error function ? ?? ?
erfc erf(3M) complementary error function ? ?? ?
exp exp(3M) exponential 2 101 230
Printed 11/19/92 Page 1
MATH(3M) RISC/os Reference Manual MATH(3M)
expm1 exp(3M) exp(x)-1 1 1
281 281
fabs floor(3M) absolute value 0 0? ?
facos sin(3M) inverse trigonometric function
fatan sin(3M) inverse trigonometric function 3 64
fatan2 sin(3M) inverse trigonometric function 3 64
fcos sin(3M) trigonometric function 1 87
fcosh sinh(3M) hyperbolic function ? 105
fexp exp(3M) exponential 1 79
finite ieee(3M) floating point arithmetic
flog exp(3M) natural logarithm 1 100
floor floor(3M) integer no greater than 0 0? ?
fsin sin(3M) trigonometric function 1 68
fsinh sinh(3M) hyperbolic function ? 44
fsqrt sqrt(3M) square root 1 95
ftan sin(3M) trigonometric function ? 61
ftanh sinh(3M) hyperbolic function ? 116
hypot hypot(3M) Euclidean distance 1 1? ?
j0 j0(3M) bessel function ? ?? ?
j1 j0(3M) bessel function ? ?? ?
jn j0(3M) bessel function ? ?? ?
lgamma lgamma(3M) log gamma function ? ?? ?
log exp(3M) natural logarithm 2 119 217
logb ieee(3M) exponent extraction 0 0? ?
log10 exp(3M) logarithm to base 10 3 3? ?
log1p exp(3M) log(1+x) 1 1
269 269
pow exp(3M) exponential xy 60-500 60-500 ??
rint floor(3M) round to nearest integer 0 0? ?
scalb ieee(3M) exponent adjustment 0 0? ?
sin sin(3M) trigonometric function 2 101 222
sinh sinh(3M) hyperbolic function ? 379 292
sqrt sqrt(3M) square root 1 133 133
tan sin(3M) trigonometric function ? 392 287
tanh sinh(3M) hyperbolic function ? 1
356 293
y0 j0(3M) bessel function ? ?? ?
y1 j0(3M) bessel function ? ?? ?
yn j0(3M) bessel function ? ?? ?
NOTES
In 4.3 BSD, distributed from the University of California in
late 1985, most of the foregoing functions come in two ver-
sions, one for the double-precision "D" format in the DEC
VAX-11 family of computers, another for double-precision
arithmetic conforming to the IEEE Standard 754 for Binary
Floating-Point Arithmetic. The two versions behave very
similarly, as should be expected from programs more accurate
and robust than was the norm when UNIX was born. For
instance, the programs are accurate to within the numbers of
s tabulated above; an is one Unit in the Last Place. And
the programs have been cured of anomalies that afflicted the
older math library libm in which incidents like the follow-
ing had been reported:
sqrt(-1.0) = 0.0 and log(-1.0) = -1.7e38.
Page 2 Printed 11/19/92
MATH(3M) RISC/os Reference Manual MATH(3M)
cos(1.0e-11) > cos(0.0) > 1.0.
pow(x,1.0) != x when x = 2.0, 3.0, 4.0, ..., 9.0.
pow(-1.0,1.0e10) trapped on Integer Overflow.
sqrt(1.0e30) and sqrt(1.0e-30) were very slow.
MIPS machines conform to the IEEE Standard 754 for Binary
Floating-Point Arithmetic, to which only the notes for IEEE
floating-point apply and are included here.
IEEE STANDARD 754 Floating-Point Arithmetic:
This standard is on its way to becoming more widely adopted
than any other design for computer arithmetic.
The main virtue of 4.3 BSD's libm codes is that they are
intended for the public domain; they may be copied freely
provided their provenance is always acknowledged, and pro-
vided users assist the authors in their researches by
reporting experience with the codes. Therefore no user of
UNIX on a machine that conforms to IEEE 754 need use any-
thing worse than the new libm.
Properties of IEEE 754 Double-Precision:
Wordsize: 64 bits, 8 bytes. Radix: Binary.
Precision: 53 sig. bits, roughly like 16 sig.
decimals.
If x and x' are consecutive positive
Double-Precision numbers (they differ by 1 ), then
1.1e-16 < 0.553 < (x'-x)/x < 0.552 < 2.3e-16.
Range: Overflow threshold = 2.01024 = 1.8e308
Underflow threshold = 0.51022 = 2.2e-308
Overflow goes by default to a signed Infinity.
Underflow is Gradual, rounding to the nearest
integer multiple of 0.51074 = 4.9e-324.
Zero is represented ambiguously as +0 or -0.
Its sign transforms correctly through multiplica-
tion or division, and is preserved by addition of
zeros with like signs; but x-x yields +0 for every
finite x. The only operations that reveal zero's
sign are division by zero and copysign(x,+0). In
particular, comparison (x > y, x > y, etc.) can-
not be affected by the sign of zero; but if finite
x = y then Infinity = 1/(x-y) != -1/(y-x) =
-Infinity.
Infinity is signed.
it persists when added to itself or to any finite
number. Its sign transforms correctly through
multiplication and division, and
(finite)/+Infinity = +0 (nonzero)/0 = +Infinity.
But Infinity-Infinity, Infinity and
Infinity/Infinity are, like 0/0 and sqrt(-3),
invalid operations that produce . ...
Reserved operands:
Printed 11/19/92 Page 3
MATH(3M) RISC/os Reference Manual MATH(3M)
there are 253-2 of them, all called (Not a
Number). Some, called Signaling s, trap any
floating-point operation performed upon them; they
could be used to mark missing or uninitialized
values, or nonexistent elements of arrays. The
rest are Quiet s; they are the default results of
Invalid Operations, and propagate through subse-
quent arithmetic operations. If x != x then x is
; every other predicate (x > y, x = y, x < y, ...)
is FALSE if is involved.
NOTE: Trichotomy is violated by .
Besides being FALSE, predicates that entail
ordered comparison, rather than mere
(in)equality, signal Invalid Operation when
is involved.
Rounding:
Every algebraic operation (+, -, /, sqrt) is
rounded by default to within half an , and when
the rounding error is exactly half an then the
rounded value's least significant bit is zero.
This kind of rounding is usually the best kind,
sometimes provably so; for instance, for every x =
1.0, 2.0, 3.0, 4.0, ..., 2.052, we find (x/3.0)0
== x and (x/10.0).0 == x and ... despite that
both the quotients and the products have been
rounded. Only rounding like IEEE 754 can do that.
But no single kind of rounding can be proved best
for every circumstance, so IEEE 754 provides
rounding towards zero or towards +Infinity or
towards -Infinity at the programmer's option. And
the same kinds of rounding are specified for
Binary-Decimal Conversions, at least for magni-
tudes between roughly 1.0e-10 and 1.0e37.
Exceptions:
IEEE 754 recognizes five kinds of floating-point
exceptions, listed below in declining order of
probable importance.
Exception Default Result
__________________________________________
Invalid Operation , or FALSE
Overflow +Infinity
Divide by Zero +Infinity
Underflow Gradual Underflow
Inexact Rounded value
NOTE: An Exception is not an Error unless handled
badly. What makes a class of exceptions excep-
tional is that no single default response can be
satisfactory in every instance. On the other
hand, if a default response will serve most
instances satisfactorily, the unsatisfactory
Page 4 Printed 11/19/92
MATH(3M) RISC/os Reference Manual MATH(3M)
instances cannot justify aborting computation
every time the exception occurs.
For each kind of floating-point exception, IEEE 754
provides a Flag that is raised each time its exception
is signaled, and stays raised until the program resets
it. Programs may also test, save and restore a flag.
Thus, IEEE 754 provides three ways by which programs
may cope with exceptions for which the default result
might be unsatisfactory:
1) Test for a condition that might cause an exception
later, and branch to avoid the exception.
2) Test a flag to see whether an exception has
occurred since the program last reset its flag.
3) Test a result to see whether it is a value that
only an exception could have produced.
CAUTION: The only reliable ways to discover whether
Underflow has occurred are to test whether products
or quotients lie closer to zero than the underflow
threshold, or to test the Underflow flag. (Sums
and differences cannot underflow in IEEE 754; if x
!= y then x-y is correct to full precision and cer-
tainly nonzero regardless of how tiny it may be.)
Products and quotients that underflow gradually can
lose accuracy gradually without vanishing, so com-
paring them with zero (as one might on a VAX) will
not reveal the loss. Fortunately, if a gradually
underflowed value is destined to be added to some-
thing bigger than the underflow threshold, as is
almost always the case, digits lost to gradual
underflow will not be missed because they would
have been rounded off anyway. So gradual under-
flows are usually provably ignorable. The same
cannot be said of underflows flushed to 0.
At the option of an implementor conforming to IEEE 754,
other ways to cope with exceptions may be provided:
4) ABORT. This mechanism classifies an exception in
advance as an incident to be handled by means trad-
itionally associated with error-handling statements
like "ON ERROR GO TO ...". Different languages
offer different forms of this statement, but most
share the following characteristics:
- No means is provided to substitute a value for the
offending operation's result and resume computation
from what may be the middle of an expression. An
exceptional result is abandoned.
Printed 11/19/92 Page 5
MATH(3M) RISC/os Reference Manual MATH(3M)
- In a subprogram that lacks an error-handling state-
ment, an exception causes the subprogram to abort
within whatever program called it, and so on back
up the chain of calling subprograms until an
error-handling statement is encountered or the
whole task is aborted and memory is dumped.
5) STOP. This mechanism, requiring an interactive
debugging environment, is more for the programmer
than the program. It classifies an exception in
advance as a symptom of a programmer's error; the
exception suspends execution as near as it can to
the offending operation so that the programmer can
look around to see how it happened. Quite often
the first several exceptions turn out to be quite
unexceptionable, so the programmer ought ideally to
be able to resume execution after each one as if
execution had not been stopped.
6) ... Other ways lie beyond the scope of this docu-
ment.
The crucial problem for exception handling is the problem of
Scope, and the problem's solution is understood, but not
enough manpower was available to implement it fully in time
to be distributed in 4.3 BSD's libm. Ideally, each elemen-
tary function should act as if it were indivisible, or
atomic, in the sense that ...
i) No exception should be signaled that is not deserved
by the data supplied to that function.
ii) Any exception signaled should be identified with that
function rather than with one of its subroutines.
iii) The internal behavior of an atomic function should not
be disrupted when a calling program changes from one
to another of the five or so ways of handling excep-
tions listed above, although the definition of the
function may be correlated intentionally with excep-
tion handling.
Ideally, every programmer should be able conveniently to
turn a debugged subprogram into one that appears atomic to
its users. But simulating all three characteristics of an
atomic function is still a tedious affair, entailing hosts
of tests and saves-restores; work is under way to ameliorate
the inconvenience.
Meanwhile, the functions in libm are only approximately
atomic. They signal no inappropriate exception except pos-
sibly ...
Page 6 Printed 11/19/92
MATH(3M) RISC/os Reference Manual MATH(3M)
Over/Underflow
when a result, if properly computed, might have
lain barely within range, and
Inexact in cabs, cbrt, hypot, log10 and pow
when it happens to be exact, thanks to fortuitous
cancellation of errors.
Otherwise, ...
Invalid Operation is signaled only when
any result but would probably be misleading.
Overflow is signaled only when
the exact result would be finite but beyond the
overflow threshold.
Divide-by-Zero is signaled only when
a function takes exactly infinite values at finite
operands.
Underflow is signaled only when
the exact result would be nonzero but tinier than
the underflow threshold.
Inexact is signaled only when
greater range or precision would be needed to
represent the exact result.
Exceptions on MIPS machines:
The exception enables and the flags that are raised
when an exception occurs (as well as the rounding mode)
are in the floating-point control and status register.
This register can be read or written by the routines
described on the man page fpc(3). This register's lay-
out is described in the file <sys/fpu.h>.
A full implementation of IEEE 754 ``user trap
handlers'' is under development at MIPS computer sys-
tems. At which time all functions in libm will appear
atomic and the full functionality of user trap handlers
will be supported in thoses language without other
floating-point error handling intrinsics (i.e. ADA,
Pl/1, etc). For a description of these trap handlers
see section 8 of the IEEE 754 standard.
What is currently available is only the raw interface
which was only intended to be used by the code to
implement IEEE user trap handlers. IEEE floating-point
exceptions are enabled by setting the enable bit for
that exception in the floating-point control and status
register. If an exception then occurs the UNIX signal
SIGFPE is sent to the process. It is up to the signal
handler to determine the instruction that caused the
exception and to take the action specified by the user.
The instruction that caused the exception is in one of
two places. If the floating-point board is used (the
floating-point implementation revision register indi-
cates this in it's implementation field) then the
instruction that caused the exception is in the
Printed 11/19/92 Page 7
MATH(3M) RISC/os Reference Manual MATH(3M)
floating-point exception instruction register. In all
other implementations the instruction that caused the
exception is at the address of the program counter as
modified by the branch delay bit in the cause register.
Both the program counter and cause register are in the
sigcontext structure passed to the signal handler (see
signal(3C)). If the program is to be continued past
the instruction that caused the exception the program
counter in the signal context must be advanced. If the
instruction is in a branch delay slot then the branch
must be emulated to determine if the branch is taken
and then the resulting program counter can be calcu-
lated (see emulate_branch(3) and the NOTES (MIPS) sec-
tion in signal(3C)).
BUGS
When signals are appropriate, they are emitted by certain
operations within the codes, so a subroutine-trace may be
needed to identify the function with its signal in case
method 5) above is in use. And the codes all take the IEEE
754 defaults for granted; this means that a decision to trap
all divisions by zero could disrupt a code that would other-
wise get correct results despite division by zero.
SEE ALSO
emulate_branch(3), fpc(3), signal(3C).
R2010 Floating Point Coprocessor Architecture
R2360 Floating Point Board Product Description
An explanation of IEEE 754 and its proposed extension p854
was published in the IEEE magazine MICRO in August 1984
under the title "A Proposed Radix- and
Word-length-independent Standard for Floating-point Arith-
metic" by W. J. Cody et al. Articles in the IEEE magazine
COMPUTER vol. 14 no. 3 (Mar. 1981), and in the ACM SIGNUM
Newsletter Special Issue of Oct. 1979, may be helpful
although they pertain to superseded drafts of the standard.
AUTHOR
W. Kahan, with the help of Z-S. Alex Liu, Stuart I.
McDonald, Dr. Kwok-Choi Ng, Peter Tang.
Page 8 Printed 11/19/92