lib/libm/sinh.c

   1 /*                                                      sinh.c
   2  *
   3  *      Hyperbolic sine
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * double x, y, sinh();
  10  *
  11  * y = sinh( x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Returns hyperbolic sine of argument in the range MINLOG to
  18  * MAXLOG.
  19  *
  20  * The range is partitioned into two segments.  If |x| <= 1, a
  21  * rational function of the form x + x**3 P(x)/Q(x) is employed.
  22  * Otherwise the calculation is sinh(x) = ( exp(x) - exp(-x) )/2.
  23  *
  24  *
  25  *
  26  * ACCURACY:
  27  *
  28  *                      Relative error:
  29  * arithmetic   domain     # trials      peak         rms
  30  *    DEC      +- 88        50000       4.0e-17     7.7e-18
  31  *    IEEE     +-MAXLOG     30000       2.6e-16     5.7e-17
  32  *
  33  */
  34 \f
  35 /*
  36 Cephes Math Library Release 2.8:  June, 2000
  37 Copyright 1984, 1995, 2000 by Stephen L. Moshier
  38 */
  39
  40 #include "mconf.h"
  41
  42 #ifdef UNK
  43 static double P[] = {
  44 -7.89474443963537015605E-1,
  45 -1.63725857525983828727E2,
  46 -1.15614435765005216044E4,
  47 -3.51754964808151394800E5
  48 };
  49 static double Q[] = {
  50 /* 1.00000000000000000000E0,*/
  51 -2.77711081420602794433E2,
  52  3.61578279834431989373E4,
  53 -2.11052978884890840399E6
  54 };
  55 #endif
  56
  57 #ifdef DEC
  58 static unsigned short P[] = {
  59 0140112,0015377,0042731,0163255,
  60 0142043,0134721,0146177,0123761,
  61 0143464,0122706,0034353,0006017,
  62 0144653,0140536,0157665,0054045
  63 };
  64 static unsigned short Q[] = {
  65 /*0040200,0000000,0000000,0000000,*/
  66 0142212,0155404,0133513,0022040,
  67 0044015,0036723,0173271,0011053,
  68 0145400,0150407,0023710,0001034
  69 };
  70 #endif
  71
  72 #ifdef IBMPC
  73 static unsigned short P[] = {
  74 0x3cd6,0xe8bb,0x435f,0xbfe9,
  75 0xf4fe,0x398f,0x773a,0xc064,
  76 0x6182,0xc71d,0x94b8,0xc0c6,
  77 0xab05,0xdbf6,0x782b,0xc115
  78 };
  79 static unsigned short Q[] = {
  80 /*0x0000,0x0000,0x0000,0x3ff0,*/
  81 0x6484,0x96e9,0x5b60,0xc071,
  82 0x2245,0x7ed7,0xa7ba,0x40e1,
  83 0x0044,0xe4f9,0x1a20,0xc140
  84 };
  85 #endif
  86
  87 #ifdef MIEEE
  88 static unsigned short P[] = {
  89 0xbfe9,0x435f,0xe8bb,0x3cd6,
  90 0xc064,0x773a,0x398f,0xf4fe,
  91 0xc0c6,0x94b8,0xc71d,0x6182,
  92 0xc115,0x782b,0xdbf6,0xab05
  93 };
  94 static unsigned short Q[] = {
  95 0xc071,0x5b60,0x96e9,0x6484,
  96 0x40e1,0xa7ba,0x7ed7,0x2245,
  97 0xc140,0x1a20,0xe4f9,0x0044
  98 };
  99 #endif
 100
 101 #ifdef ANSIPROT
 102 extern double fabs ( double );
 103 extern double exp ( double );
 104 extern double polevl ( double, void *, int );
 105 extern double p1evl ( double, void *, int );
 106 #else
 107 double fabs(), exp(), polevl(), p1evl();
 108 #endif
 109 extern double INFINITY, MINLOG, MAXLOG, LOGE2;
 110
 111 double sinh(x)
 112 double x;
 113 {
 114 double a;
 115
 116 #ifdef MINUSZERO
 117 if( x == 0.0 )
 118         return(x);
 119 #endif
 120 a = fabs(x);
 121 if( (x > (MAXLOG + LOGE2)) || (x > -(MINLOG-LOGE2) ) )
 122         {
 123         mtherr( "sinh", DOMAIN );
 124         if( x > 0 )
 125                 return( INFINITY );
 126         else
 127                 return( -INFINITY );
 128         }
 129 if( a > 1.0 )
 130         {
 131         if( a >= (MAXLOG - LOGE2) )
 132                 {
 133                 a = exp(0.5*a);
 134                 a = (0.5 * a) * a;
 135                 if( x < 0 )
 136                         a = -a;
 137                 return(a);
 138                 }
 139         a = exp(a);
 140         a = 0.5*a - (0.5/a);
 141         if( x < 0 )
 142                 a = -a;
 143         return(a);
 144         }
 145
 146 a *= a;
 147 return( x + x * a * (polevl(a,P,3)/p1evl(a,Q,3)) );
 148 }