1 /*
2 log.c -- Logarithms.
3
4 Copyright (C) 1994-2003 K. Scott Hunziker.
5 Copyright (C) 1990-1994 The Boeing Company.
6
7 See the file COPYING for license, warranty, and permission details.
8 */
9
10 static char rcsid[] =
11 "$Id: log.c,v 1.5 2003/08/01 04:57:47 ksh Exp $";
12
13 #include "log.h"
14 #include "entity.h"
15 #include "scalar.h"
16 #include "vector.h"
17 #include "matrix.h"
18 #include "apply.h"
19 #include "cast.h"
20 #include "dense.h"
21 #include "abs.h"
22
23 /* On some lame systems (AIX and HPUX), log(0) does not set errno. */
24
25 #if LAX_LOG
26 double
xlog(x)27 xlog (x)
28 double x;
29 {
30 if (!x) errno = ERANGE;
31 return log (x);
32 }
33 #endif
34
35 COMPLEX
log_complex(z)36 log_complex (z)
37 COMPLEX z;
38 {
39 /*
40 * Complex logarithm.
41 *
42 * This function was adapted from the f2c "z_log" function.
43 * Copyright (C) 1997, 1998, 2000 Lucent Technologies
44 */
45
46 COMPLEX w;
47 REAL zr = z.real;
48 REAL zi = z.imag;
49 #if 0
50 REAL s, s0, t, t2, u, v;
51 #endif
52
53 errno = 0;
54
55 w.imag = atan2 (zi, zr);
56
57 /*
58 * I've reverted to the old calculation for the real part. The new code
59 * from f2c appears to get botched with "gcc -O2". It would be nice to
60 * track this down.
61 */
62
63 w.real = xlog (abs_complex (z));
64
65 #if 0
66 if (zi < 0) zi = -zi;
67 if (zr < 0) zr = -zr;
68 if (zr < zi)
69 {
70 t = zi;
71 zi = zr;
72 zr = t;
73 }
74 t = zi/zr;
75 s = zr * sqrt (1 + t*t);
76
77 /* now s = f__cabs(zi,zr), and zr = |zr| >= |zi| = zi */
78
79 if ((t = s - 1) < 0) t = -t;
80 if (t > .01)
81 {
82 w.real = xlog (s);
83 }
84 else
85 {
86 /*
87 * log(1+x) = x - x^2/2 + x^3/3 - x^4/4 + - ...
88 * = x(1 - x/2 + x^2/3 -+...)
89 *
90 * [sqrt(y^2 + z^2) - 1] * [sqrt(y^2 + z^2) + 1] = y^2 + z^2 - 1, so
91 * sqrt(y^2 + z^2) - 1 = (y^2 + z^2 - 1) / [sqrt(y^2 + z^2) + 1]
92 */
93
94 t = ((zr*zr - 1.0) + zi*zi) / (s + 1);
95 t2 = t*t;
96 s = 1.0 - 0.5*t;
97 u = v = 1;
98 do {
99 s0 = s;
100 u *= t2;
101 v += 2;
102 s += u/v - t*u/(v+1);
103 } while (s > s0);
104
105 w.real = s*t;
106 }
107 #endif
108
109 CHECK_MATH ();
110
111 return w;
112 }
113
114 ENTITY *
bi_log(p)115 bi_log (p)
116 ENTITY *p;
117 {
118 EASSERT (p, 0, 0);
119
120 switch (p->class)
121 {
122 case scalar:
123 return (log_scalar ((SCALAR *) p));
124 case vector:
125 return (log_vector ((VECTOR *) p));
126 case matrix:
127 return (log_matrix ((MATRIX *) p));
128 default:
129 BAD_CLASS (p->class);
130 delete_entity (p);
131 raise_exception ();
132 }
133 }
134
135 ENTITY *
log_scalar(p)136 log_scalar (p)
137 SCALAR *p;
138 {
139 SCALAR * volatile s = NULL;
140
141 EASSERT (p, scalar, 0);
142
143 WITH_HANDLING
144 {
145 switch (p->type)
146 {
147 case integer:
148 if (p->v.integer < 0)
149 {
150 s = (SCALAR *) cast_scalar_integer_complex (p);
151 s->v.complex.real = xlog (-s->v.real);
152 s->v.complex.imag = atan2 (0.0, -1.0);
153 }
154 else
155 {
156 s = (SCALAR *) cast_scalar_integer_real (p);
157 errno = 0;
158 s->v.real = xlog (s->v.real);
159 CHECK_MATH ();
160 }
161 break;
162 case real:
163 if (p->v.real < 0.0)
164 {
165 s = (SCALAR *) cast_scalar_real_complex (p);
166 s->v.complex.real = xlog (-s->v.real);
167 s->v.complex.imag = atan2 (0.0, -1.0);
168 }
169 else
170 {
171 s = (SCALAR *) dup_scalar (p);
172 errno = 0;
173 s->v.real = xlog (s->v.real);
174 CHECK_MATH ();
175 }
176 break;
177
178 case complex:
179 s = (SCALAR *) dup_scalar (p);
180 errno = 0;
181 s->v.complex = log_complex (p->v.complex);
182 CHECK_MATH ();
183 break;
184
185 default:
186 BAD_TYPE (p->type);
187 delete_scalar (p);
188 raise_exception ();
189 }
190 }
191 ON_EXCEPTION
192 {
193 delete_scalar (s);
194 }
195 END_EXCEPTION;
196
197 return (ENT (s));
198 }
199
200 ENTITY *
log_vector(p)201 log_vector (p)
202 VECTOR *p;
203 {
204 VECTOR *m;
205 int i;
206 int neg = 0;
207
208 EASSERT (p, vector, 0);
209
210 errno = 0;
211
212 switch (p->type)
213 {
214 case integer:
215 for (i = 0; i < p->nn; i++)
216 if (p->a.integer[i] < 0)
217 {
218 neg++;
219 break;
220 }
221 if (neg)
222 {
223 p = (VECTOR *) cast_vector_integer_complex (p);
224 p = (VECTOR *) dense_vector (p);
225 m = (VECTOR *) apply_vector_complex_complex (log_complex, p);
226 CHECK_MATH ();
227 }
228 else
229 {
230 p = (VECTOR *) cast_vector_integer_real (p);
231 p = (VECTOR *) dense_vector (p);
232 m = (VECTOR *) apply_vector_real_real (xlog, p);
233 CHECK_MATH ();
234 }
235 break;
236
237 case real:
238 for (i = 0; i < p->nn; i++)
239 if (p->a.real[i] < 0.0)
240 {
241 neg++;
242 break;
243 }
244 if (neg)
245 {
246 p = (VECTOR *) cast_vector_real_complex (p);
247 p = (VECTOR *) dense_vector (p);
248 m = (VECTOR *) apply_vector_complex_complex (log_complex, p);
249 CHECK_MATH ();
250 }
251 else
252 {
253 p = (VECTOR *) dense_vector (p);
254 m = (VECTOR *) apply_vector_real_real (xlog, p);
255 CHECK_MATH ();
256 }
257 break;
258
259 case complex:
260 p = (VECTOR *) dense_vector (p);
261 m = (VECTOR *) apply_vector_complex_complex (log_complex, p);
262 CHECK_MATH ();
263 break;
264
265 default:
266 BAD_TYPE (p->type);
267 delete_vector (p);
268 raise_exception ();
269 }
270 return (ENT (m));
271 }
272
273 ENTITY *
log_matrix(p)274 log_matrix (p)
275 MATRIX *p;
276 {
277 MATRIX *m;
278 int i;
279 int neg = 0;
280
281 EASSERT (p, matrix, 0);
282
283 errno = 0;
284
285 switch (p->type)
286 {
287 case integer:
288 for (i = 0; i < p->nn; i++)
289 if (p->a.integer[i] < 0)
290 {
291 neg++;
292 break;
293 }
294 if (p->d.integer != NULL && !neg)
295 for (i = 0; i < p->nr; i++)
296 if (p->d.integer[i] < 0)
297 {
298 neg++;
299 break;
300 }
301 if (neg)
302 {
303 p = (MATRIX *) cast_matrix_integer_complex (p);
304 p = (MATRIX *) dense_matrix (p);
305 m = (MATRIX *) apply_matrix_complex_complex (log_complex, p);
306 CHECK_MATH ();
307 }
308 else
309 {
310 p = (MATRIX *) cast_matrix_integer_real (p);
311 p = (MATRIX *) dense_matrix (p);
312 m = (MATRIX *) apply_matrix_real_real (xlog, p);
313 CHECK_MATH ();
314 }
315 break;
316
317 case real:
318 for (i = 0; i < p->nn; i++)
319 if (p->a.real[i] < 0)
320 {
321 neg++;
322 break;
323 }
324 if (p->d.real != NULL && !neg)
325 for (i = 0; i < p->nr; i++)
326 if (p->d.real[i] < 0)
327 {
328 neg++;
329 break;
330 }
331 if (neg)
332 {
333 p = (MATRIX *) cast_matrix_real_complex (p);
334 p = (MATRIX *) dense_matrix (p);
335 m = (MATRIX *) apply_matrix_complex_complex (log_complex, p);
336 CHECK_MATH ();
337 }
338 else
339 {
340 p = (MATRIX *) dense_matrix (p);
341 m = (MATRIX *) apply_matrix_real_real (xlog, p);
342 CHECK_MATH ();
343 }
344 break;
345
346 case complex:
347 p = (MATRIX *) dense_matrix (p);
348 m = (MATRIX *) apply_matrix_complex_complex (log_complex, p);
349 CHECK_MATH ();
350 break;
351
352 default:
353 BAD_TYPE (p->type);
354 delete_matrix (p);
355 raise_exception ();
356 }
357 return (ENT (m));
358 }
359