www.pudn.com > system.rar > cosh.c

/****************************************************************************/ 
/*  cosh   v2.54                                                            */ 
/*  Copyright (c) 1995-2004 Texas Instruments Incorporated                  */ 
/****************************************************************************/ 
#include  
#include  
#include  
 
/****************************************************************************/ 
/*  COSH() - Hyperbolic Cosine						    */ 
/*               							    */ 
/*  Based on the algorithm from "Software Manual for the Elementary         */ 
/*  Functions", Cody and Waite, Prentice Hall 1980, chapter 6.              */ 
/*               							    */ 
/*  result = (exp(x) + 1 / exp(x)) / 2		  			    */ 
/****************************************************************************/ 
double cosh(double x) 
{ 
  double g, z, q, p, r, a, b; 
  int n; 
 
  /****************************************************************************/ 
  /* cosh(x) = cosh(-x)                                                       */ 
  /****************************************************************************/ 
  x = fabs(x); 
 
  /****************************************************************************/ 
  /* check to see if overflow would occur                                     */ 
  /****************************************************************************/ 
  if (x > MAXH) { errno = ERANGE; return (HUGE_VAL); } 
 
  if (x < 0) n = (int) (x * INVLOGe2 - 0.5);       /* since (int) -1.5 = -1.0 */ 
  else       n = (int) (x * INVLOGe2 + 0.5); 
 
  /****************************************************************************/ 
  /* g = x - n * ln(2) (but more mathematically stable)                       */ 
  /****************************************************************************/ 
  g  = (x - n * C3) - n * C4; 
 
  /****************************************************************************/ 
  /* determine polynomial expression                                          */ 
  /****************************************************************************/ 
  z  = g * g; 
 
#if BITS <=29 
  p = (EXP1 * z + EXP0) * g; 
  q = EXQ1 * z + EXQ0; 
#elif BITS>=30 && BITS<=42 
  p = (EXP1 * z + EXP0) * g; 
  q = (EXQ2 * z + EXQ1) * z + EXQ0; 
#elif BITS>=43 && BITS<=56 
  p = ((EXP2 * z + EXP1) * z + EXP0) * g; 
  q = (EXQ2 * z + EXQ1) * z + EXQ0; 
#else 
  p = ((EXP2 * z + EXP1) * z + EXP0) * g; 
  q = ((EXQ3 * z + EXQ2) * z + EXQ1) * z + EXQ0; 
#endif 
 
  /****************************************************************************/ 
  /* calculate exp(g)                                                         */ 
  /****************************************************************************/ 
  r  = 0.5 + p / (q - p); 
 
  /****************************************************************************/ 
  /* exp(x)/2 = exp(g) * 2 ^ (n)                                              */ 
  /****************************************************************************/ 
  a = ldexp(r, n); 
 
  /****************************************************************************/ 
  /* exp(-x)*2 = 1 / (exp(x)/2)                                               */ 
  /****************************************************************************/ 
  b  = 1.0 / a; 
 
  /****************************************************************************/ 
  /* cosh(x) = exp(x)/2 + exp(-x)*2 / 4                                       */ 
  /****************************************************************************/ 
  return (a + b * 0.25); 
}