plugins/fn-complex/gsl-complex.c

   1 /* vim: set sw=8: -*- Mode: C; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -*- */
   2 /* complex/math.c (from the GSL 1.1.1)
   3  *
   4  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Jorma Olavi Tähtinen, Brian 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 2 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, see <https://www.gnu.org/licenses/>.
  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 /*
  45  * Gnumeric specific modifications written by Jukka-Pekka Iivonen
  46  * (jiivonen@hutcs.cs.hut.fi)
  47  *
  48  * long double modifications by Morten Welinder.
  49  */
  50
  51 #include <gnumeric-config.h>
  52 #include <glib/gi18n-lib.h>
  53 #include <gnumeric.h>
  54 #include <func.h>
  55
  56 #include <complex.h>
  57 #include "gsl-complex.h"
  58 #include <parse-util.h>
  59 #include <cell.h>
  60 #include <expr.h>
  61 #include <value.h>
  62 #include <mathfunc.h>
  63 #include <sf-trig.h>
  64
  65
  66 #define GSL_REAL(x) GNM_CRE(*(x))
  67 #define GSL_IMAG(x) GNM_CIM(*(x))
  68
  69 /***********************************************************************
  70  * Complex arithmetic operators
  71  ***********************************************************************/
  72
  73 static inline void
  74 gsl_complex_mul_imag (gnm_complex const *a, gnm_float y, gnm_complex *res)
  75 {                               /* z=a*iy */
  76         *res = GNM_CMAKE (-y * GSL_IMAG (a), y * GSL_REAL (a));
  77 }
  78
  79 void
  80 gsl_complex_inverse (gnm_complex const *a, gnm_complex *res)
  81 {                               /* z=1/a */
  82         gnm_float s = 1.0 / GNM_CABS (*a);
  83
  84         *res = GNM_CMAKE ((GSL_REAL (a) * s) * s, -(GSL_IMAG (a) * s) * s);
  85 }
  86
  87 /**********************************************************************
  88  * Inverse Complex Trigonometric Functions
  89  **********************************************************************/
  90
  91 static void
  92 gsl_complex_arcsin_real (gnm_float a, gnm_complex *res)
  93 {                               /* z = arcsin(a) */
  94          if (gnm_abs (a) <= 1.0) {
  95                  *res = GNM_CMAKE (gnm_asin (a), 0.0);
  96          } else {
  97                  if (a < 0.0) {
  98                          *res = GNM_CMAKE (-M_PI_2gnum, gnm_acosh (-a));
  99                  } else {
 100                          *res = GNM_CMAKE (M_PI_2gnum, -gnm_acosh (a));
 101                  }
 102          }
 103 }
 104
 105 void
 106 gsl_complex_arcsin (gnm_complex const *a, gnm_complex *res)
 107 {                               /* z = arcsin(a) */
 108         gnm_float R = GSL_REAL (a), I = GSL_IMAG (a);
 109
 110         if (I == 0) {
 111                 gsl_complex_arcsin_real (R, res);
 112         } else {
 113                 gnm_float x = gnm_abs (R), y = gnm_abs (I);
 114                 gnm_float r = gnm_hypot (x + 1, y);
 115                 gnm_float s = gnm_hypot (x - 1, y);
 116                 gnm_float A = 0.5 * (r + s);
 117                 gnm_float B = x / A;
 118                 gnm_float y2 = y * y;
 119
 120                 gnm_float real, imag;
 121
 122                 const gnm_float A_crossover = 1.5, B_crossover = 0.6417;
 123
 124                 if (B <= B_crossover) {
 125                         real = gnm_asin (B);
 126                 } else {
 127                         if (x <= 1) {
 128                                 gnm_float D = 0.5 * (A + x) *
 129                                         (y2 / (r + x + 1) + (s + (1 - x)));
 130                                 real = gnm_atan (x / gnm_sqrt (D));
 131                         } else {
 132                                 gnm_float Apx = A + x;
 133                                 gnm_float D = 0.5 * (Apx / (r + x + 1)
 134                                                       + Apx / (s + (x - 1)));
 135                                 real = gnm_atan (x / (y * gnm_sqrt (D)));
 136                         }
 137                 }
 138
 139                 if (A <= A_crossover) {
 140                         gnm_float Am1;
 141
 142                         if (x < 1) {
 143                                 Am1 = 0.5 * (y2 / (r + (x + 1)) + y2 /
 144                                              (s + (1 - x)));
 145                         } else {
 146                                 Am1 = 0.5 * (y2 / (r + (x + 1)) +
 147                                              (s + (x - 1)));
 148                         }
 149
 150                         imag = gnm_log1p (Am1 + gnm_sqrt (Am1 * (A + 1)));
 151                 } else {
 152                         imag = gnm_log (A + gnm_sqrt (A * A - 1));
 153                 }
 154
 155                 *res = GNM_CMAKE ((R >= 0) ? real : -real, (I >= 0) ?
 156                               imag : -imag);
 157         }
 158 }
 159
 160 static void
 161 gsl_complex_arccos_real (gnm_float a, gnm_complex *res)
 162 {                               /* z = arccos(a) */
 163         if (gnm_abs (a) <= 1.0) {
 164                 *res = GNM_CMAKE (gnm_acos (a), 0);
 165         } else {
 166                 if (a < 0.0) {
 167                         *res = GNM_CMAKE (M_PIgnum, -gnm_acosh (-a));
 168                 } else {
 169                         *res = GNM_CMAKE (0, gnm_acosh (a));
 170                 }
 171         }
 172 }
 173
 174 void
 175 gsl_complex_arccos (gnm_complex const *a, gnm_complex *res)
 176 {                               /* z = arccos(a) */
 177         gnm_float R = GSL_REAL (a), I = GSL_IMAG (a);
 178
 179         if (I == 0) {
 180                 gsl_complex_arccos_real (R, res);
 181         } else {
 182                 gnm_float x = gnm_abs (R);
 183                 gnm_float y = gnm_abs (I);
 184                 gnm_float r = gnm_hypot (x + 1, y);
 185                 gnm_float s = gnm_hypot (x - 1, y);
 186                 gnm_float A = 0.5 * (r + s);
 187                 gnm_float B = x / A;
 188                 gnm_float y2 = y * y;
 189
 190                 gnm_float real, imag;
 191
 192                 const gnm_float A_crossover = 1.5;
 193                 const gnm_float B_crossover = 0.6417;
 194
 195                 if (B <= B_crossover) {
 196                         real = gnm_acos (B);
 197                 } else {
 198                         if (x <= 1) {
 199                                 gnm_float D = 0.5 * (A + x) *
 200                                         (y2 / (r + x + 1) + (s + (1 - x)));
 201                                 real = gnm_atan (gnm_sqrt (D) / x);
 202                         } else {
 203                                 gnm_float Apx = A + x;
 204                                 gnm_float D = 0.5 * (Apx / (r + x + 1) + Apx /
 205                                                       (s + (x - 1)));
 206                                 real = gnm_atan ((y * gnm_sqrt (D)) / x);
 207                         }
 208                 }
 209                 if (A <= A_crossover) {
 210                         gnm_float Am1;
 211
 212                         if (x < 1) {
 213                                 Am1 = 0.5 * (y2 / (r + (x + 1)) + y2 /
 214                                              (s + (1 - x)));
 215                         } else {
 216                                 Am1 = 0.5 * (y2 / (r + (x + 1)) +
 217                                              (s + (x - 1)));
 218                         }
 219
 220                         imag = gnm_log1p (Am1 + gnm_sqrt (Am1 * (A + 1)));
 221                 } else {
 222                         imag = gnm_log (A + gnm_sqrt (A * A - 1));
 223                 }
 224
 225                 *res = GNM_CMAKE ((R >= 0) ? real : M_PIgnum - real, (I >= 0) ?
 226                               -imag : imag);
 227         }
 228 }
 229
 230 void
 231 gsl_complex_arctan (gnm_complex const *a, gnm_complex *res)
 232 {                               /* z = arctan(a) */
 233         gnm_float R = GSL_REAL (a), I = GSL_IMAG (a);
 234
 235         if (I == 0) {
 236                 *res = GNM_CMAKE (gnm_atan (R), 0);
 237         } else {
 238                 /* FIXME: This is a naive implementation which does not fully
 239                  * take into account cancellation errors, overflow, underflow
 240                  * etc.  It would benefit from the Hull et al treatment. */
 241
 242                 gnm_float r = gnm_hypot (R, I);
 243
 244                 gnm_float imag;
 245
 246                 gnm_float u = 2 * I / (1 + r * r);
 247
 248                 /* FIXME: the following cross-over should be optimized but 0.1
 249                  * seems to work ok */
 250
 251                 if (gnm_abs (u) < 0.1) {
 252                         imag = 0.25 * (gnm_log1p (u) - gnm_log1p (-u));
 253                 } else {
 254                         gnm_float A = gnm_hypot (R, I + 1);
 255                         gnm_float B = gnm_hypot (R, I - 1);
 256                         imag = 0.5 * gnm_log (A / B);
 257                 }
 258                 if (R == 0) {
 259                         if (I > 1) {
 260                                 *res = GNM_CMAKE (M_PI_2gnum, imag);
 261                         } else if (I < -1) {
 262                                 *res = GNM_CMAKE (-M_PI_2gnum, imag);
 263                         } else {
 264                                 *res = GNM_CMAKE (0, imag);
 265                         }
 266                 } else {
 267                         *res = GNM_CMAKE (0.5 * gnm_atan2 (2 * R,
 268                                                             ((1 + r) * (1 - r))),
 269                                       imag);
 270                 }
 271         }
 272 }
 273
 274 void
 275 gsl_complex_arcsec (gnm_complex const *a, gnm_complex *res)
 276 {                               /* z = arcsec(a) */
 277         gsl_complex_inverse (a, res);
 278         gsl_complex_arccos (res, res);
 279 }
 280
 281 void
 282 gsl_complex_arccsc (gnm_complex const *a, gnm_complex *res)
 283 {                               /* z = arccsc(a) */
 284         gsl_complex_inverse (a, res);
 285         gsl_complex_arcsin (res, res);
 286 }
 287
 288 void
 289 gsl_complex_arccot (gnm_complex const *a, gnm_complex *res)
 290 {                               /* z = arccot(a) */
 291         if (GSL_REAL (a) == 0.0 && GSL_IMAG (a) == 0.0) {
 292                 *res = GNM_CMAKE (M_PI_2gnum, 0);
 293         } else {
 294                 gsl_complex_inverse (a, res);
 295                 gsl_complex_arctan (res, res);
 296         }
 297 }
 298
 299 /**********************************************************************
 300  * Complex Hyperbolic Functions
 301  **********************************************************************/
 302
 303 void
 304 gsl_complex_sinh (gnm_complex const *a, gnm_complex *res)
 305 {                               /* z = sinh(a) */
 306         gnm_float R = GSL_REAL (a), I = GSL_IMAG (a);
 307
 308         *res = GNM_CMAKE (gnm_sinh (R) * gnm_cos (I), gnm_cosh (R) * gnm_sin (I));
 309 }
 310
 311 void
 312 gsl_complex_cosh (gnm_complex const *a, gnm_complex *res)
 313 {                               /* z = cosh(a) */
 314         gnm_float R = GSL_REAL (a), I = GSL_IMAG (a);
 315
 316         *res = GNM_CMAKE (gnm_cosh (R) * gnm_cos (I), gnm_sinh (R) * gnm_sin (I));
 317 }
 318
 319 void
 320 gsl_complex_tanh (gnm_complex const *a, gnm_complex *res)
 321 {                               /* z = tanh(a) */
 322         gnm_float R = GSL_REAL (a), I = GSL_IMAG (a);
 323
 324         if (gnm_abs (R) < 1.0) {
 325                  gnm_float D =
 326                          gnm_pow (gnm_cos (I), 2.0) +
 327                          gnm_pow (gnm_sinh (R), 2.0);
 328
 329                  *res = GNM_CMAKE (gnm_sinh (R) * gnm_cosh (R) / D,
 330                                0.5 * gnm_sin (2 * I) / D);
 331         } else {
 332                  gnm_float D =
 333                          gnm_pow (gnm_cos (I), 2.0) +
 334                          gnm_pow (gnm_sinh (R), 2.0);
 335                  gnm_float F = 1 + gnm_pow (gnm_cos (I) / gnm_sinh (R), 2.0);
 336
 337                  *res = GNM_CMAKE (1.0 / (gnm_tanh (R) * F),
 338                                0.5 * gnm_sin (2 * I) / D);
 339         }
 340 }
 341
 342 void
 343 gsl_complex_sech (gnm_complex const *a, gnm_complex *res)
 344 {                               /* z = sech(a) */
 345         gsl_complex_cosh (a, res);
 346         gsl_complex_inverse (res, res);
 347 }
 348
 349 void
 350 gsl_complex_csch (gnm_complex const *a, gnm_complex *res)
 351 {                               /* z = csch(a) */
 352         gsl_complex_sinh (a, res);
 353         gsl_complex_inverse (res, res);
 354 }
 355
 356 void
 357 gsl_complex_coth (gnm_complex const *a, gnm_complex *res)
 358 {                               /* z = coth(a) */
 359         gsl_complex_tanh (a, res);
 360         gsl_complex_inverse (res, res);
 361 }
 362
 363 /**********************************************************************
 364  * Inverse Complex Hyperbolic Functions
 365  **********************************************************************/
 366
 367 void
 368 gsl_complex_arcsinh (gnm_complex const *a, gnm_complex *res)
 369 {                               /* z = arcsinh(a) */
 370         gsl_complex_mul_imag (a, 1.0, res);
 371         gsl_complex_arcsin (res, res);
 372         gsl_complex_mul_imag (res, -1.0, res);
 373 }
 374
 375 void
 376 gsl_complex_arccosh (gnm_complex const *a, gnm_complex *res)
 377 {                               /* z = arccosh(a) */
 378         if (GSL_IMAG (a) == 0 && GSL_REAL (a) == 1)
 379                 *res = GNM_C0;
 380         else {
 381                 gsl_complex_arccos (a, res);
 382                 gsl_complex_mul_imag (res, GSL_IMAG (res) > 0 ? -1.0 : 1.0, res);
 383         }
 384 }
 385
 386 static void
 387 gsl_complex_arctanh_real (gnm_float a, gnm_complex *res)
 388 {                               /* z = arctanh(a) */
 389         if (a > -1.0 && a < 1.0) {
 390                 *res = GNM_CMAKE (gnm_atanh (a), 0);
 391         } else {
 392                 *res = GNM_CMAKE (gnm_acoth (a),
 393                               (a < 0) ? M_PI_2gnum : -M_PI_2gnum);
 394         }
 395 }
 396
 397 void
 398 gsl_complex_arctanh (gnm_complex const *a, gnm_complex *res)
 399 {                               /* z = arctanh(a) */
 400         if (GSL_IMAG (a) == 0.0) {
 401                 gsl_complex_arctanh_real (GSL_REAL (a), res);
 402         } else {
 403                  gsl_complex_mul_imag (a, 1.0, res);
 404                  gsl_complex_arctan (res, res);
 405                  gsl_complex_mul_imag (res, -1.0, res);
 406         }
 407 }
 408
 409 void
 410 gsl_complex_arcsech (gnm_complex const *a, gnm_complex *res)
 411 {                               /* z = arcsech(a); */
 412         gsl_complex_inverse (a, res);
 413         gsl_complex_arccosh (res, res);
 414 }
 415
 416 void
 417 gsl_complex_arccsch (gnm_complex const *a, gnm_complex *res)
 418 {                               /* z = arccsch(a); */
 419         gsl_complex_inverse (a, res);
 420         gsl_complex_arcsinh (res, res);
 421 }
 422
 423 void
 424 gsl_complex_arccoth (gnm_complex const *a, gnm_complex *res)
 425 {                               /* z = arccoth(a); */
 426         gsl_complex_inverse (a, res);
 427         gsl_complex_arctanh (res, res);
 428 }