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