1 .\" Copyright (c) 1985, 1993
2 .\" The Regents of the University of California. 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. All advertising materials mentioning features or use of this software
13 .\" must display the following acknowledgement:
14 .\" This product includes software developed by the University of
15 .\" California, Berkeley and its contributors.
16 .\" 4. Neither the name of the University nor the names of its contributors
17 .\" may be used to endorse or promote products derived from this software
18 .\" without specific prior written permission.
20 .\" THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
21 .\" ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
22 .\" IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
23 .\" ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
24 .\" FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
25 .\" DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
26 .\" OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
27 .\" HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
28 .\" LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
29 .\" OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
32 .\" @(#)math.3 8.2 (Berkeley) 5/5/94
33 .\" $FreeBSD: src/lib/libm/common_source/math.3,v 1.4.2.1 2001/07/22 12:07:16 dd Exp $
34 .\" $DragonFly: src/lib/libm/common_source/Attic/math.3,v 1.2 2003/06/17 04:26:50 dillon Exp $
36 .TH MATH 3 "May 5, 1994"
47 math \- introduction to mathematical library functions
49 These functions constitute the C math library,
51 The link editor searches this library under the \*(lq\-lm\*(rq option.
52 Declarations for these functions may be obtained from the include file
54 The Fortran math library is described in ``man 3f intro''.
55 .SH "LIST OF FUNCTIONS"
58 .ta \w'copysign'u+2n +\w'infnan.3m'u+10n +\w'inverse trigonometric func'u
59 \fIName\fP \fIAppears on Page\fP \fIDescription\fP \fIError Bound (ULPs)\fP
60 .ta \w'copysign'u+4n +\w'infnan.3m'u+4n +\w'inverse trigonometric function'u+6nC
62 acos sin.3m inverse trigonometric function 3
63 acosh asinh.3m inverse hyperbolic function 3
64 asin sin.3m inverse trigonometric function 3
65 asinh asinh.3m inverse hyperbolic function 3
66 atan sin.3m inverse trigonometric function 1
67 atanh asinh.3m inverse hyperbolic function 3
68 atan2 sin.3m inverse trigonometric function 2
69 cabs hypot.3m complex absolute value 1
70 cbrt sqrt.3m cube root 1
71 ceil floor.3m integer no less than 0
72 copysign ieee.3m copy sign bit 0
73 cos sin.3m trigonometric function 1
74 cosh sinh.3m hyperbolic function 3
75 drem ieee.3m remainder 0
76 erf erf.3m error function ???
77 erfc erf.3m complementary error function ???
78 exp exp.3m exponential 1
79 expm1 exp.3m exp(x)\-1 1
80 fabs floor.3m absolute value 0
81 floor floor.3m integer no greater than 0
82 hypot hypot.3m Euclidean distance 1
83 infnan infnan.3m signals exceptions
84 j0 j0.3m bessel function ???
85 j1 j0.3m bessel function ???
86 jn j0.3m bessel function ???
87 lgamma lgamma.3m log gamma function; (formerly gamma.3m)
88 log exp.3m natural logarithm 1
89 logb ieee.3m exponent extraction 0
90 log10 exp.3m logarithm to base 10 3
91 log1p exp.3m log(1+x) 1
92 pow exp.3m exponential x**y 60\-500
93 rint floor.3m round to nearest integer 0
94 scalb ieee.3m exponent adjustment 0
95 sin sin.3m trigonometric function 1
96 sinh sinh.3m hyperbolic function 3
97 sqrt sqrt.3m square root 1
98 tan sin.3m trigonometric function 3
99 tanh sinh.3m hyperbolic function 3
100 y0 j0.3m bessel function ???
101 y1 j0.3m bessel function ???
102 yn j0.3m bessel function ???
106 In 4.3 BSD, distributed from the University of California
107 in late 1985, most of the foregoing functions come in two
108 versions, one for the double\-precision "D" format in the
109 DEC VAX\-11 family of computers, another for double\-precision
110 arithmetic conforming to the IEEE Standard 754 for Binary
111 Floating\-Point Arithmetic. The two versions behave very
112 similarly, as should be expected from programs more accurate
113 and robust than was the norm when UNIX was born. For
114 instance, the programs are accurate to within the numbers
115 of \*(ups tabulated above; an \*(up is one \fIU\fRnit in the \fIL\fRast
116 \fIP\fRlace. And the programs have been cured of anomalies that
117 afflicted the older math library \fIlibm\fR in which incidents like
118 the following had been reported:
120 sqrt(\-1.0) = 0.0 and log(\-1.0) = \-1.7e38.
122 cos(1.0e\-11) > cos(0.0) > 1.0.
129 x when x = 2.0, 3.0, 4.0, ..., 9.0.
131 pow(\-1.0,1.0e10) trapped on Integer Overflow.
133 sqrt(1.0e30) and sqrt(1.0e\-30) were very slow.
135 However the two versions do differ in ways that have to be
136 explained, to which end the following notes are provided.
138 \fBDEC VAX\-11 D_floating\-point:\fR
140 This is the format for which the original math library \fIlibm\fR
141 was developed, and to which this manual is still principally
142 dedicated. It is \fIthe\fR double\-precision format for the PDP\-11
143 and the earlier VAX\-11 machines; VAX\-11s after 1983 were
144 provided with an optional "G" format closer to the IEEE
145 double\-precision format. The earlier DEC MicroVAXs have no
146 D format, only G double\-precision.
149 Properties of D_floating\-point:
151 Wordsize: 64 bits, 8 bytes. Radix: Binary.
158 bits, roughly like 17
165 If x and x' are consecutive positive D_floating\-point
166 numbers (they differ by 1 \*(up), then
168 1.3e\-17 < 0.5**56 < (x'\-x)/x \(<= 0.5**55 < 2.8e\-17.
171 .ta \w'Range:'u+1n +\w'Underflow threshold'u+1n +\w'= 2.0**127'u+1n
172 Range: Overflow threshold = 2.0**127 = 1.7e38.
173 Underflow threshold = 0.5**128 = 2.9e\-39.
174 NOTE: THIS RANGE IS COMPARATIVELY NARROW.
178 Overflow customarily stops computation.
180 Underflow is customarily flushed quietly to zero.
184 It is possible to have x
190 x\-y = 0 because of underflow. Similarly
191 x > y > 0 cannot prevent either x\(**y = 0
192 or y/x = 0 from happening without warning.
195 Zero is represented ambiguously.
197 Although 2**55 different representations of zero are accepted by
198 the hardware, only the obvious representation is ever produced.
199 There is no \-0 on a VAX.
202 is not part of the VAX architecture.
206 of the 2**55 that the hardware
207 recognizes, only one of them is ever produced.
208 Any floating\-point operation upon a reserved
209 operand, even a MOVF or MOVD, customarily stops
210 computation, so they are not much used.
214 Divisions by zero and operations that
215 overflow are invalid operations that customarily
216 stop computation or, in earlier machines, produce
217 reserved operands that will stop computation.
221 Every rational operation (+, \-, \(**, /) on a
222 VAX (but not necessarily on a PDP\-11), if not an
223 over/underflow nor division by zero, is rounded to
224 within half an \*(up, and when the rounding error is
225 exactly half an \*(up then rounding is away from 0.
229 Except for its narrow range, D_floating\-point is one of the
230 better computer arithmetics designed in the 1960's.
231 Its properties are reflected fairly faithfully in the elementary
232 functions for a VAX distributed in 4.3 BSD.
233 They over/underflow only if their results have to lie out of range
234 or very nearly so, and then they behave much as any rational
235 arithmetic operation that over/underflowed would behave.
236 Similarly, expressions like log(0) and atanh(1) behave
237 like 1/0; and sqrt(\-3) and acos(3) behave like 0/0;
238 they all produce reserved operands and/or stop computation!
239 The situation is described in more detail in manual pages.
242 \fIThis response seems excessively punitive, so it is destined
243 to be replaced at some time in the foreseeable future by a
244 more flexible but still uniform scheme being developed to
245 handle all floating\-point arithmetic exceptions neatly.
246 See infnan(3M) for the present state of affairs.\fR
250 How do the functions in 4.3 BSD's new \fIlibm\fR for UNIX
251 compare with their counterparts in DEC's VAX/VMS library?
252 Some of the VMS functions are a little faster, some are
253 a little more accurate, some are more puritanical about
254 exceptions (like pow(0.0,0.0) and atan2(0.0,0.0)),
255 and most occupy much more memory than their counterparts in
257 The VMS codes interpolate in large table to achieve
258 speed and accuracy; the \fIlibm\fR codes use tricky formulas
259 compact enough that all of them may some day fit into a ROM.
261 More important, DEC regards the VMS codes as proprietary
262 and guards them zealously against unauthorized use. But the
263 \fIlibm\fR codes in 4.3 BSD are intended for the public domain;
264 they may be copied freely provided their provenance is always
265 acknowledged, and provided users assist the authors in their
266 researches by reporting experience with the codes.
267 Therefore no user of UNIX on a machine whose arithmetic resembles
268 VAX D_floating\-point need use anything worse than the new \fIlibm\fR.
270 \fBIEEE STANDARD 754 Floating\-Point Arithmetic:\fR
272 This standard is on its way to becoming more widely adopted
273 than any other design for computer arithmetic.
274 VLSI chips that conform to some version of that standard have been
275 produced by a host of manufacturers, among them ...
277 .ta 0.5i +\w'Intel i8070, i80287'u+6n
278 Intel i8087, i80287 National Semiconductor 32081
279 Motorola 68881 Weitek WTL-1032, ... , -1165
280 Zilog Z8070 Western Electric (AT&T) WE32106.
283 Other implementations range from software, done thoroughly
284 in the Apple Macintosh, through VLSI in the Hewlett\-Packard
285 9000 series, to the ELXSI 6400 running ECL at 3 Megaflops.
286 Several other companies have adopted the formats
287 of IEEE 754 without, alas, adhering to the standard's way
288 of handling rounding and exceptions like over/underflow.
289 The DEC VAX G_floating\-point format is very similar to the IEEE
290 754 Double format, so similar that the C programs for the
291 IEEE versions of most of the elementary functions listed
292 above could easily be converted to run on a MicroVAX, though
293 nobody has volunteered to do that yet.
295 The codes in 4.3 BSD's \fIlibm\fR for machines that conform to
296 IEEE 754 are intended primarily for the National Semi. 32081
297 and WTL 1164/65. To use these codes with the Intel or Zilog
298 chips, or with the Apple Macintosh or ELXSI 6400, is to
299 forego the use of better codes provided (perhaps freely) by
300 those companies and designed by some of the authors of the
302 Except for \fIatan\fR, \fIcabs\fR, \fIcbrt\fR, \fIerf\fR,
303 \fIerfc\fR, \fIhypot\fR, \fIj0\-jn\fR, \fIlgamma\fR, \fIpow\fR
305 the Motorola 68881 has all the functions in \fIlibm\fR on chip,
306 and faster and more accurate;
307 it, Apple, the i8087, Z8070 and WE32106 all use 64
313 The main virtue of 4.3 BSD's
314 \fIlibm\fR codes is that they are intended for the public domain;
315 they may be copied freely provided their provenance is always
316 acknowledged, and provided users assist the authors in their
317 researches by reporting experience with the codes.
318 Therefore no user of UNIX on a machine that conforms to
319 IEEE 754 need use anything worse than the new \fIlibm\fR.
321 Properties of IEEE 754 Double\-Precision:
323 Wordsize: 64 bits, 8 bytes. Radix: Binary.
330 bits, roughly like 16
337 If x and x' are consecutive positive Double\-Precision
338 numbers (they differ by 1 \*(up), then
340 1.1e\-16 < 0.5**53 < (x'\-x)/x \(<= 0.5**52 < 2.3e\-16.
343 .ta \w'Range:'u+1n +\w'Underflow threshold'u+1n +\w'= 2.0**1024'u+1n
344 Range: Overflow threshold = 2.0**1024 = 1.8e308
345 Underflow threshold = 0.5**1022 = 2.2e\-308
349 Overflow goes by default to a signed
352 Underflow is \fIGradual,\fR rounding to the nearest
353 integer multiple of 0.5**1074 = 4.9e\-324.
355 Zero is represented ambiguously as +0 or \-0.
357 Its sign transforms correctly through multiplication or
358 division, and is preserved by addition of zeros
359 with like signs; but x\-x yields +0 for every
360 finite x. The only operations that reveal zero's
361 sign are division by zero and copysign(x,\(+-0).
362 In particular, comparison (x > y, x \(>= y, etc.)
363 cannot be affected by the sign of zero; but if
377 it persists when added to itself
378 or to any finite number. Its sign transforms
379 correctly through multiplication and division, and
380 .If (finite)/\(+- \0=\0\(+-0
385 Infinity\-Infinity, Infinity\(**0 and Infinity/Infinity
387 \(if\-\(if, \(if\(**0 and \(if/\(if
388 are, like 0/0 and sqrt(\-3),
389 invalid operations that produce \*(nn. ...
393 there are 2**53\-2 of them, all
394 called \*(nn (\fIN\fRot \fIa N\fRumber).
395 Some, called Signaling \*(nns, trap any floating\-point operation
396 performed upon them; they are used to mark missing
397 or uninitialized values, or nonexistent elements
398 of arrays. The rest are Quiet \*(nns; they are
399 the default results of Invalid Operations, and
400 propagate through subsequent arithmetic operations.
406 x then x is \*(nn; every other predicate
407 (x > y, x = y, x < y, ...) is FALSE if \*(nn is involved.
409 NOTE: Trichotomy is violated by \*(nn.
411 Besides being FALSE, predicates that entail ordered
412 comparison, rather than mere (in)equality,
413 signal Invalid Operation when \*(nn is involved.
418 Every algebraic operation (+, \-, \(**, /,
423 is rounded by default to within half an \*(up, and
424 when the rounding error is exactly half an \*(up then
425 the rounded value's least significant bit is zero.
426 This kind of rounding is usually the best kind,
427 sometimes provably so; for instance, for every
428 x = 1.0, 2.0, 3.0, 4.0, ..., 2.0**52, we find
429 (x/3.0)\(**3.0 == x and (x/10.0)\(**10.0 == x and ...
430 despite that both the quotients and the products
431 have been rounded. Only rounding like IEEE 754
432 can do that. But no single kind of rounding can be
433 proved best for every circumstance, so IEEE 754
434 provides rounding towards zero or towards
438 at the programmer's option. And the
439 same kinds of rounding are specified for
440 Binary\-Decimal Conversions, at least for magnitudes
441 between roughly 1.0e\-10 and 1.0e37.
445 IEEE 754 recognizes five kinds of floating\-point exceptions,
446 listed below in declining order of probable importance.
449 .ta \w'Invalid Operation'u+6n +\w'Gradual Underflow'u+2n
450 Exception Default Result
452 Invalid Operation \*(nn, or FALSE
454 Overflow \(+-Infinity
455 Divide by Zero \(+-Infinity \}
458 Divide by Zero \(+-\(if \}
459 Underflow Gradual Underflow
460 Inexact Rounded value
464 NOTE: An Exception is not an Error unless handled
465 badly. What makes a class of exceptions exceptional
466 is that no single default response can be satisfactory
467 in every instance. On the other hand, if a default
468 response will serve most instances satisfactorily,
469 the unsatisfactory instances cannot justify aborting
470 computation every time the exception occurs.
473 For each kind of floating\-point exception, IEEE 754
474 provides a Flag that is raised each time its exception
475 is signaled, and stays raised until the program resets
476 it. Programs may also test, save and restore a flag.
477 Thus, IEEE 754 provides three ways by which programs
478 may cope with exceptions for which the default result
479 might be unsatisfactory:
481 Test for a condition that might cause an exception
482 later, and branch to avoid the exception.
484 Test a flag to see whether an exception has occurred
485 since the program last reset its flag.
487 Test a result to see whether it is a value that only
488 an exception could have produced.
490 CAUTION: The only reliable ways to discover
491 whether Underflow has occurred are to test whether
492 products or quotients lie closer to zero than the
493 underflow threshold, or to test the Underflow
494 flag. (Sums and differences cannot underflow in
500 y then x\-y is correct to
501 full precision and certainly nonzero regardless of
502 how tiny it may be.) Products and quotients that
503 underflow gradually can lose accuracy gradually
504 without vanishing, so comparing them with zero
505 (as one might on a VAX) will not reveal the loss.
506 Fortunately, if a gradually underflowed value is
507 destined to be added to something bigger than the
508 underflow threshold, as is almost always the case,
509 digits lost to gradual underflow will not be missed
510 because they would have been rounded off anyway.
511 So gradual underflows are usually \fIprovably\fR ignorable.
512 The same cannot be said of underflows flushed to 0.
515 At the option of an implementor conforming to IEEE 754,
516 other ways to cope with exceptions may be provided:
518 ABORT. This mechanism classifies an exception in
519 advance as an incident to be handled by means
520 traditionally associated with error\-handling
521 statements like "ON ERROR GO TO ...". Different
522 languages offer different forms of this statement,
523 but most share the following characteristics:
524 .IP \(em \w'\0\0\0\0'u
525 No means is provided to substitute a value for
526 the offending operation's result and resume
527 computation from what may be the middle of an
528 expression. An exceptional result is abandoned.
529 .IP \(em \w'\0\0\0\0'u
530 In a subprogram that lacks an error\-handling
531 statement, an exception causes the subprogram to
532 abort within whatever program called it, and so
533 on back up the chain of calling subprograms until
534 an error\-handling statement is encountered or the
535 whole task is aborted and memory is dumped.
537 STOP. This mechanism, requiring an interactive
538 debugging environment, is more for the programmer
539 than the program. It classifies an exception in
540 advance as a symptom of a programmer's error; the
541 exception suspends execution as near as it can to
542 the offending operation so that the programmer can
543 look around to see how it happened. Quite often
544 the first several exceptions turn out to be quite
545 unexceptionable, so the programmer ought ideally
546 to be able to resume execution after each one as if
547 execution had not been stopped.
549 \&... Other ways lie beyond the scope of this document.
552 The crucial problem for exception handling is the problem of
553 Scope, and the problem's solution is understood, but not
554 enough manpower was available to implement it fully in time
555 to be distributed in 4.3 BSD's \fIlibm\fR. Ideally, each
556 elementary function should act as if it were indivisible, or
557 atomic, in the sense that ...
559 No exception should be signaled that is not deserved by
560 the data supplied to that function.
562 Any exception signaled should be identified with that
563 function rather than with one of its subroutines.
564 .IP iii) \w'iii)'u+2n
565 The internal behavior of an atomic function should not
566 be disrupted when a calling program changes from
567 one to another of the five or so ways of handling
568 exceptions listed above, although the definition
569 of the function may be correlated intentionally
570 with exception handling.
572 Ideally, every programmer should be able \fIconveniently\fR to
573 turn a debugged subprogram into one that appears atomic to
574 its users. But simulating all three characteristics of an
575 atomic function is still a tedious affair, entailing hosts
576 of tests and saves\-restores; work is under way to ameliorate
579 Meanwhile, the functions in \fIlibm\fR are only approximately
580 atomic. They signal no inappropriate exception except
585 when a result, if properly computed, might have lain barely within range, and
587 Inexact in \fIcabs\fR, \fIcbrt\fR, \fIhypot\fR, \fIlog10\fR and \fIpow\fR
589 when it happens to be exact, thanks to fortuitous cancellation of errors.
594 Invalid Operation is signaled only when
596 any result but \*(nn would probably be misleading.
598 Overflow is signaled only when
600 the exact result would be finite but beyond the overflow threshold.
602 Divide\-by\-Zero is signaled only when
604 a function takes exactly infinite values at finite operands.
606 Underflow is signaled only when
608 the exact result would be nonzero but tinier than the underflow threshold.
610 Inexact is signaled only when
612 greater range or precision would be needed to represent the exact result.
616 When signals are appropriate, they are emitted by certain
617 operations within the codes, so a subroutine\-trace may be
618 needed to identify the function with its signal in case
619 method 5) above is in use. And the codes all take the
620 IEEE 754 defaults for granted; this means that a decision to
621 trap all divisions by zero could disrupt a code that would
622 otherwise get correct results despite division by zero.
624 An explanation of IEEE 754 and its proposed extension p854
625 was published in the IEEE magazine MICRO in August 1984 under
626 the title "A Proposed Radix\- and Word\-length\-independent
627 Standard for Floating\-point Arithmetic" by W. J. Cody et al.
628 The manuals for Pascal, C and BASIC on the Apple Macintosh
629 document the features of IEEE 754 pretty well.
630 Articles in the IEEE magazine COMPUTER vol. 14 no. 3 (Mar.
631 1981), and in the ACM SIGNUM Newsletter Special Issue of
632 Oct. 1979, may be helpful although they pertain to
633 superseded drafts of the standard.
projects / dragonfly.git / blob