Gecko [ 167.86.127.34 ]

Name	Size	Permission
big	[ DIR ]	drwxr-xr-x
bits	[ DIR ]	drwxr-xr-x
cmplx	[ DIR ]	drwxr-xr-x
rand	[ DIR ]	drwxr-xr-x
abs.go	363 B	-rw-r--r--
acos_s390x.s	3.73 KB	-rw-r--r--
acosh.go	1.72 KB	-rw-r--r--
acosh_s390x.s	4.32 KB	-rw-r--r--
all_test.go	84.21 KB	-rw-r--r--
arith_s390x.go	2.01 KB	-rw-r--r--
arith_s390x_test.go	10.78 KB	-rw-r--r--
asin.go	983 B	-rw-r--r--
asin_386.s	1.09 KB	-rw-r--r--
asin_s390x.s	4.16 KB	-rw-r--r--
asinh.go	1.85 KB	-rw-r--r--
asinh_s390x.s	5.74 KB	-rw-r--r--
atan.go	2.97 KB	-rw-r--r--
atan2.go	1.45 KB	-rw-r--r--
atan2_386.s	396 B	-rw-r--r--
atan2_s390x.s	6.88 KB	-rw-r--r--
atan_386.s	376 B	-rw-r--r--
atan_s390x.s	3.69 KB	-rw-r--r--
atanh.go	1.93 KB	-rw-r--r--
atanh_s390x.s	5.36 KB	-rw-r--r--
bits.go	1.87 KB	-rw-r--r--
cbrt.go	2.25 KB	-rw-r--r--
cbrt_s390x.s	4.89 KB	-rw-r--r--
const.go	2.18 KB	-rw-r--r--
copysign.go	378 B	-rw-r--r--
cosh_s390x.s	5.59 KB	-rw-r--r--
dim.go	1.55 KB	-rw-r--r--
dim_amd64.s	1.9 KB	-rw-r--r--
dim_arm64.s	947 B	-rw-r--r--
dim_s390x.s	1.96 KB	-rw-r--r--
erf.go	11.38 KB	-rw-r--r--
erf_s390x.s	8.5 KB	-rw-r--r--
erfc_s390x.s	14.4 KB	-rw-r--r--
erfinv.go	3.36 KB	-rw-r--r--
example_test.go	3.33 KB	-rw-r--r--
exp.go	5.25 KB	-rw-r--r--
exp2_386.s	1.12 KB	-rw-r--r--
exp_amd64.s	4.23 KB	-rw-r--r--
exp_arm64.s	5.35 KB	-rw-r--r--
exp_asm.go	260 B	-rw-r--r--
exp_s390x.s	4.65 KB	-rw-r--r--
expm1.go	7.84 KB	-rw-r--r--
expm1_386.s	1.84 KB	-rw-r--r--
expm1_s390x.s	5.29 KB	-rw-r--r--
export_s390x_test.go	732 B	-rw-r--r--
export_test.go	357 B	-rw-r--r--
floor.go	3.09 KB	-rw-r--r--
floor_386.s	1.45 KB	-rw-r--r--
floor_amd64.s	1.98 KB	-rw-r--r--
floor_arm64.s	549 B	-rw-r--r--
floor_ppc64x.s	483 B	-rw-r--r--
floor_s390x.s	555 B	-rw-r--r--
floor_wasm.s	447 B	-rw-r--r--
fma.go	4.49 KB	-rw-r--r--
frexp.go	860 B	-rw-r--r--
frexp_386.s	689 B	-rw-r--r--
gamma.go	5.52 KB	-rw-r--r--
huge_test.go	2.24 KB	-rw-r--r--
hypot.go	773 B	-rw-r--r--
hypot_386.s	1.8 KB	-rw-r--r--
hypot_amd64.s	1.05 KB	-rw-r--r--
j0.go	13.6 KB	-rw-r--r--
j1.go	13.3 KB	-rw-r--r--
jn.go	7.17 KB	-rw-r--r--
ldexp.go	990 B	-rw-r--r--
ldexp_386.s	441 B	-rw-r--r--
lgamma.go	11.02 KB	-rw-r--r--
log.go	3.8 KB	-rw-r--r--
log10.go	740 B	-rw-r--r--
log10_386.s	596 B	-rw-r--r--
log10_s390x.s	4.73 KB	-rw-r--r--
log1p.go	6.3 KB	-rw-r--r--
log1p_386.s	860 B	-rw-r--r--
log1p_s390x.s	5.15 KB	-rw-r--r--
log_386.s	394 B	-rw-r--r--
log_amd64.s	3.66 KB	-rw-r--r--
log_s390x.s	4.31 KB	-rw-r--r--
logb.go	1014 B	-rw-r--r--
mod.go	834 B	-rw-r--r--
mod_386.s	531 B	-rw-r--r--
modf.go	847 B	-rw-r--r--
modf_386.s	1.07 KB	-rw-r--r--
modf_arm64.s	439 B	-rw-r--r--
modf_ppc64x.s	404 B	-rw-r--r--
nextafter.go	1.2 KB	-rw-r--r--
pow.go	3.16 KB	-rw-r--r--
pow10.go	1.24 KB	-rw-r--r--
pow_s390x.s	16.27 KB	-rw-r--r--
remainder.go	1.96 KB	-rw-r--r--
remainder_386.s	543 B	-rw-r--r--
signbit.go	302 B	-rw-r--r--
sin.go	6.23 KB	-rw-r--r--
sin_s390x.s	8.34 KB	-rw-r--r--
sincos.go	1.75 KB	-rw-r--r--
sinh.go	1.56 KB	-rw-r--r--
sinh_s390x.s	5.98 KB	-rw-r--r--
sqrt.go	4.84 KB	-rw-r--r--
sqrt_386.s	296 B	-rw-r--r--
sqrt_amd64.s	326 B	-rw-r--r--
sqrt_arm.s	521 B	-rw-r--r--
sqrt_arm64.s	302 B	-rw-r--r--
sqrt_mipsx.s	375 B	-rw-r--r--
sqrt_ppc64x.s	326 B	-rw-r--r--
sqrt_s390x.s	301 B	-rw-r--r--
sqrt_wasm.s	269 B	-rw-r--r--
stubs_386.s	908 B	-rw-r--r--
stubs_amd64.s	1.41 KB	-rw-r--r--
stubs_arm.s	1.7 KB	-rw-r--r--
stubs_arm64.s	1.37 KB	-rw-r--r--
stubs_mips64x.s	1.83 KB	-rw-r--r--
stubs_mipsx.s	1.79 KB	-rw-r--r--
stubs_ppc64x.s	1.58 KB	-rw-r--r--
stubs_riscv64.s	1.75 KB	-rw-r--r--
stubs_s390x.s	13.24 KB	-rw-r--r--
stubs_wasm.s	1.63 KB	-rw-r--r--
tan.go	3.61 KB	-rw-r--r--
tan_s390x.s	2.73 KB	-rw-r--r--
tanh.go	2.59 KB	-rw-r--r--
tanh_s390x.s	4.57 KB	-rw-r--r--
trig_reduce.go	2.94 KB	-rw-r--r--
unsafe.go	1.27 KB	-rw-r--r--

