Blame view

Modif/epsilon-master/liba/src/external/openbsd/s_asinh.c 1.48 KB
6663b6c9   adorian   projet complet av...
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
  /* @(#)s_asinh.c 5.1 93/09/24 */
  /*
   * ====================================================
   * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   *
   * Developed at SunPro, a Sun Microsystems, Inc. business.
   * Permission to use, copy, modify, and distribute this
   * software is freely granted, provided that this notice 
   * is preserved.
   * ====================================================
   */
  
  /* asinh(x)
   * Method :
   *	Based on 
   *		asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ]
   *	we have
   *	asinh(x) := x  if  1+x*x=1,
   *		 := sign(x)*(log(x)+ln2)) for large |x|, else
   *		 := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else
   *		 := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))  
   */
  
  #include "math.h"
  #include "math_private.h"
  
  static const double 
  one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
  ln2 =  6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
  huge=  1.00000000000000000000e+300; 
  
  double
  asinh(double x)
  {	
  	double t,w;
  	int32_t hx,ix;
  	GET_HIGH_WORD(hx,x);
  	ix = hx&0x7fffffff;
  	if(ix>=0x7ff00000) return x+x;	/* x is inf or NaN */
  	if(ix< 0x3e300000) {	/* |x|<2**-28 */
  	    if(huge+x>one) return x;	/* return x inexact except 0 */
  	} 
  	if(ix>0x41b00000) {	/* |x| > 2**28 */
  	    w = log(fabs(x))+ln2;
  	} else if (ix>0x40000000) {	/* 2**28 > |x| > 2.0 */
  	    t = fabs(x);
  	    w = log(2.0*t+one/(sqrt(x*x+one)+t));
  	} else {		/* 2.0 > |x| > 2**-28 */
  	    t = x*x;
  	    w =log1p(fabs(x)+t/(one+sqrt(one+t)));
  	}
  	if(hx>0) return w; else return -w;
  }