2021-01-18 13:20:06 +01:00
|
|
|
// Copyright (c) 2019-2021 Alexander Medvednikov. All rights reserved.
|
2019-06-22 20:20:28 +02:00
|
|
|
// Use of this source code is governed by an MIT license
|
|
|
|
// that can be found in the LICENSE file.
|
|
|
|
module math
|
|
|
|
|
2021-02-27 10:18:26 +01:00
|
|
|
// aprox_sin returns an approximation of sin(a) made using lolremez
|
2020-05-07 07:47:24 +02:00
|
|
|
pub fn aprox_sin(a f64) f64 {
|
|
|
|
a0 := 1.91059300966915117e-31
|
|
|
|
a1 := 1.00086760103908896
|
|
|
|
a2 := -1.21276126894734565e-2
|
|
|
|
a3 := -1.38078780785773762e-1
|
|
|
|
a4 := -2.67353392911981221e-2
|
|
|
|
a5 := 2.08026600266304389e-2
|
|
|
|
a6 := -3.03996055049204407e-3
|
|
|
|
a7 := 1.38235642404333740e-4
|
|
|
|
return a0 + a * (a1 + a * (a2 + a * (a3 + a * (a4 + a * (a5 + a * (a6 + a * a7))))))
|
2019-06-22 20:20:28 +02:00
|
|
|
}
|
|
|
|
|
2021-02-27 10:18:26 +01:00
|
|
|
// aprox_cos returns an approximation of sin(a) made using lolremez
|
2020-05-07 07:47:24 +02:00
|
|
|
pub fn aprox_cos(a f64) f64 {
|
|
|
|
a0 := 9.9995999154986614e-1
|
|
|
|
a1 := 1.2548995793001028e-3
|
|
|
|
a2 := -5.0648546280678015e-1
|
|
|
|
a3 := 1.2942246466519995e-2
|
|
|
|
a4 := 2.8668384702547972e-2
|
|
|
|
a5 := 7.3726485210586547e-3
|
|
|
|
a6 := -3.8510875386947414e-3
|
|
|
|
a7 := 4.7196604604366623e-4
|
|
|
|
a8 := -1.8776444013090451e-5
|
|
|
|
return a0 + a * (a1 + a * (a2 + a * (a3 + a * (a4 + a * (a5 + a * (a6 + a * (a7 + a * a8)))))))
|
2020-02-26 13:31:54 +01:00
|
|
|
}
|
|
|
|
|
2020-05-07 07:47:24 +02:00
|
|
|
// copysign returns a value with the magnitude of x and the sign of y
|
2021-05-08 13:32:18 +02:00
|
|
|
[inline]
|
2020-10-15 14:42:16 +02:00
|
|
|
pub fn copysign(x f64, y f64) f64 {
|
2020-05-07 07:47:24 +02:00
|
|
|
return f64_from_bits((f64_bits(x) & ~sign_mask) | (f64_bits(y) & sign_mask))
|
2019-06-23 08:19:37 +02:00
|
|
|
}
|
|
|
|
|
2019-07-12 07:46:40 +02:00
|
|
|
// degrees convert from degrees to radians.
|
2021-05-08 13:32:18 +02:00
|
|
|
[inline]
|
2019-07-12 07:46:40 +02:00
|
|
|
pub fn degrees(radians f64) f64 {
|
2019-10-12 21:31:05 +02:00
|
|
|
return radians * (180.0 / pi)
|
2019-07-12 07:46:40 +02:00
|
|
|
}
|
|
|
|
|
2019-07-02 12:50:33 +02:00
|
|
|
// digits returns an array of the digits of n in the given base.
|
2020-10-15 14:42:16 +02:00
|
|
|
pub fn digits(_n int, base int) []int {
|
2019-11-11 03:21:47 +01:00
|
|
|
if base < 2 {
|
|
|
|
panic('digits: Cannot find digits of n with base $base')
|
|
|
|
}
|
2019-09-16 20:26:05 +02:00
|
|
|
mut n := _n
|
2019-07-02 12:50:33 +02:00
|
|
|
mut sign := 1
|
|
|
|
if n < 0 {
|
|
|
|
sign = -1
|
|
|
|
n = -n
|
|
|
|
}
|
2020-04-26 11:42:44 +02:00
|
|
|
mut res := []int{}
|
2019-07-02 12:50:33 +02:00
|
|
|
for n != 0 {
|
|
|
|
res << (n % base) * sign
|
|
|
|
n /= base
|
|
|
|
}
|
|
|
|
return res
|
|
|
|
}
|
|
|
|
|
2021-05-08 13:32:18 +02:00
|
|
|
[inline]
|
2020-05-14 16:58:20 +02:00
|
|
|
pub fn fabs(x f64) f64 {
|
|
|
|
if x < 0.0 {
|
|
|
|
return -x
|
|
|
|
}
|
|
|
|
return x
|
|
|
|
}
|
|
|
|
|
2019-07-02 12:50:33 +02:00
|
|
|
// gcd calculates greatest common (positive) divisor (or zero if a and b are both zero).
|
2020-10-15 14:42:16 +02:00
|
|
|
pub fn gcd(a_ i64, b_ i64) i64 {
|
2019-09-16 20:26:05 +02:00
|
|
|
mut a := a_
|
|
|
|
mut b := b_
|
2019-06-29 17:24:55 +02:00
|
|
|
if a < 0 {
|
|
|
|
a = -a
|
|
|
|
}
|
|
|
|
if b < 0 {
|
|
|
|
b = -b
|
|
|
|
}
|
|
|
|
for b != 0 {
|
|
|
|
a %= b
|
|
|
|
if a == 0 {
|
|
|
|
return b
|
|
|
|
}
|
|
|
|
b %= a
|
|
|
|
}
|
|
|
|
return a
|
|
|
|
}
|
|
|
|
|
|
|
|
// lcm calculates least common (non-negative) multiple.
|
2020-10-15 14:42:16 +02:00
|
|
|
pub fn lcm(a i64, b i64) i64 {
|
2019-06-29 17:24:55 +02:00
|
|
|
if a == 0 {
|
|
|
|
return a
|
|
|
|
}
|
|
|
|
res := a * (b / gcd(b, a))
|
|
|
|
if res < 0 {
|
|
|
|
return -res
|
|
|
|
}
|
|
|
|
return res
|
|
|
|
}
|
|
|
|
|
2019-07-02 12:50:33 +02:00
|
|
|
// max returns the maximum value of the two provided.
|
2021-05-08 13:32:18 +02:00
|
|
|
[inline]
|
2020-10-15 14:42:16 +02:00
|
|
|
pub fn max(a f64, b f64) f64 {
|
2019-06-22 20:20:28 +02:00
|
|
|
if a > b {
|
|
|
|
return a
|
|
|
|
}
|
|
|
|
return b
|
|
|
|
}
|
|
|
|
|
2019-07-12 20:45:56 +02:00
|
|
|
// min returns the minimum value of the two provided.
|
2021-05-08 13:32:18 +02:00
|
|
|
[inline]
|
2020-10-15 14:42:16 +02:00
|
|
|
pub fn min(a f64, b f64) f64 {
|
2019-06-22 20:20:28 +02:00
|
|
|
if a < b {
|
|
|
|
return a
|
|
|
|
}
|
|
|
|
return b
|
|
|
|
}
|
|
|
|
|
2021-05-08 13:32:18 +02:00
|
|
|
// sign returns the corresponding sign -1.0, 1.0 of the provided number.
|
|
|
|
// if n is not a number, its sign is nan too.
|
|
|
|
[inline]
|
|
|
|
pub fn sign(n f64) f64 {
|
|
|
|
if is_nan(n) {
|
|
|
|
return nan()
|
|
|
|
}
|
|
|
|
return copysign(1.0, n)
|
|
|
|
}
|
|
|
|
|
|
|
|
// signi returns the corresponding sign -1.0, 1.0 of the provided number.
|
|
|
|
[inline]
|
|
|
|
pub fn signi(n f64) int {
|
|
|
|
return int(copysign(1.0, n))
|
|
|
|
}
|
|
|
|
|
2019-07-02 12:50:33 +02:00
|
|
|
// radians convert from radians to degrees.
|
2021-05-08 13:32:18 +02:00
|
|
|
[inline]
|
2019-06-26 17:49:50 +02:00
|
|
|
pub fn radians(degrees f64) f64 {
|
2019-10-12 21:31:05 +02:00
|
|
|
return degrees * (pi / 180.0)
|
2019-06-22 20:20:28 +02:00
|
|
|
}
|
2021-03-13 08:05:02 +01:00
|
|
|
|
|
|
|
// signbit returns a value with the boolean representation of the sign for x
|
2021-05-08 13:32:18 +02:00
|
|
|
[inline]
|
2021-03-13 08:05:02 +01:00
|
|
|
pub fn signbit(x f64) bool {
|
|
|
|
return f64_bits(x) & sign_mask != 0
|
|
|
|
}
|
2021-08-14 07:48:25 +02:00
|
|
|
|
|
|
|
pub fn tolerance(a f64, b f64, tol f64) bool {
|
|
|
|
mut ee := tol
|
|
|
|
// Multiplying by ee here can underflow denormal values to zero.
|
|
|
|
// Check a==b so that at least if a and b are small and identical
|
|
|
|
// we say they match.
|
|
|
|
if a == b {
|
|
|
|
return true
|
|
|
|
}
|
|
|
|
mut d := a - b
|
|
|
|
if d < 0 {
|
|
|
|
d = -d
|
|
|
|
}
|
|
|
|
// note: b is correct (expected) value, a is actual value.
|
|
|
|
// make error tolerance a fraction of b, not a.
|
|
|
|
if b != 0 {
|
|
|
|
ee = ee * b
|
|
|
|
if ee < 0 {
|
|
|
|
ee = -ee
|
|
|
|
}
|
|
|
|
}
|
|
|
|
return d < ee
|
|
|
|
}
|
|
|
|
|
|
|
|
pub fn close(a f64, b f64) bool {
|
|
|
|
return tolerance(a, b, 1e-14)
|
|
|
|
}
|
|
|
|
|
|
|
|
pub fn veryclose(a f64, b f64) bool {
|
|
|
|
return tolerance(a, b, 4e-16)
|
|
|
|
}
|
|
|
|
|
|
|
|
pub fn alike(a f64, b f64) bool {
|
|
|
|
if is_nan(a) && is_nan(b) {
|
|
|
|
return true
|
|
|
|
} else if a == b {
|
|
|
|
return signbit(a) == signbit(b)
|
|
|
|
}
|
|
|
|
return false
|
|
|
|
}
|