v / vlib / math
Raw file | 113 loc (109 sloc) | 2.5 KB | Latest commit hash 1cfc4198f
1module math
2
3const (
4 tan_p = [
5 -1.30936939181383777646e+4,
6 1.15351664838587416140e+6,
7 -1.79565251976484877988e+7,
8 ]
9 tan_q = [
10 1.00000000000000000000e+0,
11 1.36812963470692954678e+4,
12 -1.32089234440210967447e+6,
13 2.50083801823357915839e+7,
14 -5.38695755929454629881e+7,
15 ]
16 tan_dp1 = 7.853981554508209228515625e-1
17 tan_dp2 = 7.94662735614792836714e-9
18 tan_dp3 = 3.06161699786838294307e-17
19 tan_lossth = 1.073741824e+9
20)
21
22// tan calculates tangent of a number
23pub fn tan(a f64) f64 {
24 mut x := a
25 if x == 0.0 || is_nan(x) {
26 return x
27 }
28 if is_inf(x, 0) {
29 return nan()
30 }
31 mut sign := 1 // make argument positive but save the sign
32 if x < 0 {
33 x = -x
34 sign = -1
35 }
36 if x > math.tan_lossth {
37 return 0.0
38 }
39 // compute x mod pi_4
40 mut y := floor(x * 4.0 / pi) // strip high bits of integer part
41 mut z := ldexp(y, -3)
42 z = floor(z) // integer part of y/8
43 z = y - ldexp(z, 3) // y - 16 * (y/16) // integer and fractional part modulo one octant
44 mut octant := int(z) // map zeros and singularities to origin
45 if (octant & 1) == 1 {
46 octant++
47 y += 1.0
48 }
49 z = ((x - y * math.tan_dp1) - y * math.tan_dp2) - y * math.tan_dp3
50 zz := z * z
51 if zz > 1.0e-14 {
52 y = z + z * (zz * (((math.tan_p[0] * zz) + math.tan_p[1]) * zz + math.tan_p[2]) / ((((zz +
53 math.tan_q[1]) * zz + math.tan_q[2]) * zz + math.tan_q[3]) * zz + math.tan_q[4]))
54 } else {
55 y = z
56 }
57 if (octant & 2) == 2 {
58 y = -1.0 / y
59 }
60 if sign < 0 {
61 y = -y
62 }
63 return y
64}
65
66// tanf calculates tangent. (float32)
67[inline]
68pub fn tanf(a f32) f32 {
69 return f32(tan(a))
70}
71
72// tan calculates cotangent of a number
73pub fn cot(a f64) f64 {
74 mut x := a
75 if x == 0.0 {
76 return inf(1)
77 }
78 mut sign := 1 // make argument positive but save the sign
79 if x < 0 {
80 x = -x
81 sign = -1
82 }
83 if x > math.tan_lossth {
84 return 0.0
85 }
86 // compute x mod pi_4
87 mut y := floor(x * 4.0 / pi) // strip high bits of integer part
88 mut z := ldexp(y, -3)
89 z = floor(z) // integer part of y/8
90 z = y - ldexp(z, 3) // y - 16 * (y/16) // integer and fractional part modulo one octant
91 mut octant := int(z) // map zeros and singularities to origin
92 if (octant & 1) == 1 {
93 octant++
94 y += 1.0
95 }
96 z = ((x - y * math.tan_dp1) - y * math.tan_dp2) - y * math.tan_dp3
97 zz := z * z
98 if zz > 1.0e-14 {
99 y = z + z * (zz * (((math.tan_p[0] * zz) + math.tan_p[1]) * zz + math.tan_p[2]) / ((((zz +
100 math.tan_q[1]) * zz + math.tan_q[2]) * zz + math.tan_q[3]) * zz + math.tan_q[4]))
101 } else {
102 y = z
103 }
104 if (octant & 2) == 2 {
105 y = -y
106 } else {
107 y = 1.0 / y
108 }
109 if sign < 0 {
110 y = -y
111 }
112 return y
113}