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