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