lib/libm/README

   1 /*
   2  * Copyright (c) 1985, 1993
   3  *      The Regents of the University of California.  All rights reserved.
   4  *
   5  * Redistribution and use in source and binary forms, with or without
   6  * modification, are permitted provided that the following conditions
   7  * are met:
   8  * 1. Redistributions of source code must retain the above copyright
   9  *    notice, this list of conditions and the following disclaimer.
  10  * 2. Redistributions in binary form must reproduce the above copyright
  11  *    notice, this list of conditions and the following disclaimer in the
  12  *    documentation and/or other materials provided with the distribution.
  13  * 3. All advertising materials mentioning features or use of this software
  14  *    must display the following acknowledgement:
  15  *      This product includes software developed by the University of
  16  *      California, Berkeley and its contributors.
  17  * 4. Neither the name of the University nor the names of its contributors
  18  *    may be used to endorse or promote products derived from this software
  19  *    without specific prior written permission.
  20  *
  21  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  22  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  23  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  24  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  25  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  26  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  27  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  28  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  29  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  30  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  31  * SUCH DAMAGE.
  32  *
  33  * K.C. Ng, with Z-S. Alex Liu, S. McDonald, P. Tang, W. Kahan.
  34  * Revised on 5/10/85, 5/13/85, 6/14/85, 8/20/85, 8/27/85, 9/11/85.
  35  *
  36  *      @(#)README      8.1 (Berkeley) 6/4/93
  37  */
  38
  39 ******************************************************************************
  40 *  This is a description of the upgraded elementary functions (listed in 1). *
  41 *  Bessel functions (j0, j1, jn, y0, y1, yn), floor, and fabs passed over    *
  42 *  from 4.2BSD without change except perhaps for the way floating point      *
  43 *  exception is signaled on a VAX.  Three lines that contain "errno" in erf.c*
  44 *  (error functions erf, erfc) have been deleted to prevent overriding the   *
  45 *  system "errno".                                                           *
  46 ******************************************************************************
  47
  48 0. Total number of files: 40
  49
  50         IEEE/Makefile   VAX/Makefile    VAX/support.s   erf.c       lgamma.c
  51         IEEE/atan2.c    VAX/argred.s    VAX/tan.s       exp.c       log.c
  52         IEEE/cabs.c     VAX/atan2.s     acosh.c         exp__E.c    log10.c
  53         IEEE/cbrt.c     VAX/cabs.s      asincos.c       expm1.c     log1p.c
  54         IEEE/support.c  VAX/cbrt.s      asinh.c         floor.c     log__L.c
  55         IEEE/trig.c     VAX/infnan.s    atan.c          j0.c        pow.c
  56         Makefile        VAX/sincos.s    atanh.c         j1.c        sinh.c
  57         README          VAX/sqrt.s      cosh.c          jn.c        tanh.c
  58
  59 1. Functions implemented :
  60     (A). Standard elementary functions (total 22) :
  61         acos(x)                 ...in file  asincos.c
  62         asin(x)                 ...in file  asincos.c
  63         atan(x)                 ...in file  atan.c
  64         atan2(x,y)              ...in files IEEE/atan2.c, VAX/atan2.s
  65         sin(x)                  ...in files IEEE/trig.c,  VAX/sincos.s
  66         cos(x)                  ...in files IEEE/trig.c,  VAX/sincos.s
  67         tan(x)                  ...in files IEEE/trig.c,  VAX/tan.s
  68         cabs(x,y)               ...in files IEEE/cabs.c,  VAX/cabs.s
  69         hypot(x,y)              ...in files IEEE/cabs.c,  VAX/cabs.s
  70         cbrt(x)                 ...in files IEEE/cbrt.c,  VAX/cbrt.s
  71         exp(x)                  ...in file  exp.c
  72         expm1(x):=exp(x)-1      ...in file  expm1.c
  73         log(x)                  ...in file  log.c
  74         log10(x)                ...in file  log10.c
  75         log1p(x):=log(1+x)      ...in file  log1p.c
  76         pow(x,y)                ...in file  pow.c
  77         sinh(x)                 ...in file  sinh.c
  78         cosh(x)                 ...in file  cosh.c
  79         tanh(x)                 ...in file  tanh.c
  80         asinh(x)                ...in file  asinh.c
  81         acosh(x)                ...in file  acosh.c
  82         atanh(x)                ...in file  atanh.c
  83
  84     (B). Kernel functions :
  85         exp__E(x,c) ...in file exp__E.c, used by expm1/exp/pow/cosh
  86         log__L(s)   ...in file log__L.c, used by log1p/log/pow
  87         libm$argred ...in file VAX/argred.s, used by VAX version of sin/cos/tan
  88
  89     (C). System supported functions :
  90         sqrt()      ...in files IEEE/support.c, VAX/sqrt.s
  91         drem()      ...in files IEEE/support.c, VAX/support.s
  92         finite()    ...in files IEEE/support.c, VAX/support.s
  93         logb()      ...in files IEEE/support.c, VAX/support.s
  94         scalb()     ...in files IEEE/support.c, VAX/support.s
  95         copysign()  ...in files IEEE/support.c, VAX/support.s
  96         rint()      ...in file  floor.c
  97
  98
  99    Notes:
 100        i. The codes in files ending with ".s" are written in VAX assembly
 101           language. They are intended for VAX computers.
 102
 103           Files that end with ".c" are written in C. They are intended
 104           for either a VAX or a machine that conforms to the IEEE
 105           standard 754 for double precision floating-point arithmetic.
 106
 107       ii. On other than VAX or IEEE machines, run the original math
 108           library, formerly "/usr/lib/libm.a", now "/usr/lib/libom.a", if
 109           nothing better is available.
 110
 111      iii. The trigonometric functions sin/cos/tan/atan2 in files "VAX/sincos.s",
 112           "VAX/tan.s" and "VAX/atan2.s" are different from those in
 113           "IEEE/trig.c" and "IEEE/atan2.c".  The VAX assembler code uses the
 114           true value of pi to perform argument reduction, while the C code uses
 115           a machine value of PI (see "IEEE/trig.c").
 116
 117
 118 2. A computer system that conforms to IEEE standard 754 should provide
 119                 sqrt(x),
 120                 drem(x,p), (double precision remainder function)
 121                 copysign(x,y),
 122                 finite(x),
 123                 scalb(x,N),
 124                 logb(x) and
 125                 rint(x).
 126    These functions are either required or recommended by the standard.
 127    For convenience, a (slow) C implementation of these functions is
 128    provided in the file "IEEE/support.c".
 129
 130    Warning: The functions in IEEE/support.c are somewhat machine dependent.
 131    Some modifications may be necessary to run them on a different machine.
 132    Currently, if compiled with a suitable flag, "IEEE/support.c" will work
 133    on a National 32000, a Zilog 8000, a VAX, and a SUN (cf. the "Makefile"
 134    in this directory). Invoke the C compiler thus:
 135
 136         cc -c -DVAX IEEE/support.c              ... on a VAX, D-format
 137         cc -c -DNATIONAL IEEE/support.c         ... on a National 32000
 138         cc -c  IEEE/support.c                   ... on other IEEE machines,
 139                                                     we hope.
 140
 141    Notes:
 142       1. Faster versions of "drem" and "sqrt" for IEEE double precision
 143          (coded in C but intended for assembly language) are given at the
 144          end of "IEEE/support.c" but commented out since they require certain
 145          machine-dependent functions.
 146
 147       2. A fast VAX assembler version of the system supported functions
 148          copysign(), logb(), scalb(), finite(), and drem() appears in file
 149          "VAX/support.s".  A fast VAX assembler version of sqrt() is in
 150          file "VAX/sqrt.s".
 151
 152 3. Two formats are supported by all the standard elementary functions:
 153    the VAX D-format (56-bit precision), and the IEEE double format
 154    (53-bit precision).  The cbrt() in "IEEE/cbrt.c" is for IEEE machines
 155    only. The functions in files that end with ".s" are for VAX computers
 156    only. The functions in files that end with ".c" (except "IEEE/cbrt.c")
 157    are for VAX and IEEE machines. To use the VAX D-format, compile the code
 158    with -DVAX; to use IEEE double format on various IEEE machines, see
 159    "Makefile" in this directory).
 160
 161     Example:
 162         cc -c -DVAX sin.c               ... for VAX D-format
 163
 164        Warning: The values of floating-point constants used in the code are
 165                 given in both hexadecimal and decimal.  The hexadecimal values
 166                 are the intended ones. The decimal values may be used provided
 167                 that the compiler converts from decimal to binary accurately
 168                 enough to produce the hexadecimal values shown. If the
 169                 conversion is inaccurate, then one must know the exact machine
 170                 representation of the constants and alter the assembly
 171                 language output from the compiler, or play tricks like
 172                 the following in a C program.
 173
 174                         Example: to store the floating-point constant
 175
 176                              p1= 2^-6 * .F83ABE67E1066A (Hexadecimal)
 177
 178                         on a VAX in C, we use two longwords to store its
 179                         machine value and define p1 to be the double constant
 180                         at the location of these two longwords:
 181
 182                         static long  p1x[] = { 0x3abe3d78, 0x066a67e1};
 183                         #define      p1      (*(double*)p1x)
 184
 185     Note:  On a VAX, some functions have two codes. For example, cabs() has
 186            one implementation in "IEEE/cabs.c", and another in "VAX/cabs.s".
 187            In this case, the assembly language version is preferred.
 188
 189
 190 4. Accuracy.
 191
 192             The errors in expm1(), log1p(), exp(), log(), cabs(), hypot()
 193             and cbrt() are below 1 ULP (Unit in the Last Place).
 194
 195             The error in pow(x,y) grows with the size of y. Nevertheless,
 196             for integers x and y, pow(x,y) returns the correct integer value
 197             on all tested machines (VAX, SUN, NATIONAL, ZILOG), provided that
 198             x to the power of y is representable exactly.
 199
 200             cosh, sinh, acosh, asinh, tanh, atanh and log10 have errors below
 201             about 3 ULPs.
 202
 203             For trigonometric and inverse trigonometric functions:
 204
 205                 Let [trig(x)] denote the value actually computed for trig(x),
 206
 207                 1) Those codes using the machine's value PI (true pi rounded):
 208                    (source codes: IEEE/{trig.c,atan2.c}, asincos.c and atan.c)
 209
 210                    The errors in [sin(x)], [cos(x)], and [atan(x)] are below
 211                    1 ULP compared with sin(x*pi/PI), cos(x*pi/PI), and
 212                    atan(x)*PI/pi respectively, where PI is the machine's
 213                    value of pi rounded. [tan(x)] returns tan(x*pi/PI) within
 214                    about 2 ULPs; [acos(x)], [asin(x)], and [atan2(y,x)]
 215                    return acos(x)*PI/pi, asin(x)*PI/pi, and atan2(y,x)*PI/pi
 216                    respectively to similar accuracy.
 217
 218
 219                 2) Those using true pi (for VAX D-format only):
 220                    (source codes: VAX/{sincos.s,tan.s,atan2.s}, asincos.c and
 221                    atan.c)
 222
 223                    The errors in [sin(x)], [cos(x)], and [atan(x)] are below
 224                    1 ULP. [tan(x)], [atan2(y,x)], [acos(x)], and [asin(x)]
 225                    have errors below about 2 ULPs.
 226
 227
 228             Here are the results of some test runs to find worst errors on
 229             the VAX :
 230
 231     tan   :  2.09 ULPs          ...1,024,000 random arguments (machine PI)
 232     sin   :  .861 ULPs          ...1,024,000 random arguments (machine PI)
 233     cos   :  .857 ULPs          ...1,024,000 random arguments (machine PI)
 234     (compared with tan, sin, cos of (x*pi/PI))
 235
 236     acos  :  2.07 ULPs          .....200,000 random arguments (machine PI)
 237     asin  :  2.06 ULPs          .....200,000 random arguments (machine PI)
 238     atan2 :  1.41 ULPs          .....356,000 random arguments (machine PI)
 239     atan  :  0.86 ULPs          ...1,536,000 random arguments (machine PI)
 240     (compared with (PI/pi)*(atan, asin, acos, atan2 of x))
 241
 242     tan   :  2.15 ULPs          ...1,024,000 random arguments (true pi)
 243     sin   :  .814 ULPs          ...1,024,000 random arguments (true pi)
 244     cos   :  .792 ULPs          ...1,024,000 random arguments (true pi)
 245     acos  :  2.15 ULPs          ...1,024,000 random arguments (true pi)
 246     asin  :  1.99 ULPs          ...1,024,000 random arguments (true pi)
 247     atan2 :  1.48 ULPs          ...1,024,000 random arguments (true pi)
 248     atan  :  .850 ULPs          ...1,024,000 random arguments (true pi)
 249
 250     acosh :  3.30 ULPs          .....512,000 random arguments
 251     asinh :  1.58 ULPs          .....512,000 random arguments
 252     atanh :  1.71 ULPs          .....512,000 random arguments
 253     cosh  :  1.23 ULPs          .....768,000 random arguments
 254     sinh  :  1.93 ULPs          ...1,024,000 random arguments
 255     tanh  :  2.22 ULPs          ...1,024,000 random arguments
 256     log10 :  1.74 ULPs          ...1,536,000 random arguments
 257     pow   :  1.79 ULPs          .....100,000 random arguments, 0 < x, y < 20.
 258
 259     exp   :  .768 ULPs          ...1,156,000 random arguments
 260     expm1 :  .844 ULPs          ...1,166,000 random arguments
 261     log1p :  .846 ULPs          ...1,536,000 random arguments
 262     log   :  .826 ULPs          ...1,536,000 random arguments
 263     cabs  :  .959 ULPs          .....500,000 random arguments
 264     cbrt  :  .666 ULPs          ...5,120,000 random arguments
 265
 266
 267 5. Speed.
 268
 269         Some functions coded in VAX assembly language (cabs(), hypot() and
 270         sqrt()) are significantly faster than the corresponding ones in 4.2BSD.
 271         In general, to improve performance, all functions in "IEEE/support.c"
 272         should be written in assembly language and, whenever possible, should
 273         be called via short subroutine calls.
 274
 275
 276 6. j0, j1, jn.
 277
 278         The modifications to these routines were only in how an invalid
 279         floating point operations is signaled.