htswanalysis/MACS/lib/gsl/gsl-1.11/sys/invhyp.c

   1 /* sys/invhyp.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 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 #include <config.h>
  21 #include <math.h>
  22 #include <gsl/gsl_math.h>
  23 #include <gsl/gsl_machine.h>
  24
  25 double
  26 gsl_acosh (const double x)
  27 {
  28   if (x > 1.0 / GSL_SQRT_DBL_EPSILON)
  29     {
  30       return log (x) + M_LN2;
  31     }
  32   else if (x > 2)
  33     {
  34       return log (2 * x - 1 / (sqrt (x * x - 1) + x));
  35     }
  36   else if (x > 1)
  37     {
  38       double t = x - 1;
  39       return log1p (t + sqrt (2 * t + t * t));
  40     }
  41   else if (x == 1)
  42     {
  43       return 0;
  44     }
  45   else
  46     {
  47       return GSL_NAN;
  48     }
  49 }
  50
  51 double
  52 gsl_asinh (const double x)
  53 {
  54   double a = fabs (x);
  55   double s = (x < 0) ? -1 : 1;
  56
  57   if (a > 1 / GSL_SQRT_DBL_EPSILON)
  58     {
  59       return s * (log (a) + M_LN2);
  60     }
  61   else if (a > 2)
  62     {
  63       return s * log (2 * a + 1 / (a + sqrt (a * a + 1)));
  64     }
  65   else if (a > GSL_SQRT_DBL_EPSILON)
  66     {
  67       double a2 = a * a;
  68       return s * log1p (a + a2 / (1 + sqrt (1 + a2)));
  69     }
  70   else
  71     {
  72       return x;
  73     }
  74 }
  75
  76 double
  77 gsl_atanh (const double x)
  78 {
  79   double a = fabs (x);
  80   double s = (x < 0) ? -1 : 1;
  81
  82   if (a > 1)
  83     {
  84       return GSL_NAN;
  85     }
  86   else if (a == 1)
  87     {
  88       return (x < 0) ? GSL_NEGINF : GSL_POSINF;
  89     }
  90   else if (a >= 0.5)
  91     {
  92       return s * 0.5 * log1p (2 * a / (1 - a));
  93     }
  94   else if (a > GSL_DBL_EPSILON)
  95     {
  96       return s * 0.5 * log1p (2 * a + 2 * a * a / (1 - a));
  97     }
  98   else
  99     {
 100       return x;
 101     }
 102 }