![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to log.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [CSRG BSD Unix] / 43BSDReno / lib / libm / common_source|
Return to log.c CVS log | Up to [CSRG BSD Unix] / 43BSDReno / lib / libm / common_source |
1.1 root 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 are permitted
6: * provided that: (1) source distributions retain this entire copyright
7: * notice and comment, and (2) distributions including binaries display
8: * the following acknowledgement: ``This product includes software
9: * developed by the University of California, Berkeley and its contributors''
10: * in the documentation or other materials provided with the distribution
11: * and in all advertising materials mentioning features or use of this
12: * software. Neither the name of the University nor the names of its
13: * contributors may be used to endorse or promote products derived
14: * from this software without specific prior written permission.
15: * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
16: * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
17: * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
18: *
19: * All recipients should regard themselves as participants in an ongoing
20: * research project and hence should feel obligated to report their
21: * experiences (good or bad) with these elementary function codes, using
22: * the sendbug(8) program, to the authors.
23: */
24:
25: #ifndef lint
26: static char sccsid[] = "@(#)log.c 5.5 (Berkeley) 6/1/90";
27: #endif /* not lint */
28:
29: /* LOG(X)
30: * RETURN THE LOGARITHM OF x
31: * DOUBLE PRECISION (VAX D FORMAT 56 bits or IEEE DOUBLE 53 BITS)
32: * CODED IN C BY K.C. NG, 1/19/85;
33: * REVISED BY K.C. NG on 2/7/85, 3/7/85, 3/24/85, 4/16/85.
34: *
35: * Required system supported functions:
36: * scalb(x,n)
37: * copysign(x,y)
38: * logb(x)
39: * finite(x)
40: *
41: * Required kernel function:
42: * log__L(z)
43: *
44: * Method :
45: * 1. Argument Reduction: find k and f such that
46: * x = 2^k * (1+f),
47: * where sqrt(2)/2 < 1+f < sqrt(2) .
48: *
49: * 2. Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
50: * = 2s + 2/3 s**3 + 2/5 s**5 + .....,
51: * log(1+f) is computed by
52: *
53: * log(1+f) = 2s + s*log__L(s*s)
54: * where
55: * log__L(z) = z*(L1 + z*(L2 + z*(... (L6 + z*L7)...)))
56: *
57: * See log__L() for the values of the coefficients.
58: *
59: * 3. Finally, log(x) = k*ln2 + log(1+f). (Here n*ln2 will be stored
60: * in two floating point number: n*ln2hi + n*ln2lo, n*ln2hi is exact
61: * since the last 20 bits of ln2hi is 0.)
62: *
63: * Special cases:
64: * log(x) is NaN with signal if x < 0 (including -INF) ;
65: * log(+INF) is +INF; log(0) is -INF with signal;
66: * log(NaN) is that NaN with no signal.
67: *
68: * Accuracy:
69: * log(x) returns the exact log(x) nearly rounded. In a test run with
70: * 1,536,000 random arguments on a VAX, the maximum observed error was
71: * .826 ulps (units in the last place).
72: *
73: * Constants:
74: * The hexadecimal values are the intended ones for the following constants.
75: * The decimal values may be used, provided that the compiler will convert
76: * from decimal to binary accurately enough to produce the hexadecimal values
77: * shown.
78: */
79:
80: #include <errno.h>
81: #include "mathimpl.h"
82:
83: vc(ln2hi, 6.9314718055829871446E-1 ,7217,4031,0000,f7d0, 0, .B17217F7D00000)
84: vc(ln2lo, 1.6465949582897081279E-12 ,bcd5,2ce7,d9cc,e4f1, -39, .E7BCD5E4F1D9CC)
85: vc(sqrt2, 1.4142135623730950622E0 ,04f3,40b5,de65,33f9, 1, .B504F333F9DE65)
86:
87: ic(ln2hi, 6.9314718036912381649E-1, -1, 1.62E42FEE00000)
88: ic(ln2lo, 1.9082149292705877000E-10, -33, 1.A39EF35793C76)
89: ic(sqrt2, 1.4142135623730951455E0, 0, 1.6A09E667F3BCD)
90:
91: #ifdef vccast
92: #define ln2hi vccast(ln2hi)
93: #define ln2lo vccast(ln2lo)
94: #define sqrt2 vccast(sqrt2)
95: #endif
96:
97:
98: double log(x)
99: double x;
100: {
101: const static double zero=0.0, negone= -1.0, half=1.0/2.0;
102: double s,z,t;
103: int k,n;
104:
105: #if !defined(vax)&&!defined(tahoe)
106: if(x!=x) return(x); /* x is NaN */
107: #endif /* !defined(vax)&&!defined(tahoe) */
108: if(finite(x)) {
109: if( x > zero ) {
110:
111: /* argument reduction */
112: k=logb(x); x=scalb(x,-k);
113: if(k == -1022) /* subnormal no. */
114: {n=logb(x); x=scalb(x,-n); k+=n;}
115: if(x >= sqrt2 ) {k += 1; x *= half;}
116: x += negone ;
117:
118: /* compute log(1+x) */
119: s=x/(2+x); t=x*x*half;
120: z=k*ln2lo+s*(t+log__L(s*s));
121: x += (z - t) ;
122:
123: return(k*ln2hi+x);
124: }
125: /* end of if (x > zero) */
126:
127: else {
128: #if defined(vax)||defined(tahoe)
129: if ( x == zero )
130: return (infnan(-ERANGE)); /* -INF */
131: else
132: return (infnan(EDOM)); /* NaN */
133: #else /* defined(vax)||defined(tahoe) */
134: /* zero argument, return -INF with signal */
135: if ( x == zero )
136: return( negone/zero );
137:
138: /* negative argument, return NaN with signal */
139: else
140: return ( zero / zero );
141: #endif /* defined(vax)||defined(tahoe) */
142: }
143: }
144: /* end of if (finite(x)) */
145: /* NOTREACHED if defined(vax)||defined(tahoe) */
146:
147: /* log(-INF) is NaN with signal */
148: else if (x<0)
149: return(zero/zero);
150:
151: /* log(+INF) is +INF */
152: else return(x);
153:
154: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.