1 /* complex/math.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Jorma Olavi T�htinen, Brian
4 * Gough
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 3 of the License, or (at
9 * your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful, but
12 * WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301,
19 * USA.
20 */
21
22 /* Basic complex arithmetic functions
23
24 * Original version by Jorma Olavi T�htinen <jotahtin@cc.hut.fi>
25 *
26 * Modified for GSL by Brian Gough, 3/2000
27 */
28
29 /* The following references describe the methods used in these
30 * functions,
31 *
32 * T. E. Hull and Thomas F. Fairgrieve and Ping Tak Peter Tang,
33 * "Implementing Complex Elementary Functions Using Exception
34 * Handling", ACM Transactions on Mathematical Software, Volume 20
35 * (1994), pp 215-244, Corrigenda, p553
36 *
37 * Hull et al, "Implementing the complex arcsin and arccosine
38 * functions using exception handling", ACM Transactions on
39 * Mathematical Software, Volume 23 (1997) pp 299-335
40 *
41 * Abramowitz and Stegun, Handbook of Mathematical Functions, "Inverse
42 * Circular Functions in Terms of Real and Imaginary Parts", Formulas
43 * 4.4.37, 4.4.38, 4.4.39
44 */
45
46 #include "config.h"
47 #include <math.h>
48 #include <gsl/gsl_math.h>
49 #include <gsl/gsl_complex.h>
50 #include <gsl/gsl_complex_math.h>
51
52 /**********************************************************************
53 * Complex numbers
54 **********************************************************************/
55
gsl_complex_polar(double r,double theta)56 gsl_complex gsl_complex_polar(double r, double theta)
57 { /* return z = r exp(i theta) */
58 gsl_complex z;
59 GSL_SET_COMPLEX(&z, r * cos(theta), r * sin(theta));
60 return z;
61 }
62
63 /**********************************************************************
64 * Properties of complex numbers
65 **********************************************************************/
66
gsl_complex_arg(gsl_complex z)67 double gsl_complex_arg(gsl_complex z)
68 { /* return arg(z), -pi < arg(z) <= +pi */
69 double x = GSL_REAL(z);
70 double y = GSL_IMAG(z);
71
72 if (x == 0.0 && y == 0.0) {
73 return 0;
74 }
75
76 return atan2(y, x);
77 }
78
gsl_complex_abs(gsl_complex z)79 double gsl_complex_abs(gsl_complex z)
80 { /* return |z| */
81 return hypot(GSL_REAL(z), GSL_IMAG(z));
82 }
83
gsl_complex_abs2(gsl_complex z)84 double gsl_complex_abs2(gsl_complex z)
85 { /* return |z|^2 */
86 double x = GSL_REAL(z);
87 double y = GSL_IMAG(z);
88
89 return (x * x + y * y);
90 }
91
gsl_complex_logabs(gsl_complex z)92 double gsl_complex_logabs(gsl_complex z)
93 { /* return log|z| */
94 double xabs = fabs(GSL_REAL(z));
95 double yabs = fabs(GSL_IMAG(z));
96 double max, u;
97
98 if (xabs >= yabs) {
99 max = xabs;
100 u = yabs / xabs;
101 } else {
102 max = yabs;
103 u = xabs / yabs;
104 }
105
106 /* Handle underflow when u is close to 0 */
107
108 return log(max) + 0.5 * log1p(u * u);
109 }
110
111 /***********************************************************************
112 * Complex arithmetic operators
113 ***********************************************************************/
114
gsl_complex_add(gsl_complex a,gsl_complex b)115 gsl_complex gsl_complex_add(gsl_complex a, gsl_complex b)
116 { /* z=a+b */
117 double ar = GSL_REAL(a), ai = GSL_IMAG(a);
118 double br = GSL_REAL(b), bi = GSL_IMAG(b);
119
120 gsl_complex z;
121 GSL_SET_COMPLEX(&z, ar + br, ai + bi);
122 return z;
123 }
124
gsl_complex_add_real(gsl_complex a,double x)125 gsl_complex gsl_complex_add_real(gsl_complex a, double x)
126 { /* z=a+x */
127 gsl_complex z;
128 GSL_SET_COMPLEX(&z, GSL_REAL(a) + x, GSL_IMAG(a));
129 return z;
130 }
131
gsl_complex_add_imag(gsl_complex a,double y)132 gsl_complex gsl_complex_add_imag(gsl_complex a, double y)
133 { /* z=a+iy */
134 gsl_complex z;
135 GSL_SET_COMPLEX(&z, GSL_REAL(a), GSL_IMAG(a) + y);
136 return z;
137 }
138
gsl_complex_sub(gsl_complex a,gsl_complex b)139 gsl_complex gsl_complex_sub(gsl_complex a, gsl_complex b)
140 { /* z=a-b */
141 double ar = GSL_REAL(a), ai = GSL_IMAG(a);
142 double br = GSL_REAL(b), bi = GSL_IMAG(b);
143
144 gsl_complex z;
145 GSL_SET_COMPLEX(&z, ar - br, ai - bi);
146 return z;
147 }
148
gsl_complex_sub_real(gsl_complex a,double x)149 gsl_complex gsl_complex_sub_real(gsl_complex a, double x)
150 { /* z=a-x */
151 gsl_complex z;
152 GSL_SET_COMPLEX(&z, GSL_REAL(a) - x, GSL_IMAG(a));
153 return z;
154 }
155
gsl_complex_sub_imag(gsl_complex a,double y)156 gsl_complex gsl_complex_sub_imag(gsl_complex a, double y)
157 { /* z=a-iy */
158 gsl_complex z;
159 GSL_SET_COMPLEX(&z, GSL_REAL(a), GSL_IMAG(a) - y);
160 return z;
161 }
162
gsl_complex_mul(gsl_complex a,gsl_complex b)163 gsl_complex gsl_complex_mul(gsl_complex a, gsl_complex b)
164 { /* z=a*b */
165 double ar = GSL_REAL(a), ai = GSL_IMAG(a);
166 double br = GSL_REAL(b), bi = GSL_IMAG(b);
167
168 gsl_complex z;
169 GSL_SET_COMPLEX(&z, ar * br - ai * bi, ar * bi + ai * br);
170 return z;
171 }
172
gsl_complex_mul_real(gsl_complex a,double x)173 gsl_complex gsl_complex_mul_real(gsl_complex a, double x)
174 { /* z=a*x */
175 gsl_complex z;
176 GSL_SET_COMPLEX(&z, x * GSL_REAL(a), x * GSL_IMAG(a));
177 return z;
178 }
179
gsl_complex_mul_imag(gsl_complex a,double y)180 gsl_complex gsl_complex_mul_imag(gsl_complex a, double y)
181 { /* z=a*iy */
182 gsl_complex z;
183 GSL_SET_COMPLEX(&z, -y * GSL_IMAG(a), y * GSL_REAL(a));
184 return z;
185 }
186
gsl_complex_div(gsl_complex a,gsl_complex b)187 gsl_complex gsl_complex_div(gsl_complex a, gsl_complex b)
188 { /* z=a/b */
189 double ar = GSL_REAL(a), ai = GSL_IMAG(a);
190 double br = GSL_REAL(b), bi = GSL_IMAG(b);
191
192 double s = 1.0 / gsl_complex_abs(b);
193
194 double sbr = s * br;
195 double sbi = s * bi;
196
197 double zr = (ar * sbr + ai * sbi) * s;
198 double zi = (ai * sbr - ar * sbi) * s;
199
200 gsl_complex z;
201 GSL_SET_COMPLEX(&z, zr, zi);
202 return z;
203 }
204
gsl_complex_div_real(gsl_complex a,double x)205 gsl_complex gsl_complex_div_real(gsl_complex a, double x)
206 { /* z=a/x */
207 gsl_complex z;
208 GSL_SET_COMPLEX(&z, GSL_REAL(a) / x, GSL_IMAG(a) / x);
209 return z;
210 }
211
gsl_complex_div_imag(gsl_complex a,double y)212 gsl_complex gsl_complex_div_imag(gsl_complex a, double y)
213 { /* z=a/(iy) */
214 gsl_complex z;
215 GSL_SET_COMPLEX(&z, GSL_IMAG(a) / y, -GSL_REAL(a) / y);
216 return z;
217 }
218
gsl_complex_conjugate(gsl_complex a)219 gsl_complex gsl_complex_conjugate(gsl_complex a)
220 { /* z=conj(a) */
221 gsl_complex z;
222 GSL_SET_COMPLEX(&z, GSL_REAL(a), -GSL_IMAG(a));
223 return z;
224 }
225
gsl_complex_negative(gsl_complex a)226 gsl_complex gsl_complex_negative(gsl_complex a)
227 { /* z=-a */
228 gsl_complex z;
229 GSL_SET_COMPLEX(&z, -GSL_REAL(a), -GSL_IMAG(a));
230 return z;
231 }
232
gsl_complex_inverse(gsl_complex a)233 gsl_complex gsl_complex_inverse(gsl_complex a)
234 { /* z=1/a */
235 double s = 1.0 / gsl_complex_abs(a);
236
237 gsl_complex z;
238 GSL_SET_COMPLEX(&z, (GSL_REAL(a) * s) * s, -(GSL_IMAG(a) * s) * s);
239 return z;
240 }
241
242 /**********************************************************************
243 * Elementary complex functions
244 **********************************************************************/
245
gsl_complex_sqrt(gsl_complex a)246 gsl_complex gsl_complex_sqrt(gsl_complex a)
247 { /* z=sqrt(a) */
248 gsl_complex z;
249
250 if (GSL_REAL(a) == 0.0 && GSL_IMAG(a) == 0.0) {
251 GSL_SET_COMPLEX(&z, 0, 0);
252 } else {
253 double x = fabs(GSL_REAL(a));
254 double y = fabs(GSL_IMAG(a));
255 double w;
256
257 if (x >= y) {
258 double t = y / x;
259 w = sqrt(x) * sqrt(0.5 * (1.0 + sqrt(1.0 + t * t)));
260 } else {
261 double t = x / y;
262 w = sqrt(y) * sqrt(0.5 * (t + sqrt(1.0 + t * t)));
263 }
264
265 if (GSL_REAL(a) >= 0.0) {
266 double ai = GSL_IMAG(a);
267 GSL_SET_COMPLEX(&z, w, ai / (2.0 * w));
268 } else {
269 double ai = GSL_IMAG(a);
270 double vi = (ai >= 0) ? w : -w;
271 GSL_SET_COMPLEX(&z, ai / (2.0 * vi), vi);
272 }
273 }
274
275 return z;
276 }
277
gsl_complex_sqrt_real(double x)278 gsl_complex gsl_complex_sqrt_real(double x)
279 { /* z=sqrt(x) */
280 gsl_complex z;
281
282 if (x >= 0) {
283 GSL_SET_COMPLEX(&z, sqrt(x), 0.0);
284 } else {
285 GSL_SET_COMPLEX(&z, 0.0, sqrt(-x));
286 }
287
288 return z;
289 }
290
gsl_complex_exp(gsl_complex a)291 gsl_complex gsl_complex_exp(gsl_complex a)
292 { /* z=exp(a) */
293 double rho = exp(GSL_REAL(a));
294 double theta = GSL_IMAG(a);
295
296 gsl_complex z;
297 GSL_SET_COMPLEX(&z, rho * cos(theta), rho * sin(theta));
298 return z;
299 }
300
gsl_complex_pow(gsl_complex a,gsl_complex b)301 gsl_complex gsl_complex_pow(gsl_complex a, gsl_complex b)
302 { /* z=a^b */
303 gsl_complex z;
304
305 if (GSL_REAL(a) == 0 && GSL_IMAG(a) == 0.0) {
306 if (GSL_REAL(b) == 0 && GSL_IMAG(b) == 0.0) {
307 GSL_SET_COMPLEX(&z, 1.0, 0.0);
308 } else {
309 GSL_SET_COMPLEX(&z, 0.0, 0.0);
310 }
311 } else if (GSL_REAL(b) == 1.0 && GSL_IMAG(b) == 0.0) {
312 return a;
313 } else if (GSL_REAL(b) == -1.0 && GSL_IMAG(b) == 0.0) {
314 return gsl_complex_inverse(a);
315 } else {
316 double logr = gsl_complex_logabs(a);
317 double theta = gsl_complex_arg(a);
318
319 double br = GSL_REAL(b), bi = GSL_IMAG(b);
320
321 double rho = exp(logr * br - bi * theta);
322 double beta = theta * br + bi * logr;
323
324 GSL_SET_COMPLEX(&z, rho * cos(beta), rho * sin(beta));
325 }
326
327 return z;
328 }
329
gsl_complex_pow_real(gsl_complex a,double b)330 gsl_complex gsl_complex_pow_real(gsl_complex a, double b)
331 { /* z=a^b */
332 gsl_complex z;
333
334 if (GSL_REAL(a) == 0 && GSL_IMAG(a) == 0) {
335 if (b == 0) {
336 GSL_SET_COMPLEX(&z, 1, 0);
337 } else {
338 GSL_SET_COMPLEX(&z, 0, 0);
339 }
340 } else {
341 double logr = gsl_complex_logabs(a);
342 double theta = gsl_complex_arg(a);
343 double rho = exp(logr * b);
344 double beta = theta * b;
345 GSL_SET_COMPLEX(&z, rho * cos(beta), rho * sin(beta));
346 }
347
348 return z;
349 }
350
gsl_complex_log(gsl_complex a)351 gsl_complex gsl_complex_log(gsl_complex a)
352 { /* z=log(a) */
353 double logr = gsl_complex_logabs(a);
354 double theta = gsl_complex_arg(a);
355
356 gsl_complex z;
357 GSL_SET_COMPLEX(&z, logr, theta);
358 return z;
359 }
360
gsl_complex_log10(gsl_complex a)361 gsl_complex gsl_complex_log10(gsl_complex a)
362 { /* z = log10(a) */
363 return gsl_complex_mul_real(gsl_complex_log(a), 1 / log(10.));
364 }
365
gsl_complex_log_b(gsl_complex a,gsl_complex b)366 gsl_complex gsl_complex_log_b(gsl_complex a, gsl_complex b)
367 {
368 return gsl_complex_div(gsl_complex_log(a), gsl_complex_log(b));
369 }
370
371 /***********************************************************************
372 * Complex trigonometric functions
373 ***********************************************************************/
374
gsl_complex_sin(gsl_complex a)375 gsl_complex gsl_complex_sin(gsl_complex a)
376 { /* z = sin(a) */
377 double R = GSL_REAL(a), I = GSL_IMAG(a);
378
379 gsl_complex z;
380
381 if (I == 0.0) {
382 /* avoid returning negative zero (-0.0) for the imaginary part */
383
384 GSL_SET_COMPLEX(&z, sin(R), 0.0);
385 } else {
386 GSL_SET_COMPLEX(&z, sin(R) * cosh(I), cos(R) * sinh(I));
387 }
388
389 return z;
390 }
391
gsl_complex_cos(gsl_complex a)392 gsl_complex gsl_complex_cos(gsl_complex a)
393 { /* z = cos(a) */
394 double R = GSL_REAL(a), I = GSL_IMAG(a);
395
396 gsl_complex z;
397
398 if (I == 0.0) {
399 /* avoid returning negative zero (-0.0) for the imaginary part */
400
401 GSL_SET_COMPLEX(&z, cos(R), 0.0);
402 } else {
403 GSL_SET_COMPLEX(&z, cos(R) * cosh(I), sin(R) * sinh(-I));
404 }
405
406 return z;
407 }
408
gsl_complex_tan(gsl_complex a)409 gsl_complex gsl_complex_tan(gsl_complex a)
410 { /* z = tan(a) */
411 double R = GSL_REAL(a), I = GSL_IMAG(a);
412
413 gsl_complex z;
414
415 if (fabs(I) < 1) {
416 double D = pow(cos(R), 2.0) + pow(sinh(I), 2.0);
417
418 GSL_SET_COMPLEX(&z, 0.5 * sin(2 * R) / D, 0.5 * sinh(2 * I) / D);
419 } else {
420 double D = pow(cos(R), 2.0) + pow(sinh(I), 2.0);
421 double F = 1 + pow(cos(R) / sinh(I), 2.0);
422
423 GSL_SET_COMPLEX(&z, 0.5 * sin(2 * R) / D, 1 / (tanh(I) * F));
424 }
425
426 return z;
427 }
428
gsl_complex_sec(gsl_complex a)429 gsl_complex gsl_complex_sec(gsl_complex a)
430 { /* z = sec(a) */
431 gsl_complex z = gsl_complex_cos(a);
432 return gsl_complex_inverse(z);
433 }
434
gsl_complex_csc(gsl_complex a)435 gsl_complex gsl_complex_csc(gsl_complex a)
436 { /* z = csc(a) */
437 gsl_complex z = gsl_complex_sin(a);
438 return gsl_complex_inverse(z);
439 }
440
gsl_complex_cot(gsl_complex a)441 gsl_complex gsl_complex_cot(gsl_complex a)
442 { /* z = cot(a) */
443 gsl_complex z = gsl_complex_tan(a);
444 return gsl_complex_inverse(z);
445 }
446
447 /**********************************************************************
448 * Inverse Complex Trigonometric Functions
449 **********************************************************************/
450
gsl_complex_arcsin(gsl_complex a)451 gsl_complex gsl_complex_arcsin(gsl_complex a)
452 { /* z = arcsin(a) */
453 double R = GSL_REAL(a), I = GSL_IMAG(a);
454 gsl_complex z;
455
456 if (I == 0) {
457 z = gsl_complex_arcsin_real(R);
458 } else {
459 double x = fabs(R), y = fabs(I);
460 double r = hypot(x + 1, y), s = hypot(x - 1, y);
461 double A = 0.5 * (r + s);
462 double B = x / A;
463 double y2 = y * y;
464
465 double real, imag;
466
467 const double A_crossover = 1.5, B_crossover = 0.6417;
468
469 if (B <= B_crossover) {
470 real = asin(B);
471 } else {
472 if (x <= 1) {
473 double D = 0.5 * (A + x) * (y2 / (r + x + 1) + (s + (1 - x)));
474 real = atan(x / sqrt(D));
475 } else {
476 double Apx = A + x;
477 double D = 0.5 * (Apx / (r + x + 1) + Apx / (s + (x - 1)));
478 real = atan(x / (y * sqrt(D)));
479 }
480 }
481
482 if (A <= A_crossover) {
483 double Am1;
484
485 if (x < 1) {
486 Am1 = 0.5 * (y2 / (r + (x + 1)) + y2 / (s + (1 - x)));
487 } else {
488 Am1 = 0.5 * (y2 / (r + (x + 1)) + (s + (x - 1)));
489 }
490
491 imag = log1p(Am1 + sqrt(Am1 * (A + 1)));
492 } else {
493 imag = log(A + sqrt(A * A - 1));
494 }
495
496 GSL_SET_COMPLEX(&z, (R >= 0) ? real : -real, (I >= 0) ? imag : -imag);
497 }
498
499 return z;
500 }
501
gsl_complex_arcsin_real(double a)502 gsl_complex gsl_complex_arcsin_real(double a)
503 { /* z = arcsin(a) */
504 gsl_complex z;
505
506 if (fabs(a) <= 1.0) {
507 GSL_SET_COMPLEX(&z, asin(a), 0.0);
508 } else {
509 if (a < 0.0) {
510 GSL_SET_COMPLEX(&z, -M_PI_2, acosh(-a));
511 } else {
512 GSL_SET_COMPLEX(&z, M_PI_2, -acosh(a));
513 }
514 }
515
516 return z;
517 }
518
gsl_complex_arccos(gsl_complex a)519 gsl_complex gsl_complex_arccos(gsl_complex a)
520 { /* z = arccos(a) */
521 double R = GSL_REAL(a), I = GSL_IMAG(a);
522 gsl_complex z;
523
524 if (I == 0) {
525 z = gsl_complex_arccos_real(R);
526 } else {
527 double x = fabs(R), y = fabs(I);
528 double r = hypot(x + 1, y), s = hypot(x - 1, y);
529 double A = 0.5 * (r + s);
530 double B = x / A;
531 double y2 = y * y;
532
533 double real, imag;
534
535 const double A_crossover = 1.5, B_crossover = 0.6417;
536
537 if (B <= B_crossover) {
538 real = acos(B);
539 } else {
540 if (x <= 1) {
541 double D = 0.5 * (A + x) * (y2 / (r + x + 1) + (s + (1 - x)));
542 real = atan(sqrt(D) / x);
543 } else {
544 double Apx = A + x;
545 double D = 0.5 * (Apx / (r + x + 1) + Apx / (s + (x - 1)));
546 real = atan((y * sqrt(D)) / x);
547 }
548 }
549
550 if (A <= A_crossover) {
551 double Am1;
552
553 if (x < 1) {
554 Am1 = 0.5 * (y2 / (r + (x + 1)) + y2 / (s + (1 - x)));
555 } else {
556 Am1 = 0.5 * (y2 / (r + (x + 1)) + (s + (x - 1)));
557 }
558
559 imag = log1p(Am1 + sqrt(Am1 * (A + 1)));
560 } else {
561 imag = log(A + sqrt(A * A - 1));
562 }
563
564 GSL_SET_COMPLEX(&z, (R >= 0) ? real : M_PI - real,
565 (I >= 0) ? -imag : imag);
566 }
567
568 return z;
569 }
570
gsl_complex_arccos_real(double a)571 gsl_complex gsl_complex_arccos_real(double a)
572 { /* z = arccos(a) */
573 gsl_complex z;
574
575 if (fabs(a) <= 1.0) {
576 GSL_SET_COMPLEX(&z, acos(a), 0);
577 } else {
578 if (a < 0.0) {
579 GSL_SET_COMPLEX(&z, M_PI, -acosh(-a));
580 } else {
581 GSL_SET_COMPLEX(&z, 0, acosh(a));
582 }
583 }
584
585 return z;
586 }
587
gsl_complex_arctan(gsl_complex a)588 gsl_complex gsl_complex_arctan(gsl_complex a)
589 { /* z = arctan(a) */
590 double R = GSL_REAL(a), I = GSL_IMAG(a);
591 gsl_complex z;
592
593 if (I == 0) {
594 GSL_SET_COMPLEX(&z, atan(R), 0);
595 } else {
596 /* FIXME: This is a naive implementation which does not fully
597 take into account cancellation errors, overflow, underflow
598 etc. It would benefit from the Hull et al treatment. */
599
600 double r = hypot(R, I);
601
602 double imag;
603
604 double u = 2 * I / (1 + r * r);
605
606 /* FIXME: the following cross-over should be optimized but 0.1
607 seems to work ok */
608
609 if (fabs(u) < 0.1) {
610 imag = 0.25 * (log1p(u) - log1p(-u));
611 } else {
612 double A = hypot(R, I + 1);
613 double B = hypot(R, I - 1);
614 imag = 0.5 * log(A / B);
615 }
616
617 if (R == 0) {
618 if (I > 1) {
619 GSL_SET_COMPLEX(&z, M_PI_2, imag);
620 } else if (I < -1) {
621 GSL_SET_COMPLEX(&z, -M_PI_2, imag);
622 } else {
623 GSL_SET_COMPLEX(&z, 0, imag);
624 };
625 } else {
626 GSL_SET_COMPLEX(&z, 0.5 * atan2(2 * R, ((1 + r) * (1 - r))), imag);
627 }
628 }
629
630 return z;
631 }
632
gsl_complex_arcsec(gsl_complex a)633 gsl_complex gsl_complex_arcsec(gsl_complex a)
634 { /* z = arcsec(a) */
635 gsl_complex z = gsl_complex_inverse(a);
636 return gsl_complex_arccos(z);
637 }
638
gsl_complex_arcsec_real(double a)639 gsl_complex gsl_complex_arcsec_real(double a)
640 { /* z = arcsec(a) */
641 gsl_complex z;
642
643 if (a <= -1.0 || a >= 1.0) {
644 GSL_SET_COMPLEX(&z, acos(1 / a), 0.0);
645 } else {
646 if (a >= 0.0) {
647 GSL_SET_COMPLEX(&z, 0, acosh(1 / a));
648 } else {
649 GSL_SET_COMPLEX(&z, M_PI, -acosh(-1 / a));
650 }
651 }
652
653 return z;
654 }
655
gsl_complex_arccsc(gsl_complex a)656 gsl_complex gsl_complex_arccsc(gsl_complex a)
657 { /* z = arccsc(a) */
658 gsl_complex z = gsl_complex_inverse(a);
659 return gsl_complex_arcsin(z);
660 }
661
gsl_complex_arccsc_real(double a)662 gsl_complex gsl_complex_arccsc_real(double a)
663 { /* z = arccsc(a) */
664 gsl_complex z;
665
666 if (a <= -1.0 || a >= 1.0) {
667 GSL_SET_COMPLEX(&z, asin(1 / a), 0.0);
668 } else {
669 if (a >= 0.0) {
670 GSL_SET_COMPLEX(&z, M_PI_2, -acosh(1 / a));
671 } else {
672 GSL_SET_COMPLEX(&z, -M_PI_2, acosh(-1 / a));
673 }
674 }
675
676 return z;
677 }
678
gsl_complex_arccot(gsl_complex a)679 gsl_complex gsl_complex_arccot(gsl_complex a)
680 { /* z = arccot(a) */
681 gsl_complex z;
682
683 if (GSL_REAL(a) == 0.0 && GSL_IMAG(a) == 0.0) {
684 GSL_SET_COMPLEX(&z, M_PI_2, 0);
685 } else {
686 z = gsl_complex_inverse(a);
687 z = gsl_complex_arctan(z);
688 }
689
690 return z;
691 }
692
693 /**********************************************************************
694 * Complex Hyperbolic Functions
695 **********************************************************************/
696
gsl_complex_sinh(gsl_complex a)697 gsl_complex gsl_complex_sinh(gsl_complex a)
698 { /* z = sinh(a) */
699 double R = GSL_REAL(a), I = GSL_IMAG(a);
700
701 gsl_complex z;
702 GSL_SET_COMPLEX(&z, sinh(R) * cos(I), cosh(R) * sin(I));
703 return z;
704 }
705
gsl_complex_cosh(gsl_complex a)706 gsl_complex gsl_complex_cosh(gsl_complex a)
707 { /* z = cosh(a) */
708 double R = GSL_REAL(a), I = GSL_IMAG(a);
709
710 gsl_complex z;
711 GSL_SET_COMPLEX(&z, cosh(R) * cos(I), sinh(R) * sin(I));
712 return z;
713 }
714
gsl_complex_tanh(gsl_complex a)715 gsl_complex gsl_complex_tanh(gsl_complex a)
716 { /* z = tanh(a) */
717 double R = GSL_REAL(a), I = GSL_IMAG(a);
718
719 gsl_complex z;
720
721 if (fabs(R) < 1.0) {
722 double D = pow(cos(I), 2.0) + pow(sinh(R), 2.0);
723
724 GSL_SET_COMPLEX(&z, sinh(R) * cosh(R) / D, 0.5 * sin(2 * I) / D);
725 } else {
726 double D = pow(cos(I), 2.0) + pow(sinh(R), 2.0);
727 double F = 1 + pow(cos(I) / sinh(R), 2.0);
728
729 GSL_SET_COMPLEX(&z, 1.0 / (tanh(R) * F), 0.5 * sin(2 * I) / D);
730 }
731
732 return z;
733 }
734
gsl_complex_sech(gsl_complex a)735 gsl_complex gsl_complex_sech(gsl_complex a)
736 { /* z = sech(a) */
737 gsl_complex z = gsl_complex_cosh(a);
738 return gsl_complex_inverse(z);
739 }
740
gsl_complex_csch(gsl_complex a)741 gsl_complex gsl_complex_csch(gsl_complex a)
742 { /* z = csch(a) */
743 gsl_complex z = gsl_complex_sinh(a);
744 return gsl_complex_inverse(z);
745 }
746
gsl_complex_coth(gsl_complex a)747 gsl_complex gsl_complex_coth(gsl_complex a)
748 { /* z = coth(a) */
749 gsl_complex z = gsl_complex_tanh(a);
750 return gsl_complex_inverse(z);
751 }
752
753 /**********************************************************************
754 * Inverse Complex Hyperbolic Functions
755 **********************************************************************/
756
gsl_complex_arcsinh(gsl_complex a)757 gsl_complex gsl_complex_arcsinh(gsl_complex a)
758 { /* z = arcsinh(a) */
759 gsl_complex z = gsl_complex_mul_imag(a, 1.0);
760 z = gsl_complex_arcsin(z);
761 z = gsl_complex_mul_imag(z, -1.0);
762 return z;
763 }
764
gsl_complex_arccosh(gsl_complex a)765 gsl_complex gsl_complex_arccosh(gsl_complex a)
766 { /* z = arccosh(a) */
767 gsl_complex z = gsl_complex_arccos(a);
768 z = gsl_complex_mul_imag(z, GSL_IMAG(z) > 0 ? -1.0 : 1.0);
769 return z;
770 }
771
gsl_complex_arccosh_real(double a)772 gsl_complex gsl_complex_arccosh_real(double a)
773 { /* z = arccosh(a) */
774 gsl_complex z;
775
776 if (a >= 1) {
777 GSL_SET_COMPLEX(&z, acosh(a), 0);
778 } else {
779 if (a >= -1.0) {
780 GSL_SET_COMPLEX(&z, 0, acos(a));
781 } else {
782 GSL_SET_COMPLEX(&z, acosh(-a), M_PI);
783 }
784 }
785
786 return z;
787 }
788
gsl_complex_arctanh(gsl_complex a)789 gsl_complex gsl_complex_arctanh(gsl_complex a)
790 { /* z = arctanh(a) */
791 if (GSL_IMAG(a) == 0.0) {
792 return gsl_complex_arctanh_real(GSL_REAL(a));
793 } else {
794 gsl_complex z = gsl_complex_mul_imag(a, 1.0);
795 z = gsl_complex_arctan(z);
796 z = gsl_complex_mul_imag(z, -1.0);
797 return z;
798 }
799 }
800
gsl_complex_arctanh_real(double a)801 gsl_complex gsl_complex_arctanh_real(double a)
802 { /* z = arctanh(a) */
803 gsl_complex z;
804
805 if (a > -1.0 && a < 1.0) {
806 GSL_SET_COMPLEX(&z, atanh(a), 0);
807 } else {
808 GSL_SET_COMPLEX(&z, atanh(1 / a), (a < 0) ? M_PI_2 : -M_PI_2);
809 }
810
811 return z;
812 }
813
gsl_complex_arcsech(gsl_complex a)814 gsl_complex gsl_complex_arcsech(gsl_complex a)
815 { /* z = arcsech(a); */
816 gsl_complex t = gsl_complex_inverse(a);
817 return gsl_complex_arccosh(t);
818 }
819
gsl_complex_arccsch(gsl_complex a)820 gsl_complex gsl_complex_arccsch(gsl_complex a)
821 { /* z = arccsch(a) */
822 gsl_complex t = gsl_complex_inverse(a);
823 return gsl_complex_arcsinh(t);
824 }
825
gsl_complex_arccoth(gsl_complex a)826 gsl_complex gsl_complex_arccoth(gsl_complex a)
827 { /* z = arccoth(a) */
828 gsl_complex t = gsl_complex_inverse(a);
829 return gsl_complex_arctanh(t);
830 }
831