Home | History | Annotate | Line # | Download | only in noieee_src
n_cosh.c revision 1.7.22.1
      1  1.7.22.1    matt /*      n_cosh.c,v 1.7 2003/08/07 16:44:50 agc Exp */
      2       1.1   ragge /*
      3       1.1   ragge  * Copyright (c) 1985, 1993
      4       1.1   ragge  *	The Regents of the University of California.  All rights reserved.
      5       1.1   ragge  *
      6       1.1   ragge  * Redistribution and use in source and binary forms, with or without
      7       1.1   ragge  * modification, are permitted provided that the following conditions
      8       1.1   ragge  * are met:
      9       1.1   ragge  * 1. Redistributions of source code must retain the above copyright
     10       1.1   ragge  *    notice, this list of conditions and the following disclaimer.
     11       1.1   ragge  * 2. Redistributions in binary form must reproduce the above copyright
     12       1.1   ragge  *    notice, this list of conditions and the following disclaimer in the
     13       1.1   ragge  *    documentation and/or other materials provided with the distribution.
     14       1.7     agc  * 3. Neither the name of the University nor the names of its contributors
     15       1.1   ragge  *    may be used to endorse or promote products derived from this software
     16       1.1   ragge  *    without specific prior written permission.
     17       1.1   ragge  *
     18       1.1   ragge  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
     19       1.1   ragge  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
     20       1.1   ragge  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
     21       1.1   ragge  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
     22       1.1   ragge  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
     23       1.1   ragge  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
     24       1.1   ragge  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
     25       1.1   ragge  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
     26       1.1   ragge  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
     27       1.1   ragge  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
     28       1.1   ragge  * SUCH DAMAGE.
     29       1.1   ragge  */
     30       1.1   ragge 
     31       1.1   ragge #ifndef lint
     32       1.2   ragge #if 0
     33       1.1   ragge static char sccsid[] = "@(#)cosh.c	8.1 (Berkeley) 6/4/93";
     34       1.2   ragge #endif
     35       1.1   ragge #endif /* not lint */
     36       1.1   ragge 
     37       1.1   ragge /* COSH(X)
     38       1.1   ragge  * RETURN THE HYPERBOLIC COSINE OF X
     39       1.1   ragge  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
     40       1.5  simonb  * CODED IN C BY K.C. NG, 1/8/85;
     41       1.1   ragge  * REVISED BY K.C. NG on 2/8/85, 2/23/85, 3/7/85, 3/29/85, 4/16/85.
     42       1.1   ragge  *
     43       1.1   ragge  * Required system supported functions :
     44       1.1   ragge  *	copysign(x,y)
     45       1.1   ragge  *	scalb(x,N)
     46       1.1   ragge  *
     47       1.1   ragge  * Required kernel function:
     48       1.5  simonb  *	exp(x)
     49       1.1   ragge  *	exp__E(x,c)	...return exp(x+c)-1-x for |x|<0.3465
     50       1.1   ragge  *
     51       1.1   ragge  * Method :
     52       1.5  simonb  *	1. Replace x by |x|.
     53       1.5  simonb  *	2.
     54       1.5  simonb  *		                                        [ exp(x) - 1 ]^2
     55       1.1   ragge  *	    0        <= x <= 0.3465  :  cosh(x) := 1 + -------------------
     56       1.1   ragge  *			       			           2*exp(x)
     57       1.1   ragge  *
     58       1.1   ragge  *		                                   exp(x) +  1/exp(x)
     59       1.1   ragge  *	    0.3465   <= x <= 22      :  cosh(x) := -------------------
     60       1.1   ragge  *			       			           2
     61       1.5  simonb  *	    22       <= x <= lnovfl  :  cosh(x) := exp(x)/2
     62       1.1   ragge  *	    lnovfl   <= x <= lnovfl+log(2)
     63       1.1   ragge  *				     :  cosh(x) := exp(x)/2 (avoid overflow)
     64       1.1   ragge  *	    log(2)+lnovfl <  x <  INF:  overflow to INF
     65       1.1   ragge  *
     66       1.1   ragge  *	Note: .3465 is a number near one half of ln2.
     67       1.1   ragge  *
     68       1.1   ragge  * Special cases:
     69       1.1   ragge  *	cosh(x) is x if x is +INF, -INF, or NaN.
     70       1.1   ragge  *	only cosh(0)=1 is exact for finite x.
     71       1.1   ragge  *
     72       1.1   ragge  * Accuracy:
     73       1.1   ragge  *	cosh(x) returns the exact hyperbolic cosine of x nearly rounded.
     74       1.1   ragge  *	In a test run with 768,000 random arguments on a VAX, the maximum
     75       1.1   ragge  *	observed error was 1.23 ulps (units in the last place).
     76       1.1   ragge  *
     77       1.1   ragge  * Constants:
     78       1.1   ragge  * The hexadecimal values are the intended ones for the following constants.
     79       1.1   ragge  * The decimal values may be used, provided that the compiler will convert
     80       1.1   ragge  * from decimal to binary accurately enough to produce the hexadecimal values
     81       1.1   ragge  * shown.
     82       1.1   ragge  */
     83       1.1   ragge 
     84       1.6    matt #define _LIBM_STATIC
     85  1.7.22.1    matt #include "../src/namespace.h"
     86       1.1   ragge #include "mathimpl.h"
     87       1.1   ragge 
     88  1.7.22.1    matt #ifdef __weak_alias
     89  1.7.22.1    matt __weak_alias(cosh, _cosh);
     90  1.7.22.1    matt __weak_alias(coshf, _coshf);
     91  1.7.22.1    matt #endif
     92  1.7.22.1    matt 
     93       1.1   ragge vc(mln2hi, 8.8029691931113054792E1   ,0f33,43b0,2bdb,c7e2,   7, .B00F33C7E22BDB)
     94       1.1   ragge vc(mln2lo,-4.9650192275318476525E-16 ,1b60,a70f,582a,279e, -50,-.8F1B60279E582A)
     95       1.1   ragge vc(lnovfl, 8.8029691931113053016E1   ,0f33,43b0,2bda,c7e2,   7, .B00F33C7E22BDA)
     96       1.1   ragge 
     97       1.1   ragge ic(mln2hi, 7.0978271289338397310E2,    10, 1.62E42FEFA39EF)
     98       1.1   ragge ic(mln2lo, 2.3747039373786107478E-14, -45, 1.ABC9E3B39803F)
     99       1.1   ragge ic(lnovfl, 7.0978271289338397310E2,     9, 1.62E42FEFA39EF)
    100       1.1   ragge 
    101       1.1   ragge #ifdef vccast
    102       1.1   ragge #define   mln2hi    vccast(mln2hi)
    103       1.1   ragge #define   mln2lo    vccast(mln2lo)
    104       1.1   ragge #define   lnovfl    vccast(lnovfl)
    105       1.1   ragge #endif
    106       1.1   ragge 
    107       1.3    matt #if defined(__vax__)||defined(tahoe)
    108       1.6    matt #define EXPMAX 126
    109       1.6    matt #else
    110       1.6    matt #define EXPMAX 1023
    111       1.3    matt #endif	/* defined(__vax__)||defined(tahoe) */
    112       1.1   ragge 
    113       1.6    matt double
    114       1.6    matt cosh(double x)
    115       1.5  simonb {
    116       1.1   ragge 	static const double half=1.0/2.0,
    117       1.1   ragge 		one=1.0, small=1.0E-18; /* fl(1+small)==1 */
    118       1.1   ragge 	double t;
    119       1.1   ragge 
    120       1.4   ragge #if !defined(__vax__)&&!defined(tahoe)
    121       1.1   ragge 	if(x!=x) return(x);	/* x is NaN */
    122       1.4   ragge #endif	/* !defined(__vax__)&&!defined(tahoe) */
    123       1.3    matt 	if((x=copysign(x,one)) <= 22) {
    124       1.3    matt 	    if(x<0.3465) {
    125       1.3    matt 		if(x<small) { return(one+x); }
    126       1.1   ragge 		else {t=x+__exp__E(x,0.0);x=t+t; return(one+t*t/(2.0+x)); }
    127       1.1   ragge 
    128       1.3    matt 	    } else /* for x lies in [0.3465,22] */
    129       1.1   ragge 	        { t=exp(x); return((t+one/t)*half); }
    130       1.3    matt 	}
    131       1.1   ragge 
    132       1.5  simonb 	if( lnovfl <= x && x <= (lnovfl+0.7))
    133       1.6    matt         /* for x lies in [lnovfl, lnovfl+ln2], decrease x by ln(2^(EXPMAX+1))
    134       1.6    matt          * and return 2^EXPMAX*exp(x) to avoid unnecessary overflow
    135       1.1   ragge          */
    136       1.6    matt 	    return(scalb(exp((x-mln2hi)-mln2lo), EXPMAX));
    137       1.1   ragge 
    138       1.5  simonb 	else
    139       1.1   ragge 	    return(exp(x)*half);	/* for large x,  cosh(x)=exp(x)/2 */
    140       1.1   ragge }
    141  1.7.22.1    matt 
    142  1.7.22.1    matt float
    143  1.7.22.1    matt coshf(float x)
    144  1.7.22.1    matt {
    145  1.7.22.1    matt 	return(cosh((double)x));
    146  1.7.22.1    matt }
    147