1*84d9c625SLionel Sambuc.\" $NetBSD: math.3,v 1.26 2012/11/10 15:59:58 njoly Exp $ 22fe8fb19SBen Gras.\" 32fe8fb19SBen Gras.\" Copyright (c) 1985 Regents of the University of California. 42fe8fb19SBen Gras.\" All rights reserved. 52fe8fb19SBen Gras.\" 62fe8fb19SBen Gras.\" Redistribution and use in source and binary forms, with or without 72fe8fb19SBen Gras.\" modification, are permitted provided that the following conditions 82fe8fb19SBen Gras.\" are met: 92fe8fb19SBen Gras.\" 1. Redistributions of source code must retain the above copyright 102fe8fb19SBen Gras.\" notice, this list of conditions and the following disclaimer. 112fe8fb19SBen Gras.\" 2. Redistributions in binary form must reproduce the above copyright 122fe8fb19SBen Gras.\" notice, this list of conditions and the following disclaimer in the 132fe8fb19SBen Gras.\" documentation and/or other materials provided with the distribution. 142fe8fb19SBen Gras.\" 3. Neither the name of the University nor the names of its contributors 152fe8fb19SBen Gras.\" may be used to endorse or promote products derived from this software 162fe8fb19SBen Gras.\" without specific prior written permission. 172fe8fb19SBen Gras.\" 182fe8fb19SBen Gras.\" THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 192fe8fb19SBen Gras.\" ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 202fe8fb19SBen Gras.\" IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 212fe8fb19SBen Gras.\" ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 222fe8fb19SBen Gras.\" FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 232fe8fb19SBen Gras.\" DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 242fe8fb19SBen Gras.\" OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 252fe8fb19SBen Gras.\" HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 262fe8fb19SBen Gras.\" LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 272fe8fb19SBen Gras.\" OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 282fe8fb19SBen Gras.\" SUCH DAMAGE. 292fe8fb19SBen Gras.\" 302fe8fb19SBen Gras.\" from: @(#)math.3 6.10 (Berkeley) 5/6/91 312fe8fb19SBen Gras.\" 32*84d9c625SLionel Sambuc.Dd February 23, 2007 33*84d9c625SLionel Sambuc.Dt MATH 3 34*84d9c625SLionel Sambuc.Os 35*84d9c625SLionel Sambuc.Sh NAME 36*84d9c625SLionel Sambuc.Nm math 37*84d9c625SLionel Sambuc.Nd introduction to mathematical library functions 38*84d9c625SLionel Sambuc.Sh LIBRARY 39*84d9c625SLionel Sambuc.Lb libm 40*84d9c625SLionel Sambuc.Sh SYNOPSIS 41*84d9c625SLionel Sambuc.In math.h 42*84d9c625SLionel Sambuc.Sh DESCRIPTION 43*84d9c625SLionel SambucThese functions constitute the C 44*84d9c625SLionel Sambuc.Lb libm . 452fe8fb19SBen GrasDeclarations for these functions may be obtained from the include file 46*84d9c625SLionel Sambuc.In math.h . 472fe8fb19SBen Gras.\" The Fortran math library is described in ``man 3f intro''. 48*84d9c625SLionel Sambuc.Ss List of Functions 49*84d9c625SLionel Sambuc.Bl -column "copysignX" "gammaX3XX" "inverse trigonometric funcX" 50*84d9c625SLionel Sambuc.It Sy Name Ta Sy Man page Ta Sy Description Ta Sy Error Bound Dv ( ULP Ns No s) 51*84d9c625SLionel Sambuc.It acos Ta Xr acos 3 Ta inverse trigonometric function Ta 3 52*84d9c625SLionel Sambuc.It acosh Ta Xr acosh 3 Ta inverse hyperbolic function Ta 3 53*84d9c625SLionel Sambuc.It asin Ta Xr asin 3 Ta inverse trigonometric function Ta 3 54*84d9c625SLionel Sambuc.It asinh Ta Xr asinh 3 Ta inverse hyperbolic function Ta 3 55*84d9c625SLionel Sambuc.It atan Ta Xr atan 3 Ta inverse trigonometric function Ta 1 56*84d9c625SLionel Sambuc.It atanh Ta Xr atanh 3 Ta inverse hyperbolic function Ta 3 57*84d9c625SLionel Sambuc.It atan2 Ta Xr atan2 3 Ta inverse trigonometric function Ta 2 58*84d9c625SLionel Sambuc.It cbrt Ta Xr sqrt 3 Ta cube root Ta 1 59*84d9c625SLionel Sambuc.It ceil Ta Xr ceil 3 Ta integer no less than Ta 0 60*84d9c625SLionel Sambuc.It copysign Ta Xr copysign 3 Ta copy sign bit Ta 0 61*84d9c625SLionel Sambuc.It cos Ta Xr cos 3 Ta trigonometric function Ta 1 62*84d9c625SLionel Sambuc.It cosh Ta Xr cosh 3 Ta hyperbolic function Ta 3 63*84d9c625SLionel Sambuc.It erf Ta Xr erf 3 Ta error function Ta ??? 64*84d9c625SLionel Sambuc.It erfc Ta Xr erf 3 Ta complementary error function Ta ??? 65*84d9c625SLionel Sambuc.It exp Ta Xr exp 3 Ta exponential Ta 1 66*84d9c625SLionel Sambuc.It expm1 Ta Xr exp 3 Ta exp(x)\-1 Ta 1 67*84d9c625SLionel Sambuc.It fabs Ta Xr fabs 3 Ta absolute value Ta 0 68*84d9c625SLionel Sambuc.It finite Ta Xr finite 3 Ta test for finity Ta 0 69*84d9c625SLionel Sambuc.It floor Ta Xr floor 3 Ta integer no greater than Ta 0 70*84d9c625SLionel Sambuc.It fmod Ta Xr fmod 3 Ta remainder Ta ??? 71*84d9c625SLionel Sambuc.It hypot Ta Xr hypot 3 Ta Euclidean distance Ta 1 72*84d9c625SLionel Sambuc.It ilogb Ta Xr ilogb 3 Ta exponent extraction Ta 0 73*84d9c625SLionel Sambuc.It isinf Ta Xr isinf 3 Ta test for infinity Ta 0 74*84d9c625SLionel Sambuc.It isnan Ta Xr isnan 3 Ta test for not-a-number Ta 0 75*84d9c625SLionel Sambuc.It j0 Ta Xr j0 3 Ta Bessel function Ta ??? 76*84d9c625SLionel Sambuc.It j1 Ta Xr j0 3 Ta Bessel function Ta ??? 77*84d9c625SLionel Sambuc.It jn Ta Xr j0 3 Ta Bessel function Ta ??? 78*84d9c625SLionel Sambuc.It lgamma Ta Xr lgamma 3 Ta log gamma function Ta ??? 79*84d9c625SLionel Sambuc.It log Ta Xr log 3 Ta natural logarithm Ta 1 80*84d9c625SLionel Sambuc.It log10 Ta Xr log 3 Ta logarithm to base 10 Ta 3 81*84d9c625SLionel Sambuc.It log1p Ta Xr log 3 Ta log(1+x) Ta 1 82*84d9c625SLionel Sambuc.It nan Ta Xr nan 3 Ta return quiet \*(Na Ta 0 83*84d9c625SLionel Sambuc.It nextafter Ta Xr nextafter 3 Ta next representable number Ta 0 84*84d9c625SLionel Sambuc.It pow Ta Xr pow 3 Ta exponential x**y Ta 60\-500 85*84d9c625SLionel Sambuc.It remainder Ta Xr remainder 3 Ta remainder Ta 0 86*84d9c625SLionel Sambuc.It rint Ta Xr rint 3 Ta round to nearest integer Ta 0 87*84d9c625SLionel Sambuc.It scalbn Ta Xr scalbn 3 Ta exponent adjustment Ta 0 88*84d9c625SLionel Sambuc.It sin Ta Xr sin 3 Ta trigonometric function Ta 1 89*84d9c625SLionel Sambuc.It sinh Ta Xr sinh 3 Ta hyperbolic function Ta 3 90*84d9c625SLionel Sambuc.It sqrt Ta Xr sqrt 3 Ta square root Ta 1 91*84d9c625SLionel Sambuc.It tan Ta Xr tan 3 Ta trigonometric function Ta 3 92*84d9c625SLionel Sambuc.It tanh Ta Xr tanh 3 Ta hyperbolic function Ta 3 93*84d9c625SLionel Sambuc.It trunc Ta Xr trunc 3 Ta nearest integral value Ta 3 94*84d9c625SLionel Sambuc.It y0 Ta Xr j0 3 Ta Bessel function Ta ??? 95*84d9c625SLionel Sambuc.It y1 Ta Xr j0 3 Ta Bessel function Ta ??? 96*84d9c625SLionel Sambuc.It yn Ta Xr j0 3 Ta Bessel function Ta ??? 97*84d9c625SLionel Sambuc.El 98*84d9c625SLionel Sambuc.Ss List of Defined Values 99*84d9c625SLionel Sambuc.Bl -column "M_2_SQRTPIXX" "1.12837916709551257390XX" "2/sqrt(pi)XXX" 100*84d9c625SLionel Sambuc.It Sy Name Ta Sy Value Ta Sy Description 101*84d9c625SLionel Sambuc.It M_E 2.7182818284590452354 e 102*84d9c625SLionel Sambuc.It M_LOG2E 1.4426950408889634074 log 2e 103*84d9c625SLionel Sambuc.It M_LOG10E 0.43429448190325182765 log 10e 104*84d9c625SLionel Sambuc.It M_LN2 0.69314718055994530942 log e2 105*84d9c625SLionel Sambuc.It M_LN10 2.30258509299404568402 log e10 106*84d9c625SLionel Sambuc.It M_PI 3.14159265358979323846 pi 107*84d9c625SLionel Sambuc.It M_PI_2 1.57079632679489661923 pi/2 108*84d9c625SLionel Sambuc.It M_PI_4 0.78539816339744830962 pi/4 109*84d9c625SLionel Sambuc.It M_1_PI 0.31830988618379067154 1/pi 110*84d9c625SLionel Sambuc.It M_2_PI 0.63661977236758134308 2/pi 111*84d9c625SLionel Sambuc.It M_2_SQRTPI 1.12837916709551257390 2/sqrt(pi) 112*84d9c625SLionel Sambuc.It M_SQRT2 1.41421356237309504880 sqrt(2) 113*84d9c625SLionel Sambuc.It M_SQRT1_2 0.70710678118654752440 1/sqrt(2) 114*84d9c625SLionel Sambuc.El 115*84d9c625SLionel Sambuc.Sh NOTES 1162fe8fb19SBen GrasIn 4.3 BSD, distributed from the University of California 1172fe8fb19SBen Grasin late 1985, most of the foregoing functions come in two 1182fe8fb19SBen Grasversions, one for the double\-precision "D" format in the 1192fe8fb19SBen GrasDEC VAX\-11 family of computers, another for double\-precision 1202fe8fb19SBen Grasarithmetic conforming to the IEEE Standard 754 for Binary 1212fe8fb19SBen GrasFloating\-Point Arithmetic. 1222fe8fb19SBen GrasThe two versions behave very 1232fe8fb19SBen Grassimilarly, as should be expected from programs more accurate 1242fe8fb19SBen Grasand robust than was the norm when UNIX was born. 1252fe8fb19SBen GrasFor instance, the programs are accurate to within the numbers 126*84d9c625SLionel Sambucof 127*84d9c625SLionel Sambuc.Dv ULPs 128*84d9c625SLionel Sambuctabulated above; an 129*84d9c625SLionel Sambuc.Dv ULP 130*84d9c625SLionel Sambucis one Unit in the Last Place. 1312fe8fb19SBen GrasAnd the programs have been cured of anomalies that 132*84d9c625SLionel Sambucafflicted the older math library 133*84d9c625SLionel Sambucin which incidents like 1342fe8fb19SBen Grasthe following had been reported: 135*84d9c625SLionel Sambuc.Bd -literal -offset indent 1362fe8fb19SBen Grassqrt(\-1.0) = 0.0 and log(\-1.0) = \-1.7e38. 1372fe8fb19SBen Grascos(1.0e\-11) \*[Gt] cos(0.0) \*[Gt] 1.0. 138*84d9c625SLionel Sambucpow(x,1.0) \(!= x when x = 2.0, 3.0, 4.0, ..., 9.0. 1392fe8fb19SBen Graspow(\-1.0,1.0e10) trapped on Integer Overflow. 1402fe8fb19SBen Grassqrt(1.0e30) and sqrt(1.0e\-30) were very slow. 141*84d9c625SLionel Sambuc.Ed 1422fe8fb19SBen GrasHowever the two versions do differ in ways that have to be 1432fe8fb19SBen Grasexplained, to which end the following notes are provided. 144*84d9c625SLionel Sambuc.Ss DEC VAX\-11 D_floating\-point 145*84d9c625SLionel SambucThis is the format for which the original math library 1462fe8fb19SBen Graswas developed, and to which this manual is still principally dedicated. 147*84d9c625SLionel SambucIt is 148*84d9c625SLionel Sambuc.Em the 149*84d9c625SLionel Sambucdouble\-precision format for the PDP\-11 1502fe8fb19SBen Grasand the earlier VAX\-11 machines; VAX\-11s after 1983 were 1512fe8fb19SBen Grasprovided with an optional "G" format closer to the IEEE 1522fe8fb19SBen Grasdouble\-precision format. 1532fe8fb19SBen GrasThe earlier DEC MicroVAXs have no D format, only G double\-precision. 1542fe8fb19SBen Gras(Why? 1552fe8fb19SBen GrasWhy not?) 156*84d9c625SLionel Sambuc.Pp 1572fe8fb19SBen GrasProperties of D_floating\-point: 158*84d9c625SLionel Sambuc.Bl -hang -offset indent 159*84d9c625SLionel Sambuc.It Wordsize : 160*84d9c625SLionel Sambuc64 bits, 8 bytes. 161*84d9c625SLionel Sambuc.It Radix : 162*84d9c625SLionel SambucBinary. 163*84d9c625SLionel Sambuc.It Precision : 164*84d9c625SLionel Sambuc56 significant bits, roughly like 17 significant decimals. 1652fe8fb19SBen GrasIf x and x' are consecutive positive D_floating\-point 166*84d9c625SLionel Sambucnumbers (they differ by 1 167*84d9c625SLionel Sambuc.Dv ULP ) , 168*84d9c625SLionel Sambucthen 169*84d9c625SLionel Sambuc.Dl 1.3e\-17 \*[Lt] 0.5**56 \*[Lt] (x'\-x)/x \*[Le] 0.5**55 \*[Lt] 2.8e\-17. 170*84d9c625SLionel Sambuc.It Range : 171*84d9c625SLionel Sambuc.Bl -column "Underflow thresholdX" "2.0**127X" 172*84d9c625SLionel Sambuc.It Overflow threshold = 2.0**127 = 1.7e38. 173*84d9c625SLionel Sambuc.It Underflow threshold = 0.5**128 = 2.9e\-39. 174*84d9c625SLionel Sambuc.El 175*84d9c625SLionel Sambuc.Em NOTE: THIS RANGE IS COMPARATIVELY NARROW. 176*84d9c625SLionel Sambuc.Pp 1772fe8fb19SBen GrasOverflow customarily stops computation. 1782fe8fb19SBen GrasUnderflow is customarily flushed quietly to zero. 179*84d9c625SLionel Sambuc.Em CAUTION : 1802fe8fb19SBen GrasIt is possible to have x 1812fe8fb19SBen Gras\(!= 182*84d9c625SLionel Sambucy and yet x\-y = 0 because of underflow. 1832fe8fb19SBen GrasSimilarly x \*[Gt] y \*[Gt] 0 cannot prevent either x\(**y = 0 1842fe8fb19SBen Grasor y/x = 0 from happening without warning. 185*84d9c625SLionel Sambuc.It Zero is represented ambiguously : 1862fe8fb19SBen GrasAlthough 2**55 different representations of zero are accepted by 1872fe8fb19SBen Grasthe hardware, only the obvious representation is ever produced. 1882fe8fb19SBen GrasThere is no \-0 on a VAX. 189*84d9c625SLionel Sambuc.It \*(If is not part of the VAX architecture . 190*84d9c625SLionel Sambuc.It Reserved operands : 1912fe8fb19SBen Grasof the 2**55 that the hardware 1922fe8fb19SBen Grasrecognizes, only one of them is ever produced. 1932fe8fb19SBen GrasAny floating\-point operation upon a reserved 1942fe8fb19SBen Grasoperand, even a MOVF or MOVD, customarily stops 1952fe8fb19SBen Grascomputation, so they are not much used. 196*84d9c625SLionel Sambuc.It Exceptions : 1972fe8fb19SBen GrasDivisions by zero and operations that 1982fe8fb19SBen Grasoverflow are invalid operations that customarily 1992fe8fb19SBen Grasstop computation or, in earlier machines, produce 2002fe8fb19SBen Grasreserved operands that will stop computation. 201*84d9c625SLionel Sambuc.It Rounding : 2022fe8fb19SBen GrasEvery rational operation (+, \-, \(**, /) on a 2032fe8fb19SBen GrasVAX (but not necessarily on a PDP\-11), if not an 2042fe8fb19SBen Grasover/underflow nor division by zero, is rounded to 205*84d9c625SLionel Sambucwithin half an 206*84d9c625SLionel Sambuc.Dv ULP , 207*84d9c625SLionel Sambucand when the rounding error is 208*84d9c625SLionel Sambucexactly half an 209*84d9c625SLionel Sambuc.Dv ULP 210*84d9c625SLionel Sambucthen rounding is away from 0. 211*84d9c625SLionel Sambuc.El 212*84d9c625SLionel Sambuc.Pp 2132fe8fb19SBen GrasExcept for its narrow range, D_floating\-point is one of the 2142fe8fb19SBen Grasbetter computer arithmetics designed in the 1960's. 2152fe8fb19SBen GrasIts properties are reflected fairly faithfully in the elementary 2162fe8fb19SBen Grasfunctions for a VAX distributed in 4.3 BSD. 2172fe8fb19SBen GrasThey over/underflow only if their results have to lie out of range 2182fe8fb19SBen Grasor very nearly so, and then they behave much as any rational 2192fe8fb19SBen Grasarithmetic operation that over/underflowed would behave. 2202fe8fb19SBen GrasSimilarly, expressions like log(0) and atanh(1) behave 2212fe8fb19SBen Graslike 1/0; and sqrt(\-3) and acos(3) behave like 0/0; 2222fe8fb19SBen Grasthey all produce reserved operands and/or stop computation! 2232fe8fb19SBen GrasThe situation is described in more detail in manual pages. 224*84d9c625SLionel Sambuc.Pp 225*84d9c625SLionel Sambuc.Em This response seems excessively punitive, so it is destined 226*84d9c625SLionel Sambuc.Em to be replaced at some time in the foreseeable future by a 227*84d9c625SLionel Sambuc.Em more flexible but still uniform scheme being developed to 228*84d9c625SLionel Sambuc.Em handle all floating\-point arithmetic exceptions neatly. 229*84d9c625SLionel Sambuc.Pp 230*84d9c625SLionel SambucHow do the functions in 4.3 BSD's new math library for UNIX 2312fe8fb19SBen Grascompare with their counterparts in DEC's VAX/VMS library? 2322fe8fb19SBen GrasSome of the VMS functions are a little faster, some are 2332fe8fb19SBen Grasa little more accurate, some are more puritanical about 2342fe8fb19SBen Grasexceptions (like pow(0.0,0.0) and atan2(0.0,0.0)), 2352fe8fb19SBen Grasand most occupy much more memory than their counterparts in 236*84d9c625SLionel Sambuclibm. 2372fe8fb19SBen GrasThe VMS codes interpolate in large table to achieve 238*84d9c625SLionel Sambucspeed and accuracy; the libm codes use tricky formulas 2392fe8fb19SBen Grascompact enough that all of them may some day fit into a ROM. 240*84d9c625SLionel Sambuc.Pp 2412fe8fb19SBen GrasMore important, DEC regards the VMS codes as proprietary 2422fe8fb19SBen Grasand guards them zealously against unauthorized use. 243*84d9c625SLionel SambucBut the libm codes in 4.3 BSD are intended for the public domain; 2442fe8fb19SBen Grasthey may be copied freely provided their provenance is always 2452fe8fb19SBen Grasacknowledged, and provided users assist the authors in their 2462fe8fb19SBen Grasresearches by reporting experience with the codes. 2472fe8fb19SBen GrasTherefore no user of UNIX on a machine whose arithmetic resembles 248*84d9c625SLionel SambucVAX D_floating\-point need use anything worse than the new libm. 249*84d9c625SLionel Sambuc.Ss IEEE STANDARD 754 Floating\-Point Arithmetic 2502fe8fb19SBen GrasThis standard is on its way to becoming more widely adopted 2512fe8fb19SBen Grasthan any other design for computer arithmetic. 2522fe8fb19SBen GrasVLSI chips that conform to some version of that standard have been 2532fe8fb19SBen Grasproduced by a host of manufacturers, among them ... 254*84d9c625SLionel Sambuc.Bl -column "Intel i8070, i80287XX" 255*84d9c625SLionel Sambuc.It Intel i8087, i80287 National Semiconductor 32081 256*84d9c625SLionel Sambuc.It 68881 Weitek WTL-1032, ... , -1165 257*84d9c625SLionel Sambuc.It Zilog Z8070 Western Electric (AT\*[Am]T) WE32106. 258*84d9c625SLionel Sambuc.El 2592fe8fb19SBen GrasOther implementations range from software, done thoroughly 2602fe8fb19SBen Grasin the Apple Macintosh, through VLSI in the Hewlett\-Packard 2612fe8fb19SBen Gras9000 series, to the ELXSI 6400 running ECL at 3 Megaflops. 2622fe8fb19SBen GrasSeveral other companies have adopted the formats 2632fe8fb19SBen Grasof IEEE 754 without, alas, adhering to the standard's way 2642fe8fb19SBen Grasof handling rounding and exceptions like over/underflow. 2652fe8fb19SBen GrasThe DEC VAX G_floating\-point format is very similar to the IEEE 2662fe8fb19SBen Gras754 Double format, so similar that the C programs for the 2672fe8fb19SBen GrasIEEE versions of most of the elementary functions listed 2682fe8fb19SBen Grasabove could easily be converted to run on a MicroVAX, though 2692fe8fb19SBen Grasnobody has volunteered to do that yet. 270*84d9c625SLionel Sambuc.Pp 271*84d9c625SLionel SambucThe codes in 4.3 BSD's libm for machines that conform to 272*84d9c625SLionel SambucIEEE 754 are intended primarily for the National Semiconductor 32081 2732fe8fb19SBen Grasand WTL 1164/65. 2742fe8fb19SBen GrasTo use these codes with the Intel or Zilog 2752fe8fb19SBen Graschips, or with the Apple Macintosh or ELXSI 6400, is to 2762fe8fb19SBen Grasforego the use of better codes provided (perhaps freely) by 2772fe8fb19SBen Grasthose companies and designed by some of the authors of the 2782fe8fb19SBen Grascodes above. 279*84d9c625SLionel SambucExcept for 280*84d9c625SLionel Sambuc.Fn atan , 281*84d9c625SLionel Sambuc.Fn cbrt , 282*84d9c625SLionel Sambuc.Fn erf , 283*84d9c625SLionel Sambuc.Fn erfc , 284*84d9c625SLionel Sambuc.Fn hypot , 285*84d9c625SLionel Sambuc.Fn j0-jn , 286*84d9c625SLionel Sambuc.Fn lgamma , 287*84d9c625SLionel Sambuc.Fn pow , 288*84d9c625SLionel Sambucand 289*84d9c625SLionel Sambuc.Fn y0\-yn , 290*84d9c625SLionel Sambucthe Motorola 68881 has all the functions in libm on chip, 2912fe8fb19SBen Grasand faster and more accurate; 292*84d9c625SLionel Sambucit, Apple, the i8087, Z8070 and WE32106 all use 64 significant bits. 2932fe8fb19SBen GrasThe main virtue of 4.3 BSD's 294*84d9c625SLionel Sambuclibm codes is that they are intended for the public domain; 2952fe8fb19SBen Grasthey may be copied freely provided their provenance is always 2962fe8fb19SBen Grasacknowledged, and provided users assist the authors in their 2972fe8fb19SBen Grasresearches by reporting experience with the codes. 2982fe8fb19SBen GrasTherefore no user of UNIX on a machine that conforms to 299*84d9c625SLionel SambucIEEE 754 need use anything worse than the new libm. 300*84d9c625SLionel Sambuc.Pp 3012fe8fb19SBen GrasProperties of IEEE 754 Double\-Precision: 302*84d9c625SLionel Sambuc.Bl -hang -offset indent 303*84d9c625SLionel Sambuc.It Wordsize : 304*84d9c625SLionel Sambuc64 bits, 8 bytes. 305*84d9c625SLionel Sambuc.It Radix : 306*84d9c625SLionel SambucBinary. 307*84d9c625SLionel Sambuc.It Precision : 308*84d9c625SLionel Sambuc53 significant bits, roughly like 16 significant decimals. 3092fe8fb19SBen GrasIf x and x' are consecutive positive Double\-Precision 310*84d9c625SLionel Sambucnumbers (they differ by 1 311*84d9c625SLionel Sambuc.Dv ULP ) , 312*84d9c625SLionel Sambucthen 313*84d9c625SLionel Sambuc.Dl 1.1e\-16 \*[Lt] 0.5**53 \*[Lt] (x'\-x)/x \*[Le] 0.5**52 \*[Lt] 2.3e\-16. 314*84d9c625SLionel Sambuc.It Range : 315*84d9c625SLionel Sambuc.Bl -column "Underflow thresholdX" "2.0**1024X" 316*84d9c625SLionel Sambuc.It Overflow threshold = 2.0**1024 = 1.8e308 317*84d9c625SLionel Sambuc.It Underflow threshold = 0.5**1022 = 2.2e\-308 318*84d9c625SLionel Sambuc.El 319*84d9c625SLionel SambucOverflow goes by default to a signed \*(If. 320*84d9c625SLionel SambucUnderflow is 321*84d9c625SLionel Sambuc.Sy Gradual , 322*84d9c625SLionel Sambucrounding to the nearest 3232fe8fb19SBen Grasinteger multiple of 0.5**1074 = 4.9e\-324. 324*84d9c625SLionel Sambuc.It Zero is represented ambiguously as +0 or \-0: 3252fe8fb19SBen GrasIts sign transforms correctly through multiplication or 3262fe8fb19SBen Grasdivision, and is preserved by addition of zeros 3272fe8fb19SBen Graswith like signs; but x\-x yields +0 for every 3282fe8fb19SBen Grasfinite x. 3292fe8fb19SBen GrasThe only operations that reveal zero's 3302fe8fb19SBen Grassign are division by zero and copysign(x,\(+-0). 3312fe8fb19SBen GrasIn particular, comparison (x \*[Gt] y, x \*[Ge] y, etc.) 3322fe8fb19SBen Grascannot be affected by the sign of zero; but if 333*84d9c625SLionel Sambucfinite x = y then \*(If 3342fe8fb19SBen Gras\&= 1/(x\-y) 3352fe8fb19SBen Gras\(!= 3362fe8fb19SBen Gras\-1/(y\-x) = 337*84d9c625SLionel Sambuc\- \*(If . 338*84d9c625SLionel Sambuc.It \*(If is signed : 3392fe8fb19SBen Grasit persists when added to itself 3402fe8fb19SBen Grasor to any finite number. 3412fe8fb19SBen GrasIts sign transforms 3422fe8fb19SBen Grascorrectly through multiplication and division, and 343*84d9c625SLionel Sambuc\*(If (finite)/\(+- \0=\0\(+-0 3442fe8fb19SBen Gras(nonzero)/0 = 345*84d9c625SLionel Sambuc\(+- \*(If. 3462fe8fb19SBen GrasBut 3472fe8fb19SBen Gras\(if\-\(if, \(if\(**0 and \(if/\(if 3482fe8fb19SBen Grasare, like 0/0 and sqrt(\-3), 349*84d9c625SLionel Sambucinvalid operations that produce \*(Na. 350*84d9c625SLionel Sambuc.It Reserved operands : 3512fe8fb19SBen Grasthere are 2**53\-2 of them, all 352*84d9c625SLionel Sambuccalled \*(Na (Not A Number). 353*84d9c625SLionel SambucSome, called Signaling \*[Na]s, trap any floating\-point operation 3542fe8fb19SBen Grasperformed upon them; they are used to mark missing 3552fe8fb19SBen Grasor uninitialized values, or nonexistent elements of arrays. 356*84d9c625SLionel SambucThe rest are Quiet \*[Na]s; they are 3572fe8fb19SBen Grasthe default results of Invalid Operations, and 3582fe8fb19SBen Graspropagate through subsequent arithmetic operations. 3592fe8fb19SBen GrasIf x 3602fe8fb19SBen Gras\(!= 361*84d9c625SLionel Sambucx then x is \*(Na; every other predicate 362*84d9c625SLionel Sambuc(x \*[Gt] y, x = y, x \*[Lt] y, ...) is FALSE if \*(Na is involved. 363*84d9c625SLionel Sambuc.Pp 364*84d9c625SLionel Sambuc.Em NOTE : 365*84d9c625SLionel SambucTrichotomy is violated by \*(Na. 3662fe8fb19SBen GrasBesides being FALSE, predicates that entail ordered 3672fe8fb19SBen Grascomparison, rather than mere (in)equality, 368*84d9c625SLionel Sambucsignal Invalid Operation when \*(Na is involved. 369*84d9c625SLionel Sambuc.It Rounding : 3702fe8fb19SBen GrasEvery algebraic operation (+, \-, \(**, /, 3712fe8fb19SBen Gras\(sr) 372*84d9c625SLionel Sambucis rounded by default to within half an 373*84d9c625SLionel Sambuc.Dv ULP , 374*84d9c625SLionel Sambucand when the rounding error is exactly half an 375*84d9c625SLionel Sambuc.Dv ULP 376*84d9c625SLionel Sambucthen the rounded value's least significant bit is zero. 3772fe8fb19SBen GrasThis kind of rounding is usually the best kind, 3782fe8fb19SBen Grassometimes provably so; for instance, for every 3792fe8fb19SBen Grasx = 1.0, 2.0, 3.0, 4.0, ..., 2.0**52, we find 3802fe8fb19SBen Gras(x/3.0)\(**3.0 == x and (x/10.0)\(**10.0 == x and ... 3812fe8fb19SBen Grasdespite that both the quotients and the products 3822fe8fb19SBen Grashave been rounded. 3832fe8fb19SBen GrasOnly rounding like IEEE 754 can do that. 3842fe8fb19SBen GrasBut no single kind of rounding can be 3852fe8fb19SBen Grasproved best for every circumstance, so IEEE 754 3862fe8fb19SBen Grasprovides rounding towards zero or towards 387*84d9c625SLionel Sambuc+\*(If 3882fe8fb19SBen Grasor towards 389*84d9c625SLionel Sambuc\-\*(If 3902fe8fb19SBen Grasat the programmer's option. 3912fe8fb19SBen GrasAnd the same kinds of rounding are specified for 3922fe8fb19SBen GrasBinary\-Decimal Conversions, at least for magnitudes 3932fe8fb19SBen Grasbetween roughly 1.0e\-10 and 1.0e37. 394*84d9c625SLionel Sambuc.It Exceptions : 3952fe8fb19SBen GrasIEEE 754 recognizes five kinds of floating\-point exceptions, 3962fe8fb19SBen Graslisted below in declining order of probable importance. 397*84d9c625SLionel Sambuc.Bl -column "Invalid OperationX" "Gradual OverflowX" 398*84d9c625SLionel Sambuc.It Sy Exception Ta Sy Default Result 399*84d9c625SLionel Sambuc.It Invalid Operation \*(Na, or FALSE 400*84d9c625SLionel Sambuc.It Overflow \(+-\(if 401*84d9c625SLionel Sambuc.It Divide by Zero \(+-\(if \} 402*84d9c625SLionel Sambuc.It Underflow Gradual Underflow 403*84d9c625SLionel Sambuc.It Inexact Rounded value 404*84d9c625SLionel Sambuc.El 405*84d9c625SLionel Sambuc.Pp 406*84d9c625SLionel Sambuc.Em NOTE : 407*84d9c625SLionel SambucAn Exception is not an Error unless handled badly. 4082fe8fb19SBen GrasWhat makes a class of exceptions exceptional 4092fe8fb19SBen Grasis that no single default response can be satisfactory 4102fe8fb19SBen Grasin every instance. 4112fe8fb19SBen GrasOn the other hand, if a default 4122fe8fb19SBen Grasresponse will serve most instances satisfactorily, 4132fe8fb19SBen Grasthe unsatisfactory instances cannot justify aborting 4142fe8fb19SBen Grascomputation every time the exception occurs. 415*84d9c625SLionel Sambuc.El 416*84d9c625SLionel Sambuc.Pp 4172fe8fb19SBen GrasFor each kind of floating\-point exception, IEEE 754 4182fe8fb19SBen Grasprovides a Flag that is raised each time its exception 4192fe8fb19SBen Grasis signaled, and stays raised until the program resets it. 4202fe8fb19SBen GrasPrograms may also test, save and restore a flag. 4212fe8fb19SBen GrasThus, IEEE 754 provides three ways by which programs 4222fe8fb19SBen Grasmay cope with exceptions for which the default result 4232fe8fb19SBen Grasmight be unsatisfactory: 424*84d9c625SLionel Sambuc.Bl -enum 425*84d9c625SLionel Sambuc.It 4262fe8fb19SBen GrasTest for a condition that might cause an exception 4272fe8fb19SBen Graslater, and branch to avoid the exception. 428*84d9c625SLionel Sambuc.It 4292fe8fb19SBen GrasTest a flag to see whether an exception has occurred 4302fe8fb19SBen Grassince the program last reset its flag. 431*84d9c625SLionel Sambuc.It 4322fe8fb19SBen GrasTest a result to see whether it is a value that only 4332fe8fb19SBen Grasan exception could have produced. 434*84d9c625SLionel Sambuc.Em CAUTION : 435*84d9c625SLionel SambucThe only reliable ways to discover 4362fe8fb19SBen Graswhether Underflow has occurred are to test whether 4372fe8fb19SBen Grasproducts or quotients lie closer to zero than the 4382fe8fb19SBen Grasunderflow threshold, or to test the Underflow flag. 4392fe8fb19SBen Gras(Sums and differences cannot underflow in 4402fe8fb19SBen GrasIEEE 754; if x 4412fe8fb19SBen Gras\(!= 4422fe8fb19SBen Grasy then x\-y is correct to 4432fe8fb19SBen Grasfull precision and certainly nonzero regardless of 4442fe8fb19SBen Grashow tiny it may be.) 4452fe8fb19SBen GrasProducts and quotients that 4462fe8fb19SBen Grasunderflow gradually can lose accuracy gradually 4472fe8fb19SBen Graswithout vanishing, so comparing them with zero 4482fe8fb19SBen Gras(as one might on a VAX) will not reveal the loss. 4492fe8fb19SBen GrasFortunately, if a gradually underflowed value is 4502fe8fb19SBen Grasdestined to be added to something bigger than the 4512fe8fb19SBen Grasunderflow threshold, as is almost always the case, 4522fe8fb19SBen Grasdigits lost to gradual underflow will not be missed 4532fe8fb19SBen Grasbecause they would have been rounded off anyway. 454*84d9c625SLionel SambucSo gradual underflows are usually 455*84d9c625SLionel Sambuc.Em provably 456*84d9c625SLionel Sambucignorable. 4572fe8fb19SBen GrasThe same cannot be said of underflows flushed to 0. 458*84d9c625SLionel Sambuc.Pp 4592fe8fb19SBen GrasAt the option of an implementor conforming to IEEE 754, 4602fe8fb19SBen Grasother ways to cope with exceptions may be provided: 461*84d9c625SLionel Sambuc.It 4622fe8fb19SBen GrasABORT. 4632fe8fb19SBen GrasThis mechanism classifies an exception in 4642fe8fb19SBen Grasadvance as an incident to be handled by means 4652fe8fb19SBen Grastraditionally associated with error\-handling 4662fe8fb19SBen Grasstatements like "ON ERROR GO TO ...". 4672fe8fb19SBen GrasDifferent languages offer different forms of this statement, 4682fe8fb19SBen Grasbut most share the following characteristics: 469*84d9c625SLionel Sambuc.Bl -dash 470*84d9c625SLionel Sambuc.It 4712fe8fb19SBen GrasNo means is provided to substitute a value for 4722fe8fb19SBen Grasthe offending operation's result and resume 4732fe8fb19SBen Grascomputation from what may be the middle of an expression. 4742fe8fb19SBen GrasAn exceptional result is abandoned. 475*84d9c625SLionel Sambuc.It 4762fe8fb19SBen GrasIn a subprogram that lacks an error\-handling 4772fe8fb19SBen Grasstatement, an exception causes the subprogram to 4782fe8fb19SBen Grasabort within whatever program called it, and so 4792fe8fb19SBen Grason back up the chain of calling subprograms until 4802fe8fb19SBen Grasan error\-handling statement is encountered or the 4812fe8fb19SBen Graswhole task is aborted and memory is dumped. 482*84d9c625SLionel Sambuc.El 483*84d9c625SLionel Sambuc.It 4842fe8fb19SBen GrasSTOP. 4852fe8fb19SBen GrasThis mechanism, requiring an interactive 4862fe8fb19SBen Grasdebugging environment, is more for the programmer 4872fe8fb19SBen Grasthan the program. 4882fe8fb19SBen GrasIt classifies an exception in 4892fe8fb19SBen Grasadvance as a symptom of a programmer's error; the 4902fe8fb19SBen Grasexception suspends execution as near as it can to 4912fe8fb19SBen Grasthe offending operation so that the programmer can 4922fe8fb19SBen Graslook around to see how it happened. 4932fe8fb19SBen GrasQuite often 4942fe8fb19SBen Grasthe first several exceptions turn out to be quite 4952fe8fb19SBen Grasunexceptionable, so the programmer ought ideally 4962fe8fb19SBen Grasto be able to resume execution after each one as if 4972fe8fb19SBen Grasexecution had not been stopped. 498*84d9c625SLionel Sambuc.It 4992fe8fb19SBen Gras\&... Other ways lie beyond the scope of this document. 500*84d9c625SLionel Sambuc.El 501*84d9c625SLionel Sambuc.Pp 5022fe8fb19SBen GrasThe crucial problem for exception handling is the problem of 5032fe8fb19SBen GrasScope, and the problem's solution is understood, but not 5042fe8fb19SBen Grasenough manpower was available to implement it fully in time 505*84d9c625SLionel Sambucto be distributed in 4.3 BSD's libm. 5062fe8fb19SBen GrasIdeally, each elementary function should act 5072fe8fb19SBen Grasas if it were indivisible, or atomic, in the sense that ... 508*84d9c625SLionel Sambuc.Bl -enum 509*84d9c625SLionel Sambuc.It 5102fe8fb19SBen GrasNo exception should be signaled that is not deserved by 5112fe8fb19SBen Grasthe data supplied to that function. 512*84d9c625SLionel Sambuc.It 5132fe8fb19SBen GrasAny exception signaled should be identified with that 5142fe8fb19SBen Grasfunction rather than with one of its subroutines. 515*84d9c625SLionel Sambuc.It 5162fe8fb19SBen GrasThe internal behavior of an atomic function should not 5172fe8fb19SBen Grasbe disrupted when a calling program changes from 5182fe8fb19SBen Grasone to another of the five or so ways of handling 5192fe8fb19SBen Grasexceptions listed above, although the definition 5202fe8fb19SBen Grasof the function may be correlated intentionally 5212fe8fb19SBen Graswith exception handling. 522*84d9c625SLionel Sambuc.El 523*84d9c625SLionel Sambuc.Pp 524*84d9c625SLionel SambucIdeally, every programmer should be able 525*84d9c625SLionel Sambuc.Em conveniently 526*84d9c625SLionel Sambucto turn a debugged subprogram into one that appears atomic to 5272fe8fb19SBen Grasits users. 5282fe8fb19SBen GrasBut simulating all three characteristics of an 5292fe8fb19SBen Grasatomic function is still a tedious affair, entailing hosts 5302fe8fb19SBen Grasof tests and saves\-restores; work is under way to ameliorate 5312fe8fb19SBen Grasthe inconvenience. 532*84d9c625SLionel Sambuc.Pp 533*84d9c625SLionel SambucMeanwhile, the functions in libm are only approximately atomic. 5342fe8fb19SBen GrasThey signal no inappropriate exception except possibly ... 535*84d9c625SLionel Sambuc.Bl -ohang -offset indent 536*84d9c625SLionel Sambuc.It Over/Underflow 5372fe8fb19SBen Graswhen a result, if properly computed, might have lain barely within range, and 538*84d9c625SLionel Sambuc.It Inexact in Fn cbrt , Fn hypot , Fn log10 No and Fn pow 5392fe8fb19SBen Graswhen it happens to be exact, thanks to fortuitous cancellation of errors. 540*84d9c625SLionel Sambuc.El 5412fe8fb19SBen GrasOtherwise, ... 542*84d9c625SLionel Sambuc.Bl -ohang -offset indent 543*84d9c625SLionel Sambuc.It Invalid Operation is signaled only when 544*84d9c625SLionel Sambucany result but \*(Na would probably be misleading. 545*84d9c625SLionel Sambuc.It Overflow is signaled only when 5462fe8fb19SBen Grasthe exact result would be finite but beyond the overflow threshold. 547*84d9c625SLionel Sambuc.It Divide\-by\-Zero is signaled only when 5482fe8fb19SBen Grasa function takes exactly infinite values at finite operands. 549*84d9c625SLionel Sambuc.It Underflow is signaled only when 5502fe8fb19SBen Grasthe exact result would be nonzero but tinier than the underflow threshold. 551*84d9c625SLionel Sambuc.It Inexact is signaled only when 5522fe8fb19SBen Grasgreater range or precision would be needed to represent the exact result. 553*84d9c625SLionel Sambuc.El 5542fe8fb19SBen Gras.\" .Sh FILES 5552fe8fb19SBen Gras.\" .Bl -tag -width /usr/lib/libm_p.a -compact 5562fe8fb19SBen Gras.\" .It Pa /usr/lib/libm.a 5572fe8fb19SBen Gras.\" the static math library 5582fe8fb19SBen Gras.\" .It Pa /usr/lib/libm.so 5592fe8fb19SBen Gras.\" the dynamic math library 5602fe8fb19SBen Gras.\" .It Pa /usr/lib/libm_p.a 5612fe8fb19SBen Gras.\" the static math library compiled for profiling 5622fe8fb19SBen Gras.\" .El 563*84d9c625SLionel Sambuc.Sh SEE ALSO 5642fe8fb19SBen GrasAn explanation of IEEE 754 and its proposed extension p854 5652fe8fb19SBen Graswas published in the IEEE magazine MICRO in August 1984 under 5662fe8fb19SBen Grasthe title "A Proposed Radix\- and Word\-length\-independent 5672fe8fb19SBen GrasStandard for Floating\-point Arithmetic" by W. J. Cody et al. 5682fe8fb19SBen GrasThe manuals for Pascal, C and BASIC on the Apple Macintosh 5692fe8fb19SBen Grasdocument the features of IEEE 754 pretty well. 5702fe8fb19SBen GrasArticles in the IEEE magazine COMPUTER vol. 14 no. 3 (Mar. 1981), 5712fe8fb19SBen Grasand in the ACM SIGNUM Newsletter Special Issue of 5722fe8fb19SBen GrasOct. 1979, may be helpful although they pertain to 5732fe8fb19SBen Grassuperseded drafts of the standard. 574*84d9c625SLionel Sambuc.Sh BUGS 5752fe8fb19SBen GrasWhen signals are appropriate, they are emitted by certain 5762fe8fb19SBen Grasoperations within the codes, so a subroutine\-trace may be 5772fe8fb19SBen Grasneeded to identify the function with its signal in case 5782fe8fb19SBen Grasmethod 5) above is in use. 5792fe8fb19SBen GrasAnd the codes all take the 5802fe8fb19SBen GrasIEEE 754 defaults for granted; this means that a decision to 5812fe8fb19SBen Grastrap all divisions by zero could disrupt a code that would 5822fe8fb19SBen Grasotherwise get correct results despite division by zero. 583