Switch to SpeakEasy.net DSL

The Modular Manual Browser

Home Page
Manual: (4.4BSD-Lite2)
Apropos / Subsearch:
optional field

MATH(3)              BSD Programmer's Manual              MATH(3)

       math - introduction to mathematical library functions

       These  functions constitute the C math library, libm.  The
       link editor searches this library under the "-lm"  option.
       Declarations  for these functions may be obtained from the
       include  file  <math.h>.   The  Fortran  math  library  is
       described in ``man 3f intro''.

       Name      Appears on Page    Description               Error Bound (ULPs)
       acos        sin.3m       inverse trigonometric function      3
       acosh       asinh.3m     inverse hyperbolic function         3
       asin        sin.3m       inverse trigonometric function      3
       asinh       asinh.3m     inverse hyperbolic function         3
       atan        sin.3m       inverse trigonometric function      1
       atanh       asinh.3m     inverse hyperbolic function         3
       atan2       sin.3m       inverse trigonometric function      2
       cabs        hypot.3m     complex absolute value              1
       cbrt        sqrt.3m      cube root                           1
       ceil        floor.3m     integer no less than                0
       copysign    ieee.3m      copy sign bit                       0
       cos         sin.3m       trigonometric function              1
       cosh        sinh.3m      hyperbolic function                 3
       drem        ieee.3m      remainder                           0
       erf         erf.3m       error function                     ???
       erfc        erf.3m       complementary error function       ???
       exp         exp.3m       exponential                         1
       expm1       exp.3m       exp(x)-1                            1
       fabs        floor.3m     absolute value                      0
       floor       floor.3m     integer no greater than             0
       hypot       hypot.3m     Euclidean distance                  1
       infnan      infnan.3m    signals exceptions
       j0          j0.3m        bessel function                    ???
       j1          j0.3m        bessel function                    ???
       jn          j0.3m        bessel function                    ???
       lgamma      lgamma.3m    log gamma function; (formerly gamma.3m)
       log         exp.3m       natural logarithm                   1
       logb        ieee.3m      exponent extraction                 0
       log10       exp.3m       logarithm to base 10                3
       log1p       exp.3m       log(1+x)                            1
       pow         exp.3m       exponential x**y                 60-500
       rint        floor.3m     round to nearest integer            0
       scalb       ieee.3m      exponent adjustment                 0
       sin         sin.3m       trigonometric function              1
       sinh        sinh.3m      hyperbolic function                 3
       sqrt        sqrt.3m      square root                         1
       tan         sin.3m       trigonometric function              3
       tanh        sinh.3m      hyperbolic function                 3
       y0          j0.3m        bessel function                    ???

4th Berkeley Distribution  May 5, 1994                          1