Code Editor : expm1.go

// Copyright 2010 The Go Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.

package math

// The original C code, the long comment, and the constants
// below are from FreeBSD's /usr/src/lib/msun/src/s_expm1.c
// and came with this notice. The go code is a simplified
// version of the original C.
//
// ====================================================
// 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.
// ====================================================
//
// expm1(x)
// Returns exp(x)-1, the exponential of x minus 1.
//
// Method
//   1. Argument reduction:
//      Given x, find r and integer k such that
//
//               x = k*ln2 + r,  |r| <= 0.5*ln2 ~ 0.34658
//
//      Here a correction term c will be computed to compensate
//      the error in r when rounded to a floating-point number.
//
//   2. Approximating expm1(r) by a special rational function on
//      the interval [0,0.34658]:
//      Since
//          r*(exp(r)+1)/(exp(r)-1) = 2+ r**2/6 - r**4/360 + ...
//      we define R1(r*r) by
//          r*(exp(r)+1)/(exp(r)-1) = 2+ r**2/6 * R1(r*r)
//      That is,
//          R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r)
//                   = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r))
//                   = 1 - r**2/60 + r**4/2520 - r**6/100800 + ...
//      We use a special Reme algorithm on [0,0.347] to generate
//      a polynomial of degree 5 in r*r to approximate R1. The
//      maximum error of this polynomial approximation is bounded
//      by 2**-61. In other words,
//          R1(z) ~ 1.0 + Q1*z + Q2*z**2 + Q3*z**3 + Q4*z**4 + Q5*z**5
//      where   Q1  =  -1.6666666666666567384E-2,
//              Q2  =   3.9682539681370365873E-4,
//              Q3  =  -9.9206344733435987357E-6,
//              Q4  =   2.5051361420808517002E-7,
//              Q5  =  -6.2843505682382617102E-9;
//      (where z=r*r, and the values of Q1 to Q5 are listed below)
//      with error bounded by
//          |                  5           |     -61
//          | 1.0+Q1*z+...+Q5*z   -  R1(z) | <= 2
//          |                              |
//
//      expm1(r) = exp(r)-1 is then computed by the following
//      specific way which minimize the accumulation rounding error:
//                             2     3
//                            r     r    [ 3 - (R1 + R1*r/2)  ]
//            expm1(r) = r + --- + --- * [--------------------]
//                            2     2    [ 6 - r*(3 - R1*r/2) ]
//
//      To compensate the error in the argument reduction, we use
//              expm1(r+c) = expm1(r) + c + expm1(r)*c
//                         ~ expm1(r) + c + r*c
//      Thus c+r*c will be added in as the correction terms for
//      expm1(r+c). Now rearrange the term to avoid optimization
//      screw up:
//                      (      2                                    2 )
//                      ({  ( r    [ R1 -  (3 - R1*r/2) ]  )  }    r  )
//       expm1(r+c)~r - ({r*(--- * [--------------------]-c)-c} - --- )
//                      ({  ( 2    [ 6 - r*(3 - R1*r/2) ]  )  }    2  )
//                      (                                             )
//
//                 = r - E
//   3. Scale back to obtain expm1(x):
//      From step 1, we have
//         expm1(x) = either 2**k*[expm1(r)+1] - 1
//                  = or     2**k*[expm1(r) + (1-2**-k)]
//   4. Implementation notes:
//      (A). To save one multiplication, we scale the coefficient Qi
//           to Qi*2**i, and replace z by (x**2)/2.
//      (B). To achieve maximum accuracy, we compute expm1(x) by
//        (i)   if x < -56*ln2, return -1.0, (raise inexact if x!=inf)
//        (ii)  if k=0, return r-E
//        (iii) if k=-1, return 0.5*(r-E)-0.5
//        (iv)  if k=1 if r < -0.25, return 2*((r+0.5)- E)
//                     else          return  1.0+2.0*(r-E);
//        (v)   if (k<-2||k>56) return 2**k(1-(E-r)) - 1 (or exp(x)-1)
//        (vi)  if k <= 20, return 2**k((1-2**-k)-(E-r)), else
//        (vii) return 2**k(1-((E+2**-k)-r))
//
// Special cases:
//      expm1(INF) is INF, expm1(NaN) is NaN;
//      expm1(-INF) is -1, and
//      for finite argument, only expm1(0)=0 is exact.
//
// Accuracy:
//      according to an error analysis, the error is always less than
//      1 ulp (unit in the last place).
//
// Misc. info.
//      For IEEE double
//          if x >  7.09782712893383973096e+02 then expm1(x) overflow
//
// Constants:
// The hexadecimal values are the intended ones for the following
// constants. The decimal values may be used, provided that the
// compiler will convert from decimal to binary accurately enough
// to produce the hexadecimal values shown.
//

