1 .\" Copyright (c) 1985 Regents of the University of California.
2 .\" All rights reserved.
4 .\" Redistribution and use in source and binary forms, with or without
5 .\" modification, are permitted provided that the following conditions
7 .\" 1. Redistributions of source code must retain the above copyright
8 .\" notice, this list of conditions and the following disclaimer.
9 .\" 2. Redistributions in binary form must reproduce the above copyright
10 .\" notice, this list of conditions and the following disclaimer in the
11 .\" documentation and/or other materials provided with the distribution.
12 .\" 3. Neither the name of the University nor the names of its contributors
13 .\" may be used to endorse or promote products derived from this software
14 .\" without specific prior written permission.
16 .\" THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
17 .\" ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
18 .\" IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
19 .\" ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
20 .\" FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
21 .\" DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
22 .\" OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
23 .\" HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
24 .\" LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
25 .\" OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
28 .\" from: @(#)math.3 6.10 (Berkeley) 5/6/91
29 .\" $NetBSD: math.3,v 1.18 2003/12/03 23:31:21 jschauma Exp $
30 .\" $DragonFly: src/lib/libm/man/math.3,v 1.3 2007/06/17 17:46:01 pavalos Exp $
32 .TH MATH 3 "July 12, 2009"
43 math \- introduction to mathematical library functions
45 These functions constitute the C math library,
47 The link editor searches this library under the \*(lq\-lm\*(rq option.
48 Declarations for these functions may be obtained from the include file
49 .RI \*[Lt] math.h \*[Gt].
50 .\" The Fortran math library is described in ``man 3f intro''.
51 .SH "LIST OF FUNCTIONS"
54 .ta \w'copysign'u+2n +\w'lgamma.3'u+10n +\w'inverse trigonometric func'u
55 \fIName\fP \fIAppears on Page\fP \fIDescription\fP \fIError Bound (ULPs)\fP
56 .ta \w'copysign'u+4n +\w'lgamma.3'u+4n +\w'inverse trigonometric function'u+6nC
58 acos acos.3 inverse trigonometric function 3
59 acosh acosh.3 inverse hyperbolic function 3
60 asin asin.3 inverse trigonometric function 3
61 asinh asinh.3 inverse hyperbolic function 3
62 atan atan.3 inverse trigonometric function 1
63 atanh atanh.3 inverse hyperbolic function 3
64 atan2 atan2.3 inverse trigonometric function 2
65 cabs hypot.3 complex absolute value 1
66 cbrt sqrt.3 cube root 1
67 ceil ceil.3 integer no less than 0
68 copysign ieee.3 copy sign bit 0
69 cos cos.3 trigonometric function 1
70 cosh cosh.3 hyperbolic function 3
71 erf erf.3 error function ???
72 erfc erf.3 complementary error function ???
73 exp exp.3 exponential 1
74 expm1 exp.3 exp(x)\-1 1
75 fabs fabs.3 absolute value 0
76 fdim fdim.3 positive difference ???
77 finite ieee.3 test for finity 0
78 floor floor.3 integer no greater than 0
79 fmax fmax.3 maximum function ???
80 fmin fmin.3 minimum function ???
81 fmod fmod.3 remainder ???
82 hypot hypot.3 Euclidean distance 1
83 ilogb ieee.3 exponent extraction 0
84 isinf isinf.3 test for infinity 0
85 isnan isnan.3 test for not-a-number 0
86 j0 j0.3 Bessel function ???
87 j1 j0.3 Bessel function ???
88 jn j0.3 Bessel function ???
89 lgamma lgamma.3 log gamma function ???
90 log exp.3 natural logarithm 1
91 log10 exp.3 logarithm to base 10 3
92 log1p exp.3 log(1+x) 1
93 nan nan.3 return quiet \*(nn 0
94 nextafter ieee.3 next representable number 0
95 pow exp.3 exponential x**y 60\-500
96 remainder ieee.3 remainder 0
97 rint rint.3 round to nearest integer 0
98 scalbn ieee.3 exponent adjustment 0
99 sin sin.3 trigonometric function 1
100 sinh sinh.3 hyperbolic function 3
101 sqrt sqrt.3 square root 1
102 tan tan.3 trigonometric function 3
103 tanh tanh.3 hyperbolic function 3
104 trunc trunc.3 nearest integral value 3
105 y0 j0.3 Bessel function ???
106 y1 j0.3 Bessel function ???
107 yn j0.3 Bessel function ???
110 .SH "LIST OF DEFINED VALUES"
113 .ta \w'M_2_SQRTPI'u+2n +\w'1.12837916709551257390'u+4n +\w'2/sqrt(pi)'u+6nC
114 \fIName\fP \fIValue\fP \fIDescription\fP
115 .ta \w'M_2_SQRTPI'u+2n +\w'1.12837916709551257390'u+4n +\w'2/sqrt(pi)'u+6nC
117 M_E 2.7182818284590452354 e
118 M_LOG2E 1.4426950408889634074 log 2e
119 M_LOG10E 0.43429448190325182765 log 10e
120 M_LN2 0.69314718055994530942 log e2
121 M_LN10 2.30258509299404568402 log e10
122 M_PI 3.14159265358979323846 pi
123 M_PI_2 1.57079632679489661923 pi/2
124 M_PI_4 0.78539816339744830962 pi/4
125 M_1_PI 0.31830988618379067154 1/pi
126 M_2_PI 0.63661977236758134308 2/pi
127 M_2_SQRTPI 1.12837916709551257390 2/sqrt(pi)
128 M_SQRT2 1.41421356237309504880 sqrt(2)
129 M_SQRT1_2 0.70710678118654752440 1/sqrt(2)
133 In 4.3 BSD, distributed from the University of California
134 in late 1985, most of the foregoing functions come in two
135 versions, one for the double\-precision "D" format in the
136 DEC VAX\-11 family of computers, another for double\-precision
137 arithmetic conforming to the IEEE Standard 754 for Binary
138 Floating\-Point Arithmetic.
139 The two versions behave very
140 similarly, as should be expected from programs more accurate
141 and robust than was the norm when UNIX was born.
142 For instance, the programs are accurate to within the numbers
143 of \*(ups tabulated above; an \*(up is one \fIU\fRnit in the \fIL\fRast
145 And the programs have been cured of anomalies that
146 afflicted the older math library \fIlibm\fR in which incidents like
147 the following had been reported:
149 sqrt(\-1.0) = 0.0 and log(\-1.0) = \-1.7e38.
151 cos(1.0e\-11) \*[Gt] cos(0.0) \*[Gt] 1.0.
158 x when x = 2.0, 3.0, 4.0, ..., 9.0.
160 pow(\-1.0,1.0e10) trapped on Integer Overflow.
162 sqrt(1.0e30) and sqrt(1.0e\-30) were very slow.
164 However the two versions do differ in ways that have to be
165 explained, to which end the following notes are provided.
167 \fBDEC VAX\-11 D_floating\-point:\fR
169 This is the format for which the original math library \fIlibm\fR
170 was developed, and to which this manual is still principally dedicated.
171 It is \fIthe\fR double\-precision format for the PDP\-11
172 and the earlier VAX\-11 machines; VAX\-11s after 1983 were
173 provided with an optional "G" format closer to the IEEE
174 double\-precision format.
175 The earlier DEC MicroVAXs have no D format, only G double\-precision.
179 Properties of D_floating\-point:
181 Wordsize: 64 bits, 8 bytes.
189 bits, roughly like 17
196 If x and x' are consecutive positive D_floating\-point
197 numbers (they differ by 1 \*(up), then
199 1.3e\-17 \*[Lt] 0.5**56 \*[Lt] (x'\-x)/x \*[Le] 0.5**55 \*[Lt] 2.8e\-17.
202 .ta \w'Range:'u+1n +\w'Underflow threshold'u+1n +\w'= 2.0**127'u+1n
203 Range: Overflow threshold = 2.0**127 = 1.7e38.
204 Underflow threshold = 0.5**128 = 2.9e\-39.
205 NOTE: THIS RANGE IS COMPARATIVELY NARROW.
209 Overflow customarily stops computation.
211 Underflow is customarily flushed quietly to zero.
215 It is possible to have x
221 x\-y = 0 because of underflow.
222 Similarly x \*[Gt] y \*[Gt] 0 cannot prevent either x\(**y = 0
223 or y/x = 0 from happening without warning.
226 Zero is represented ambiguously.
228 Although 2**55 different representations of zero are accepted by
229 the hardware, only the obvious representation is ever produced.
230 There is no \-0 on a VAX.
233 is not part of the VAX architecture.
237 of the 2**55 that the hardware
238 recognizes, only one of them is ever produced.
239 Any floating\-point operation upon a reserved
240 operand, even a MOVF or MOVD, customarily stops
241 computation, so they are not much used.
245 Divisions by zero and operations that
246 overflow are invalid operations that customarily
247 stop computation or, in earlier machines, produce
248 reserved operands that will stop computation.
252 Every rational operation (+, \-, \(**, /) on a
253 VAX (but not necessarily on a PDP\-11), if not an
254 over/underflow nor division by zero, is rounded to
255 within half an \*(up, and when the rounding error is
256 exactly half an \*(up then rounding is away from 0.
260 Except for its narrow range, D_floating\-point is one of the
261 better computer arithmetics designed in the 1960's.
262 Its properties are reflected fairly faithfully in the elementary
263 functions for a VAX distributed in 4.3 BSD.
264 They over/underflow only if their results have to lie out of range
265 or very nearly so, and then they behave much as any rational
266 arithmetic operation that over/underflowed would behave.
267 Similarly, expressions like log(0) and atanh(1) behave
268 like 1/0; and sqrt(\-3) and acos(3) behave like 0/0;
269 they all produce reserved operands and/or stop computation!
270 The situation is described in more detail in manual pages.
273 \fIThis response seems excessively punitive, so it is destined
274 to be replaced at some time in the foreseeable future by a
275 more flexible but still uniform scheme being developed to
276 handle all floating\-point arithmetic exceptions neatly.\fR
280 How do the functions in 4.3 BSD's new \fIlibm\fR for UNIX
281 compare with their counterparts in DEC's VAX/VMS library?
282 Some of the VMS functions are a little faster, some are
283 a little more accurate, some are more puritanical about
284 exceptions (like pow(0.0,0.0) and atan2(0.0,0.0)),
285 and most occupy much more memory than their counterparts in
287 The VMS codes interpolate in large table to achieve
288 speed and accuracy; the \fIlibm\fR codes use tricky formulas
289 compact enough that all of them may some day fit into a ROM.
291 More important, DEC regards the VMS codes as proprietary
292 and guards them zealously against unauthorized use.
293 But the \fIlibm\fR codes in 4.3 BSD are intended for the public domain;
294 they may be copied freely provided their provenance is always
295 acknowledged, and provided users assist the authors in their
296 researches by reporting experience with the codes.
297 Therefore no user of UNIX on a machine whose arithmetic resembles
298 VAX D_floating\-point need use anything worse than the new \fIlibm\fR.
300 \fBIEEE STANDARD 754 Floating\-Point Arithmetic:\fR
302 This standard is on its way to becoming more widely adopted
303 than any other design for computer arithmetic.
304 VLSI chips that conform to some version of that standard have been
305 produced by a host of manufacturers, among them ...
307 .ta 0.5i +\w'Intel i8070, i80287'u+6n
308 Intel i8087, i80287 National Semiconductor 32081
309 Motorola 68881 Weitek WTL-1032, ... , -1165
310 Zilog Z8070 Western Electric (AT\*[Am]T) WE32106.
313 Other implementations range from software, done thoroughly
314 in the Apple Macintosh, through VLSI in the Hewlett\-Packard
315 9000 series, to the ELXSI 6400 running ECL at 3 Megaflops.
316 Several other companies have adopted the formats
317 of IEEE 754 without, alas, adhering to the standard's way
318 of handling rounding and exceptions like over/underflow.
319 The DEC VAX G_floating\-point format is very similar to the IEEE
320 754 Double format, so similar that the C programs for the
321 IEEE versions of most of the elementary functions listed
322 above could easily be converted to run on a MicroVAX, though
323 nobody has volunteered to do that yet.
325 The codes in 4.3 BSD's \fIlibm\fR for machines that conform to
326 IEEE 754 are intended primarily for the National Semi. 32081
328 To use these codes with the Intel or Zilog
329 chips, or with the Apple Macintosh or ELXSI 6400, is to
330 forego the use of better codes provided (perhaps freely) by
331 those companies and designed by some of the authors of the
333 Except for \fIatan\fR, \fIcabs\fR, \fIcbrt\fR, \fIerf\fR,
334 \fIerfc\fR, \fIhypot\fR, \fIj0\-jn\fR, \fIlgamma\fR, \fIpow\fR
336 the Motorola 68881 has all the functions in \fIlibm\fR on chip,
337 and faster and more accurate;
338 it, Apple, the i8087, Z8070 and WE32106 all use 64
344 The main virtue of 4.3 BSD's
345 \fIlibm\fR codes is that they are intended for the public domain;
346 they may be copied freely provided their provenance is always
347 acknowledged, and provided users assist the authors in their
348 researches by reporting experience with the codes.
349 Therefore no user of UNIX on a machine that conforms to
350 IEEE 754 need use anything worse than the new \fIlibm\fR.
352 Properties of IEEE 754 Double\-Precision:
354 Wordsize: 64 bits, 8 bytes.
362 bits, roughly like 16
369 If x and x' are consecutive positive Double\-Precision
370 numbers (they differ by 1 \*(up), then
372 1.1e\-16 \*[Lt] 0.5**53 \*[Lt] (x'\-x)/x \*[Le] 0.5**52 \*[Lt] 2.3e\-16.
375 .ta \w'Range:'u+1n +\w'Underflow threshold'u+1n +\w'= 2.0**1024'u+1n
376 Range: Overflow threshold = 2.0**1024 = 1.8e308
377 Underflow threshold = 0.5**1022 = 2.2e\-308
381 Overflow goes by default to a signed
384 Underflow is \fIGradual,\fR rounding to the nearest
385 integer multiple of 0.5**1074 = 4.9e\-324.
387 Zero is represented ambiguously as +0 or \-0.
389 Its sign transforms correctly through multiplication or
390 division, and is preserved by addition of zeros
391 with like signs; but x\-x yields +0 for every
393 The only operations that reveal zero's
394 sign are division by zero and copysign(x,\(+-0).
395 In particular, comparison (x \*[Gt] y, x \*[Ge] y, etc.)
396 cannot be affected by the sign of zero; but if
410 it persists when added to itself
411 or to any finite number.
413 correctly through multiplication and division, and
414 .If (finite)/\(+- \0=\0\(+-0
419 Infinity\-Infinity, Infinity\(**0 and Infinity/Infinity
421 \(if\-\(if, \(if\(**0 and \(if/\(if
422 are, like 0/0 and sqrt(\-3),
423 invalid operations that produce \*(nn. ...
427 there are 2**53\-2 of them, all
428 called \*(nn (\fIN\fRot \fIa N\fRumber).
429 Some, called Signaling \*(nns, trap any floating\-point operation
430 performed upon them; they are used to mark missing
431 or uninitialized values, or nonexistent elements of arrays.
432 The rest are Quiet \*(nns; they are
433 the default results of Invalid Operations, and
434 propagate through subsequent arithmetic operations.
440 x then x is \*(nn; every other predicate
441 (x \*[Gt] y, x = y, x \*[Lt] y, ...) is FALSE if \*(nn is involved.
443 NOTE: Trichotomy is violated by \*(nn.
445 Besides being FALSE, predicates that entail ordered
446 comparison, rather than mere (in)equality,
447 signal Invalid Operation when \*(nn is involved.
452 Every algebraic operation (+, \-, \(**, /,
457 is rounded by default to within half an \*(up, and
458 when the rounding error is exactly half an \*(up then
459 the rounded value's least significant bit is zero.
460 This kind of rounding is usually the best kind,
461 sometimes provably so; for instance, for every
462 x = 1.0, 2.0, 3.0, 4.0, ..., 2.0**52, we find
463 (x/3.0)\(**3.0 == x and (x/10.0)\(**10.0 == x and ...
464 despite that both the quotients and the products
466 Only rounding like IEEE 754 can do that.
467 But no single kind of rounding can be
468 proved best for every circumstance, so IEEE 754
469 provides rounding towards zero or towards
473 at the programmer's option.
474 And the same kinds of rounding are specified for
475 Binary\-Decimal Conversions, at least for magnitudes
476 between roughly 1.0e\-10 and 1.0e37.
480 IEEE 754 recognizes five kinds of floating\-point exceptions,
481 listed below in declining order of probable importance.
484 .ta \w'Invalid Operation'u+6n +\w'Gradual Underflow'u+2n
485 Exception Default Result
489 Invalid Operation \*(nn, or FALSE
491 Overflow \(+-Infinity
492 Divide by Zero \(+-Infinity \}
495 Divide by Zero \(+-\(if \}
496 Underflow Gradual Underflow
497 Inexact Rounded value
501 NOTE: An Exception is not an Error unless handled badly.
502 What makes a class of exceptions exceptional
503 is that no single default response can be satisfactory
505 On the other hand, if a default
506 response will serve most instances satisfactorily,
507 the unsatisfactory instances cannot justify aborting
508 computation every time the exception occurs.
511 For each kind of floating\-point exception, IEEE 754
512 provides a Flag that is raised each time its exception
513 is signaled, and stays raised until the program resets it.
514 Programs may also test, save and restore a flag.
515 Thus, IEEE 754 provides three ways by which programs
516 may cope with exceptions for which the default result
517 might be unsatisfactory:
519 Test for a condition that might cause an exception
520 later, and branch to avoid the exception.
522 Test a flag to see whether an exception has occurred
523 since the program last reset its flag.
525 Test a result to see whether it is a value that only
526 an exception could have produced.
528 CAUTION: The only reliable ways to discover
529 whether Underflow has occurred are to test whether
530 products or quotients lie closer to zero than the
531 underflow threshold, or to test the Underflow flag.
532 (Sums and differences cannot underflow in
538 y then x\-y is correct to
539 full precision and certainly nonzero regardless of
541 Products and quotients that
542 underflow gradually can lose accuracy gradually
543 without vanishing, so comparing them with zero
544 (as one might on a VAX) will not reveal the loss.
545 Fortunately, if a gradually underflowed value is
546 destined to be added to something bigger than the
547 underflow threshold, as is almost always the case,
548 digits lost to gradual underflow will not be missed
549 because they would have been rounded off anyway.
550 So gradual underflows are usually \fIprovably\fR ignorable.
551 The same cannot be said of underflows flushed to 0.
554 At the option of an implementor conforming to IEEE 754,
555 other ways to cope with exceptions may be provided:
558 This mechanism classifies an exception in
559 advance as an incident to be handled by means
560 traditionally associated with error\-handling
561 statements like "ON ERROR GO TO ...".
562 Different languages offer different forms of this statement,
563 but most share the following characteristics:
564 .IP \(em \w'\0\0\0\0'u
565 No means is provided to substitute a value for
566 the offending operation's result and resume
567 computation from what may be the middle of an expression.
568 An exceptional result is abandoned.
569 .IP \(em \w'\0\0\0\0'u
570 In a subprogram that lacks an error\-handling
571 statement, an exception causes the subprogram to
572 abort within whatever program called it, and so
573 on back up the chain of calling subprograms until
574 an error\-handling statement is encountered or the
575 whole task is aborted and memory is dumped.
578 This mechanism, requiring an interactive
579 debugging environment, is more for the programmer
581 It classifies an exception in
582 advance as a symptom of a programmer's error; the
583 exception suspends execution as near as it can to
584 the offending operation so that the programmer can
585 look around to see how it happened.
587 the first several exceptions turn out to be quite
588 unexceptionable, so the programmer ought ideally
589 to be able to resume execution after each one as if
590 execution had not been stopped.
592 \&... Other ways lie beyond the scope of this document.
595 The crucial problem for exception handling is the problem of
596 Scope, and the problem's solution is understood, but not
597 enough manpower was available to implement it fully in time
598 to be distributed in 4.3 BSD's \fIlibm\fR.
599 Ideally, each elementary function should act
600 as if it were indivisible, or atomic, in the sense that ...
602 No exception should be signaled that is not deserved by
603 the data supplied to that function.
605 Any exception signaled should be identified with that
606 function rather than with one of its subroutines.
607 .IP iii) \w'iii)'u+2n
608 The internal behavior of an atomic function should not
609 be disrupted when a calling program changes from
610 one to another of the five or so ways of handling
611 exceptions listed above, although the definition
612 of the function may be correlated intentionally
613 with exception handling.
615 Ideally, every programmer should be able \fIconveniently\fR to
616 turn a debugged subprogram into one that appears atomic to
618 But simulating all three characteristics of an
619 atomic function is still a tedious affair, entailing hosts
620 of tests and saves\-restores; work is under way to ameliorate
623 Meanwhile, the functions in \fIlibm\fR are only approximately atomic.
624 They signal no inappropriate exception except possibly ...
628 when a result, if properly computed, might have lain barely within range, and
630 Inexact in \fIcabs\fR, \fIcbrt\fR, \fIhypot\fR, \fIlog10\fR and \fIpow\fR
632 when it happens to be exact, thanks to fortuitous cancellation of errors.
637 Invalid Operation is signaled only when
639 any result but \*(nn would probably be misleading.
641 Overflow is signaled only when
643 the exact result would be finite but beyond the overflow threshold.
645 Divide\-by\-Zero is signaled only when
647 a function takes exactly infinite values at finite operands.
649 Underflow is signaled only when
651 the exact result would be nonzero but tinier than the underflow threshold.
653 Inexact is signaled only when
655 greater range or precision would be needed to represent the exact result.
659 .\" .Bl -tag -width /usr/lib/libm_p.a -compact
660 .\" .It Pa /usr/lib/libm.a
661 .\" the static math library
662 .\" .It Pa /usr/lib/libm.so
663 .\" the dynamic math library
664 .\" .It Pa /usr/lib/libm_p.a
665 .\" the static math library compiled for profiling
668 An explanation of IEEE 754 and its proposed extension p854
669 was published in the IEEE magazine MICRO in August 1984 under
670 the title "A Proposed Radix\- and Word\-length\-independent
671 Standard for Floating\-point Arithmetic" by W. J. Cody et al.
672 The manuals for Pascal, C and BASIC on the Apple Macintosh
673 document the features of IEEE 754 pretty well.
674 Articles in the IEEE magazine COMPUTER vol. 14 no. 3 (Mar. 1981),
675 and in the ACM SIGNUM Newsletter Special Issue of
676 Oct. 1979, may be helpful although they pertain to
677 superseded drafts of the standard.
679 When signals are appropriate, they are emitted by certain
680 operations within the codes, so a subroutine\-trace may be
681 needed to identify the function with its signal in case
682 method 5) above is in use.
683 And the codes all take the
684 IEEE 754 defaults for granted; this means that a decision to
685 trap all divisions by zero could disrupt a code that would
686 otherwise get correct results despite division by zero.