cosh.c

OMAP1030 处理器的ARM 侧硬件测试代码 OMAP1030 是TI的双核处理器

/****************************************************************************/
/*  cosh   v2.24                                                           */
/*  Copyright (c) 1995-2002 Texas Instruments Incorporated                  */
/****************************************************************************/
#include <math.h>
#include <values.h>
#include <errno.h>

/****************************************************************************/
/*  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);
}