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