1***************************************************************************
2* *
3* Copyright (c) 1985 Regents of the University of California. *
4* *
5* Use and reproduction of this software are granted in accordance with *
6* the terms and conditions specified in the Berkeley Software License *
7* Agreement (in particular, this entails acknowledgement of the programs' *
8* source, and inclusion of this notice) with the additional understanding *
9* that all recipients should regard themselves as participants in an *
10* ongoing research project and hence should feel obligated to report *
11* their experiences (good or bad) with these elementary function codes, *
12* using "sendbug 4bsd-bugs@BERKELEY", to the authors. *
13* *
14* K.C. Ng, with Z-S. Alex Liu, S. McDonald, P. Tang, W. Kahan. *
15* Revised on 5/10/85, 5/13/85, 6/14/85, 8/20/85, 8/27/85, 9/11/85. *
16* *
17***************************************************************************
18
19/* @(#)README 1.6 (Berkeley) 09/12/85 */
20
21NB. The machine-independent Version 7 math library found in 4.2BSD
22 is now /usr/lib/libom.a. To compile with those routines use -lom.
23
24******************************************************************************
25* This is a description of the upgraded elementary functions (listed in 1). *
26* Bessel functions (j0, j1, jn, y0, y1, yn), floor, and fabs passed over *
27* from 4.2BSD without change except perhaps for the way floating point *
28* exception is signaled on a VAX. Three lines that contain "errno" in erf.c*
29* (error functions erf, erfc) have been deleted to prevent overriding the *
30* system "errno". *
31******************************************************************************
32
330. Total number of files: 40
34
35 IEEE/Makefile VAX/Makefile VAX/support.s erf.c lgamma.c
36 IEEE/atan2.c VAX/argred.s VAX/tan.s exp.c log.c
37 IEEE/cabs.c VAX/atan2.s acosh.c exp__E.c log10.c
38 IEEE/cbrt.c VAX/cabs.s asincos.c expm1.c log1p.c
39 IEEE/support.c VAX/cbrt.s asinh.c floor.c log__L.c
40 IEEE/trig.c VAX/infnan.s atan.c j0.c pow.c
41 Makefile VAX/sincos.s atanh.c j1.c sinh.c
42 README VAX/sqrt.s cosh.c jn.c tanh.c
43
441. Functions implemented :
45 (A). Standard elementary functions (total 22) :
46 acos(x) ...in file asincos.c
47 asin(x) ...in file asincos.c
48 atan(x) ...in file atan.c
49 atan2(x,y) ...in files IEEE/atan2.c, VAX/atan2.s
50 sin(x) ...in files IEEE/trig.c, VAX/sincos.s
51 cos(x) ...in files IEEE/trig.c, VAX/sincos.s
52 tan(x) ...in files IEEE/trig.c, VAX/tan.s
53 cabs(x,y) ...in files IEEE/cabs.c, VAX/cabs.s
54 hypot(x,y) ...in files IEEE/cabs.c, VAX/cabs.s
55 cbrt(x) ...in files IEEE/cbrt.c, VAX/cbrt.s
56 exp(x) ...in file exp.c
57 expm1(x):=exp(x)-1 ...in file expm1.c
58 log(x) ...in file log.c
59 log10(x) ...in file log10.c
60 log1p(x):=log(1+x) ...in file log1p.c
61 pow(x,y) ...in file pow.c
62 sinh(x) ...in file sinh.c
63 cosh(x) ...in file cosh.c
64 tanh(x) ...in file tanh.c
65 asinh(x) ...in file asinh.c
66 acosh(x) ...in file acosh.c
67 atanh(x) ...in file atanh.c
68
69 (B). Kernel functions :
70 exp__E(x,c) ...in file exp__E.c, used by expm1/exp/pow/cosh
71 log__L(s) ...in file log__L.c, used by log1p/log/pow
72 libm$argred ...in file VAX/argred.s, used by VAX version of sin/cos/tan
73
74 (C). System supported functions :
75 sqrt() ...in files IEEE/support.c, VAX/sqrt.s
76 drem() ...in files IEEE/support.c, VAX/support.s
77 finite() ...in files IEEE/support.c, VAX/support.s
78 logb() ...in files IEEE/support.c, VAX/support.s
79 scalb() ...in files IEEE/support.c, VAX/support.s
80 copysign() ...in files IEEE/support.c, VAX/support.s
81 rint() ...in file floor.c
82
83
84 Notes:
85 i. The codes in files ending with ".s" are written in VAX assembly
86 language. They are intended for VAX computers.
87
88 Files that end with ".c" are written in C. They are intended
89 for either a VAX or a machine that conforms to the IEEE
90 standard 754 for double precision floating-point arithmetic.
91
92 ii. On other than VAX or IEEE machines, run the original math
93 library, formerly "/usr/lib/libm.a", now "/usr/lib/libom.a", if
94 nothing better is available.
95
96 iii. The trigonometric functions sin/cos/tan/atan2 in files "VAX/sincos.s",
97 "VAX/tan.s" and "VAX/atan2.s" are different from those in
98 "IEEE/trig.c" and "IEEE/atan2.c". The VAX assembler code uses the
99 true value of pi to perform argument reduction, while the C code uses
100 a machine value of PI (see "IEEE/trig.c").
101
102
1032. A computer system that conforms to IEEE standard 754 should provide
104 sqrt(x),
105 drem(x,p), (double precision remainder function)
106 copysign(x,y),
107 finite(x),
108 scalb(x,N),
109 logb(x) and
110 rint(x).
111 These functions are either required or recommended by the standard.
112 For convenience, a (slow) C implementation of these functions is
113 provided in the file "IEEE/support.c".
114
115 Warning: The functions in IEEE/support.c are somewhat machine dependent.
116 Some modifications may be necessary to run them on a different machine.
117 Currently, if compiled with a suitable flag, "IEEE/support.c" will work
118 on a National 32000, a Zilog 8000, a VAX, and a SUN (cf. the "Makefile"
119 in this directory). Invoke the C compiler thus:
120
121 cc -c -DVAX IEEE/support.c ... on a VAX, D-format
122 cc -c -DNATIONAL IEEE/support.c ... on a National 32000
123 cc -c IEEE/support.c ... on other IEEE machines,
124 we hope.
125
126 Notes:
127 1. Faster versions of "drem" and "sqrt" for IEEE double precision
128 (coded in C but intended for assembly language) are given at the
129 end of "IEEE/support.c" but commented out since they require certain
130 machine-dependent functions.
131
132 2. A fast VAX assembler version of the system supported functions
133 copysign(), logb(), scalb(), finite(), and drem() appears in file
134 "VAX/support.s". A fast VAX assembler version of sqrt() is in
135 file "VAX/sqrt.s".
136
1373. Two formats are supported by all the standard elementary functions:
138 the VAX D-format (56-bit precision), and the IEEE double format
139 (53-bit precision). The cbrt() in "IEEE/cbrt.c" is for IEEE machines
140 only. The functions in files that end with ".s" are for VAX computers
141 only. The functions in files that end with ".c" (except "IEEE/cbrt.c")
142 are for VAX and IEEE machines. To use the VAX D-format, compile the code
143 with -DVAX; to use IEEE double format on various IEEE machines, see
144 "Makefile" in this directory).
145
146 Example:
147 cc -c -DVAX sin.c ... for VAX D-format
148
149 Warning: The values of floating-point constants used in the code are
150 given in both hexadecimal and decimal. The hexadecimal values
151 are the intended ones. The decimal values may be used provided
152 that the compiler converts from decimal to binary accurately
153 enough to produce the hexadecimal values shown. If the
154 conversion is inaccurate, then one must know the exact machine
155 representation of the constants and alter the assembly
156 language output from the compiler, or play tricks like
157 the following in a C program.
158
159 Example: to store the floating-point constant
160
161 p1= 2^-6 * .F83ABE67E1066A (Hexadecimal)
162
163 on a VAX in C, we use two longwords to store its
164 machine value and define p1 to be the double constant
165 at the location of these two longwords:
166
167 static long p1x[] = { 0x3abe3d78, 0x066a67e1};
168 #define p1 (*(double*)p1x)
169
170 Note: On a VAX, some functions have two codes. For example, cabs() has
171 one implementation in "IEEE/cabs.c", and another in "VAX/cabs.s".
172 In this case, the assembly language version is preferred.
173
174
1754. Accuracy.
176
177 The errors in expm1(), log1p(), exp(), log(), cabs(), hypot()
178 and cbrt() are below 1 ULP (Unit in the Last Place).
179
180 The error in pow(x,y) grows with the size of y. Nevertheless,
181 for integers x and y, pow(x,y) returns the correct integer value
182 on all tested machines (VAX, SUN, NATIONAL, ZILOG), provided that
183 x to the power of y is representable exactly.
184
185 cosh, sinh, acosh, asinh, tanh, atanh and log10 have errors below
186 about 3 ULPs.
187
188 For trigonometric and inverse trigonometric functions:
189
190 Let [trig(x)] denote the value actually computed for trig(x),
191
192 1) Those codes using the machine's value PI (true pi rounded):
193 (source codes: IEEE/{trig.c,atan2.c}, asincos.c and atan.c)
194
195 The errors in [sin(x)], [cos(x)], and [atan(x)] are below
196 1 ULP compared with sin(x*pi/PI), cos(x*pi/PI), and
197 atan(x)*PI/pi respectively, where PI is the machine's
198 value of pi rounded. [tan(x)] returns tan(x*pi/PI) within
199 about 2 ULPs; [acos(x)], [asin(x)], and [atan2(y,x)]
200 return acos(x)*PI/pi, asin(x)*PI/pi, and atan2(y,x)*PI/pi
201 respectively to similar accuracy.
202
203
204 2) Those using true pi (for VAX D-format only):
205 (source codes: VAX/{sincos.s,tan.s,atan2.s}, asincos.c and
206 atan.c)
207
208 The errors in [sin(x)], [cos(x)], and [atan(x)] are below
209 1 ULP. [tan(x)], [atan2(y,x)], [acos(x)], and [asin(x)]
210 have errors below about 2 ULPs.
211
212
213 Here are the results of some test runs to find worst errors on
214 the VAX :
215
216 tan : 2.09 ULPs ...1,024,000 random arguments (machine PI)
217 sin : .861 ULPs ...1,024,000 random arguments (machine PI)
218 cos : .857 ULPs ...1,024,000 random arguments (machine PI)
219 (compared with tan, sin, cos of (x*pi/PI))
220
221 acos : 2.07 ULPs .....200,000 random arguments (machine PI)
222 asin : 2.06 ULPs .....200,000 random arguments (machine PI)
223 atan2 : 1.41 ULPs .....356,000 random arguments (machine PI)
224 atan : 0.86 ULPs ...1,536,000 random arguments (machine PI)
225 (compared with (PI/pi)*(atan, asin, acos, atan2 of x))
226
227 tan : 2.15 ULPs ...1,024,000 random arguments (true pi)
228 sin : .814 ULPs ...1,024,000 random arguments (true pi)
229 cos : .792 ULPs ...1,024,000 random arguments (true pi)
230 acos : 2.15 ULPs ...1,024,000 random arguments (true pi)
231 asin : 1.99 ULPs ...1,024,000 random arguments (true pi)
232 atan2 : 1.48 ULPs ...1,024,000 random arguments (true pi)
233 atan : .850 ULPs ...1,024,000 random arguments (true pi)
234
235 acosh : 3.30 ULPs .....512,000 random arguments
236 asinh : 1.58 ULPs .....512,000 random arguments
237 atanh : 1.71 ULPs .....512,000 random arguments
238 cosh : 1.23 ULPs .....768,000 random arguments
239 sinh : 1.93 ULPs ...1,024,000 random arguments
240 tanh : 2.22 ULPs ...1,024,000 random arguments
241 log10 : 1.74 ULPs ...1,536,000 random arguments
242 pow : 1.79 ULPs .....100,000 random arguments, 0 < x, y < 20.
243
244 exp : .768 ULPs ...1,156,000 random arguments
245 expm1 : .844 ULPs ...1,166,000 random arguments
246 log1p : .846 ULPs ...1,536,000 random arguments
247 log : .826 ULPs ...1,536,000 random arguments
248 cabs : .959 ULPs .....500,000 random arguments
249 cbrt : .666 ULPs ...5,120,000 random arguments
250
251
2525. Speed.
253
254 Some functions coded in VAX assembly language (cabs(), hypot() and
255 sqrt()) are significantly faster than the corresponding ones in 4.2BSD.
256 In general, to improve performance, all functions in "IEEE/support.c"
257 should be written in assembly language and, whenever possible, should
258 be called via short subroutine calls.
259
260
2616. j0, j1, jn.
262
263 The modifications to these routines were only in how an invalid
264 floating point operations is signaled.
265