Commit | Line | Data |
---|---|---|
a62df508 | 1 | .\" @(#)README.txt 1.2 (ucb.elefunt) %G% |
079afbd6 ZAL |
2 | .\" troff -ms README.txt |
3 | .de Pi | |
4 | .if n \ | |
5 | pi | |
6 | .if t \ | |
7 | \\(*p | |
8 | .. | |
9 | .ND | |
10 | .ds CH | |
11 | .if t \ | |
12 | .KS | |
13 | .IP \fB...\fR | |
14 | The machine-independent Version 7 math library found in 4.2BSD | |
15 | is now \*Q/usr/lib/libom.a\*U. To compile with those routines use \-lom. | |
16 | .LP | |
17 | .nf | |
18 | K.C. Ng, March 7, 1985, with Z-S. Alex Liu, S. McDonald, P. Tang, W. Kahan. | |
a62df508 | 19 | Revised on 5/10/85, 5/13/85, 6/14/85, 8/20/85, 8/27/85, 9/11/85. |
079afbd6 ZAL |
20 | .fi |
21 | .if n \{\ | |
22 | .LP | |
23 | .nf | |
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 function erf, erfc) have been deleted to prevent overriding the * | |
30 | * system "errno". * | |
31 | ******************************************************************************* | |
32 | .fi \} | |
33 | .if t \{\ | |
34 | .sp 0.5 | |
35 | \fB\l'6i'\fR | |
36 | .QP | |
37 | \fIThis is a description of the upgraded elementary functions (listed in \fB\(sc1\fP). | |
38 | Bessel functions (j0, j1, jn, y0, y1, yn), floor, and fabs passed over | |
39 | from 4.2BSD without change except perhaps for the way floating point | |
40 | exception is signaled on a VAX. Three lines that contain \*Qerrno\*U in erf.c | |
41 | (the error functions erf, erfc) have been deleted to prevent overriding the | |
42 | system \*Qerrno\*U. | |
43 | .LP | |
44 | .sp -0.5 | |
45 | \fB\l'6i'\fR | |
46 | .sp 0.5 \} | |
47 | .IP \fB\(sc0.\fR | |
48 | Total number of files: 40 | |
49 | .sp 0.5 | |
50 | .nf | |
51 | .ta +\w'IEEE/support.c'u+2n +\w'VAX/Makefile'u+2n +\w'VAX/support.s'u+2n +\w'exp__E.c'u+2n \w'log__L.c'u+2n | |
a62df508 | 52 | IEEE/Makefile VAX/Makefile VAX/support.s erf.c lgamma.c |
079afbd6 ZAL |
53 | IEEE/atan2.c VAX/argred.s VAX/tan.s exp.c log.c |
54 | IEEE/cabs.c VAX/atan2.s acosh.c exp__E.c log10.c | |
55 | IEEE/cbrt.c VAX/cabs.s asincos.c expm1.c log1p.c | |
56 | IEEE/support.c VAX/cbrt.s asinh.c floor.c log__L.c | |
57 | IEEE/trig.c VAX/infnan.s atan.c j0.c pow.c | |
58 | Makefile VAX/sincos.s atanh.c j1.c sinh.c | |
59 | README VAX/sqrt.s cosh.c jn.c tanh.c | |
60 | .ta | |
61 | .fi | |
62 | .sp 0.5 | |
63 | .IP \fB\(sc1.\fR | |
64 | Functions implemented: | |
65 | .RS | |
66 | .IP (A). | |
67 | Standard elementary functions (total 22): | |
68 | .nf | |
69 | .ta +\w'expm1(x):=exp(x)\-1'u+4n +\w'... in files'u+1n +\w'\0IEEE/atan2.c\0,'u+1n +\w'\0VAX/sincos.s\0'u+1n | |
70 | acos(x) ... in file \*Qasincos.c\*U | |
71 | asin(x) ... in file \*Qasincos.c\*U | |
72 | atan(x) ... in file \*Qatan.c\*U | |
73 | atan2(x,y) ... in files \*QIEEE/atan2.c\*U, \*QVAX/atan2.s\*U | |
74 | sin(x) ... in files \*QIEEE/trig.c\*U, \*QVAX/sincos.s\*U | |
75 | cos(x) ... in files \*QIEEE/trig.c\*U, \*QVAX/sincos.s\*U | |
76 | tan(x) ... in files \*QIEEE/trig.c\*U, \*QVAX/tan.s\*U | |
77 | cabs(x,y) ... in files \*QIEEE/cabs.c\*U, \*QVAX/cabs.s\*U | |
78 | hypot(x,y) ... in files \*QIEEE/cabs.c\*U, \*QVAX/cabs.s\*U | |
79 | cbrt(x) ... in files \*QIEEE/cbrt.c\*U, \*QVAX/cbrt.s\*U | |
80 | exp(x) ... in file \*Qexp.c\*U | |
81 | expm1(x):=exp(x)\-1 ... in file \*Qexpm1.c\*U | |
82 | log(x) ... in file \*Qlog.c\*U | |
83 | log10(x) ... in file \*Qlog10.c\*U | |
84 | log1p(x):=log(1+x) ... in file \*Qlog1p.c\*U | |
85 | pow(x,y) ... in file \*Qpow.c\*U | |
86 | sinh(x) ... in file \*Qsinh.c\*U | |
87 | cosh(x) ... in file \*Qcosh.c\*U | |
88 | tanh(x) ... in file \*Qtanh.c\*U | |
89 | asinh(x) ... in file \*Qasinh.c\*U | |
90 | acosh(x) ... in file \*Qacosh.c\*U | |
91 | atanh(x) ... in file \*Qatanh.c\*U | |
92 | .ta | |
93 | .fi | |
94 | .sp 0.25 | |
95 | .IP (B). | |
96 | Kernel functions: | |
97 | .nf | |
98 | .ta +\w'libm$argred 'u+2n | |
99 | exp__E(x,c) ... in file \*Qexp__E.c\*U, used by expm1(), exp(), pow() and cosh() | |
100 | log__L(s) ... in file \*Qlog__L.c\*U, used by log1p(), log() and pow() | |
101 | libm$argred ... in file \*QVAX/argred.s\*U, used by VAX version of sin(), cos() and tan() | |
102 | .ta | |
103 | .fi | |
104 | .if t \{\ | |
105 | .RE | |
106 | .KE | |
107 | .bp | |
108 | .KS | |
109 | .RS \} | |
110 | .if n \ | |
111 | .sp 0.25 | |
112 | .IP (C). | |
113 | System supported functions: | |
114 | .nf | |
115 | .ta +\w'copysign()'u+4n +\w'... in files'u+1n | |
116 | sqrt() ... in files \*QIEEE/support.c\*U, \*QVAX/sqrt.s\*U | |
117 | drem() ... in files \*QIEEE/support.c\*U, \*QVAX/support.s\*U | |
118 | finite() ... in files \*QIEEE/support.c\*U, \*QVAX/support.s\*U | |
119 | logb() ... in files \*QIEEE/support.c\*U, \*QVAX/support.s\*U | |
120 | scalb() ... in files \*QIEEE/support.c\*U, \*QVAX/support.s\*U | |
121 | copysign() ... in files \*QIEEE/support.c\*U, \*QVAX/support.s\*U | |
122 | rint() ... in file \*Qfloor.c\*U | |
123 | .ta | |
124 | .fi | |
125 | .sp 0.25 | |
126 | .LP | |
127 | Notes: | |
128 | .IP \fBi\fR. | |
129 | The codes in files ending with \*Q.s\*U are written in VAX assembly | |
130 | language. They are intended for VAX computers. | |
131 | .IP | |
132 | Files that end with \*Q.c\*U are written in C. They are intended | |
133 | for either a VAX or a machine that conforms to the IEEE | |
134 | standard 754 for double precision floating-point arithmetic. | |
135 | .IP \fBii\fR. | |
136 | On other than VAX or IEEE machines, run the original math | |
137 | library, formerly \*Q/usr/lib/libm.a\*U, now \*Q/usr/lib/libom.a\*U, | |
138 | if nothing better is available. | |
139 | .IP \fBiii\fR. | |
140 | The trigonometric functions sin(), cos(), tan() and atan2() in files | |
141 | \*QVAX/sincos.s\*U, \*QVAX/tan.s\*U and \*QVAX/atan2.s\*U are different | |
142 | from those in \*QIEEE/trig.c\*U and \*QIEEE/atan2.c\*U. | |
143 | The VAX assembler code uses the true value of | |
144 | .Pi | |
145 | to perform argument reduction, while the C code uses | |
146 | the machine's value of PI rounded (see \*QIEEE/trig.c\*U). | |
147 | .RE | |
148 | .sp 0.5 | |
149 | .IP \fB\(sc2.\fR | |
150 | A computer system that conforms to IEEE standard 754 should provide | |
151 | .LP | |
152 | .RS | |
153 | .RS | |
154 | .nf | |
155 | sqrt(x), | |
156 | drem(x,p), (double precision remainder function) | |
157 | copysign(x,y), | |
158 | finite(x), | |
159 | scalb(x,N), | |
160 | logb(x) and | |
161 | rint(x). | |
162 | .fi | |
163 | .RE | |
164 | .LP | |
165 | These functions are either required or recommended by the standard. | |
166 | .LP | |
167 | For convenience, a (slow) C implementation of these functions is | |
168 | provided in the file \*QIEEE/support.c\*U. | |
169 | .LP | |
170 | \fBWarning\fR: The functions in \*QIEEE/support.c\*U are somewhat machine | |
171 | dependent. | |
172 | Some modifications may be necessary to run them on a different machine. | |
173 | Currently, if compiled with a suitable flag, \*QIEEE/support.c\*U will work | |
174 | on a National 32000, a Zilog 8000, a VAX, and a SUN (cf. the \*QMakefile\*U in | |
175 | this directory). Invoke the C compiler thus: | |
176 | .RS | |
177 | .nf | |
178 | .ta +\w'cc \-c \-DNATIONAL IEEE/support.c'u+4n +\w'...'u+1n | |
179 | cc \-c \-DVAX IEEE/support.c ... on a VAX, D-format | |
180 | cc \-c \-DNATIONAL IEEE/support.c ... on a National 32000 | |
181 | cc \-c IEEE/support.c ... on other IEEE machines, we hope. | |
182 | .ta | |
183 | .fi | |
184 | .RE | |
185 | .LP | |
186 | Notes: | |
187 | .IP \fBi\fR. | |
188 | Faster versions of drem() and sqrt() for IEEE double precision | |
189 | (coded in C but intended for assembly language) are given at the | |
190 | end of \*QIEEE/support.c\*U but commented out since they require certain | |
191 | machine-dependent functions. | |
192 | .IP \fBii\fR. | |
193 | A fast VAX assembler version of the system supported functions | |
194 | copysign(), logb(), scalb(), finite(), and drem() appears in file | |
195 | \*QVAX/support.s\*U. A fast VAX assembler version of sqrt() is in | |
196 | file \*QVAX/sqrt.s\*U. | |
197 | .RE | |
198 | .if n \ | |
199 | .sp | |
200 | .if t \{\ | |
201 | .RE | |
202 | .KE | |
203 | .bp | |
204 | .KS \} | |
205 | .IP \fB\(sc3.\fR | |
206 | Two formats are supported by all the standard elementary functions: | |
207 | .RS | |
208 | .LP | |
209 | the VAX D-format (56-bit precision), and the IEEE double format | |
210 | (53-bit precision). The cbrt() in \*QIEEE/cbrt.c\*U is for IEEE machines | |
211 | only. The functions in files that end with \*Q.s\*U are for VAX computers | |
212 | only. The functions in files that end with \*Q.c\*U (except \*QIEEE/cbrt.c\*U) | |
213 | are for VAX and IEEE machines. To use the VAX D-format, compile the code | |
214 | with \-DVAX; to use IEEE double format on various IEEE machines, see | |
215 | \*QMakefile\*U in this directory). | |
216 | .LP | |
217 | Example: | |
218 | .RS | |
219 | .nf | |
220 | cc \-c \-DVAX sin.c ... for VAX D-format | |
221 | .fi | |
222 | .RE | |
223 | .sp 0.25 | |
224 | .LP | |
225 | \fBWarning\fR: | |
226 | The values of floating-point constants used in the code are | |
227 | given in both hexadecimal and decimal. The hexadecimal values | |
228 | are the intended ones. The decimal values may be used provided | |
229 | that the compiler converts from decimal to binary accurately | |
230 | enough to produce the hexadecimal values shown. If the | |
231 | conversion is inaccurate, then one must know the exact machine | |
232 | representation of the constants and alter the | |
233 | assembly-language output from the compiler, or play tricks like | |
234 | the following in a C program. | |
235 | .sp 0.25 | |
236 | .RS | |
237 | Example: to store the floating-point constant | |
238 | .sp 0.25 | |
239 | .RS | |
240 | .if n \ | |
241 | p1 = 2**\-6 \(** .F83ABE67E1066A (hexadecimal) | |
242 | .if t \ | |
243 | p1 = 2\u\s-2\-6\s+2\d \(** .F83ABE67E1066A (hexadecimal) | |
244 | .RE | |
245 | on a VAX in C, we use two longwords to store its | |
246 | machine value and define p1 to be the double constant | |
247 | at the location of these two longwords: | |
248 | .sp 0.25 | |
249 | .nf | |
250 | .ta +\w'static long'u+1n +\w'p1x[] ='u+1n | |
251 | static long p1x[] = {0x3abe3d78, 0x066a67e1}; | |
252 | #define p1 (\(**(double\(**)p1x) | |
253 | .ta | |
254 | .fi | |
255 | .RE | |
256 | .IP Note: \w'Note:'u | |
257 | \0On a VAX, some functions have two codes. For example, cabs() | |
258 | has one implementation in \*QIEEE/cabs.c\*U, and the other in | |
259 | \*QVAX/cabs.s\*U. | |
260 | In this case, the assembly language version is preferred. | |
261 | .RE | |
262 | .sp 0.5 | |
263 | .IP \fB\(sc4.\fR | |
264 | Accuracy. | |
265 | .RS | |
266 | .LP | |
267 | The errors in expm1(), log1p(), exp(), log(), cabs(), hypot() | |
268 | and cbrt() are below 1 ULP (Unit in the Last Place). | |
269 | .LP | |
270 | The error in pow(x,y) grows with the size of y. Nevertheless, | |
271 | for integers x and y, pow(x,y) returns the correct integer value | |
272 | on all tested machines (VAX, SUN, National 32000, Zilog 8000), provided that | |
273 | x to the power of y is representable exactly. | |
274 | .LP | |
275 | cosh(), sinh(), acosh(), asinh(), tanh(), atanh() and log10() have errors | |
276 | below about 3 ULPs. | |
277 | .LP | |
278 | For trigonometric and inverse trigonometric functions, | |
279 | let [trig(x)] denote the value actually computed for trig(x). | |
280 | .IP 1) | |
281 | Those codes using PI, the machine's value of | |
282 | .Pi | |
283 | rounded): | |
284 | .nf | |
285 | (in files \*QIEEE/trig.c\*U, \*QIEEE/atan2.c\*U, \*Qasincos.c\*U and \*Qatan.c\*U.) | |
286 | .fi | |
287 | .IP | |
288 | The errors in [sin(x)], [cos(x)], and [atan(x)] are below | |
289 | 1 ULP compared with | |
290 | .if n \ | |
291 | sin(x\(**pi/PI), cos(x\(**pi/PI), and atan(x)\(**PI/pi | |
292 | .if t \ | |
293 | sin(x\(**\(*p/PI), cos(x\(**\(*p/PI), and atan(x)\(**PI/\(*p | |
294 | respectively, where PI is the machine's | |
295 | value of | |
296 | .Pi | |
297 | rounded. [tan(x)] returns | |
298 | .if n \ | |
299 | tan(x\(**pi/PI) | |
300 | .if t \ | |
301 | tan(x\(**\(*p/PI) | |
302 | within about 2 ULPs; [acos(x)], [asin(x)], and [atan2(y,x)] return | |
303 | .if n \ | |
304 | acos(x)\(**PI/pi, asin(x)\(**PI/pi, and atan2(y,x)\(**PI/pi | |
305 | .if t \ | |
306 | acos(x)\(**PI/\(*p, asin(x)\(**PI/\(*p, and atan2(y,x)\(**PI/\(*p | |
307 | respectively to similar accuracy. | |
308 | .IP 2) | |
309 | Those using true | |
310 | .Pi | |
311 | (for VAX D-format only): | |
312 | .br | |
313 | .nf | |
314 | (in files \*QVAX/sincos.s\*U, \*QVAX/tan.s\*U, \*QVAX/atan2.s\*U, \*Qasincos.c\*U and \*Qatan.c\*U.) | |
315 | .fi | |
316 | .IP | |
317 | The errors in [sin(x)], [cos(x)], and [atan(x)] are below | |
318 | 1 ULP. [tan(x)], [atan2(y,x)], [acos(x)] and [asin(x)] | |
319 | have errors below about 2 ULPs. | |
320 | .RE | |
321 | .if t \{\ | |
322 | .KE | |
323 | .bp | |
324 | .KS \} | |
325 | .IP | |
326 | Here are the results of some test runs to find worst errors on a VAX: | |
327 | .nf | |
328 | .ta +\w'expm1'u+1n +\w':'u+2n +\w'2.09'u+1n +\w'ULPs'u+4n +\w'... 1,024,000'uR +1n | |
329 | .sp 0.25 | |
330 | tan : 2.09 ULPs ... 1,024,000 random arguments (machine PI) | |
331 | sin : .861 ULPs ... 1,024,000 random arguments (machine PI) | |
332 | cos : .857 ULPs ... 1,024,000 random arguments (machine PI) | |
333 | .if n \ | |
334 | (compared with tan, sin, cos of (x\(**pi/PI)) | |
335 | .if t \ | |
336 | (compared with tan, sin, cos of (x\(**\(*p/PI)) | |
337 | .sp 0.25 | |
338 | atan : 0.86 ULPs ... 1,536,000 random arguments (machine PI) | |
339 | asin : 2.06 ULPs ... 200,000 random arguments (machine PI) | |
340 | acos : 2.07 ULPs ... 200,000 random arguments (machine PI) | |
341 | atan2 : 1.41 ULPs ... 356,000 random arguments (machine PI) | |
342 | .if n \ | |
343 | (compared with (PI/pi)\(**(atan, asin, acos, atan2 of x)) | |
344 | .if t \ | |
345 | (compared with (PI/\(*p)\(**(atan, asin, acos, atan2 of x)) | |
346 | .sp 0.25 | |
347 | .if n \ | |
348 | tan : 2.15 ULPs ... 1,024,000 random arguments (true pi) | |
349 | .if t \ | |
350 | tan : 2.15 ULPs ... 1,024,000 random arguments (true \(*p) | |
351 | .if n \ | |
352 | sin : .814 ULPs ... 1,024,000 random arguments (true pi) | |
353 | .if t \ | |
354 | sin : .814 ULPs ... 1,024,000 random arguments (true \(*p) | |
355 | .if n \ | |
356 | cos : .792 ULPs ... 1,024,000 random arguments (true pi) | |
357 | .if t \ | |
358 | cos : .792 ULPs ... 1,024,000 random arguments (true \(*p) | |
359 | .if n \ | |
360 | acos : 2.15 ULPs ... 1,024,000 random arguments (true pi) | |
361 | .if t \ | |
362 | acos : 2.15 ULPs ... 1,024,000 random arguments (true \(*p) | |
363 | .if n \ | |
364 | asin : 1.99 ULPs ... 1,024,000 random arguments (true pi) | |
365 | .if t \ | |
366 | asin : 1.99 ULPs ... 1,024,000 random arguments (true \(*p) | |
367 | .if n \ | |
368 | atan2 : 1.48 ULPs ... 1,024,000 random arguments (true pi) | |
369 | .if t \ | |
370 | atan2 : 1.48 ULPs ... 1,024,000 random arguments (true \(*p) | |
371 | .if n \ | |
372 | atan : .850 ULPs ... 1,024,000 random arguments (true pi) | |
373 | .if n \ | |
374 | atan : .850 ULPs ... 1,024,000 random arguments (true \(*p) | |
375 | .sp 0.25 | |
376 | acosh : 3.30 ULPs ... 512,000 random arguments | |
377 | asinh : 1.58 ULPs ... 512,000 random arguments | |
378 | atanh : 1.71 ULPs ... 512,000 random arguments | |
379 | cosh : 1.23 ULPs ... 768,000 random arguments | |
380 | sinh : 1.93 ULPs ... 1,024,000 random arguments | |
381 | tanh : 2.22 ULPs ... 1,024,000 random arguments | |
382 | log10 : 1.74 ULPs ... 1,536,000 random arguments | |
383 | pow : 1.79 ULPs ... 100,000 random arguments, 0 < x, y < 20. | |
384 | .sp 0.25 | |
385 | exp : .768 ULPs ... 1,156,000 random arguments | |
386 | expm1 : .844 ULPs ... 1,166,000 random arguments | |
387 | log1p : .846 ULPs ... 1,536,000 random arguments | |
388 | log : .826 ULPs ... 1,536,000 random arguments | |
389 | cabs : .959 ULPs ... 500,000 random arguments | |
390 | cbrt : .666 ULPs ... 5,120,000 random arguments | |
391 | .ta | |
392 | .fi | |
393 | .sp 0.5 | |
394 | .IP \fB\(sc5.\fR | |
395 | Speed. | |
396 | .IP | |
397 | Some functions coded in VAX assembly language (cabs(), hypot() and sqrt()) | |
398 | are significantly faster than the corresponding ones in 4.2BSD. | |
399 | In general, to improve performance, all functions in \*QIEEE/support.c\*U | |
400 | should be written in assembly language and, whenever possible, should be | |
401 | called via short subroutine calls. | |
402 | .sp 0.5 | |
403 | .IP \fB\(sc6.\fR | |
404 | j0, j1, jn. | |
405 | .IP | |
406 | The modifications to these routines were only in how an invalid | |
407 | floating point operation is signaled on a VAX. | |
408 | .sp 0.5 | |
409 | .if n \ | |
410 | .KS | |
411 | .LP | |
412 | .nf | |
413 | \fB\(sc7.\fR Copyright notice, and Disclaimer: | |
414 | .fi | |
415 | .if n \{\ | |
416 | .LP | |
417 | .nf | |
418 | *************************************************************************** | |
419 | * * | |
420 | * Copyright (c) 1985 Regents of the University of California. * | |
421 | * * | |
422 | * Use and reproduction of this software are granted in accordance with * | |
423 | * the terms and conditions specified in the Berkeley Software License * | |
424 | * Agreement (in particular, this entails acknowledgement of the programs' * | |
425 | * source, and inclusion of this notice) with the additional understanding * | |
426 | * that all recipients should regard themselves as participants in an * | |
427 | * ongoing research project and hence should feel obligated to report * | |
428 | * their experiences (good or bad) with these elementary function codes, * | |
429 | * using "sendbug 4bsd-bugs@BERKELEY", to the authors. * | |
430 | * * | |
431 | *************************************************************************** | |
432 | .fi \} | |
433 | .if t \{\ | |
434 | .sp 0.25 | |
435 | \fB\l'6i'\fR | |
436 | .QP | |
437 | \fICopyright (c) 1985 Regents of the University of California.\fR | |
438 | .QP | |
439 | \fIUse and reproduction of this software are granted in accordance with | |
440 | the terms and conditions specified in the Berkeley Software License | |
441 | Agreement (in particular, this entails acknowledgement of the programs' | |
442 | source, and inclusion of this notice) with the additional understanding | |
443 | that all recipients should regard themselves as participants in an | |
444 | ongoing research project and hence should feel obligated to report | |
445 | their experiences (good or bad) with these elementary function codes, | |
446 | using \*Qsendbug 4bsd-bugs@BERKELEY\*U, to the authors.\fR | |
447 | .LP | |
448 | .sp -0.5 | |
449 | \fB\l'6i'\fR \} | |
450 | .KE |