![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to math.3 CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [CSRG BSD Unix] / 43BSDReno / lib / libm / common_source|
Return to math.3 CVS log | Up to [CSRG BSD Unix] / 43BSDReno / lib / libm / common_source |
1.1 root 1: .\" Copyright (c) 1985 Regents of the University of California.
2: .\" All rights reserved. The Berkeley software License Agreement
3: .\" specifies the terms and conditions for redistribution.
4: .\"
5: .\" @(#)math.3 6.8 (Berkeley) 5/27/86
6: .\"
7: .TH MATH 3M "May 27, 1986"
8: .UC 4
9: .ds up \fIulp\fR
10: .ds nn \fINaN\fR
11: .de If
12: .if n \\
13: \\$1Infinity\\$2
14: .if t \\
15: \\$1\\(if\\$2
16: ..
17: .SH NAME
18: math \- introduction to mathematical library functions
19: .SH DESCRIPTION
20: These functions constitute the C math library,
21: .I libm.
22: The link editor searches this library under the \*(lq\-lm\*(rq option.
23: Declarations for these functions may be obtained from the include file
24: .RI < math.h >.
25: The Fortran math library is described in ``man 3f intro''.
26: .SH "LIST OF FUNCTIONS"
27: .sp 2
28: .nf
29: .ta \w'copysign'u+2n +\w'infnan.3m'u+10n +\w'inverse trigonometric func'u
30: \fIName\fP \fIAppears on Page\fP \fIDescription\fP \fIError Bound (ULPs)\fP
31: .ta \w'copysign'u+4n +\w'infnan.3m'u+4n +\w'inverse trigonometric function'u+6nC
32: .sp 5p
33: acos sin.3m inverse trigonometric function 3
34: acosh asinh.3m inverse hyperbolic function 3
35: asin sin.3m inverse trigonometric function 3
36: asinh asinh.3m inverse hyperbolic function 3
37: atan sin.3m inverse trigonometric function 1
38: atanh asinh.3m inverse hyperbolic function 3
39: atan2 sin.3m inverse trigonometric function 2
40: cabs hypot.3m complex absolute value 1
41: cbrt sqrt.3m cube root 1
42: ceil floor.3m integer no less than 0
43: copysign ieee.3m copy sign bit 0
44: cos sin.3m trigonometric function 1
45: cosh sinh.3m hyperbolic function 3
46: drem ieee.3m remainder 0
47: erf erf.3m error function ???
48: erfc erf.3m complementary error function ???
49: exp exp.3m exponential 1
50: expm1 exp.3m exp(x)\-1 1
51: fabs floor.3m absolute value 0
52: floor floor.3m integer no greater than 0
53: hypot hypot.3m Euclidean distance 1
54: infnan infnan.3m signals exceptions
55: j0 j0.3m bessel function ???
56: j1 j0.3m bessel function ???
57: jn j0.3m bessel function ???
58: lgamma lgamma.3m log gamma function; (formerly gamma.3m)
59: log exp.3m natural logarithm 1
60: logb ieee.3m exponent extraction 0
61: log10 exp.3m logarithm to base 10 3
62: log1p exp.3m log(1+x) 1
63: pow exp.3m exponential x**y 60\-500
64: rint floor.3m round to nearest integer 0
65: scalb ieee.3m exponent adjustment 0
66: sin sin.3m trigonometric function 1
67: sinh sinh.3m hyperbolic function 3
68: sqrt sqrt.3m square root 1
69: tan sin.3m trigonometric function 3
70: tanh sinh.3m hyperbolic function 3
71: y0 j0.3m bessel function ???
72: y1 j0.3m bessel function ???
73: yn j0.3m bessel function ???
74: .ta
75: .fi
76: .SH NOTES
77: In 4.3 BSD, distributed from the University of California
78: in late 1985, most of the foregoing functions come in two
79: versions, one for the double\-precision "D" format in the
80: DEC VAX\-11 family of computers, another for double\-precision
81: arithmetic conforming to the IEEE Standard 754 for Binary
82: Floating\-Point Arithmetic. The two versions behave very
83: similarly, as should be expected from programs more accurate
84: and robust than was the norm when UNIX was born. For
85: instance, the programs are accurate to within the numbers
86: of \*(ups tabulated above; an \*(up is one \fIU\fRnit in the \fIL\fRast
87: \fIP\fRlace. And the programs have been cured of anomalies that
88: afflicted the older math library \fIlibm\fR in which incidents like
89: the following had been reported:
90: .RS
91: sqrt(\-1.0) = 0.0 and log(\-1.0) = \-1.7e38.
92: .br
93: cos(1.0e\-11) > cos(0.0) > 1.0.
94: .br
95: pow(x,1.0)
96: .if n \
97: !=
98: .if t \
99: \(!=
100: x when x = 2.0, 3.0, 4.0, ..., 9.0.
101: .br
102: pow(\-1.0,1.0e10) trapped on Integer Overflow.
103: .br
104: sqrt(1.0e30) and sqrt(1.0e\-30) were very slow.
105: .RE
106: However the two versions do differ in ways that have to be
107: explained, to which end the following notes are provided.
108: .PP
109: \fBDEC VAX\-11 D_floating\-point:\fR
110: .PP
111: This is the format for which the original math library \fIlibm\fR
112: was developed, and to which this manual is still principally
113: dedicated. It is \fIthe\fR double\-precision format for the PDP\-11
114: and the earlier VAX\-11 machines; VAX\-11s after 1983 were
115: provided with an optional "G" format closer to the IEEE
116: double\-precision format. The earlier DEC MicroVAXs have no
117: D format, only G double\-precision. (Why? Why not?)
118: .PP
119: Properties of D_floating\-point:
120: .RS
121: Wordsize: 64 bits, 8 bytes. Radix: Binary.
122: .br
123: Precision: 56
124: .if n \
125: sig.
126: .if t \
127: significant
128: bits, roughly like 17
129: .if n \
130: sig.
131: .if t \
132: significant
133: decimals.
134: .RS
135: If x and x' are consecutive positive D_floating\-point
136: numbers (they differ by 1 \*(up), then
137: .br
138: 1.3e\-17 < 0.5**56 < (x'\-x)/x \(<= 0.5**55 < 2.8e\-17.
139: .RE
140: .nf
141: .ta \w'Range:'u+1n +\w'Underflow threshold'u+1n +\w'= 2.0**127'u+1n
142: Range: Overflow threshold = 2.0**127 = 1.7e38.
143: Underflow threshold = 0.5**128 = 2.9e\-39.
144: NOTE: THIS RANGE IS COMPARATIVELY NARROW.
145: .ta
146: .fi
147: .RS
148: Overflow customarily stops computation.
149: .br
150: Underflow is customarily flushed quietly to zero.
151: .br
152: CAUTION:
153: .RS
154: It is possible to have x
155: .if n \
156: !=
157: .if t \
158: \(!=
159: y and yet
160: x\-y = 0 because of underflow. Similarly
161: x > y > 0 cannot prevent either x\(**y = 0
162: or y/x = 0 from happening without warning.
163: .RE
164: .RE
165: Zero is represented ambiguously.
166: .RS
167: Although 2**55 different representations of zero are accepted by
168: the hardware, only the obvious representation is ever produced.
169: There is no \-0 on a VAX.
170: .RE
171: .If
172: is not part of the VAX architecture.
173: .br
174: Reserved operands:
175: .RS
176: of the 2**55 that the hardware
177: recognizes, only one of them is ever produced.
178: Any floating\-point operation upon a reserved
179: operand, even a MOVF or MOVD, customarily stops
180: computation, so they are not much used.
181: .RE
182: Exceptions:
183: .RS
184: Divisions by zero and operations that
185: overflow are invalid operations that customarily
186: stop computation or, in earlier machines, produce
187: reserved operands that will stop computation.
188: .RE
189: Rounding:
190: .RS
191: Every rational operation (+, \-, \(**, /) on a
192: VAX (but not necessarily on a PDP\-11), if not an
193: over/underflow nor division by zero, is rounded to
194: within half an \*(up, and when the rounding error is
195: exactly half an \*(up then rounding is away from 0.
196: .RE
197: .RE
198: .PP
199: Except for its narrow range, D_floating\-point is one of the
200: better computer arithmetics designed in the 1960's.
201: Its properties are reflected fairly faithfully in the elementary
202: functions for a VAX distributed in 4.3 BSD.
203: They over/underflow only if their results have to lie out of range
204: or very nearly so, and then they behave much as any rational
205: arithmetic operation that over/underflowed would behave.
206: Similarly, expressions like log(0) and atanh(1) behave
207: like 1/0; and sqrt(\-3) and acos(3) behave like 0/0;
208: they all produce reserved operands and/or stop computation!
209: The situation is described in more detail in manual pages.
210: .RS
211: .ll -0.5i
212: \fIThis response seems excessively punitive, so it is destined
213: to be replaced at some time in the foreseeable future by a
214: more flexible but still uniform scheme being developed to
215: handle all floating\-point arithmetic exceptions neatly.
216: See infnan(3M) for the present state of affairs.\fR
217: .ll +0.5i
218: .RE
219: .PP
220: How do the functions in 4.3 BSD's new \fIlibm\fR for UNIX
221: compare with their counterparts in DEC's VAX/VMS library?
222: Some of the VMS functions are a little faster, some are
223: a little more accurate, some are more puritanical about
224: exceptions (like pow(0.0,0.0) and atan2(0.0,0.0)),
225: and most occupy much more memory than their counterparts in
226: \fIlibm\fR.
227: The VMS codes interpolate in large table to achieve
228: speed and accuracy; the \fIlibm\fR codes use tricky formulas
229: compact enough that all of them may some day fit into a ROM.
230: .PP
231: More important, DEC regards the VMS codes as proprietary
232: and guards them zealously against unauthorized use. But the
233: \fIlibm\fR codes in 4.3 BSD are intended for the public domain;
234: they may be copied freely provided their provenance is always
235: acknowledged, and provided users assist the authors in their
236: researches by reporting experience with the codes.
237: Therefore no user of UNIX on a machine whose arithmetic resembles
238: VAX D_floating\-point need use anything worse than the new \fIlibm\fR.
239: .PP
240: \fBIEEE STANDARD 754 Floating\-Point Arithmetic:\fR
241: .PP
242: This standard is on its way to becoming more widely adopted
243: than any other design for computer arithmetic.
244: VLSI chips that conform to some version of that standard have been
245: produced by a host of manufacturers, among them ...
246: .nf
247: .ta 0.5i +\w'Intel i8070, i80287'u+6n
248: Intel i8087, i80287 National Semiconductor 32081
249: Motorola 68881 Weitek WTL-1032, ... , -1165
250: Zilog Z8070 Western Electric (AT&T) WE32106.
251: .ta
252: .fi
253: Other implementations range from software, done thoroughly
254: in the Apple Macintosh, through VLSI in the Hewlett\-Packard
255: 9000 series, to the ELXSI 6400 running ECL at 3 Megaflops.
256: Several other companies have adopted the formats
257: of IEEE 754 without, alas, adhering to the standard's way
258: of handling rounding and exceptions like over/underflow.
259: The DEC VAX G_floating\-point format is very similar to the IEEE
260: 754 Double format, so similar that the C programs for the
261: IEEE versions of most of the elementary functions listed
262: above could easily be converted to run on a MicroVAX, though
263: nobody has volunteered to do that yet.
264: .PP
265: The codes in 4.3 BSD's \fIlibm\fR for machines that conform to
266: IEEE 754 are intended primarily for the National Semi. 32081
267: and WTL 1164/65. To use these codes with the Intel or Zilog
268: chips, or with the Apple Macintosh or ELXSI 6400, is to
269: forego the use of better codes provided (perhaps freely) by
270: those companies and designed by some of the authors of the
271: codes above.
272: Except for \fIatan\fR, \fIcabs\fR, \fIcbrt\fR, \fIerf\fR,
273: \fIerfc\fR, \fIhypot\fR, \fIj0\-jn\fR, \fIlgamma\fR, \fIpow\fR
274: and \fIy0\-yn\fR,
275: the Motorola 68881 has all the functions in \fIlibm\fR on chip,
276: and faster and more accurate;
277: it, Apple, the i8087, Z8070 and WE32106 all use 64
278: .if n \
279: sig.
280: .if t \
281: significant
282: bits.
283: The main virtue of 4.3 BSD's
284: \fIlibm\fR codes is that they are intended for the public domain;
285: they may be copied freely provided their provenance is always
286: acknowledged, and provided users assist the authors in their
287: researches by reporting experience with the codes.
288: Therefore no user of UNIX on a machine that conforms to
289: IEEE 754 need use anything worse than the new \fIlibm\fR.
290: .PP
291: Properties of IEEE 754 Double\-Precision:
292: .RS
293: Wordsize: 64 bits, 8 bytes. Radix: Binary.
294: .br
295: Precision: 53
296: .if n \
297: sig.
298: .if t \
299: significant
300: bits, roughly like 16
301: .if n \
302: sig.
303: .if t \
304: significant
305: decimals.
306: .RS
307: If x and x' are consecutive positive Double\-Precision
308: numbers (they differ by 1 \*(up), then
309: .br
310: 1.1e\-16 < 0.5**53 < (x'\-x)/x \(<= 0.5**52 < 2.3e\-16.
311: .RE
312: .nf
313: .ta \w'Range:'u+1n +\w'Underflow threshold'u+1n +\w'= 2.0**1024'u+1n
314: Range: Overflow threshold = 2.0**1024 = 1.8e308
315: Underflow threshold = 0.5**1022 = 2.2e\-308
316: .ta
317: .fi
318: .RS
319: Overflow goes by default to a signed
320: .If "" .
321: .br
322: Underflow is \fIGradual,\fR rounding to the nearest
323: integer multiple of 0.5**1074 = 4.9e\-324.
324: .RE
325: Zero is represented ambiguously as +0 or \-0.
326: .RS
327: Its sign transforms correctly through multiplication or
328: division, and is preserved by addition of zeros
329: with like signs; but x\-x yields +0 for every
330: finite x. The only operations that reveal zero's
331: sign are division by zero and copysign(x,\(+-0).
332: In particular, comparison (x > y, x \(>= y, etc.)
333: cannot be affected by the sign of zero; but if
334: finite x = y then
335: .If
336: \&= 1/(x\-y)
337: .if n \
338: !=
339: .if t \
340: \(!=
341: \-1/(y\-x) =
342: .If \- .
343: .RE
344: .If
345: is signed.
346: .RS
347: it persists when added to itself
348: or to any finite number. Its sign transforms
349: correctly through multiplication and division, and
350: .If (finite)/\(+- \0=\0\(+-0
351: (nonzero)/0 =
352: .If \(+- .
353: But
354: .if n \
355: Infinity\-Infinity, Infinity\(**0 and Infinity/Infinity
356: .if t \
357: \(if\-\(if, \(if\(**0 and \(if/\(if
358: are, like 0/0 and sqrt(\-3),
359: invalid operations that produce \*(nn. ...
360: .RE
361: Reserved operands:
362: .RS
363: there are 2**53\-2 of them, all
364: called \*(nn (\fIN\fRot \fIa N\fRumber).
365: Some, called Signaling \*(nns, trap any floating\-point operation
366: performed upon them; they are used to mark missing
367: or uninitialized values, or nonexistent elements
368: of arrays. The rest are Quiet \*(nns; they are
369: the default results of Invalid Operations, and
370: propagate through subsequent arithmetic operations.
371: If x
372: .if n \
373: !=
374: .if t \
375: \(!=
376: x then x is \*(nn; every other predicate
377: (x > y, x = y, x < y, ...) is FALSE if \*(nn is involved.
378: .br
379: NOTE: Trichotomy is violated by \*(nn.
380: .RS
381: Besides being FALSE, predicates that entail ordered
382: comparison, rather than mere (in)equality,
383: signal Invalid Operation when \*(nn is involved.
384: .RE
385: .RE
386: Rounding:
387: .RS
388: Every algebraic operation (+, \-, \(**, /,
389: .if n \
390: sqrt)
391: .if t \
392: \(sr)
393: is rounded by default to within half an \*(up, and
394: when the rounding error is exactly half an \*(up then
395: the rounded value's least significant bit is zero.
396: This kind of rounding is usually the best kind,
397: sometimes provably so; for instance, for every
398: x = 1.0, 2.0, 3.0, 4.0, ..., 2.0**52, we find
399: (x/3.0)\(**3.0 == x and (x/10.0)\(**10.0 == x and ...
400: despite that both the quotients and the products
401: have been rounded. Only rounding like IEEE 754
402: can do that. But no single kind of rounding can be
403: proved best for every circumstance, so IEEE 754
404: provides rounding towards zero or towards
405: .If +
406: or towards
407: .If \-
408: at the programmer's option. And the
409: same kinds of rounding are specified for
410: Binary\-Decimal Conversions, at least for magnitudes
411: between roughly 1.0e\-10 and 1.0e37.
412: .RE
413: Exceptions:
414: .RS
415: IEEE 754 recognizes five kinds of floating\-point exceptions,
416: listed below in declining order of probable importance.
417: .RS
418: .nf
419: .ta \w'Invalid Operation'u+6n +\w'Gradual Underflow'u+2n
420: Exception Default Result
421: .tc \(ru
422:
423: .tc
424: Invalid Operation \*(nn, or FALSE
425: .if n \{\
426: Overflow \(+-Infinity
427: Divide by Zero \(+-Infinity \}
428: .if t \{\
429: Overflow \(+-\(if
430: Divide by Zero \(+-\(if \}
431: Underflow Gradual Underflow
432: Inexact Rounded value
433: .ta
434: .fi
435: .RE
436: NOTE: An Exception is not an Error unless handled
437: badly. What makes a class of exceptions exceptional
438: is that no single default response can be satisfactory
439: in every instance. On the other hand, if a default
440: response will serve most instances satisfactorily,
441: the unsatisfactory instances cannot justify aborting
442: computation every time the exception occurs.
443: .RE
444: .PP
445: For each kind of floating\-point exception, IEEE 754
446: provides a Flag that is raised each time its exception
447: is signaled, and stays raised until the program resets
448: it. Programs may also test, save and restore a flag.
449: Thus, IEEE 754 provides three ways by which programs
450: may cope with exceptions for which the default result
451: might be unsatisfactory:
452: .IP 1) \w'\0\0\0\0'u
453: Test for a condition that might cause an exception
454: later, and branch to avoid the exception.
455: .IP 2) \w'\0\0\0\0'u
456: Test a flag to see whether an exception has occurred
457: since the program last reset its flag.
458: .IP 3) \w'\0\0\0\0'u
459: Test a result to see whether it is a value that only
460: an exception could have produced.
461: .RS
462: CAUTION: The only reliable ways to discover
463: whether Underflow has occurred are to test whether
464: products or quotients lie closer to zero than the
465: underflow threshold, or to test the Underflow
466: flag. (Sums and differences cannot underflow in
467: IEEE 754; if x
468: .if n \
469: !=
470: .if t \
471: \(!=
472: y then x\-y is correct to
473: full precision and certainly nonzero regardless of
474: how tiny it may be.) Products and quotients that
475: underflow gradually can lose accuracy gradually
476: without vanishing, so comparing them with zero
477: (as one might on a VAX) will not reveal the loss.
478: Fortunately, if a gradually underflowed value is
479: destined to be added to something bigger than the
480: underflow threshold, as is almost always the case,
481: digits lost to gradual underflow will not be missed
482: because they would have been rounded off anyway.
483: So gradual underflows are usually \fIprovably\fR ignorable.
484: The same cannot be said of underflows flushed to 0.
485: .RE
486: .PP
487: At the option of an implementor conforming to IEEE 754,
488: other ways to cope with exceptions may be provided:
489: .IP 4) \w'\0\0\0\0'u
490: ABORT. This mechanism classifies an exception in
491: advance as an incident to be handled by means
492: traditionally associated with error\-handling
493: statements like "ON ERROR GO TO ...". Different
494: languages offer different forms of this statement,
495: but most share the following characteristics:
496: .IP \(em \w'\0\0\0\0'u
497: No means is provided to substitute a value for
498: the offending operation's result and resume
499: computation from what may be the middle of an
500: expression. An exceptional result is abandoned.
501: .IP \(em \w'\0\0\0\0'u
502: In a subprogram that lacks an error\-handling
503: statement, an exception causes the subprogram to
504: abort within whatever program called it, and so
505: on back up the chain of calling subprograms until
506: an error\-handling statement is encountered or the
507: whole task is aborted and memory is dumped.
508: .IP 5) \w'\0\0\0\0'u
509: STOP. This mechanism, requiring an interactive
510: debugging environment, is more for the programmer
511: than the program. It classifies an exception in
512: advance as a symptom of a programmer's error; the
513: exception suspends execution as near as it can to
514: the offending operation so that the programmer can
515: look around to see how it happened. Quite often
516: the first several exceptions turn out to be quite
517: unexceptionable, so the programmer ought ideally
518: to be able to resume execution after each one as if
519: execution had not been stopped.
520: .IP 6) \w'\0\0\0\0'u
521: \&... Other ways lie beyond the scope of this document.
522: .RE
523: .PP
524: The crucial problem for exception handling is the problem of
525: Scope, and the problem's solution is understood, but not
526: enough manpower was available to implement it fully in time
527: to be distributed in 4.3 BSD's \fIlibm\fR. Ideally, each
528: elementary function should act as if it were indivisible, or
529: atomic, in the sense that ...
530: .IP i) \w'iii)'u+2n
531: No exception should be signaled that is not deserved by
532: the data supplied to that function.
533: .IP ii) \w'iii)'u+2n
534: Any exception signaled should be identified with that
535: function rather than with one of its subroutines.
536: .IP iii) \w'iii)'u+2n
537: The internal behavior of an atomic function should not
538: be disrupted when a calling program changes from
539: one to another of the five or so ways of handling
540: exceptions listed above, although the definition
541: of the function may be correlated intentionally
542: with exception handling.
543: .PP
544: Ideally, every programmer should be able \fIconveniently\fR to
545: turn a debugged subprogram into one that appears atomic to
546: its users. But simulating all three characteristics of an
547: atomic function is still a tedious affair, entailing hosts
548: of tests and saves\-restores; work is under way to ameliorate
549: the inconvenience.
550: .PP
551: Meanwhile, the functions in \fIlibm\fR are only approximately
552: atomic. They signal no inappropriate exception except
553: possibly ...
554: .RS
555: Over/Underflow
556: .RS
557: when a result, if properly computed, might have lain barely within range, and
558: .RE
559: Inexact in \fIcabs\fR, \fIcbrt\fR, \fIhypot\fR, \fIlog10\fR and \fIpow\fR
560: .RS
561: when it happens to be exact, thanks to fortuitous cancellation of errors.
562: .RE
563: .RE
564: Otherwise, ...
565: .RS
566: Invalid Operation is signaled only when
567: .RS
568: any result but \*(nn would probably be misleading.
569: .RE
570: Overflow is signaled only when
571: .RS
572: the exact result would be finite but beyond the overflow threshold.
573: .RE
574: Divide\-by\-Zero is signaled only when
575: .RS
576: a function takes exactly infinite values at finite operands.
577: .RE
578: Underflow is signaled only when
579: .RS
580: the exact result would be nonzero but tinier than the underflow threshold.
581: .RE
582: Inexact is signaled only when
583: .RS
584: greater range or precision would be needed to represent the exact result.
585: .RE
586: .RE
587: .SH BUGS
588: When signals are appropriate, they are emitted by certain
589: operations within the codes, so a subroutine\-trace may be
590: needed to identify the function with its signal in case
591: method 5) above is in use. And the codes all take the
592: IEEE 754 defaults for granted; this means that a decision to
593: trap all divisions by zero could disrupt a code that would
594: otherwise get correct results despite division by zero.
595: .SH SEE ALSO
596: An explanation of IEEE 754 and its proposed extension p854
597: was published in the IEEE magazine MICRO in August 1984 under
598: the title "A Proposed Radix\- and Word\-length\-independent
599: Standard for Floating\-point Arithmetic" by W. J. Cody et al.
600: The manuals for Pascal, C and BASIC on the Apple Macintosh
601: document the features of IEEE 754 pretty well.
602: Articles in the IEEE magazine COMPUTER vol. 14 no. 3 (Mar.
603: 1981), and in the ACM SIGNUM Newsletter Special Issue of
604: Oct. 1979, may be helpful although they pertain to
605: superseded drafts of the standard.
606: .SH AUTHOR
607: W. Kahan, with the help of Z\-S. Alex Liu, Stuart I. McDonald,
608: Dr. Kwok\-Choi Ng, Peter Tang.
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.