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