1 .\" $NetBSD: math.3,v 1.26 2012/11/10 15:59:58 njoly Exp $
3 .\" Copyright (c) 1985 Regents of the University of California.
4 .\" All rights reserved.
6 .\" Redistribution and use in source and binary forms, with or without
7 .\" modification, are permitted provided that the following conditions
9 .\" 1. Redistributions of source code must retain the above copyright
10 .\" notice, this list of conditions and the following disclaimer.
11 .\" 2. Redistributions in binary form must reproduce the above copyright
12 .\" notice, this list of conditions and the following disclaimer in the
13 .\" documentation and/or other materials provided with the distribution.
14 .\" 3. Neither the name of the University nor the names of its contributors
15 .\" may be used to endorse or promote products derived from this software
16 .\" without specific prior written permission.
18 .\" THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
19 .\" ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 .\" IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21 .\" ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
22 .\" FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 .\" DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
24 .\" OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
25 .\" HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
26 .\" LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
27 .\" OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
30 .\" from: @(#)math.3 6.10 (Berkeley) 5/6/91
37 .Nd introduction to mathematical library functions
43 These functions constitute the C
45 Declarations for these functions may be obtained from the include file
47 .\" The Fortran math library is described in ``man 3f intro''.
49 .Bl -column "copysignX" "gammaX3XX" "inverse trigonometric funcX"
50 .It Sy Name Ta Sy Man page Ta Sy Description Ta Sy Error Bound Dv ( ULP Ns No s)
51 .It acos Ta Xr acos 3 Ta inverse trigonometric function Ta 3
52 .It acosh Ta Xr acosh 3 Ta inverse hyperbolic function Ta 3
53 .It asin Ta Xr asin 3 Ta inverse trigonometric function Ta 3
54 .It asinh Ta Xr asinh 3 Ta inverse hyperbolic function Ta 3
55 .It atan Ta Xr atan 3 Ta inverse trigonometric function Ta 1
56 .It atanh Ta Xr atanh 3 Ta inverse hyperbolic function Ta 3
57 .It atan2 Ta Xr atan2 3 Ta inverse trigonometric function Ta 2
58 .It cbrt Ta Xr sqrt 3 Ta cube root Ta 1
59 .It ceil Ta Xr ceil 3 Ta integer no less than Ta 0
60 .It copysign Ta Xr copysign 3 Ta copy sign bit Ta 0
61 .It cos Ta Xr cos 3 Ta trigonometric function Ta 1
62 .It cosh Ta Xr cosh 3 Ta hyperbolic function Ta 3
63 .It erf Ta Xr erf 3 Ta error function Ta ???
64 .It erfc Ta Xr erf 3 Ta complementary error function Ta ???
65 .It exp Ta Xr exp 3 Ta exponential Ta 1
66 .It expm1 Ta Xr exp 3 Ta exp(x)\-1 Ta 1
67 .It fabs Ta Xr fabs 3 Ta absolute value Ta 0
68 .It finite Ta Xr finite 3 Ta test for finity Ta 0
69 .It floor Ta Xr floor 3 Ta integer no greater than Ta 0
70 .It fmod Ta Xr fmod 3 Ta remainder Ta ???
71 .It hypot Ta Xr hypot 3 Ta Euclidean distance Ta 1
72 .It ilogb Ta Xr ilogb 3 Ta exponent extraction Ta 0
73 .It isinf Ta Xr isinf 3 Ta test for infinity Ta 0
74 .It isnan Ta Xr isnan 3 Ta test for not-a-number Ta 0
75 .It j0 Ta Xr j0 3 Ta Bessel function Ta ???
76 .It j1 Ta Xr j0 3 Ta Bessel function Ta ???
77 .It jn Ta Xr j0 3 Ta Bessel function Ta ???
78 .It lgamma Ta Xr lgamma 3 Ta log gamma function Ta ???
79 .It log Ta Xr log 3 Ta natural logarithm Ta 1
80 .It log10 Ta Xr log 3 Ta logarithm to base 10 Ta 3
81 .It log1p Ta Xr log 3 Ta log(1+x) Ta 1
82 .It nan Ta Xr nan 3 Ta return quiet \*(Na Ta 0
83 .It nextafter Ta Xr nextafter 3 Ta next representable number Ta 0
84 .It pow Ta Xr pow 3 Ta exponential x**y Ta 60\-500
85 .It remainder Ta Xr remainder 3 Ta remainder Ta 0
86 .It rint Ta Xr rint 3 Ta round to nearest integer Ta 0
87 .It scalbn Ta Xr scalbn 3 Ta exponent adjustment Ta 0
88 .It sin Ta Xr sin 3 Ta trigonometric function Ta 1
89 .It sinh Ta Xr sinh 3 Ta hyperbolic function Ta 3
90 .It sqrt Ta Xr sqrt 3 Ta square root Ta 1
91 .It tan Ta Xr tan 3 Ta trigonometric function Ta 3
92 .It tanh Ta Xr tanh 3 Ta hyperbolic function Ta 3
93 .It trunc Ta Xr trunc 3 Ta nearest integral value Ta 3
94 .It y0 Ta Xr j0 3 Ta Bessel function Ta ???
95 .It y1 Ta Xr j0 3 Ta Bessel function Ta ???
96 .It yn Ta Xr j0 3 Ta Bessel function Ta ???
98 .Ss List of Defined Values
99 .Bl -column "M_2_SQRTPIXX" "1.12837916709551257390XX" "2/sqrt(pi)XXX"
100 .It Sy Name Ta Sy Value Ta Sy Description
101 .It M_E 2.7182818284590452354 e
102 .It M_LOG2E 1.4426950408889634074 log 2e
103 .It M_LOG10E 0.43429448190325182765 log 10e
104 .It M_LN2 0.69314718055994530942 log e2
105 .It M_LN10 2.30258509299404568402 log e10
106 .It M_PI 3.14159265358979323846 pi
107 .It M_PI_2 1.57079632679489661923 pi/2
108 .It M_PI_4 0.78539816339744830962 pi/4
109 .It M_1_PI 0.31830988618379067154 1/pi
110 .It M_2_PI 0.63661977236758134308 2/pi
111 .It M_2_SQRTPI 1.12837916709551257390 2/sqrt(pi)
112 .It M_SQRT2 1.41421356237309504880 sqrt(2)
113 .It M_SQRT1_2 0.70710678118654752440 1/sqrt(2)
116 In 4.3 BSD, distributed from the University of California
117 in late 1985, most of the foregoing functions come in two
118 versions, one for the double\-precision "D" format in the
119 DEC VAX\-11 family of computers, another for double\-precision
120 arithmetic conforming to the IEEE Standard 754 for Binary
121 Floating\-Point Arithmetic.
122 The two versions behave very
123 similarly, as should be expected from programs more accurate
124 and robust than was the norm when UNIX was born.
125 For instance, the programs are accurate to within the numbers
130 is one Unit in the Last Place.
131 And the programs have been cured of anomalies that
132 afflicted the older math library
133 in which incidents like
134 the following had been reported:
135 .Bd -literal -offset indent
136 sqrt(\-1.0) = 0.0 and log(\-1.0) = \-1.7e38.
137 cos(1.0e\-11) \*[Gt] cos(0.0) \*[Gt] 1.0.
138 pow(x,1.0) \(!= x when x = 2.0, 3.0, 4.0, ..., 9.0.
139 pow(\-1.0,1.0e10) trapped on Integer Overflow.
140 sqrt(1.0e30) and sqrt(1.0e\-30) were very slow.
142 However the two versions do differ in ways that have to be
143 explained, to which end the following notes are provided.
144 .Ss DEC VAX\-11 D_floating\-point
145 This is the format for which the original math library
146 was developed, and to which this manual is still principally dedicated.
149 double\-precision format for the PDP\-11
150 and the earlier VAX\-11 machines; VAX\-11s after 1983 were
151 provided with an optional "G" format closer to the IEEE
152 double\-precision format.
153 The earlier DEC MicroVAXs have no D format, only G double\-precision.
157 Properties of D_floating\-point:
158 .Bl -hang -offset indent
164 56 significant bits, roughly like 17 significant decimals.
165 If x and x' are consecutive positive D_floating\-point
166 numbers (they differ by 1
169 .Dl 1.3e\-17 \*[Lt] 0.5**56 \*[Lt] (x'\-x)/x \*[Le] 0.5**55 \*[Lt] 2.8e\-17.
171 .Bl -column "Underflow thresholdX" "2.0**127X"
172 .It Overflow threshold = 2.0**127 = 1.7e38.
173 .It Underflow threshold = 0.5**128 = 2.9e\-39.
175 .Em NOTE: THIS RANGE IS COMPARATIVELY NARROW.
177 Overflow customarily stops computation.
178 Underflow is customarily flushed quietly to zero.
180 It is possible to have x
182 y and yet x\-y = 0 because of underflow.
183 Similarly x \*[Gt] y \*[Gt] 0 cannot prevent either x\(**y = 0
184 or y/x = 0 from happening without warning.
185 .It Zero is represented ambiguously :
186 Although 2**55 different representations of zero are accepted by
187 the hardware, only the obvious representation is ever produced.
188 There is no \-0 on a VAX.
189 .It \*(If is not part of the VAX architecture .
190 .It Reserved operands :
191 of the 2**55 that the hardware
192 recognizes, only one of them is ever produced.
193 Any floating\-point operation upon a reserved
194 operand, even a MOVF or MOVD, customarily stops
195 computation, so they are not much used.
197 Divisions by zero and operations that
198 overflow are invalid operations that customarily
199 stop computation or, in earlier machines, produce
200 reserved operands that will stop computation.
202 Every rational operation (+, \-, \(**, /) on a
203 VAX (but not necessarily on a PDP\-11), if not an
204 over/underflow nor division by zero, is rounded to
207 and when the rounding error is
210 then rounding is away from 0.
213 Except for its narrow range, D_floating\-point is one of the
214 better computer arithmetics designed in the 1960's.
215 Its properties are reflected fairly faithfully in the elementary
216 functions for a VAX distributed in 4.3 BSD.
217 They over/underflow only if their results have to lie out of range
218 or very nearly so, and then they behave much as any rational
219 arithmetic operation that over/underflowed would behave.
220 Similarly, expressions like log(0) and atanh(1) behave
221 like 1/0; and sqrt(\-3) and acos(3) behave like 0/0;
222 they all produce reserved operands and/or stop computation!
223 The situation is described in more detail in manual pages.
225 .Em This response seems excessively punitive, so it is destined
226 .Em to be replaced at some time in the foreseeable future by a
227 .Em more flexible but still uniform scheme being developed to
228 .Em handle all floating\-point arithmetic exceptions neatly.
230 How do the functions in 4.3 BSD's new math library for UNIX
231 compare with their counterparts in DEC's VAX/VMS library?
232 Some of the VMS functions are a little faster, some are
233 a little more accurate, some are more puritanical about
234 exceptions (like pow(0.0,0.0) and atan2(0.0,0.0)),
235 and most occupy much more memory than their counterparts in
237 The VMS codes interpolate in large table to achieve
238 speed and accuracy; the libm codes use tricky formulas
239 compact enough that all of them may some day fit into a ROM.
241 More important, DEC regards the VMS codes as proprietary
242 and guards them zealously against unauthorized use.
243 But the libm codes in 4.3 BSD are intended for the public domain;
244 they may be copied freely provided their provenance is always
245 acknowledged, and provided users assist the authors in their
246 researches by reporting experience with the codes.
247 Therefore no user of UNIX on a machine whose arithmetic resembles
248 VAX D_floating\-point need use anything worse than the new libm.
249 .Ss IEEE STANDARD 754 Floating\-Point Arithmetic
250 This standard is on its way to becoming more widely adopted
251 than any other design for computer arithmetic.
252 VLSI chips that conform to some version of that standard have been
253 produced by a host of manufacturers, among them ...
254 .Bl -column "Intel i8070, i80287XX"
255 .It Intel i8087, i80287 National Semiconductor 32081
256 .It 68881 Weitek WTL-1032, ... , -1165
257 .It Zilog Z8070 Western Electric (AT\*[Am]T) WE32106.
259 Other implementations range from software, done thoroughly
260 in the Apple Macintosh, through VLSI in the Hewlett\-Packard
261 9000 series, to the ELXSI 6400 running ECL at 3 Megaflops.
262 Several other companies have adopted the formats
263 of IEEE 754 without, alas, adhering to the standard's way
264 of handling rounding and exceptions like over/underflow.
265 The DEC VAX G_floating\-point format is very similar to the IEEE
266 754 Double format, so similar that the C programs for the
267 IEEE versions of most of the elementary functions listed
268 above could easily be converted to run on a MicroVAX, though
269 nobody has volunteered to do that yet.
271 The codes in 4.3 BSD's libm for machines that conform to
272 IEEE 754 are intended primarily for the National Semiconductor 32081
274 To use these codes with the Intel or Zilog
275 chips, or with the Apple Macintosh or ELXSI 6400, is to
276 forego the use of better codes provided (perhaps freely) by
277 those companies and designed by some of the authors of the
290 the Motorola 68881 has all the functions in libm on chip,
291 and faster and more accurate;
292 it, Apple, the i8087, Z8070 and WE32106 all use 64 significant bits.
293 The main virtue of 4.3 BSD's
294 libm codes is that they are intended for the public domain;
295 they may be copied freely provided their provenance is always
296 acknowledged, and provided users assist the authors in their
297 researches by reporting experience with the codes.
298 Therefore no user of UNIX on a machine that conforms to
299 IEEE 754 need use anything worse than the new libm.
301 Properties of IEEE 754 Double\-Precision:
302 .Bl -hang -offset indent
308 53 significant bits, roughly like 16 significant decimals.
309 If x and x' are consecutive positive Double\-Precision
310 numbers (they differ by 1
313 .Dl 1.1e\-16 \*[Lt] 0.5**53 \*[Lt] (x'\-x)/x \*[Le] 0.5**52 \*[Lt] 2.3e\-16.
315 .Bl -column "Underflow thresholdX" "2.0**1024X"
316 .It Overflow threshold = 2.0**1024 = 1.8e308
317 .It Underflow threshold = 0.5**1022 = 2.2e\-308
319 Overflow goes by default to a signed \*(If.
322 rounding to the nearest
323 integer multiple of 0.5**1074 = 4.9e\-324.
324 .It Zero is represented ambiguously as +0 or \-0:
325 Its sign transforms correctly through multiplication or
326 division, and is preserved by addition of zeros
327 with like signs; but x\-x yields +0 for every
329 The only operations that reveal zero's
330 sign are division by zero and copysign(x,\(+-0).
331 In particular, comparison (x \*[Gt] y, x \*[Ge] y, etc.)
332 cannot be affected by the sign of zero; but if
333 finite x = y then \*(If
338 .It \*(If is signed :
339 it persists when added to itself
340 or to any finite number.
342 correctly through multiplication and division, and
343 \*(If (finite)/\(+- \0=\0\(+-0
347 \(if\-\(if, \(if\(**0 and \(if/\(if
348 are, like 0/0 and sqrt(\-3),
349 invalid operations that produce \*(Na.
350 .It Reserved operands :
351 there are 2**53\-2 of them, all
352 called \*(Na (Not A Number).
353 Some, called Signaling \*[Na]s, trap any floating\-point operation
354 performed upon them; they are used to mark missing
355 or uninitialized values, or nonexistent elements of arrays.
356 The rest are Quiet \*[Na]s; they are
357 the default results of Invalid Operations, and
358 propagate through subsequent arithmetic operations.
361 x then x is \*(Na; every other predicate
362 (x \*[Gt] y, x = y, x \*[Lt] y, ...) is FALSE if \*(Na is involved.
365 Trichotomy is violated by \*(Na.
366 Besides being FALSE, predicates that entail ordered
367 comparison, rather than mere (in)equality,
368 signal Invalid Operation when \*(Na is involved.
370 Every algebraic operation (+, \-, \(**, /,
372 is rounded by default to within half an
374 and when the rounding error is exactly half an
376 then the rounded value's least significant bit is zero.
377 This kind of rounding is usually the best kind,
378 sometimes provably so; for instance, for every
379 x = 1.0, 2.0, 3.0, 4.0, ..., 2.0**52, we find
380 (x/3.0)\(**3.0 == x and (x/10.0)\(**10.0 == x and ...
381 despite that both the quotients and the products
383 Only rounding like IEEE 754 can do that.
384 But no single kind of rounding can be
385 proved best for every circumstance, so IEEE 754
386 provides rounding towards zero or towards
390 at the programmer's option.
391 And the same kinds of rounding are specified for
392 Binary\-Decimal Conversions, at least for magnitudes
393 between roughly 1.0e\-10 and 1.0e37.
395 IEEE 754 recognizes five kinds of floating\-point exceptions,
396 listed below in declining order of probable importance.
397 .Bl -column "Invalid OperationX" "Gradual OverflowX"
398 .It Sy Exception Ta Sy Default Result
399 .It Invalid Operation \*(Na, or FALSE
400 .It Overflow \(+-\(if
401 .It Divide by Zero \(+-\(if \}
402 .It Underflow Gradual Underflow
403 .It Inexact Rounded value
407 An Exception is not an Error unless handled badly.
408 What makes a class of exceptions exceptional
409 is that no single default response can be satisfactory
411 On the other hand, if a default
412 response will serve most instances satisfactorily,
413 the unsatisfactory instances cannot justify aborting
414 computation every time the exception occurs.
417 For each kind of floating\-point exception, IEEE 754
418 provides a Flag that is raised each time its exception
419 is signaled, and stays raised until the program resets it.
420 Programs may also test, save and restore a flag.
421 Thus, IEEE 754 provides three ways by which programs
422 may cope with exceptions for which the default result
423 might be unsatisfactory:
426 Test for a condition that might cause an exception
427 later, and branch to avoid the exception.
429 Test a flag to see whether an exception has occurred
430 since the program last reset its flag.
432 Test a result to see whether it is a value that only
433 an exception could have produced.
435 The only reliable ways to discover
436 whether Underflow has occurred are to test whether
437 products or quotients lie closer to zero than the
438 underflow threshold, or to test the Underflow flag.
439 (Sums and differences cannot underflow in
442 y then x\-y is correct to
443 full precision and certainly nonzero regardless of
445 Products and quotients that
446 underflow gradually can lose accuracy gradually
447 without vanishing, so comparing them with zero
448 (as one might on a VAX) will not reveal the loss.
449 Fortunately, if a gradually underflowed value is
450 destined to be added to something bigger than the
451 underflow threshold, as is almost always the case,
452 digits lost to gradual underflow will not be missed
453 because they would have been rounded off anyway.
454 So gradual underflows are usually
457 The same cannot be said of underflows flushed to 0.
459 At the option of an implementor conforming to IEEE 754,
460 other ways to cope with exceptions may be provided:
463 This mechanism classifies an exception in
464 advance as an incident to be handled by means
465 traditionally associated with error\-handling
466 statements like "ON ERROR GO TO ...".
467 Different languages offer different forms of this statement,
468 but most share the following characteristics:
471 No means is provided to substitute a value for
472 the offending operation's result and resume
473 computation from what may be the middle of an expression.
474 An exceptional result is abandoned.
476 In a subprogram that lacks an error\-handling
477 statement, an exception causes the subprogram to
478 abort within whatever program called it, and so
479 on back up the chain of calling subprograms until
480 an error\-handling statement is encountered or the
481 whole task is aborted and memory is dumped.
485 This mechanism, requiring an interactive
486 debugging environment, is more for the programmer
488 It classifies an exception in
489 advance as a symptom of a programmer's error; the
490 exception suspends execution as near as it can to
491 the offending operation so that the programmer can
492 look around to see how it happened.
494 the first several exceptions turn out to be quite
495 unexceptionable, so the programmer ought ideally
496 to be able to resume execution after each one as if
497 execution had not been stopped.
499 \&... Other ways lie beyond the scope of this document.
502 The crucial problem for exception handling is the problem of
503 Scope, and the problem's solution is understood, but not
504 enough manpower was available to implement it fully in time
505 to be distributed in 4.3 BSD's libm.
506 Ideally, each elementary function should act
507 as if it were indivisible, or atomic, in the sense that ...
510 No exception should be signaled that is not deserved by
511 the data supplied to that function.
513 Any exception signaled should be identified with that
514 function rather than with one of its subroutines.
516 The internal behavior of an atomic function should not
517 be disrupted when a calling program changes from
518 one to another of the five or so ways of handling
519 exceptions listed above, although the definition
520 of the function may be correlated intentionally
521 with exception handling.
524 Ideally, every programmer should be able
526 to turn a debugged subprogram into one that appears atomic to
528 But simulating all three characteristics of an
529 atomic function is still a tedious affair, entailing hosts
530 of tests and saves\-restores; work is under way to ameliorate
533 Meanwhile, the functions in libm are only approximately atomic.
534 They signal no inappropriate exception except possibly ...
535 .Bl -ohang -offset indent
537 when a result, if properly computed, might have lain barely within range, and
538 .It Inexact in Fn cbrt , Fn hypot , Fn log10 No and Fn pow
539 when it happens to be exact, thanks to fortuitous cancellation of errors.
542 .Bl -ohang -offset indent
543 .It Invalid Operation is signaled only when
544 any result but \*(Na would probably be misleading.
545 .It Overflow is signaled only when
546 the exact result would be finite but beyond the overflow threshold.
547 .It Divide\-by\-Zero is signaled only when
548 a function takes exactly infinite values at finite operands.
549 .It Underflow is signaled only when
550 the exact result would be nonzero but tinier than the underflow threshold.
551 .It Inexact is signaled only when
552 greater range or precision would be needed to represent the exact result.
555 .\" .Bl -tag -width /usr/lib/libm_p.a -compact
556 .\" .It Pa /usr/lib/libm.a
557 .\" the static math library
558 .\" .It Pa /usr/lib/libm.so
559 .\" the dynamic math library
560 .\" .It Pa /usr/lib/libm_p.a
561 .\" the static math library compiled for profiling
564 An explanation of IEEE 754 and its proposed extension p854
565 was published in the IEEE magazine MICRO in August 1984 under
566 the title "A Proposed Radix\- and Word\-length\-independent
567 Standard for Floating\-point Arithmetic" by W. J. Cody et al.
568 The manuals for Pascal, C and BASIC on the Apple Macintosh
569 document the features of IEEE 754 pretty well.
570 Articles in the IEEE magazine COMPUTER vol. 14 no. 3 (Mar. 1981),
571 and in the ACM SIGNUM Newsletter Special Issue of
572 Oct. 1979, may be helpful although they pertain to
573 superseded drafts of the standard.
575 When signals are appropriate, they are emitted by certain
576 operations within the codes, so a subroutine\-trace may be
577 needed to identify the function with its signal in case
578 method 5) above is in use.
579 And the codes all take the
580 IEEE 754 defaults for granted; this means that a decision to
581 trap all divisions by zero could disrupt a code that would
582 otherwise get correct results despite division by zero.