+K.C. Ng, March 7, 1985
+Revised on 5/10/85, 5/13/85.
+
+******************************************************************************
+* This is a description of the upgraded elementary functions (listed in 1). *
+* Bessel functions (j0, j1, jn, y0, y1, yn), gamma function (gamma), error *
+* functions (erf, erfc), floor, and fabs are passed over from 4.2BSD *
+* without change. They will not be discussed here. *
+******************************************************************************
+
+0. Total number of files: 34
+
+ ( unchanges files from 4.2BSD, total 6 )
+ erf.c floor.c gamma.c j0.c j1.c jn.c
+
+ ( new files for 4.3BSD, total 27 )
+ E.c Read_me atanh.c pow.c
+ IEEE/atan2.c VAX/atan2.s atrig.c sinh.c
+ IEEE/cabs.c VAX/cabs.s cosh.c tanh.c
+ IEEE/cbrt.c VAX/cbrt.s exp.c
+ IEEE/support.c VAX/support.s exp__E.c
+ IEEE/trig.c VAX/trig.s log.c
+ L.c acosh.c log10.c
+ Makefile asinh.c log__L.c
+
+
+1. Functions implemented (total 22):
+ (A). Standard elementary functions:
+ acos(x) ...in file atrig.c
+ asin(x) ...in file atrig.c
+ atan(x) ...in file atrig.c
+ atan2(x,y) ...in files IEEE/atan2.c, VAX/atan2.s
+ sin(x) ...in files IEEE/trig.c , VAX/trig.s
+ cos(x) ...in files IEEE/trig.c , VAX/trig.s
+ tan(x) ...in files IEEE/trig.c , VAX/trig.s
+ cabs(x,y) ...in files IEEE/cabs.c , VAX/cabs.s
+ hypot(x,y) ...in file IEEE/cabs.c , VAX/cabs.s
+ cbrt(x) ...in files IEEE/cbrt.c , VAX/cbrt.s
+ exp(x) ...in file exp.c
+ E(x):=exp(x)-1 ...in file E.c
+ log(x) ...in file log.c
+ L(x):=log(1+x) ...in file L.c
+ log10(x) ...in file log10.c
+ pow(x,y) ...in file pow.c
+ sinh(x) ...in file sinh.c
+ cosh(x) ...in file cosh.c
+ tanh(x) ...in file cosh.c
+ asinh(x) ...in file asinh.c
+ acosh(x) ...in file acosh.c
+ atanh(x) ...in file atanh.c
+
+ (B). Kernel functions :
+ sin__S(y) ...in file IEEE/trig.c, used by sin/cos/tan in IEEE/trig.c
+ cos__C(y) ...in file IEEE/trig.c, used by sin/cos/tan in IEEE/trig.c
+ exp__E(x,c) ...in file exp__E.c, used by E/exp/pow/cosh
+ log__L(s) ...in file log__L.c, used by L/log/pow
+
+ (C). System supported functions :
+ drem() ...in files IEEE/support.c , VAX/support.s
+ sqrt() ...in files IEEE/support.c , VAX/support.s
+ finite() ...in files IEEE/support.c , VAX/support.s
+ logb() ...in files IEEE/support.c , VAX/support.s
+ scalb() ...in files IEEE/support.c , VAX/support.s
+ copysign() ...in files IEEE/support.c , VAX/support.s
+
+
+ Notes:
+ i. The codes in files ending with .s are written in VAX assembly
+ language. They are intended for VAX computers.
+
+ Files that end with .c are written in C. They are intended
+ for either a VAX or a machine that conforms to the IEEE
+ standard 754 for (double-precision) floating-point arithmetic.
+
+ ii. On other than VAX or IEEE machines, run the original math
+ library, formerly libm.a, now olibm.a, if nothing better
+ is available.
+
+
+ iii. The trigonometric functions sin/cos/tan/atan2 in files "VAX/trig.s"
+ and "VAX/atan2.s" are different from those in "IEEE/trig.c" and
+ "IEEE/atan2.c". The VAX assembler code uses the true value of pi to
+ perform argument reduction, while the C code uses a machine value of
+ PI (see "IEEE/trig.c").
+
+
+2. A computer system that conforms to IEEE standard 754 should provide
+ sqrt(x),
+ drem(x,p), (double precision remainder function)
+ copysign(x,y),
+ finite(x),
+ scalb(x,N),
+ logb(x) .
+ These functions are required or recommended by the standard.
+ For convenience, a (slow) C implementation of these functions is
+ provided in the file IEEE/support.c.
+
+ Warning: The functions in IEEE/support.c are somewhat machine dependent.
+ Some modifications may be necessary to run them on a different machine.
+ Currently, if compiled with a suitable flag, IEEE/support.c will work on a
+ National 32000, a Zilog 8000, a VAX, and a SUN (cf. the "Makefile" in
+ this directory). Invoke the C compiler thus:
+
+ cc -c -DVAX IEEE/support.c ... on a VAX, D-format
+ cc -c -DNATIONAL IEEE/support.c ... on a National 32081
+ cc -c IEEE/support.c ... on other IEEE machines,
+ we hope.
+
+ Notes:
+ 1. Faster versions of "drem" and "sqrt" for IEEE double precision
+ (coded in C but intended for assembly language) are given at the
+ end of support.c but commented out since they require certain
+ machine-dependent functions.
+
+ 2. A fast VAX assembler version of IEEE/support.c (copysign(), logb(),
+ scalb(), finite(), drem() and sqrt()) appears in file VAX/support.s.
+
+3. Two formats are supported by all the standard elementary functions:
+ the VAX D-format (56 bits' precision), and the IEEE double format
+ (53 bits' precision). The cbrt() in IEEE/cbrt.c is for IEEE machines
+ only. The functions in files that end with ".s" are for VAX computers
+ only. The functions in files that end with ".c" (except IEEE/cbrt.c) are
+ for VAX and IEEE machines. To use the VAX D-format, compile the code
+ with -DVAX; to use IEEE double format on various IEEE machines, see
+ Makefile in this directory).
+
+ Example:
+ cc -c -DVAX sin.c ... for VAX D-format
+
+ Warning: The values of floating constants used in the code are given
+ in both hexadecimal and decimal. The hexadecimal values are
+ the intended ones. The decimal values may be use provided
+ that the compiler converts from decimal to binary accurately
+ enough to produce the hexadecimal values shown. If the
+ conversion is inaccurate, then one must know the exact machine
+ representation of the constants and alter the assembly-
+ language output from the compiler, or apply tricks like
+ the following in a C program.
+
+ Example: to store the floating point constant
+
+ p1= 2^-6 * .F83ABE67E1066A (Hexadecimal)
+
+ on a VAX in C, we use two long word to store its
+ machine value and define p1 to be the double constant
+ at the location of these two long words:
+
+ static long p1x[] = { 0x3abe3d78, 0x066a67e1};
+ #define p1 (*(double*)p1x)
+
+ Note: On a VAX, some functions have two set of codes. For example,
+ cabs() has one implementation in cabs.c, and another in
+ VAX/cabs.s. In this case, the assembly version is preferred.
+
+
+4. Accuracy.
+
+ The errors in E(), L(), exp(), log(), cabs(), hypot(), and cbrt()
+ are below 1 ULP (Unit in the Last Place).
+
+ The error in pow(x,y) grows with the size of y. Nevertheless,
+ for integers x and y, pow(x,y) returns the correct integer value
+ on all tested machines (VAX, SUN, NATIONAL, ZILOG), provided that
+ x to the power of y is representable exactly.
+
+ cosh, sinh, acosh, asinh, tanh, atanh and log10 have errors below
+ about 3 ULPs.
+
+ For trig and atrig functions:
+
+ Let [trig(x)] denote the value actually computed for trig(x),
+
+ 1) Those using machine value PI (true pi rounded):
+ (source codes: IEEE/trig.c, IEEE/atan2.c and atrig.c)
+
+ The errors in [sin(x)], [cos(x)], and [atan(x)] are below
+ 1 ULP compared with sin(x*pi/PI), cos(x*pi/PI), and
+ atan(x)*PI/pi respectively, where PI is the machine's
+ value of pi rounded. [Tan(x)] returns tan(x*pi/PI) within
+ about 2 ULPs; [acos(x)], [asin(x)], and [atan2(y,x)]
+ return acos(x)*PI/pi, asin(x)*PI/pi, and atan2(y,x)*PI/pi
+ respectively to similar accuracy.
+
+
+ 2) Those using true pi (for VAX D-format only)
+ (source codes: VAX/trig.s, VAX/atan2.s and atrig.c)
+
+ The errors in [sin(x)], [cos(x)], and [atan(x)] are below
+ 1 ULP. [Tan(x)], [atan2(y,x)], [acos(x)], and [asin(x)]
+ have errors below about 2 ULPs.
+
+
+ Here are the results of some test runs to find worst errors on
+ the VAX :
+
+ tan : 2.08 ULPs ...1,024,000 random arguments (machine PI)
+ sin : .861 ULPs ...1,024,000 random arguments (machine PI)
+ cos : .857 ULPs ...1,024,000 random arguments (machine PI)
+ (compared with tan, sin, cos of (x*pi/PI))
+
+ acos : 2.07 ULPs .....200,000 random arguments (machine PI)
+ asin : 2.06 ULPs .....200,000 random arguments (machine PI)
+ atan2 : 1.41 ULPs .....356,000 random arguments (machine PI)
+ atan : 0.86 ULPs ...1,536,000 random arguments (machine PI)
+ (compared with (PI/pi)*(atan, asin, acos, atan2 of x))
+
+ tan : 2.15 ULPs ...1,024,000 random arguments (true pi)
+ sin : .814 ULPs ...1,024,000 random arguments (true pi)
+ cos : .792 ULPs ...1,024,000 random arguments (true pi)
+ acos : 2.15 ULPs ...1,024,000 random arguments (true pi)
+ asin : 1.99 ULPs ...1,024,000 random arguments (true pi)
+ atan2 : 1.48 ULPs ...1,024,000 random arguments (true pi)
+ atan : .850 ULPs ...1,024,000 random arguments (true pi)
+
+ acosh : 3.20 ULPs .....500,000 random arguments
+ asinh : 1.58 ULPs ..... 52,000 random arguments
+ atanh : 1.45 ULPs .....200,000 random arguments
+ cosh : 1.23 ULPs .....768,000 random arguments
+ sinh : 1.93 ULPs ...1,024,000 random arguments
+ tanh : 2.22 ULPs ...1,024,000 random arguments
+ log10 : 1.74 ULPs ...1,536,000 random arguments
+ pow : 1.79 ULPs .....100,000 random arguments, 0 < x, y < 20.
+
+ exp : .768 ULPs ...1,156,000 random arguments
+ E : .844 ULPs ...1,166,000 random arguments
+ L : .846 ULPs ...1,536,000 random arguments
+ log : .826 ULPs ...1,536,000 random arguments
+ cabs : .959 ULPs .....500,000 random arguments
+ cbrt : .666 ULPs ...5,120,000 random arguments
+
+
+5. Speed.
+
+ The functions coded in VAX assembly language (sin, cos, tan, atan2,
+ cabs, hypot, sqrt, cbrt) are significantly faster than the
+ corresponding ones in 4.2BSD. In general, to improve performance,
+ all functions in IEEE/support.c should be written in assembler and,
+ whenever possible, should be called via short subroutine calls.
+
+
+6. Copyright notice and Miscellaneous.
+
+***************************************************************************
+* *
+* Copyright (c) 1985 Regents of the University of California. *
+* *
+* Use and reproduction of this software are granted in accordance with *
+* the terms and conditions specified in the Berkeley Software License *
+* Agreement (in particular, this entails acknowledgement of the programs' *
+* source, and inclusion of this notice) with the additional understanding *
+* that all recipients should regard themselves as participants in an *
+* ongoing research project and hence should feel obligated to report *
+* their experiences (good or bad) with these elementary function codes, *
+* using "sendbug 4bsd-bugs@BERKELEY", to the authors. *
+* *
+***************************************************************************
+
projects / unix-history / commitdiff