v/vlib/math/tan.v

module math

const (
	tan_p      = [
		-1.30936939181383777646e+4,
		1.15351664838587416140e+6,
		-1.79565251976484877988e+7,
	]
	tan_q      = [
		1.00000000000000000000e+0,
		1.36812963470692954678e+4,
		-1.32089234440210967447e+6,
		2.50083801823357915839e+7,
		-5.38695755929454629881e+7,
	]
	tan_dp1    = 7.853981554508209228515625e-1
	tan_dp2    = 7.94662735614792836714e-9
	tan_dp3    = 3.06161699786838294307e-17
	tan_lossth = 1.073741824e+9
)

// tan calculates tangent of a number
pub fn tan(a f64) f64 {
	mut x := a
	if x == 0.0 || is_nan(x) {
		return x
	}
	if is_inf(x, 0) {
		return nan()
	}
	mut sign := 1 // make argument positive but save the sign
	if x < 0 {
		x = -x
		sign = -1
	}
	if x > math.tan_lossth {
		return 0.0
	}
	// compute x mod pi_4
	mut y := floor(x * 4.0 / pi) // strip high bits of integer part
	mut z := ldexp(y, -3)
	z = floor(z) // integer part of y/8
	z = y - ldexp(z, 3) // y - 16 * (y/16) // integer and fractional part modulo one octant
	mut octant := int(z) // map zeros and singularities to origin
	if (octant & 1) == 1 {
		octant++
		y += 1.0
	}
	z = ((x - y * math.tan_dp1) - y * math.tan_dp2) - y * math.tan_dp3
	zz := z * z
	if zz > 1.0e-14 {
		y = z + z * (zz * (((math.tan_p[0] * zz) + math.tan_p[1]) * zz + math.tan_p[2]) / ((((zz +
			math.tan_q[1]) * zz + math.tan_q[2]) * zz + math.tan_q[3]) * zz + math.tan_q[4]))
	} else {
		y = z
	}
	if (octant & 2) == 2 {
		y = -1.0 / y
	}
	if sign < 0 {
		y = -y
	}
	return y
}

// tanf calculates tangent. (float32)
[inline]
pub fn tanf(a f32) f32 {
	return f32(tan(a))
}

// tan calculates cotangent of a number
pub fn cot(a f64) f64 {
	mut x := a
	if x == 0.0 {
		return inf(1)
	}
	mut sign := 1 // make argument positive but save the sign
	if x < 0 {
		x = -x
		sign = -1
	}
	if x > math.tan_lossth {
		return 0.0
	}
	// compute x mod pi_4
	mut y := floor(x * 4.0 / pi) // strip high bits of integer part
	mut z := ldexp(y, -3)
	z = floor(z) // integer part of y/8
	z = y - ldexp(z, 3) // y - 16 * (y/16) // integer and fractional part modulo one octant
	mut octant := int(z) // map zeros and singularities to origin
	if (octant & 1) == 1 {
		octant++
		y += 1.0
	}
	z = ((x - y * math.tan_dp1) - y * math.tan_dp2) - y * math.tan_dp3
	zz := z * z
	if zz > 1.0e-14 {
		y = z + z * (zz * (((math.tan_p[0] * zz) + math.tan_p[1]) * zz + math.tan_p[2]) / ((((zz +
			math.tan_q[1]) * zz + math.tan_q[2]) * zz + math.tan_q[3]) * zz + math.tan_q[4]))
	} else {
		y = z
	}
	if (octant & 2) == 2 {
		y = -y
	} else {
		y = 1.0 / y
	}
	if sign < 0 {
		y = -y
	}
	return y
}