MATH(3)              BSD Programmer's Manual              MATH(3)

       y1          j0.3m        bessel function                    ???
       yn          j0.3m        bessel function                    ???

       In  4.3 BSD, distributed from the University of California
       in late 1985, most of the foregoing functions come in  two
       versions,  one  for the double-precision "D" format in the
       DEC  VAX-11  family  of  computers,   another   for   dou-
       ble-precision  arithmetic  conforming to the IEEE Standard
       754 for Binary Floating-Point Arithmetic.   The  two  ver-
       sions  behave  very  similarly, as should be expected from
       programs more accurate and robust than was the  norm  when
       UNIX was born.  For instance, the programs are accurate to
       within the numbers of ulps tabulated above; an ulp is  one
       Unit  in the Last Place.  And the programs have been cured
       of anomalies that afflicted the older math library libm in
       which incidents like the following had been reported:
              sqrt(-1.0) = 0.0 and log(-1.0) = -1.7e38.
              cos(1.0e-11) > cos(0.0) > 1.0.
              pow(x,1.0) != x when x = 2.0, 3.0, 4.0, ..., 9.0.
              pow(-1.0,1.0e10) trapped on Integer Overflow.
              sqrt(1.0e30) and sqrt(1.0e-30) were very slow.
       However the two versions do differ in ways that have to be
       explained, to which end the following notes are  provided.

       DEC VAX-11 D_floating-point:

       This  is  the  format  for which the original math library
       libm was developed, and to  which  this  manual  is  still
       principally  dedicated.  It is the double-precision format
       for the PDP-11 and the earlier  VAX-11  machines;  VAX-11s
       after  1983  were  provided  with  an  optional "G" format
       closer to the IEEE double-precision format.   The  earlier
       DEC  MicroVAXs  have no D format, only G double-precision.
       (Why?  Why not?)

       Properties of D_floating-point:
              Wordsize: 64 bits, 8 bytes.  Radix: Binary.
              Precision: 56 sig.   bits,  roughly  like  17  sig.
                     If   x   and  x'  are  consecutive  positive
                     D_floating-point numbers (they differ  by  1
                     ulp), then
                     1.3e-17  <  0.5**56  < (x'-x)/x <= 0.5**55 <
              Range: Overflow threshold  = 2.0**127 = 1.7e38.
                     Underflow threshold = 0.5**128 = 2.9e-39.
                     Overflow customarily stops computation.
                     Underflow is customarily flushed quietly  to

4th Berkeley Distribution  May 5, 1994                          2

MATH(3)              BSD Programmer's Manual              MATH(3)

                             It  is  possible  to have x != y and
                             yet x-y = 0  because  of  underflow.
                             Similarly  x  > y > 0 cannot prevent
                             either x*y = 0 or  y/x = 0 from hap-
                             pening without warning.
              Zero is represented ambiguously.
                     Although  2**55 different representations of
                     zero are accepted by the hardware, only  the
                     obvious  representation  is  ever  produced.
                     There is no -0 on a VAX.
              Infinity is not part of the VAX architecture.
              Reserved operands:
                     of the 2**55 that the  hardware  recognizes,
                     only  one  of  them  is  ever produced.  Any
                     floating-point  operation  upon  a  reserved
                     operand,  even  a  MOVF or MOVD, customarily
                     stops computation,  so  they  are  not  much
                     Divisions  by zero and operations that over-
                     flow are invalid operations that customarily
                     stop  computation  or,  in earlier machines,
                     produce reserved  operands  that  will  stop
                     Every  rational operation  (+, -, *, /) on a
                     VAX (but not necessarily on  a  PDP-11),  if
                     not  an over/underflow nor division by zero,
                     is rounded to within half an ulp,  and  when
                     the  rounding  error  is exactly half an ulp
                     then rounding is away from 0.

       Except for its narrow range, D_floating-point  is  one  of
       the  better  computer  arithmetics designed in the 1960's.
       Its properties are reflected fairly faithfully in the ele-
       mentary  functions for a VAX distributed in 4.3 BSD.  They
       over/underflow only if their results have to  lie  out  of
       range  or very nearly so, and then they behave much as any
       rational arithmetic operation that over/underflowed  would
       behave.   Similarly,  expressions like log(0) and atanh(1)
       behave like 1/0; and sqrt(-3) and acos(3) behave like 0/0;
       they  all  produce  reserved operands and/or stop computa-
       tion!  The situation is described in more detail in manual
              This  response  seems excessively punitive, so
              it is destined to be replaced at some time  in
              the  foreseeable future by a more flexible but
              still uniform scheme being developed to handle
              all   floating-point   arithmetic   exceptions
              neatly.  See infnan(3M) for the present  state

4th Berkeley Distribution  May 5, 1994                          3

MATH(3)              BSD Programmer's Manual              MATH(3)

              of affairs.

       How  do  the functions in 4.3 BSD's new libm for UNIX com-
       pare with their counterparts  in  DEC's  VAX/VMS  library?
       Some  of the VMS functions are a little faster, some are a
       little more accurate,  some  are  more  puritanical  about
       exceptions  (like  pow(0.0,0.0)  and  atan2(0.0,0.0)), and
       most occupy much more memory than  their  counterparts  in
       libm.  The VMS codes interpolate in large table to achieve
       speed and accuracy; the libm  codes  use  tricky  formulas
       compact  enough  that  all of them may some day fit into a

       More important, DEC regards the VMS codes  as  proprietary
       and  guards  them zealously against unauthorized use.  But
       the libm codes in 4.3 BSD  are  intended  for  the  public
       domain;  they  may  be copied freely provided their prove-
       nance is always acknowledged, and  provided  users  assist
       the  authors  in  their researches by reporting experience
       with the codes.  Therefore no user of UNIX  on  a  machine
       whose  arithmetic  resembles VAX D_floating-point need use
       anything worse than the new libm.

       IEEE STANDARD 754 Floating-Point Arithmetic:

       This standard is  on  its  way  to  becoming  more  widely
       adopted  than  any  other  design for computer arithmetic.
       VLSI chips that conform to some version of  that  standard
       have  been produced by a host of manufacturers, among them
            Intel i8087, i80287      National Semiconductor  32081
            Motorola 68881           Weitek WTL-1032, ... , -1165
            Zilog Z8070              Western Electric (AT&T) WE32106.
       Other implementations range from software, done thoroughly
       in    the   Apple   Macintosh,   through   VLSI   in   the
       Hewlett-Packard 9000 series, to the ELXSI 6400 running ECL
       at  3 Megaflops.  Several other companies have adopted the
       formats of IEEE 754 without, alas, adhering to  the  stan-
       dard's  way  of  handling  rounding  and  exceptions  like
       over/underflow.  The DEC VAX  G_floating-point  format  is
       very  similar  to  the  IEEE 754 Double format, so similar
       that the C programs for the IEEE versions of most  of  the
       elementary  functions  listed  above  could easily be con-
       verted to run on a MicroVAX, though nobody has volunteered
       to do that yet.

       The  codes  in 4.3 BSD's libm for machines that conform to
       IEEE 754 are intended primarily  for  the  National  Semi.
       32081  and WTL 1164/65.  To use these codes with the Intel
       or Zilog chips, or with the Apple Macintosh or ELXSI 6400,
       is  to  forego  the  use of better codes provided (perhaps

4th Berkeley Distribution  May 5, 1994                          4

MATH(3)              BSD Programmer's Manual              MATH(3)

       freely) by those companies and designed  by  some  of  the
       authors  of the codes above.  Except for atan, cabs, cbrt,
       erf, erfc,  hypot,  j0-jn,  lgamma,  pow  and  y0-yn,  the
       Motorola  68881 has all the functions in libm on chip, and
       faster and more accurate; it, Apple, the i8087, Z8070  and
       WE32106  all  use  64  sig.  bits.  The main virtue of 4.3
       BSD's libm codes is that they are intended for the  public
       domain;  they  may  be copied freely provided their prove-
       nance is always acknowledged, and  provided  users  assist
       the  authors  in  their researches by reporting experience
       with the codes.  Therefore no user of UNIX  on  a  machine
       that conforms to IEEE 754 need use anything worse than the
       new libm.

       Properties of IEEE 754 Double-Precision:
              Wordsize: 64 bits, 8 bytes.  Radix: Binary.
              Precision: 53 sig.   bits,  roughly  like  16  sig.
                     If  x  and  x' are consecutive positive Dou-
                     ble-Precision  numbers  (they  differ  by  1
                     ulp), then
                     1.1e-16  <  0.5**53  < (x'-x)/x <= 0.5**52 <
              Range: Overflow threshold  = 2.0**1024 = 1.8e308
                     Underflow threshold = 0.5**1022 = 2.2e-308
                     Overflow goes by default to a signed  Infin-
                     Underflow  is Gradual, rounding to the near-
                     est  integer   multiple   of   0.5**1074   =
              Zero is represented ambiguously as +0 or -0.
                     Its sign transforms correctly through multi-
                     plication or division, and is  preserved  by
                     addition  of  zeros with like signs; but x-x
                     yields +0 for  every  finite  x.   The  only
                     operations that reveal zero's sign are divi-
                     sion by zero and copysign(x,+-0).   In  par-
                     ticular,  comparison  (x  > y, x >= y, etc.)
                     cannot be affected by the sign of zero;  but
                     if  finite  x = y then Infinity = 1/(x-y) !=
                     -1/(y-x) = -Infinity.
              Infinity is signed.
                     it persists when added to itself or  to  any
                     finite  number.   Its  sign  transforms cor-
                     rectly through multiplication and  division,
                     and  (finite)/+-Infinity = +-0 (nonzero)/0 =
                     +-Infinity.  But  Infinity-Infinity,  Infin-
                     ity*0  and  Infinity/Infinity  are, like 0/0
                     and sqrt(-3), invalid operations  that  pro-
                     duce NaN. ...
              Reserved operands:

4th Berkeley Distribution  May 5, 1994                          5

MATH(3)              BSD Programmer's Manual              MATH(3)

                     there  are  2**53-2  of them, all called NaN
                     (Not  a  Number).   Some,  called  Signaling
                     NaNs, trap any floating-point operation per-
                     formed upon them;  they  are  used  to  mark
                     missing or uninitialized values, or nonexis-
                     tent elements of arrays.  The rest are Quiet
                     NaNs;   they  are  the  default  results  of
                     Invalid Operations,  and  propagate  through
                     subsequent arithmetic operations.  If x != x
                     then x is NaN; every other predicate (x > y,
                     x  =  y,  x  <  y,  ...)  is FALSE if NaN is
                     NOTE: Trichotomy is violated by NaN.
                             Besides being FALSE, predicates that
                             entail  ordered  comparison,  rather
                             than   mere   (in)equality,   signal
                             Invalid   Operation   when   NaN  is
                     Every algebraic operation (+, -, *, /, sqrt)
                     is rounded by default to within half an ulp,
                     and when the rounding error is exactly  half
                     an  ulp  then the rounded value's least sig-
                     nificant bit is zero.  This kind of rounding
                     is usually the best kind, sometimes provably
                     so; for instance, for every x  =  1.0,  2.0,
                     3.0,  4.0, ..., 2.0**52, we find (x/3.0)*3.0
                     == x and (x/10.0)*10.0 == x and ...  despite
                     that  both  the  quotients  and the products
                     have been rounded.  Only rounding like  IEEE
                     754  can  do  that.   But  no single kind of
                     rounding can be proved best for  every  cir-
                     cumstance,  so  IEEE  754  provides rounding
                     towards zero or towards +Infinity or towards
                     -Infinity  at  the programmer's option.  And
                     the same kinds of rounding are specified for
                     Binary-Decimal  Conversions,  at  least  for
                     magnitudes  between  roughly   1.0e-10   and
                     IEEE  754  recognizes  five  kinds of float-
                     ing-point  exceptions,   listed   below   in
                     declining order of probable importance.
                             Exception              Default Result
                             Invalid Operation      NaN, or FALSE
                             Overflow               +-Infinity
                             Divide by Zero         +-Infinity
                             Underflow              Gradual Underflow
                             Inexact                Rounded value
                     NOTE:   An  Exception is not an Error unless
                     handled  badly.   What  makes  a  class   of

4th Berkeley Distribution  May 5, 1994                          6

MATH(3)              BSD Programmer's Manual              MATH(3)

                     exceptions  exceptional  is  that  no single
                     default  response  can  be  satisfactory  in
                     every  instance.   On  the  other hand, if a
                     default response will serve  most  instances
                     satisfactorily, the unsatisfactory instances
                     cannot justify  aborting  computation  every
                     time the exception occurs.

              For each kind of floating-point exception, IEEE 754
              provides a Flag that is raised each time its excep-
              tion  is  signaled, and stays raised until the pro-
              gram resets it.  Programs may also test,  save  and
              restore a flag.  Thus, IEEE 754 provides three ways
              by which programs  may  cope  with  exceptions  for
              which the default result might be unsatisfactory:

              1)  Test for a condition that might cause an excep-
                  tion later, and branch to avoid the  exception.

              2)  Test  a  flag  to  see whether an exception has
                  occurred since the program last reset its flag.

              3)  Test a result to see whether it is a value that
                  only an exception could have produced.
                  CAUTION: The only  reliable  ways  to  discover
                  whether  Underflow  has  occurred  are  to test
                  whether products or  quotients  lie  closer  to
                  zero  than  the underflow threshold, or to test
                  the Underflow flag.  (Sums and differences can-
                  not  underflow  in IEEE 754; if x != y then x-y
                  is correct  to  full  precision  and  certainly
                  nonzero  regardless  of  how  tiny  it may be.)
                  Products and quotients that underflow gradually
                  can  lose accuracy gradually without vanishing,
                  so comparing them with zero (as one might on  a
                  VAX) will not reveal the loss.  Fortunately, if
                  a gradually underflowed value is destined to be
                  added  to  something  bigger than the underflow
                  threshold, as is almost always the case, digits
                  lost  to  gradual  underflow will not be missed
                  because they would have been rounded  off  any-
                  way.   So  gradual underflows are usually prov-
                  ably ignorable.  The same  cannot  be  said  of
                  underflows flushed to 0.

              At  the option of an implementor conforming to IEEE
              754, other ways to cope with exceptions may be pro-

              4)  ABORT.   This mechanism classifies an exception
                  in advance as an  incident  to  be  handled  by

