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.1 2005/07/26 21:15:20 joerg Exp $
32 .TH MATH 3 "Dec 3, 2003"
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 finite ieee.3 test for finity 0
77 floor floor.3 integer no greater than 0
78 fmod fmod.3 remainder ???
79 hypot hypot.3 Euclidean distance 1
80 ilogb ieee.3 exponent extraction 0
81 isinf isinf.3 test for infinity 0
82 isnan isnan.3 test for not-a-number 0
83 j0 j0.3 Bessel function ???
84 j1 j0.3 Bessel function ???
85 jn j0.3 Bessel function ???
86 lgamma lgamma.3 log gamma function ???
87 log exp.3 natural logarithm 1
88 log10 exp.3 logarithm to base 10 3
89 log1p exp.3 log(1+x) 1
90 nextafter ieee.3 next representable number 0
91 pow exp.3 exponential x**y 60\-500
92 remainder ieee.3 remainder 0
93 rint rint.3 round to nearest integer 0
94 scalbn ieee.3 exponent adjustment 0
95 sin sin.3 trigonometric function 1
96 sinh sinh.3 hyperbolic function 3
97 sqrt sqrt.3 square root 1
98 tan tan.3 trigonometric function 3
99 tanh tanh.3 hyperbolic function 3
100 y0 j0.3 Bessel function ???
101 y1 j0.3 Bessel function ???
102 yn j0.3 Bessel function ???
105 .SH "LIST OF DEFINED VALUES"
108 .ta \w'M_2_SQRTPI'u+2n +\w'1.12837916709551257390'u+4n +\w'2/sqrt(pi)'u+6nC
109 \fIName\fP \fIValue\fP \fIDescription\fP
110 .ta \w'M_2_SQRTPI'u+2n +\w'1.12837916709551257390'u+4n +\w'2/sqrt(pi)'u+6nC
112 M_E 2.7182818284590452354 e
113 M_LOG2E 1.4426950408889634074 log 2e
114 M_LOG10E 0.43429448190325182765 log 10e
115 M_LN2 0.69314718055994530942 log e2
116 M_LN10 2.30258509299404568402 log e10
117 M_PI 3.14159265358979323846 pi
118 M_PI_2 1.57079632679489661923 pi/2
119 M_PI_4 0.78539816339744830962 pi/4
120 M_1_PI 0.31830988618379067154 1/pi
121 M_2_PI 0.63661977236758134308 2/pi
122 M_2_SQRTPI 1.12837916709551257390 2/sqrt(pi)
123 M_SQRT2 1.41421356237309504880 sqrt(2)
124 M_SQRT1_2 0.70710678118654752440 1/sqrt(2)
128 In 4.3 BSD, distributed from the University of California
129 in late 1985, most of the foregoing functions come in two
130 versions, one for the double\-precision "D" format in the
131 DEC VAX\-11 family of computers, another for double\-precision
132 arithmetic conforming to the IEEE Standard 754 for Binary
133 Floating\-Point Arithmetic.
134 The two versions behave very
135 similarly, as should be expected from programs more accurate
136 and robust than was the norm when UNIX was born.
137 For instance, the programs are accurate to within the numbers
138 of \*(ups tabulated above; an \*(up is one \fIU\fRnit in the \fIL\fRast
140 And the programs have been cured of anomalies that
141 afflicted the older math library \fIlibm\fR in which incidents like
142 the following had been reported:
144 sqrt(\-1.0) = 0.0 and log(\-1.0) = \-1.7e38.
146 cos(1.0e\-11) \*[Gt] cos(0.0) \*[Gt] 1.0.
153 x when x = 2.0, 3.0, 4.0, ..., 9.0.
155 pow(\-1.0,1.0e10) trapped on Integer Overflow.
157 sqrt(1.0e30) and sqrt(1.0e\-30) were very slow.
159 However the two versions do differ in ways that have to be
160 explained, to which end the following notes are provided.
162 \fBDEC VAX\-11 D_floating\-point:\fR
164 This is the format for which the original math library \fIlibm\fR
165 was developed, and to which this manual is still principally dedicated.
166 It is \fIthe\fR double\-precision format for the PDP\-11
167 and the earlier VAX\-11 machines; VAX\-11s after 1983 were
168 provided with an optional "G" format closer to the IEEE
169 double\-precision format.
170 The earlier DEC MicroVAXs have no D format, only G double\-precision.
174 Properties of D_floating\-point:
176 Wordsize: 64 bits, 8 bytes.
184 bits, roughly like 17
191 If x and x' are consecutive positive D_floating\-point
192 numbers (they differ by 1 \*(up), then
194 1.3e\-17 \*[Lt] 0.5**56 \*[Lt] (x'\-x)/x \*[Le] 0.5**55 \*[Lt] 2.8e\-17.
197 .ta \w'Range:'u+1n +\w'Underflow threshold'u+1n +\w'= 2.0**127'u+1n
198 Range: Overflow threshold = 2.0**127 = 1.7e38.
199 Underflow threshold = 0.5**128 = 2.9e\-39.
200 NOTE: THIS RANGE IS COMPARATIVELY NARROW.
204 Overflow customarily stops computation.
206 Underflow is customarily flushed quietly to zero.
210 It is possible to have x
216 x\-y = 0 because of underflow.
217 Similarly x \*[Gt] y \*[Gt] 0 cannot prevent either x\(**y = 0
218 or y/x = 0 from happening without warning.
221 Zero is represented ambiguously.
223 Although 2**55 different representations of zero are accepted by
224 the hardware, only the obvious representation is ever produced.
225 There is no \-0 on a VAX.
228 is not part of the VAX architecture.
232 of the 2**55 that the hardware
233 recognizes, only one of them is ever produced.
234 Any floating\-point operation upon a reserved
235 operand, even a MOVF or MOVD, customarily stops
236 computation, so they are not much used.
240 Divisions by zero and operations that
241 overflow are invalid operations that customarily
242 stop computation or, in earlier machines, produce
243 reserved operands that will stop computation.
247 Every rational operation (+, \-, \(**, /) on a
248 VAX (but not necessarily on a PDP\-11), if not an
249 over/underflow nor division by zero, is rounded to
250 within half an \*(up, and when the rounding error is
251 exactly half an \*(up then rounding is away from 0.
255 Except for its narrow range, D_floating\-point is one of the
256 better computer arithmetics designed in the 1960's.
257 Its properties are reflected fairly faithfully in the elementary
258 functions for a VAX distributed in 4.3 BSD.
259 They over/underflow only if their results have to lie out of range
260 or very nearly so, and then they behave much as any rational
261 arithmetic operation that over/underflowed would behave.
262 Similarly, expressions like log(0) and atanh(1) behave
263 like 1/0; and sqrt(\-3) and acos(3) behave like 0/0;
264 they all produce reserved operands and/or stop computation!
265 The situation is described in more detail in manual pages.
268 \fIThis response seems excessively punitive, so it is destined
269 to be replaced at some time in the foreseeable future by a
270 more flexible but still uniform scheme being developed to
271 handle all floating\-point arithmetic exceptions neatly.\fR
275 How do the functions in 4.3 BSD's new \fIlibm\fR for UNIX
276 compare with their counterparts in DEC's VAX/VMS library?
277 Some of the VMS functions are a little faster, some are
278 a little more accurate, some are more puritanical about
279 exceptions (like pow(0.0,0.0) and atan2(0.0,0.0)),
280 and most occupy much more memory than their counterparts in
282 The VMS codes interpolate in large table to achieve
283 speed and accuracy; the \fIlibm\fR codes use tricky formulas
284 compact enough that all of them may some day fit into a ROM.
286 More important, DEC regards the VMS codes as proprietary
287 and guards them zealously against unauthorized use.
288 But the \fIlibm\fR codes in 4.3 BSD are intended for the public domain;
289 they may be copied freely provided their provenance is always
290 acknowledged, and provided users assist the authors in their
291 researches by reporting experience with the codes.
292 Therefore no user of UNIX on a machine whose arithmetic resembles
293 VAX D_floating\-point need use anything worse than the new \fIlibm\fR.
295 \fBIEEE STANDARD 754 Floating\-Point Arithmetic:\fR
297 This standard is on its way to becoming more widely adopted
298 than any other design for computer arithmetic.
299 VLSI chips that conform to some version of that standard have been
300 produced by a host of manufacturers, among them ...
302 .ta 0.5i +\w'Intel i8070, i80287'u+6n
303 Intel i8087, i80287 National Semiconductor 32081
304 Motorola 68881 Weitek WTL-1032, ... , -1165
305 Zilog Z8070 Western Electric (AT\*[Am]T) WE32106.
308 Other implementations range from software, done thoroughly
309 in the Apple Macintosh, through VLSI in the Hewlett\-Packard
310 9000 series, to the ELXSI 6400 running ECL at 3 Megaflops.
311 Several other companies have adopted the formats
312 of IEEE 754 without, alas, adhering to the standard's way
313 of handling rounding and exceptions like over/underflow.
314 The DEC VAX G_floating\-point format is very similar to the IEEE
315 754 Double format, so similar that the C programs for the
316 IEEE versions of most of the elementary functions listed
317 above could easily be converted to run on a MicroVAX, though
318 nobody has volunteered to do that yet.
320 The codes in 4.3 BSD's \fIlibm\fR for machines that conform to
321 IEEE 754 are intended primarily for the National Semi. 32081
323 To use these codes with the Intel or Zilog
324 chips, or with the Apple Macintosh or ELXSI 6400, is to
325 forego the use of better codes provided (perhaps freely) by
326 those companies and designed by some of the authors of the
328 Except for \fIatan\fR, \fIcabs\fR, \fIcbrt\fR, \fIerf\fR,
329 \fIerfc\fR, \fIhypot\fR, \fIj0\-jn\fR, \fIlgamma\fR, \fIpow\fR
331 the Motorola 68881 has all the functions in \fIlibm\fR on chip,
332 and faster and more accurate;
333 it, Apple, the i8087, Z8070 and WE32106 all use 64
339 The main virtue of 4.3 BSD's
340 \fIlibm\fR codes is that they are intended for the public domain;
341 they may be copied freely provided their provenance is always
342 acknowledged, and provided users assist the authors in their
343 researches by reporting experience with the codes.
344 Therefore no user of UNIX on a machine that conforms to
345 IEEE 754 need use anything worse than the new \fIlibm\fR.
347 Properties of IEEE 754 Double\-Precision:
349 Wordsize: 64 bits, 8 bytes.
357 bits, roughly like 16
364 If x and x' are consecutive positive Double\-Precision
365 numbers (they differ by 1 \*(up), then
367 1.1e\-16 \*[Lt] 0.5**53 \*[Lt] (x'\-x)/x \*[Le] 0.5**52 \*[Lt] 2.3e\-16.
370 .ta \w'Range:'u+1n +\w'Underflow threshold'u+1n +\w'= 2.0**1024'u+1n
371 Range: Overflow threshold = 2.0**1024 = 1.8e308
372 Underflow threshold = 0.5**1022 = 2.2e\-308
376 Overflow goes by default to a signed
379 Underflow is \fIGradual,\fR rounding to the nearest
380 integer multiple of 0.5**1074 = 4.9e\-324.
382 Zero is represented ambiguously as +0 or \-0.
384 Its sign transforms correctly through multiplication or
385 division, and is preserved by addition of zeros
386 with like signs; but x\-x yields +0 for every
388 The only operations that reveal zero's
389 sign are division by zero and copysign(x,\(+-0).
390 In particular, comparison (x \*[Gt] y, x \*[Ge] y, etc.)
391 cannot be affected by the sign of zero; but if
405 it persists when added to itself
406 or to any finite number.
408 correctly through multiplication and division, and
409 .If (finite)/\(+- \0=\0\(+-0
414 Infinity\-Infinity, Infinity\(**0 and Infinity/Infinity
416 \(if\-\(if, \(if\(**0 and \(if/\(if
417 are, like 0/0 and sqrt(\-3),
418 invalid operations that produce \*(nn. ...
422 there are 2**53\-2 of them, all
423 called \*(nn (\fIN\fRot \fIa N\fRumber).
424 Some, called Signaling \*(nns, trap any floating\-point operation
425 performed upon them; they are used to mark missing
426 or uninitialized values, or nonexistent elements of arrays.
427 The rest are Quiet \*(nns; they are
428 the default results of Invalid Operations, and
429 propagate through subsequent arithmetic operations.
435 x then x is \*(nn; every other predicate
436 (x \*[Gt] y, x = y, x \*[Lt] y, ...) is FALSE if \*(nn is involved.
438 NOTE: Trichotomy is violated by \*(nn.
440 Besides being FALSE, predicates that entail ordered
441 comparison, rather than mere (in)equality,
442 signal Invalid Operation when \*(nn is involved.
447 Every algebraic operation (+, \-, \(**, /,
452 is rounded by default to within half an \*(up, and
453 when the rounding error is exactly half an \*(up then
454 the rounded value's least significant bit is zero.
455 This kind of rounding is usually the best kind,
456 sometimes provably so; for instance, for every
457 x = 1.0, 2.0, 3.0, 4.0, ..., 2.0**52, we find
458 (x/3.0)\(**3.0 == x and (x/10.0)\(**10.0 == x and ...
459 despite that both the quotients and the products
461 Only rounding like IEEE 754 can do that.
462 But no single kind of rounding can be
463 proved best for every circumstance, so IEEE 754
464 provides rounding towards zero or towards
468 at the programmer's option.
469 And the same kinds of rounding are specified for
470 Binary\-Decimal Conversions, at least for magnitudes
471 between roughly 1.0e\-10 and 1.0e37.
475 IEEE 754 recognizes five kinds of floating\-point exceptions,
476 listed below in declining order of probable importance.
479 .ta \w'Invalid Operation'u+6n +\w'Gradual Underflow'u+2n
480 Exception Default Result
484 Invalid Operation \*(nn, or FALSE
486 Overflow \(+-Infinity
487 Divide by Zero \(+-Infinity \}
490 Divide by Zero \(+-\(if \}
491 Underflow Gradual Underflow
492 Inexact Rounded value
496 NOTE: An Exception is not an Error unless handled badly.
497 What makes a class of exceptions exceptional
498 is that no single default response can be satisfactory
500 On the other hand, if a default
501 response will serve most instances satisfactorily,
502 the unsatisfactory instances cannot justify aborting
503 computation every time the exception occurs.
506 For each kind of floating\-point exception, IEEE 754
507 provides a Flag that is raised each time its exception
508 is signaled, and stays raised until the program resets it.
509 Programs may also test, save and restore a flag.
510 Thus, IEEE 754 provides three ways by which programs
511 may cope with exceptions for which the default result
512 might be unsatisfactory:
514 Test for a condition that might cause an exception
515 later, and branch to avoid the exception.
517 Test a flag to see whether an exception has occurred
518 since the program last reset its flag.
520 Test a result to see whether it is a value that only
521 an exception could have produced.
523 CAUTION: The only reliable ways to discover
524 whether Underflow has occurred are to test whether
525 products or quotients lie closer to zero than the
526 underflow threshold, or to test the Underflow flag.
527 (Sums and differences cannot underflow in
533 y then x\-y is correct to
534 full precision and certainly nonzero regardless of
536 Products and quotients that
537 underflow gradually can lose accuracy gradually
538 without vanishing, so comparing them with zero
539 (as one might on a VAX) will not reveal the loss.
540 Fortunately, if a gradually underflowed value is
541 destined to be added to something bigger than the
542 underflow threshold, as is almost always the case,
543 digits lost to gradual underflow will not be missed
544 because they would have been rounded off anyway.
545 So gradual underflows are usually \fIprovably\fR ignorable.
546 The same cannot be said of underflows flushed to 0.
549 At the option of an implementor conforming to IEEE 754,
550 other ways to cope with exceptions may be provided:
553 This mechanism classifies an exception in
554 advance as an incident to be handled by means
555 traditionally associated with error\-handling
556 statements like "ON ERROR GO TO ...".
557 Different languages offer different forms of this statement,
558 but most share the following characteristics:
559 .IP \(em \w'\0\0\0\0'u
560 No means is provided to substitute a value for
561 the offending operation's result and resume
562 computation from what may be the middle of an expression.
563 An exceptional result is abandoned.
564 .IP \(em \w'\0\0\0\0'u
565 In a subprogram that lacks an error\-handling
566 statement, an exception causes the subprogram to
567 abort within whatever program called it, and so
568 on back up the chain of calling subprograms until
569 an error\-handling statement is encountered or the
570 whole task is aborted and memory is dumped.
573 This mechanism, requiring an interactive
574 debugging environment, is more for the programmer
576 It classifies an exception in
577 advance as a symptom of a programmer's error; the
578 exception suspends execution as near as it can to
579 the offending operation so that the programmer can
580 look around to see how it happened.
582 the first several exceptions turn out to be quite
583 unexceptionable, so the programmer ought ideally
584 to be able to resume execution after each one as if
585 execution had not been stopped.
587 \&... Other ways lie beyond the scope of this document.
590 The crucial problem for exception handling is the problem of
591 Scope, and the problem's solution is understood, but not
592 enough manpower was available to implement it fully in time
593 to be distributed in 4.3 BSD's \fIlibm\fR.
594 Ideally, each elementary function should act
595 as if it were indivisible, or atomic, in the sense that ...
597 No exception should be signaled that is not deserved by
598 the data supplied to that function.
600 Any exception signaled should be identified with that
601 function rather than with one of its subroutines.
602 .IP iii) \w'iii)'u+2n
603 The internal behavior of an atomic function should not
604 be disrupted when a calling program changes from
605 one to another of the five or so ways of handling
606 exceptions listed above, although the definition
607 of the function may be correlated intentionally
608 with exception handling.
610 Ideally, every programmer should be able \fIconveniently\fR to
611 turn a debugged subprogram into one that appears atomic to
613 But simulating all three characteristics of an
614 atomic function is still a tedious affair, entailing hosts
615 of tests and saves\-restores; work is under way to ameliorate
618 Meanwhile, the functions in \fIlibm\fR are only approximately atomic.
619 They signal no inappropriate exception except possibly ...
623 when a result, if properly computed, might have lain barely within range, and
625 Inexact in \fIcabs\fR, \fIcbrt\fR, \fIhypot\fR, \fIlog10\fR and \fIpow\fR
627 when it happens to be exact, thanks to fortuitous cancellation of errors.
632 Invalid Operation is signaled only when
634 any result but \*(nn would probably be misleading.
636 Overflow is signaled only when
638 the exact result would be finite but beyond the overflow threshold.
640 Divide\-by\-Zero is signaled only when
642 a function takes exactly infinite values at finite operands.
644 Underflow is signaled only when
646 the exact result would be nonzero but tinier than the underflow threshold.
648 Inexact is signaled only when
650 greater range or precision would be needed to represent the exact result.
654 .\" .Bl -tag -width /usr/lib/libm_p.a -compact
655 .\" .It Pa /usr/lib/libm.a
656 .\" the static math library
657 .\" .It Pa /usr/lib/libm.so
658 .\" the dynamic math library
659 .\" .It Pa /usr/lib/libm_p.a
660 .\" the static math library compiled for profiling
663 An explanation of IEEE 754 and its proposed extension p854
664 was published in the IEEE magazine MICRO in August 1984 under
665 the title "A Proposed Radix\- and Word\-length\-independent
666 Standard for Floating\-point Arithmetic" by W. J. Cody et al.
667 The manuals for Pascal, C and BASIC on the Apple Macintosh
668 document the features of IEEE 754 pretty well.
669 Articles in the IEEE magazine COMPUTER vol. 14 no. 3 (Mar. 1981),
670 and in the ACM SIGNUM Newsletter Special Issue of
671 Oct. 1979, may be helpful although they pertain to
672 superseded drafts of the standard.
674 When signals are appropriate, they are emitted by certain
675 operations within the codes, so a subroutine\-trace may be
676 needed to identify the function with its signal in case
677 method 5) above is in use.
678 And the codes all take the
679 IEEE 754 defaults for granted; this means that a decision to
680 trap all divisions by zero could disrupt a code that would
681 otherwise get correct results despite division by zero.