// Expm1 returns e**x - 1, the base-e exponential of x minus 1.
// It is more accurate than Exp(x) - 1 when x is near zero.
//
// Special cases are:
//	Expm1(+Inf) = +Inf
//	Expm1(-Inf) = -1
//	Expm1(NaN) = NaN
// Very large values overflow to -1 or +Inf.
func Expm1(x float64) float64

func expm1(x float64) float64 {
	const (
		Othreshold = 7.09782712893383973096e+02 // 0x40862E42FEFA39EF
		Ln2X56     = 3.88162421113569373274e+01 // 0x4043687a9f1af2b1
		Ln2HalfX3  = 1.03972077083991796413e+00 // 0x3ff0a2b23f3bab73
		Ln2Half    = 3.46573590279972654709e-01 // 0x3fd62e42fefa39ef
		Ln2Hi      = 6.93147180369123816490e-01 // 0x3fe62e42fee00000
		Ln2Lo      = 1.90821492927058770002e-10 // 0x3dea39ef35793c76
		InvLn2     = 1.44269504088896338700e+00 // 0x3ff71547652b82fe
		Tiny       = 1.0 / (1 << 54)            // 2**-54 = 0x3c90000000000000
		// scaled coefficients related to expm1
		Q1 = -3.33333333333331316428e-02 // 0xBFA11111111110F4
		Q2 = 1.58730158725481460165e-03  // 0x3F5A01A019FE5585
		Q3 = -7.93650757867487942473e-05 // 0xBF14CE199EAADBB7
		Q4 = 4.00821782732936239552e-06  // 0x3ED0CFCA86E65239
		Q5 = -2.01099218183624371326e-07 // 0xBE8AFDB76E09C32D
	)

// special cases
	switch {
	case IsInf(x, 1) || IsNaN(x):
		return x
	case IsInf(x, -1):
		return -1
	}

absx := x
	sign := false
	if x < 0 {
		absx = -absx
		sign = true
	}

// filter out huge argument
	if absx >= Ln2X56 { // if |x| >= 56 * ln2
		if sign {
			return -1 // x < -56*ln2, return -1
		}
		if absx >= Othreshold { // if |x| >= 709.78...
			return Inf(1)
		}
	}

// argument reduction
	var c float64
	var k int
	if absx > Ln2Half { // if  |x| > 0.5 * ln2
		var hi, lo float64
		if absx < Ln2HalfX3 { // and |x| < 1.5 * ln2
			if !sign {
				hi = x - Ln2Hi
				lo = Ln2Lo
				k = 1
			} else {
				hi = x + Ln2Hi
				lo = -Ln2Lo
				k = -1
			}
		} else {
			if !sign {
				k = int(InvLn2*x + 0.5)
			} else {
				k = int(InvLn2*x - 0.5)
			}
			t := float64(k)
			hi = x - t*Ln2Hi // t * Ln2Hi is exact here
			lo = t * Ln2Lo
		}
		x = hi - lo
		c = (hi - x) - lo
	} else if absx < Tiny { // when |x| < 2**-54, return x
		return x
	} else {
		k = 0
	}

// x is now in primary range
	hfx := 0.5 * x
	hxs := x * hfx
	r1 := 1 + hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5))))
	t := 3 - r1*hfx
	e := hxs * ((r1 - t) / (6.0 - x*t))
	if k == 0 {
		return x - (x*e - hxs) // c is 0
	}
	e = (x*(e-c) - c)
	e -= hxs
	switch {
	case k == -1:
		return 0.5*(x-e) - 0.5
	case k == 1:
		if x < -0.25 {
			return -2 * (e - (x + 0.5))
		}
		return 1 + 2*(x-e)
	case k <= -2 || k > 56: // suffice to return exp(x)-1
		y := 1 - (e - x)
		y = Float64frombits(Float64bits(y) + uint64(k)<<52) // add k to y's exponent
		return y - 1
	}
	if k < 20 {
		t := Float64frombits(0x3ff0000000000000 - (0x20000000000000 >> uint(k))) // t=1-2**-k
		y := t - (e - x)
		y = Float64frombits(Float64bits(y) + uint64(k)<<52) // add k to y's exponent
		return y
	}
	t = Float64frombits(uint64(0x3ff-k) << 52) // 2**-k
	y := x - (e + t)
	y++
	y = Float64frombits(Float64bits(y) + uint64(k)<<52) // add k to y's exponent
	return y
}

${this.title}

Code Editor : expm1.go