4th Berkeley Distribution  May 5, 1994                          7

MATH(3)              BSD Programmer's Manual              MATH(3)

                  means     traditionally     associated     with
                  error-handling statements like "ON ERROR GO  TO
                  ...".    Different  languages  offer  different
                  forms of this statement,  but  most  share  the
                  following characteristics:

              --  No  means is provided to substitute a value for
                  the offending  operation's  result  and  resume
                  computation  from  what may be the middle of an
                  expression.  An  exceptional  result  is  aban-

              --  In  a  subprogram  that lacks an error-handling
                  statement, an exception causes  the  subprogram
                  to abort within whatever program called it, and
                  so on back up the chain of calling  subprograms
                  until  an  error-handling  statement is encoun-
                  tered or the whole task is aborted  and  memory
                  is dumped.

              5)  STOP.  This mechanism, requiring an interactive
                  debugging environment, is more for the program-
                  mer  than the program.  It classifies an excep-
                  tion in advance as a symptom of a  programmer's
                  error; the exception suspends execution as near
                  as it can to the offending  operation  so  that
                  the  programmer  can  look around to see how it
                  happened.  Quite often the first several excep-
                  tions  turn out to be quite unexceptionable, so
                  the programmer ought  ideally  to  be  able  to
                  resume execution after each one as if execution
                  had not been stopped.

              6)  ... Other ways lie beyond  the  scope  of  this

       The  crucial problem for exception handling is the problem
       of Scope, and the problem's solution  is  understood,  but
       not enough manpower was available to implement it fully in
       time to be distributed in 4.3 BSD's libm.   Ideally,  each
       elementary  function should act as if it were indivisible,
       or atomic, in the sense that ...

       i)    No exception should be signaled that is not deserved
             by the data supplied to that function.

       ii)   Any  exception  signaled  should  be identified with
             that function rather than with one  of  its  subrou-

       iii)  The  internal  behavior of an atomic function should

