1 /*
2 * Copyright (c) 1985 Regents of the University of California.
3 *
4 * Use and reproduction of this software are granted in accordance with
5 * the terms and conditions specified in the Berkeley Software License
6 * Agreement (in particular, this entails acknowledgement of the programs'
7 * source, and inclusion of this notice) with the additional understanding
8 * that all recipients should regard themselves as participants in an
9 * ongoing research project and hence should feel obligated to report
10 * their experiences (good or bad) with these elementary function codes,
11 * using "sendbug 4bsd-bugs@BERKELEY", to the authors.
12 */
13
14 #ifndef lint
15 static char sccsid[] = "@(#)log__L.c 1.2 (Berkeley) 08/21/85";
16 #endif not lint
17
18 /* log__L(Z)
19 * LOG(1+X) - 2S X
20 * RETURN --------------- WHERE Z = S*S, S = ------- , 0 <= Z <= .0294...
21 * S 2 + X
22 *
23 * DOUBLE PRECISION (VAX D FORMAT 56 bits or IEEE DOUBLE 53 BITS)
24 * KERNEL FUNCTION FOR LOG; TO BE USED IN LOG1P, LOG, AND POW FUNCTIONS
25 * CODED IN C BY K.C. NG, 1/19/85;
26 * REVISED BY K.C. Ng, 2/3/85, 4/16/85.
27 *
28 * Method :
29 * 1. Polynomial approximation: let s = x/(2+x).
30 * Based on log(1+x) = log(1+s) - log(1-s)
31 * = 2s + 2/3 s**3 + 2/5 s**5 + .....,
32 *
33 * (log(1+x) - 2s)/s is computed by
34 *
35 * z*(L1 + z*(L2 + z*(... (L7 + z*L8)...)))
36 *
37 * where z=s*s. (See the listing below for Lk's values.) The
38 * coefficients are obtained by a special Remez algorithm.
39 *
40 * Accuracy:
41 * Assuming no rounding error, the maximum magnitude of the approximation
42 * error (absolute) is 2**(-58.49) for IEEE double, and 2**(-63.63)
43 * for VAX D format.
44 *
45 * Constants:
46 * The hexadecimal values are the intended ones for the following constants.
47 * The decimal values may be used, provided that the compiler will convert
48 * from decimal to binary accurately enough to produce the hexadecimal values
49 * shown.
50 */
51
52 #ifdef VAX /* VAX D format (56 bits) */
53 /* static double */
54 /* L1 = 6.6666666666666703212E-1 , Hex 2^ 0 * .AAAAAAAAAAAAC5 */
55 /* L2 = 3.9999999999970461961E-1 , Hex 2^ -1 * .CCCCCCCCCC2684 */
56 /* L3 = 2.8571428579395698188E-1 , Hex 2^ -1 * .92492492F85782 */
57 /* L4 = 2.2222221233634724402E-1 , Hex 2^ -2 * .E38E3839B7AF2C */
58 /* L5 = 1.8181879517064680057E-1 , Hex 2^ -2 * .BA2EB4CC39655E */
59 /* L6 = 1.5382888777946145467E-1 , Hex 2^ -2 * .9D8551E8C5781D */
60 /* L7 = 1.3338356561139403517E-1 , Hex 2^ -2 * .8895B3907FCD92 */
61 /* L8 = 1.2500000000000000000E-1 , Hex 2^ -2 * .80000000000000 */
62 static long L1x[] = { 0xaaaa402a, 0xaac5aaaa};
63 static long L2x[] = { 0xcccc3fcc, 0x2684cccc};
64 static long L3x[] = { 0x49243f92, 0x578292f8};
65 static long L4x[] = { 0x8e383f63, 0xaf2c39b7};
66 static long L5x[] = { 0x2eb43f3a, 0x655ecc39};
67 static long L6x[] = { 0x85513f1d, 0x781de8c5};
68 static long L7x[] = { 0x95b33f08, 0xcd92907f};
69 static long L8x[] = { 0x00003f00, 0x00000000};
70 #define L1 (*(double*)L1x)
71 #define L2 (*(double*)L2x)
72 #define L3 (*(double*)L3x)
73 #define L4 (*(double*)L4x)
74 #define L5 (*(double*)L5x)
75 #define L6 (*(double*)L6x)
76 #define L7 (*(double*)L7x)
77 #define L8 (*(double*)L8x)
78 #else /* IEEE double */
79 static double
80 L1 = 6.6666666666667340202E-1 , /*Hex 2^ -1 * 1.5555555555592 */
81 L2 = 3.9999999999416702146E-1 , /*Hex 2^ -2 * 1.999999997FF24 */
82 L3 = 2.8571428742008753154E-1 , /*Hex 2^ -2 * 1.24924941E07B4 */
83 L4 = 2.2222198607186277597E-1 , /*Hex 2^ -3 * 1.C71C52150BEA6 */
84 L5 = 1.8183562745289935658E-1 , /*Hex 2^ -3 * 1.74663CC94342F */
85 L6 = 1.5314087275331442206E-1 , /*Hex 2^ -3 * 1.39A1EC014045B */
86 L7 = 1.4795612545334174692E-1 ; /*Hex 2^ -3 * 1.2F039F0085122 */
87 #endif
88
log__L(z)89 double log__L(z)
90 double z;
91 {
92 #ifdef VAX
93 return(z*(L1+z*(L2+z*(L3+z*(L4+z*(L5+z*(L6+z*(L7+z*L8))))))));
94 #else /* IEEE double */
95 return(z*(L1+z*(L2+z*(L3+z*(L4+z*(L5+z*(L6+z*L7)))))));
96 #endif
97 }
98