htswanalysis/MACS/lib/gsl/gsl-1.11/complex/math.c

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