4th Berkeley Distribution  May 5, 1994                          8

MATH(3)              BSD Programmer's Manual              MATH(3)

             not be disrupted when a calling program changes from
             one  to  another  of the five or so ways of handling
             exceptions listed above, although the definition  of
             the  function  may  be correlated intentionally with
             exception handling.

       Ideally, every programmer should be able  conveniently  to
       turn a debugged subprogram into one that appears atomic to
       its users.  But simulating all three characteristics of an
       atomic function is still a tedious affair, entailing hosts
       of tests and saves-restores; work is under way to  amelio-
       rate the inconvenience.

       Meanwhile,  the  functions  in libm are only approximately
       atomic.  They signal  no  inappropriate  exception  except
       possibly ...
                     when  a  result, if properly computed, might
                     have lain barely within range, and
              Inexact in cabs, cbrt, hypot, log10 and pow
                     when it happens to be exact, thanks to  for-
                     tuitous cancellation of errors.
       Otherwise, ...
              Invalid Operation is signaled only when
                     any  result  but  NaN would probably be mis-
              Overflow is signaled only when
                     the exact result would be finite but  beyond
                     the overflow threshold.
              Divide-by-Zero is signaled only when
                     a  function takes exactly infinite values at
                     finite operands.
              Underflow is signaled only when
                     the exact result would be nonzero but tinier
                     than the underflow threshold.
              Inexact is signaled only when
                     greater  range  or precision would be needed
                     to represent the exact result.

       When signals are appropriate, they are emitted by  certain
       operations  within the codes, so a subroutine-trace may be
       needed to identify the function with its  signal  in  case
       method  5)  above  is  in use.  And the codes all take the
       IEEE 754 defaults for granted; this means that a  decision
       to  trap  all  divisions by zero could disrupt a code that
       would otherwise get correct results  despite  division  by

       An explanation of IEEE 754 and its proposed extension p854

4th Berkeley Distribution  May 5, 1994                          9

MATH(3)              BSD Programmer's Manual              MATH(3)

       was published in the IEEE magazine MICRO  in  August  1984
       under     the     title    "A    Proposed    Radix-    and
       Word-length-independent Standard for Floating-point Arith-
       metic" by W. J. Cody et al.  The manuals for Pascal, C and
       BASIC on the Apple Macintosh document the features of IEEE
       754  pretty  well.  Articles in the IEEE magazine COMPUTER
       vol. 14 no. 3 (Mar.  1981), and in the ACM SIGNUM Newslet-
       ter  Special  Issue  of Oct. 1979, may be helpful although
       they pertain to superseded drafts of the standard.

4th Berkeley Distribution  May 5, 1994                         10