| 1 | // Copyright (c) 2019-2024 Alexander Medvednikov. All rights reserved. |
| 2 | // Use of this source code is governed by an MIT license |
| 3 | // that can be found in the LICENSE file. |
| 4 | module math |
| 5 | |
| 6 | // aprox_sin returns an approximation of sin(a) made using lolremez |
| 7 | pub fn aprox_sin(a f64) f64 { |
| 8 | a0 := 1.91059300966915117e-31 |
| 9 | a1 := 1.00086760103908896 |
| 10 | a2 := -1.21276126894734565e-2 |
| 11 | a3 := -1.38078780785773762e-1 |
| 12 | a4 := -2.67353392911981221e-2 |
| 13 | a5 := 2.08026600266304389e-2 |
| 14 | a6 := -3.03996055049204407e-3 |
| 15 | a7 := 1.38235642404333740e-4 |
| 16 | tmp := a4 + a * (a5 + a * (a6 + a * a7)) |
| 17 | return a0 + a * (a1 + a * (a2 + a * (a3 + a * tmp))) |
| 18 | } |
| 19 | |
| 20 | // aprox_cos returns an approximation of cos(a) made using lolremez |
| 21 | pub fn aprox_cos(a f64) f64 { |
| 22 | a0 := 9.9995999154986614e-1 |
| 23 | a1 := 1.2548995793001028e-3 |
| 24 | a2 := -5.0648546280678015e-1 |
| 25 | a3 := 1.2942246466519995e-2 |
| 26 | a4 := 2.8668384702547972e-2 |
| 27 | a5 := 7.3726485210586547e-3 |
| 28 | a6 := -3.8510875386947414e-3 |
| 29 | a7 := 4.7196604604366623e-4 |
| 30 | a8 := -1.8776444013090451e-5 |
| 31 | tmp := a4 + a * (a5 + a * (a6 + a * (a7 + a * a8))) |
| 32 | return a0 + a * (a1 + a * (a2 + a * (a3 + a * tmp))) |
| 33 | } |
| 34 | |
| 35 | // copysign returns a value with the magnitude of x and the sign of y |
| 36 | @[inline] |
| 37 | pub fn copysign(x f64, y f64) f64 { |
| 38 | return f64_from_bits((f64_bits(x) & ~sign_mask) | (f64_bits(y) & sign_mask)) |
| 39 | } |
| 40 | |
| 41 | // degrees converts an angle in radians to a corresponding angle in degrees. |
| 42 | @[inline] |
| 43 | pub fn degrees(radians f64) f64 { |
| 44 | return radians * (180.0 / pi) |
| 45 | } |
| 46 | |
| 47 | // radians converts an angle in degrees to a corresponding angle in radians. |
| 48 | @[inline] |
| 49 | pub fn radians(degrees f64) f64 { |
| 50 | return degrees * (pi / 180.0) |
| 51 | } |
| 52 | |
| 53 | // angle_diff calculates the difference between angles in radians |
| 54 | @[inline] |
| 55 | pub fn angle_diff(radian_a f64, radian_b f64) f64 { |
| 56 | mut delta := fmod(radian_b - radian_a, tau) |
| 57 | delta = fmod(delta + 1.5 * tau, tau) |
| 58 | delta -= .5 * tau |
| 59 | return delta |
| 60 | } |
| 61 | |
| 62 | @[params] |
| 63 | pub struct DigitParams { |
| 64 | pub: |
| 65 | base int = 10 |
| 66 | reverse bool |
| 67 | } |
| 68 | |
| 69 | // digits returns an array of the digits of `num` in the given optional `base`. |
| 70 | // The `num` argument accepts any integer type (i8|i16|int|isize|i64), and will be cast to i64 |
| 71 | // The `base:` argument is optional, it will default to base: 10. |
| 72 | // This function returns an array of the digits in reverse order i.e.: |
| 73 | // Example: assert math.digits(12345, base: 10) == [5,4,3,2,1] |
| 74 | // You can also use it, with an explicit `reverse: true` parameter, |
| 75 | // (it will do a reverse of the result array internally => slower): |
| 76 | // Example: assert math.digits(12345, reverse: true) == [1,2,3,4,5] |
| 77 | pub fn digits(num i64, params DigitParams) []int { |
| 78 | // set base to 10 initially and change only if base is explicitly set. |
| 79 | mut b := params.base |
| 80 | if b < 2 { |
| 81 | panic_n('digits: Cannot find digits of n with base:', b) |
| 82 | } |
| 83 | mut n := num |
| 84 | mut sgn := 1 |
| 85 | if n < 0 { |
| 86 | sgn = -1 |
| 87 | n = -n |
| 88 | } |
| 89 | |
| 90 | mut res := []int{} |
| 91 | if n == 0 { |
| 92 | // short-circuit and return 0 |
| 93 | res << 0 |
| 94 | return res |
| 95 | } |
| 96 | for n != 0 { |
| 97 | next_n := n / b |
| 98 | res << int(n - next_n * b) |
| 99 | n = next_n |
| 100 | } |
| 101 | |
| 102 | if sgn == -1 { |
| 103 | res[res.len - 1] *= sgn |
| 104 | } |
| 105 | |
| 106 | if params.reverse { |
| 107 | res = res.reverse() |
| 108 | } |
| 109 | |
| 110 | return res |
| 111 | } |
| 112 | |
| 113 | // count_digits return the number of digits in the number passed. |
| 114 | // Number argument accepts any integer type (i8|i16|int|isize|i64) and will be cast to i64 |
| 115 | pub fn count_digits(number i64) int { |
| 116 | mut n := number |
| 117 | if n == 0 { |
| 118 | return 1 |
| 119 | } |
| 120 | mut c := 0 |
| 121 | for n != 0 { |
| 122 | n = n / 10 |
| 123 | c++ |
| 124 | } |
| 125 | return c |
| 126 | } |
| 127 | |
| 128 | // minmax returns the minimum and maximum value of the two provided. |
| 129 | pub fn minmax(a f64, b f64) (f64, f64) { |
| 130 | if a < b { |
| 131 | return a, b |
| 132 | } |
| 133 | return b, a |
| 134 | } |
| 135 | |
| 136 | // clamp returns x constrained between a and b |
| 137 | @[inline] |
| 138 | pub fn clamp(x f64, a f64, b f64) f64 { |
| 139 | if x < a { |
| 140 | return a |
| 141 | } |
| 142 | if x > b { |
| 143 | return b |
| 144 | } |
| 145 | return x |
| 146 | } |
| 147 | |
| 148 | // sign returns the corresponding sign -1.0, 1.0 of the provided number. |
| 149 | // if n is not a number, its sign is nan too. |
| 150 | @[inline] |
| 151 | pub fn sign(n f64) f64 { |
| 152 | // dump(n) |
| 153 | if is_nan(n) { |
| 154 | return nan() |
| 155 | } |
| 156 | return copysign(1.0, n) |
| 157 | } |
| 158 | |
| 159 | // signi returns the corresponding sign -1, 1 of the provided number. |
| 160 | @[inline] |
| 161 | pub fn signi(n f64) int { |
| 162 | return int(copysign(1.0, n)) |
| 163 | } |
| 164 | |
| 165 | // signbit returns a value with the boolean representation of the sign for x |
| 166 | @[inline] |
| 167 | pub fn signbit(x f64) bool { |
| 168 | return f64_bits(x) & sign_mask != 0 |
| 169 | } |
| 170 | |
| 171 | // tolerance checks if the difference between `actual` and `expected` numbers, is less than or equal to the tolerance value `tol`. |
| 172 | // Note: the `actual` and `expected` parameters are NOT symmetrical. |
| 173 | // If you compare against a known expected number, put it in the *second* argument, especially if it is 0. |
| 174 | pub fn tolerance(actual f64, expected f64, tol f64) bool { |
| 175 | // Check actual==expected so that at least if they both are small and identical we say they match. |
| 176 | if actual == expected { |
| 177 | return true |
| 178 | } |
| 179 | mut d := actual - expected |
| 180 | if d < 0 { |
| 181 | d = -d |
| 182 | } |
| 183 | mut ee := tol |
| 184 | // make error tolerance a fraction of expected, not actual. |
| 185 | if expected != 0 { |
| 186 | // Multiplying by ee here can underflow denormal values to zero. |
| 187 | ee *= expected |
| 188 | if ee < 0 { |
| 189 | ee = -ee |
| 190 | } |
| 191 | } |
| 192 | return d < ee |
| 193 | } |
| 194 | |
| 195 | // close checks if `actual` and `expected` are within 1e-14 of each other |
| 196 | // Note: the `actual` and `expected` parameters are NOT symmetrical. |
| 197 | // If you compare against a known expected number, put it in the *second* argument, especially if it is 0. |
| 198 | pub fn close(actual f64, expected f64) bool { |
| 199 | return tolerance(actual, expected, 1e-14) |
| 200 | } |
| 201 | |
| 202 | // veryclose checks if a and b are within 4e-16 of each other |
| 203 | // Note: the `actual` and `expected` parameters are NOT symmetrical. |
| 204 | // If you compare against a known expected number, put it in the *second* argument, especially if it is 0. |
| 205 | pub fn veryclose(actual f64, expected f64) bool { |
| 206 | return tolerance(actual, expected, 4e-16) |
| 207 | } |
| 208 | |
| 209 | // alike checks if a and b are equal |
| 210 | pub fn alike(a f64, b f64) bool { |
| 211 | // eprintln('>>> a: ${f64_bits(a):20} | b: ${f64_bits(b):20} | a==b: ${a == b} | a: ${a:10} | b: ${b:10}') |
| 212 | // compare a and b, ignoring their last 2 bits: |
| 213 | if f64_bits(a) & 0xFFFF_FFFF_FFFF_FFFC == f64_bits(b) & 0xFFFF_FFFF_FFFF_FFFC { |
| 214 | return true |
| 215 | } |
| 216 | if a == -0 && b == 0 { |
| 217 | return true |
| 218 | } |
| 219 | if a == 0 && b == -0 { |
| 220 | return true |
| 221 | } |
| 222 | if is_nan(a) && is_nan(b) { |
| 223 | return true |
| 224 | } |
| 225 | if a == b { |
| 226 | return signbit(a) == signbit(b) |
| 227 | } |
| 228 | return false |
| 229 | } |
| 230 | |