1 /* complex/math.c (from the GSL 1.1.1)
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Jorma Olavi T�htinen, Brian Gough
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or (at
8 * your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, see <https://www.gnu.org/licenses/>.
17 */
18
19 /* Basic complex arithmetic functions
20
21 * Original version by Jorma Olavi T�htinen <jotahtin@cc.hut.fi>
22 *
23 * Modified for GSL by Brian Gough, 3/2000
24 */
25
26 /* The following references describe the methods used in these
27 * functions,
28 *
29 * T. E. Hull and Thomas F. Fairgrieve and Ping Tak Peter Tang,
30 * "Implementing Complex Elementary Functions Using Exception
31 * Handling", ACM Transactions on Mathematical Software, Volume 20
32 * (1994), pp 215-244, Corrigenda, p553
33 *
34 * Hull et al, "Implementing the complex arcsin and arccosine
35 * functions using exception handling", ACM Transactions on
36 * Mathematical Software, Volume 23 (1997) pp 299-335
37 *
38 * Abramowitz and Stegun, Handbook of Mathematical Functions, "Inverse
39 * Circular Functions in Terms of Real and Imaginary Parts", Formulas
40 * 4.4.37, 4.4.38, 4.4.39
41 */
42
43 /*
44 * Gnumeric specific modifications written by Jukka-Pekka Iivonen
45 * (jiivonen@hutcs.cs.hut.fi)
46 *
47 * long double modifications by Morten Welinder.
48 */
49
50 #include <gnumeric-config.h>
51 #include <glib/gi18n-lib.h>
52 #include <gnumeric.h>
53 #include <func.h>
54
55 #include <complex.h>
56 #include "gsl-complex.h"
57 #include <parse-util.h>
58 #include <cell.h>
59 #include <expr.h>
60 #include <value.h>
61 #include <mathfunc.h>
62 #include <sf-trig.h>
63
64
65 #define GSL_REAL(x) GNM_CRE(*(x))
66 #define GSL_IMAG(x) GNM_CIM(*(x))
67
68 /***********************************************************************
69 * Complex arithmetic operators
70 ***********************************************************************/
71
72 static inline void
gsl_complex_mul_imag(gnm_complex const * a,gnm_float y,gnm_complex * res)73 gsl_complex_mul_imag (gnm_complex const *a, gnm_float y, gnm_complex *res)
74 { /* z=a*iy */
75 *res = GNM_CMAKE (-y * GSL_IMAG (a), y * GSL_REAL (a));
76 }
77
78 void
gsl_complex_inverse(gnm_complex const * a,gnm_complex * res)79 gsl_complex_inverse (gnm_complex const *a, gnm_complex *res)
80 { /* z=1/a */
81 gnm_float s = 1.0 / GNM_CABS (*a);
82
83 *res = GNM_CMAKE ((GSL_REAL (a) * s) * s, -(GSL_IMAG (a) * s) * s);
84 }
85
86 /**********************************************************************
87 * Inverse Complex Trigonometric Functions
88 **********************************************************************/
89
90 static void
gsl_complex_arcsin_real(gnm_float a,gnm_complex * res)91 gsl_complex_arcsin_real (gnm_float a, gnm_complex *res)
92 { /* z = arcsin(a) */
93 if (gnm_abs (a) <= 1.0) {
94 *res = GNM_CMAKE (gnm_asin (a), 0.0);
95 } else {
96 if (a < 0.0) {
97 *res = GNM_CMAKE (-M_PI_2gnum, gnm_acosh (-a));
98 } else {
99 *res = GNM_CMAKE (M_PI_2gnum, -gnm_acosh (a));
100 }
101 }
102 }
103
104 void
gsl_complex_arcsin(gnm_complex const * a,gnm_complex * res)105 gsl_complex_arcsin (gnm_complex const *a, gnm_complex *res)
106 { /* z = arcsin(a) */
107 gnm_float R = GSL_REAL (a), I = GSL_IMAG (a);
108
109 if (I == 0) {
110 gsl_complex_arcsin_real (R, res);
111 } else {
112 gnm_float x = gnm_abs (R), y = gnm_abs (I);
113 gnm_float r = gnm_hypot (x + 1, y);
114 gnm_float s = gnm_hypot (x - 1, y);
115 gnm_float A = 0.5 * (r + s);
116 gnm_float B = x / A;
117 gnm_float y2 = y * y;
118
119 gnm_float real, imag;
120
121 const gnm_float A_crossover = 1.5, B_crossover = 0.6417;
122
123 if (B <= B_crossover) {
124 real = gnm_asin (B);
125 } else {
126 if (x <= 1) {
127 gnm_float D = 0.5 * (A + x) *
128 (y2 / (r + x + 1) + (s + (1 - x)));
129 real = gnm_atan (x / gnm_sqrt (D));
130 } else {
131 gnm_float Apx = A + x;
132 gnm_float D = 0.5 * (Apx / (r + x + 1)
133 + Apx / (s + (x - 1)));
134 real = gnm_atan (x / (y * gnm_sqrt (D)));
135 }
136 }
137
138 if (A <= A_crossover) {
139 gnm_float Am1;
140
141 if (x < 1) {
142 Am1 = 0.5 * (y2 / (r + (x + 1)) + y2 /
143 (s + (1 - x)));
144 } else {
145 Am1 = 0.5 * (y2 / (r + (x + 1)) +
146 (s + (x - 1)));
147 }
148
149 imag = gnm_log1p (Am1 + gnm_sqrt (Am1 * (A + 1)));
150 } else {
151 imag = gnm_log (A + gnm_sqrt (A * A - 1));
152 }
153
154 *res = GNM_CMAKE ((R >= 0) ? real : -real, (I >= 0) ?
155 imag : -imag);
156 }
157 }
158
159 static void
gsl_complex_arccos_real(gnm_float a,gnm_complex * res)160 gsl_complex_arccos_real (gnm_float a, gnm_complex *res)
161 { /* z = arccos(a) */
162 if (gnm_abs (a) <= 1.0) {
163 *res = GNM_CMAKE (gnm_acos (a), 0);
164 } else {
165 if (a < 0.0) {
166 *res = GNM_CMAKE (M_PIgnum, -gnm_acosh (-a));
167 } else {
168 *res = GNM_CMAKE (0, gnm_acosh (a));
169 }
170 }
171 }
172
173 void
gsl_complex_arccos(gnm_complex const * a,gnm_complex * res)174 gsl_complex_arccos (gnm_complex const *a, gnm_complex *res)
175 { /* z = arccos(a) */
176 gnm_float R = GSL_REAL (a), I = GSL_IMAG (a);
177
178 if (I == 0) {
179 gsl_complex_arccos_real (R, res);
180 } else {
181 gnm_float x = gnm_abs (R);
182 gnm_float y = gnm_abs (I);
183 gnm_float r = gnm_hypot (x + 1, y);
184 gnm_float s = gnm_hypot (x - 1, y);
185 gnm_float A = 0.5 * (r + s);
186 gnm_float B = x / A;
187 gnm_float y2 = y * y;
188
189 gnm_float real, imag;
190
191 const gnm_float A_crossover = 1.5;
192 const gnm_float B_crossover = 0.6417;
193
194 if (B <= B_crossover) {
195 real = gnm_acos (B);
196 } else {
197 if (x <= 1) {
198 gnm_float D = 0.5 * (A + x) *
199 (y2 / (r + x + 1) + (s + (1 - x)));
200 real = gnm_atan (gnm_sqrt (D) / x);
201 } else {
202 gnm_float Apx = A + x;
203 gnm_float D = 0.5 * (Apx / (r + x + 1) + Apx /
204 (s + (x - 1)));
205 real = gnm_atan ((y * gnm_sqrt (D)) / x);
206 }
207 }
208 if (A <= A_crossover) {
209 gnm_float Am1;
210
211 if (x < 1) {
212 Am1 = 0.5 * (y2 / (r + (x + 1)) + y2 /
213 (s + (1 - x)));
214 } else {
215 Am1 = 0.5 * (y2 / (r + (x + 1)) +
216 (s + (x - 1)));
217 }
218
219 imag = gnm_log1p (Am1 + gnm_sqrt (Am1 * (A + 1)));
220 } else {
221 imag = gnm_log (A + gnm_sqrt (A * A - 1));
222 }
223
224 *res = GNM_CMAKE ((R >= 0) ? real : M_PIgnum - real, (I >= 0) ?
225 -imag : imag);
226 }
227 }
228
229 void
gsl_complex_arctan(gnm_complex const * a,gnm_complex * res)230 gsl_complex_arctan (gnm_complex const *a, gnm_complex *res)
231 { /* z = arctan(a) */
232 gnm_float R = GSL_REAL (a), I = GSL_IMAG (a);
233
234 if (I == 0) {
235 *res = GNM_CMAKE (gnm_atan (R), 0);
236 } else {
237 /* FIXME: This is a naive implementation which does not fully
238 * take into account cancellation errors, overflow, underflow
239 * etc. It would benefit from the Hull et al treatment. */
240
241 gnm_float r = gnm_hypot (R, I);
242
243 gnm_float imag;
244
245 gnm_float u = 2 * I / (1 + r * r);
246
247 /* FIXME: the following cross-over should be optimized but 0.1
248 * seems to work ok */
249
250 if (gnm_abs (u) < 0.1) {
251 imag = 0.25 * (gnm_log1p (u) - gnm_log1p (-u));
252 } else {
253 gnm_float A = gnm_hypot (R, I + 1);
254 gnm_float B = gnm_hypot (R, I - 1);
255 imag = 0.5 * gnm_log (A / B);
256 }
257 if (R == 0) {
258 if (I > 1) {
259 *res = GNM_CMAKE (M_PI_2gnum, imag);
260 } else if (I < -1) {
261 *res = GNM_CMAKE (-M_PI_2gnum, imag);
262 } else {
263 *res = GNM_CMAKE (0, imag);
264 }
265 } else {
266 *res = GNM_CMAKE (0.5 * gnm_atan2 (2 * R,
267 ((1 + r) * (1 - r))),
268 imag);
269 }
270 }
271 }
272
273 void
gsl_complex_arcsec(gnm_complex const * a,gnm_complex * res)274 gsl_complex_arcsec (gnm_complex const *a, gnm_complex *res)
275 { /* z = arcsec(a) */
276 gsl_complex_inverse (a, res);
277 gsl_complex_arccos (res, res);
278 }
279
280 void
gsl_complex_arccsc(gnm_complex const * a,gnm_complex * res)281 gsl_complex_arccsc (gnm_complex const *a, gnm_complex *res)
282 { /* z = arccsc(a) */
283 gsl_complex_inverse (a, res);
284 gsl_complex_arcsin (res, res);
285 }
286
287 void
gsl_complex_arccot(gnm_complex const * a,gnm_complex * res)288 gsl_complex_arccot (gnm_complex const *a, gnm_complex *res)
289 { /* z = arccot(a) */
290 if (GSL_REAL (a) == 0.0 && GSL_IMAG (a) == 0.0) {
291 *res = GNM_CMAKE (M_PI_2gnum, 0);
292 } else {
293 gsl_complex_inverse (a, res);
294 gsl_complex_arctan (res, res);
295 }
296 }
297
298 /**********************************************************************
299 * Complex Hyperbolic Functions
300 **********************************************************************/
301
302 void
gsl_complex_sinh(gnm_complex const * a,gnm_complex * res)303 gsl_complex_sinh (gnm_complex const *a, gnm_complex *res)
304 { /* z = sinh(a) */
305 gnm_float R = GSL_REAL (a), I = GSL_IMAG (a);
306
307 *res = GNM_CMAKE (gnm_sinh (R) * gnm_cos (I), gnm_cosh (R) * gnm_sin (I));
308 }
309
310 void
gsl_complex_cosh(gnm_complex const * a,gnm_complex * res)311 gsl_complex_cosh (gnm_complex const *a, gnm_complex *res)
312 { /* z = cosh(a) */
313 gnm_float R = GSL_REAL (a), I = GSL_IMAG (a);
314
315 *res = GNM_CMAKE (gnm_cosh (R) * gnm_cos (I), gnm_sinh (R) * gnm_sin (I));
316 }
317
318 void
gsl_complex_tanh(gnm_complex const * a,gnm_complex * res)319 gsl_complex_tanh (gnm_complex const *a, gnm_complex *res)
320 { /* z = tanh(a) */
321 gnm_float R = GSL_REAL (a), I = GSL_IMAG (a);
322
323 if (gnm_abs (R) < 1.0) {
324 gnm_float D =
325 gnm_pow (gnm_cos (I), 2.0) +
326 gnm_pow (gnm_sinh (R), 2.0);
327
328 *res = GNM_CMAKE (gnm_sinh (R) * gnm_cosh (R) / D,
329 0.5 * gnm_sin (2 * I) / D);
330 } else {
331 gnm_float D =
332 gnm_pow (gnm_cos (I), 2.0) +
333 gnm_pow (gnm_sinh (R), 2.0);
334 gnm_float F = 1 + gnm_pow (gnm_cos (I) / gnm_sinh (R), 2.0);
335
336 *res = GNM_CMAKE (1.0 / (gnm_tanh (R) * F),
337 0.5 * gnm_sin (2 * I) / D);
338 }
339 }
340
341 void
gsl_complex_sech(gnm_complex const * a,gnm_complex * res)342 gsl_complex_sech (gnm_complex const *a, gnm_complex *res)
343 { /* z = sech(a) */
344 gsl_complex_cosh (a, res);
345 gsl_complex_inverse (res, res);
346 }
347
348 void
gsl_complex_csch(gnm_complex const * a,gnm_complex * res)349 gsl_complex_csch (gnm_complex const *a, gnm_complex *res)
350 { /* z = csch(a) */
351 gsl_complex_sinh (a, res);
352 gsl_complex_inverse (res, res);
353 }
354
355 void
gsl_complex_coth(gnm_complex const * a,gnm_complex * res)356 gsl_complex_coth (gnm_complex const *a, gnm_complex *res)
357 { /* z = coth(a) */
358 gsl_complex_tanh (a, res);
359 gsl_complex_inverse (res, res);
360 }
361
362 /**********************************************************************
363 * Inverse Complex Hyperbolic Functions
364 **********************************************************************/
365
366 void
gsl_complex_arcsinh(gnm_complex const * a,gnm_complex * res)367 gsl_complex_arcsinh (gnm_complex const *a, gnm_complex *res)
368 { /* z = arcsinh(a) */
369 gsl_complex_mul_imag (a, 1.0, res);
370 gsl_complex_arcsin (res, res);
371 gsl_complex_mul_imag (res, -1.0, res);
372 }
373
374 void
gsl_complex_arccosh(gnm_complex const * a,gnm_complex * res)375 gsl_complex_arccosh (gnm_complex const *a, gnm_complex *res)
376 { /* z = arccosh(a) */
377 if (GSL_IMAG (a) == 0 && GSL_REAL (a) == 1)
378 *res = GNM_C0;
379 else {
380 gsl_complex_arccos (a, res);
381 gsl_complex_mul_imag (res, GSL_IMAG (res) > 0 ? -1.0 : 1.0, res);
382 }
383 }
384
385 static void
gsl_complex_arctanh_real(gnm_float a,gnm_complex * res)386 gsl_complex_arctanh_real (gnm_float a, gnm_complex *res)
387 { /* z = arctanh(a) */
388 if (a > -1.0 && a < 1.0) {
389 *res = GNM_CMAKE (gnm_atanh (a), 0);
390 } else {
391 *res = GNM_CMAKE (gnm_acoth (a),
392 (a < 0) ? M_PI_2gnum : -M_PI_2gnum);
393 }
394 }
395
396 void
gsl_complex_arctanh(gnm_complex const * a,gnm_complex * res)397 gsl_complex_arctanh (gnm_complex const *a, gnm_complex *res)
398 { /* z = arctanh(a) */
399 if (GSL_IMAG (a) == 0.0) {
400 gsl_complex_arctanh_real (GSL_REAL (a), res);
401 } else {
402 gsl_complex_mul_imag (a, 1.0, res);
403 gsl_complex_arctan (res, res);
404 gsl_complex_mul_imag (res, -1.0, res);
405 }
406 }
407
408 void
gsl_complex_arcsech(gnm_complex const * a,gnm_complex * res)409 gsl_complex_arcsech (gnm_complex const *a, gnm_complex *res)
410 { /* z = arcsech(a); */
411 gsl_complex_inverse (a, res);
412 gsl_complex_arccosh (res, res);
413 }
414
415 void
gsl_complex_arccsch(gnm_complex const * a,gnm_complex * res)416 gsl_complex_arccsch (gnm_complex const *a, gnm_complex *res)
417 { /* z = arccsch(a); */
418 gsl_complex_inverse (a, res);
419 gsl_complex_arcsinh (res, res);
420 }
421
422 void
gsl_complex_arccoth(gnm_complex const * a,gnm_complex * res)423 gsl_complex_arccoth (gnm_complex const *a, gnm_complex *res)
424 { /* z = arccoth(a); */
425 gsl_complex_inverse (a, res);
426 gsl_complex_arctanh (res, res);
427 }
428