sysdeps/ieee754/flt-32/e_atanhf.c

   1 /* Copyright (C) 2011 Free Software Foundation, Inc.
   2    This file is part of the GNU C Library.
   3    Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
   4
   5    The GNU C Library is free software; you can redistribute it and/or
   6    modify it under the terms of the GNU Lesser General Public
   7    License as published by the Free Software Foundation; either
   8    version 2.1 of the License, or (at your option) any later version.
   9
  10    The GNU C Library is distributed in the hope that it will be useful,
  11    but WITHOUT ANY WARRANTY; without even the implied warranty of
  12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13    Lesser General Public License for more details.
  14
  15    You should have received a copy of the GNU Lesser General Public
  16    License along with the GNU C Library; if not, write to the Free
  17    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
  18    02111-1307 USA.  */
  19
  20
  21 /* __ieee754_atanh(x)
  22    Method :
  23       1.Reduced x to positive by atanh(-x) = -atanh(x)
  24       2.For x>=0.5
  25                     1              2x                          x
  26         atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
  27                     2             1 - x                      1 - x
  28
  29         For x<0.5
  30         atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
  31
  32    Special cases:
  33         atanh(x) is NaN if |x| > 1 with signal;
  34         atanh(NaN) is that NaN with no signal;
  35         atanh(+-1) is +-INF with signal.
  36
  37  */
  38
  39 #include <inttypes.h>
  40 #include "math.h"
  41 #include "math_private.h"
  42
  43 static const float huge = 1e30;
  44
  45 float
  46 __ieee754_atanhf (float x)
  47 {
  48   float xa = fabsf (x);
  49   float t;
  50   if (xa < 0.5f)
  51     {
  52       if (__builtin_expect (xa < 0x1.0p-28f, 0))
  53         {
  54           math_force_eval (huge + x);
  55           return x;
  56         }
  57
  58       t = xa + xa;
  59       t = 0.5f * __log1pf (t + t * xa / (1.0f - xa));
  60     }
  61   else if (__builtin_expect (xa < 1.0f, 1))
  62     t = 0.5f * __log1pf ((xa + xa) / (1.0f - xa));
  63   else
  64     {
  65       if (xa > 1.0f)
  66         return (x - x) / (x - x);
  67
  68       return x / 0.0f;
  69     }
  70
  71   return __copysignf (t, x);
  72 }
  73 strong_alias (__ieee754_atanhf, __atanhf_finite)