4.3BSD - /usr/src/usr.lib/libm/cosh.c

   1: /*
   2:  * Copyright (c) 1985 Regents of the University of California.
   3:  *
   4:  * Use and reproduction of this software are granted  in  accordance  with
   5:  * the terms and conditions specified in  the  Berkeley  Software  License
   6:  * Agreement (in particular, this entails acknowledgement of the programs'
   7:  * source, and inclusion of this notice) with the additional understanding
   8:  * that  all  recipients  should regard themselves as participants  in  an
   9:  * ongoing  research  project and hence should  feel  obligated  to report
  10:  * their  experiences (good or bad) with these elementary function  codes,
  11:  * using "sendbug 4bsd-bugs@BERKELEY", to the authors.
  12:  */
  13: 
  14: #ifndef lint
  15: static char sccsid[] = "@(#)cosh.c	1.2 (Berkeley) 8/21/85";
  16: #endif not lint
  17: 
  18: /* COSH(X)
  19:  * RETURN THE HYPERBOLIC COSINE OF X
  20:  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
  21:  * CODED IN C BY K.C. NG, 1/8/85;
  22:  * REVISED BY K.C. NG on 2/8/85, 2/23/85, 3/7/85, 3/29/85, 4/16/85.
  23:  *
  24:  * Required system supported functions :
  25:  *	copysign(x,y)
  26:  *	scalb(x,N)
  27:  *
  28:  * Required kernel function:
  29:  *	exp(x)
  30:  *	exp__E(x,c)	...return exp(x+c)-1-x for |x|<0.3465
  31:  *
  32:  * Method :
  33:  *	1. Replace x by |x|.
  34:  *	2.
  35:  *		                                        [ exp(x) - 1 ]^2
  36:  *	    0        <= x <= 0.3465  :  cosh(x) := 1 + -------------------
  37:  *			       			           2*exp(x)
  38:  *
  39:  *		                                   exp(x) +  1/exp(x)
  40:  *	    0.3465   <= x <= 22      :  cosh(x) := -------------------
  41:  *			       			           2
  42:  *	    22       <= x <= lnovfl  :  cosh(x) := exp(x)/2
  43:  *	    lnovfl   <= x <= lnovfl+log(2)
  44:  *				     :  cosh(x) := exp(x)/2 (avoid overflow)
  45:  *	    log(2)+lnovfl <  x <  INF:  overflow to INF
  46:  *
  47:  *	Note: .3465 is a number near one half of ln2.
  48:  *
  49:  * Special cases:
  50:  *	cosh(x) is x if x is +INF, -INF, or NaN.
  51:  *	only cosh(0)=1 is exact for finite x.
  52:  *
  53:  * Accuracy:
  54:  *	cosh(x) returns the exact hyperbolic cosine of x nearly rounded.
  55:  *	In a test run with 768,000 random arguments on a VAX, the maximum
  56:  *	observed error was 1.23 ulps (units in the last place).
  57:  *
  58:  * Constants:
  59:  * The hexadecimal values are the intended ones for the following constants.
  60:  * The decimal values may be used, provided that the compiler will convert
  61:  * from decimal to binary accurately enough to produce the hexadecimal values
  62:  * shown.
  63:  */
  64: 
  65: #ifdef VAX
  66: /* double static  */
  67: /* mln2hi =  8.8029691931113054792E1     , Hex  2^  7   *  .B00F33C7E22BDB */
  68: /* mln2lo = -4.9650192275318476525E-16   , Hex  2^-50   * -.8F1B60279E582A */
  69: /* lnovfl =  8.8029691931113053016E1     ; Hex  2^  7   *  .B00F33C7E22BDA */
  70: static long    mln2hix[] = { 0x0f3343b0, 0x2bdbc7e2};
  71: static long    mln2lox[] = { 0x1b60a70f, 0x582a279e};
  72: static long    lnovflx[] = { 0x0f3343b0, 0x2bdac7e2};
  73: #define   mln2hi    (*(double*)mln2hix)
  74: #define   mln2lo    (*(double*)mln2lox)
  75: #define   lnovfl    (*(double*)lnovflx)
  76: #else   /* IEEE double */
  77: double static
  78: mln2hi =  7.0978271289338397310E2     , /*Hex  2^ 10   *  1.62E42FEFA39EF */
  79: mln2lo =  2.3747039373786107478E-14   , /*Hex  2^-45   *  1.ABC9E3B39803F */
  80: lnovfl =  7.0978271289338397310E2     ; /*Hex  2^  9   *  1.62E42FEFA39EF */
  81: #endif
  82: 
  83: #ifdef VAX
  84: static max = 126                      ;
  85: #else   /* IEEE double */
  86: static max = 1023                     ;
  87: #endif
  88: 
  89: double cosh(x)
  90: double x;
  91: {
  92:     static double half=1.0/2.0,one=1.0, small=1.0E-18; /* fl(1+small)==1 */
  93:     double scalb(),copysign(),exp(),exp__E(),t;
  94: 
  95: #ifndef VAX
  96:     if(x!=x) return(x); /* x is NaN */
  97: #endif
  98:     if((x=copysign(x,one)) <= 22)
  99:         if(x<0.3465)
 100:         if(x<small) return(one+x);
 101:         else {t=x+exp__E(x,0.0);x=t+t; return(one+t*t/(2.0+x)); }
 102: 
 103:         else /* for x lies in [0.3465,22] */
 104:             { t=exp(x); return((t+one/t)*half); }
 105: 
 106:     if( lnovfl <= x && x <= (lnovfl+0.7))
 107:         /* for x lies in [lnovfl, lnovfl+ln2], decrease x by ln(2^(max+1))
 108:          * and return 2^max*exp(x) to avoid unnecessary overflow
 109:          */
 110:         return(scalb(exp((x-mln2hi)-mln2lo), max));
 111: 
 112:     else
 113:         return(exp(x)*half);    /* for large x,  cosh(x)=exp(x)/2 */
 114: }

Defined functions

Defined variables

Defined macros