1*24602Szliu /* 2*24602Szliu * Copyright (c) 1985 Regents of the University of California. 3*24602Szliu * 4*24602Szliu * Use and reproduction of this software are granted in accordance with 5*24602Szliu * the terms and conditions specified in the Berkeley Software License 6*24602Szliu * Agreement (in particular, this entails acknowledgement of the programs' 7*24602Szliu * source, and inclusion of this notice) with the additional understanding 8*24602Szliu * that all recipients should regard themselves as participants in an 9*24602Szliu * ongoing research project and hence should feel obligated to report 10*24602Szliu * their experiences (good or bad) with these elementary function codes, 11*24602Szliu * using "sendbug 4bsd-bugs@BERKELEY", to the authors. 12*24602Szliu */ 13*24602Szliu 14*24602Szliu #ifndef lint 15*24602Szliu static char sccsid[] = "@(#)log10.c 1.1 (ELEFUNT) 09/06/85"; 16*24602Szliu #endif not lint 17*24602Szliu 18*24602Szliu /* LOG10(X) 19*24602Szliu * RETURN THE BASE 10 LOGARITHM OF x 20*24602Szliu * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) 21*24602Szliu * CODED IN C BY K.C. NG, 1/20/85; 22*24602Szliu * REVISED BY K.C. NG on 1/23/85, 3/7/85, 4/16/85. 23*24602Szliu * 24*24602Szliu * Required kernel function: 25*24602Szliu * log(x) 26*24602Szliu * 27*24602Szliu * Method : 28*24602Szliu * log(x) 29*24602Szliu * log10(x) = --------- or [1/log(10)]*log(x) 30*24602Szliu * log(10) 31*24602Szliu * 32*24602Szliu * Note: 33*24602Szliu * [log(10)] rounded to 56 bits has error .0895 ulps, 34*24602Szliu * [1/log(10)] rounded to 53 bits has error .198 ulps; 35*24602Szliu * therefore, for better accuracy, in VAX D format, we divide 36*24602Szliu * log(x) by log(10), but in IEEE Double format, we multiply 37*24602Szliu * log(x) by [1/log(10)]. 38*24602Szliu * 39*24602Szliu * Special cases: 40*24602Szliu * log10(x) is NaN with signal if x < 0; 41*24602Szliu * log10(+INF) is +INF with no signal; log10(0) is -INF with signal; 42*24602Szliu * log10(NaN) is that NaN with no signal. 43*24602Szliu * 44*24602Szliu * Accuracy: 45*24602Szliu * log10(X) returns the exact log10(x) nearly rounded. In a test run 46*24602Szliu * with 1,536,000 random arguments on a VAX, the maximum observed 47*24602Szliu * error was 1.74 ulps (units in the last place). 48*24602Szliu * 49*24602Szliu * Constants: 50*24602Szliu * The hexadecimal values are the intended ones for the following constants. 51*24602Szliu * The decimal values may be used, provided that the compiler will convert 52*24602Szliu * from decimal to binary accurately enough to produce the hexadecimal values 53*24602Szliu * shown. 54*24602Szliu */ 55*24602Szliu 56*24602Szliu #ifdef VAX /* VAX D format (56 bits) */ 57*24602Szliu /* static double */ 58*24602Szliu /* ln10hi = 2.3025850929940456790E0 ; Hex 2^ 2 * .935D8DDDAAA8AC */ 59*24602Szliu static long ln10hix[] = { 0x5d8d4113, 0xa8acddaa}; 60*24602Szliu #define ln10hi (*(double*)ln10hix) 61*24602Szliu #else /* IEEE double */ 62*24602Szliu static double 63*24602Szliu ivln10 = 4.3429448190325181667E-1 ; /*Hex 2^ -2 * 1.BCB7B1526E50E */ 64*24602Szliu #endif 65*24602Szliu 66*24602Szliu double log10(x) 67*24602Szliu double x; 68*24602Szliu { 69*24602Szliu double log(); 70*24602Szliu 71*24602Szliu #ifdef VAX 72*24602Szliu return(log(x)/ln10hi); 73*24602Szliu #else /* IEEE double */ 74*24602Szliu return(ivln10*log(x)); 75*24602Szliu #endif 76*24602Szliu } 77