troff -ms README.txt
pi \\(*p .. .ND .KS
K.C. Ng, March 7, 1985, with Z-S. Alex Liu, S. McDonald, P. Tang, W. Kahan. Revised on 5/10/85, 5/13/85, 6/14/85, 8/20/85, 8/27/85, 9/11/85.
******************************************************************************* * This is a description of the upgraded elementary functions (listed in 1). * * Bessel functions (j0, j1, jn, y0, y1, yn), floor, and fabs passed over * * from 4.2BSD without change except perhaps for the way floating point * * exception is signaled on a VAX. Three lines that contain "errno" in erf.c * * (error function erf, erfc) have been deleted to prevent overriding the * * system "errno". * *******************************************************************************\} \l'6i' .QP This is a description of the upgraded elementary functions (listed in \(sc1). Bessel functions (j0, j1, jn, y0, y1, yn), floor, and fabs passed over from 4.2BSD without change except perhaps for the way floating point exception is signaled on a VAX. Three lines that contain \*Qerrno\*U in erf.c (the error functions erf, erfc) have been deleted to prevent overriding the system \*Qerrno\*U.
\l'6i'
IEEE/Makefile VAX/Makefile VAX/support.s erf.c lgamma.c IEEE/atan2.c VAX/argred.s VAX/tan.s exp.c log.c IEEE/cabs.c VAX/atan2.s acosh.c exp__E.c log10.c IEEE/cbrt.c VAX/cabs.s asincos.c expm1.c log1p.c IEEE/support.c VAX/cbrt.s asinh.c floor.c log__L.c IEEE/trig.c VAX/infnan.s atan.c j0.c pow.c Makefile VAX/sincos.s atanh.c j1.c sinh.c README VAX/sqrt.s cosh.c jn.c tanh.c
acos(x) ... in file \*Qasincos.c\*U asin(x) ... in file \*Qasincos.c\*U atan(x) ... in file \*Qatan.c\*U atan2(x,y) ... in files \*QIEEE/atan2.c\*U, \*QVAX/atan2.s\*U sin(x) ... in files \*QIEEE/trig.c\*U, \*QVAX/sincos.s\*U cos(x) ... in files \*QIEEE/trig.c\*U, \*QVAX/sincos.s\*U tan(x) ... in files \*QIEEE/trig.c\*U, \*QVAX/tan.s\*U cabs(x,y) ... in files \*QIEEE/cabs.c\*U, \*QVAX/cabs.s\*U hypot(x,y) ... in files \*QIEEE/cabs.c\*U, \*QVAX/cabs.s\*U cbrt(x) ... in files \*QIEEE/cbrt.c\*U, \*QVAX/cbrt.s\*U exp(x) ... in file \*Qexp.c\*U expm1(x):=exp(x)-1 ... in file \*Qexpm1.c\*U log(x) ... in file \*Qlog.c\*U log10(x) ... in file \*Qlog10.c\*U log1p(x):=log(1+x) ... in file \*Qlog1p.c\*U pow(x,y) ... in file \*Qpow.c\*U sinh(x) ... in file \*Qsinh.c\*U cosh(x) ... in file \*Qcosh.c\*U tanh(x) ... in file \*Qtanh.c\*U asinh(x) ... in file \*Qasinh.c\*U acosh(x) ... in file \*Qacosh.c\*U atanh(x) ... in file \*Qatanh.c\*U
exp__E(x,c) ... in file \*Qexp__E.c\*U, used by expm1(), exp(), pow() and cosh() log__L(s) ... in file \*Qlog__L.c\*U, used by log1p(), log() and pow() libm$argred ... in file \*QVAX/argred.s\*U, used by VAX version of sin(), cos() and tan()
sqrt() ... in files \*QIEEE/support.c\*U, \*QVAX/sqrt.s\*U drem() ... in files \*QIEEE/support.c\*U, \*QVAX/support.s\*U finite() ... in files \*QIEEE/support.c\*U, \*QVAX/support.s\*U logb() ... in files \*QIEEE/support.c\*U, \*QVAX/support.s\*U scalb() ... in files \*QIEEE/support.c\*U, \*QVAX/support.s\*U copysign() ... in files \*QIEEE/support.c\*U, \*QVAX/support.s\*U rint() ... in file \*Qfloor.c\*U
Notes:
i to perform argument reduction, while the C code uses the machine's value of PI rounded (see \*QIEEE/trig.c\*U).
sqrt(x), drem(x,p), (double precision remainder function) copysign(x,y), finite(x), scalb(x,N), logb(x) and rint(x).
These functions are either required or recommended by the standard.
For convenience, a (slow) C implementation of these functions is provided in the file \*QIEEE/support.c\*U.
Warning: The functions in \*QIEEE/support.c\*U are somewhat machine dependent. Some modifications may be necessary to run them on a different machine. Currently, if compiled with a suitable flag, \*QIEEE/support.c\*U will work on a National 32000, a Zilog 8000, a VAX, and a SUN (cf. the \*QMakefile\*U 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 32000 cc -c IEEE/support.c ... on other IEEE machines, we hope.
Notes:
the VAX D-format (56-bit precision), and the IEEE double format (53-bit precision). The cbrt() in \*QIEEE/cbrt.c\*U is for IEEE machines only. The functions in files that end with \*Q.s\*U are for VAX computers only. The functions in files that end with \*Q.c\*U (except \*QIEEE/cbrt.c\*U) 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 \*QMakefile\*U in this directory).
Example:
cc -c -DVAX sin.c ... for VAX D-format
Warning: The values of floating-point constants used in the code are given in both hexadecimal and decimal. The hexadecimal values are the intended ones. The decimal values may be used 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 play tricks like the following in a C program.
Example: to store the floating-point constant
p1 = 2**-6 \(** .F83ABE67E1066A (hexadecimal) p1 = 2\u\s-2-6\s+2\d \(** .F83ABE67E1066A (hexadecimal)
static long p1x[] = {0x3abe3d78, 0x066a67e1}; #define p1 (\(**(double\(**)p1x)
The errors in expm1(), log1p(), 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 32000, Zilog 8000), 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 trigonometric and inverse trigonometric functions, let [trig(x)] denote the value actually computed for trig(x).
i rounded):
(in files \*QIEEE/trig.c\*U, \*QIEEE/atan2.c\*U, \*Qasincos.c\*U and \*Qatan.c\*U.)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 sin(x\(**\(*p/PI), cos(x\(**\(*p/PI), and atan(x)\(**PI/\(*p respectively, where PI is the machine's value of
i rounded. [tan(x)] returns tan(x\(**pi/PI) tan(x\(**\(*p/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 acos(x)\(**PI/\(*p, asin(x)\(**PI/\(*p, and atan2(y,x)\(**PI/\(*p respectively to similar accuracy.
i (for VAX D-format only):
(in files \*QVAX/sincos.s\*U, \*QVAX/tan.s\*U, \*QVAX/atan2.s\*U, \*Qasincos.c\*U and \*Qatan.c\*U.)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.
tan : 2.09 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)) (compared with tan, sin, cos of (x\(**\(*p/PI)) atan : 0.86 ULPs ... 1,536,000 random arguments (machine PI) asin : 2.06 ULPs ... 200,000 random arguments (machine PI) acos : 2.07 ULPs ... 200,000 random arguments (machine PI) atan2 : 1.41 ULPs ... 356,000 random arguments (machine PI) (compared with (PI/pi)\(**(atan, asin, acos, atan2 of x)) (compared with (PI/\(*p)\(**(atan, asin, acos, atan2 of x)) tan : 2.15 ULPs ... 1,024,000 random arguments (true pi) tan : 2.15 ULPs ... 1,024,000 random arguments (true \(*p) sin : .814 ULPs ... 1,024,000 random arguments (true pi) sin : .814 ULPs ... 1,024,000 random arguments (true \(*p) cos : .792 ULPs ... 1,024,000 random arguments (true pi) cos : .792 ULPs ... 1,024,000 random arguments (true \(*p) acos : 2.15 ULPs ... 1,024,000 random arguments (true pi) acos : 2.15 ULPs ... 1,024,000 random arguments (true \(*p) asin : 1.99 ULPs ... 1,024,000 random arguments (true pi) asin : 1.99 ULPs ... 1,024,000 random arguments (true \(*p) atan2 : 1.48 ULPs ... 1,024,000 random arguments (true pi) atan2 : 1.48 ULPs ... 1,024,000 random arguments (true \(*p) atan : .850 ULPs ... 1,024,000 random arguments (true pi) atan : .850 ULPs ... 1,024,000 random arguments (true \(*p) acosh : 3.30 ULPs ... 512,000 random arguments asinh : 1.58 ULPs ... 512,000 random arguments atanh : 1.71 ULPs ... 512,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 expm1 : .844 ULPs ... 1,166,000 random arguments log1p : .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
\(sc7. Copyright notice, and Disclaimer:
*************************************************************************** * * * 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. * * * ***************************************************************************\} \l'6i' .QP Copyright (c) 1985 Regents of the University of California. .QP 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 \*Qsendbug 4bsd-bugs@BERKELEY\*U, to the authors.
\l'6i' \} .KE