trigonometry.rs
4.1 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
//
// Some trigonometic functions with Fractions results.
// Currently only sin, cos and tan are implemented.
// As I was unable to find a really good integral approximation for them I
// implement them as a table which is predefined using the floating point
// function f64::sin and then transformed into a fraction of a given
// PRECISION.
// These approximations are quite good and for a few edge cases
// even better than the floating point implementations.
//
// Georg Hopp <georg@steffers.org>
//
// Copyright © 2019 Georg Hopp
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
use std::cmp::Ordering;
use crate::{Fractional};
pub const PI :Fractional = Fractional(355, 113); // This is a really close
// fractional approximation
// for pi
// Try to keep precision as high as possible while having a denominator
// as small as possible.
// The values are taken by triel and error.
const PRECISION :i64 = 1000000;
const MAX_DENOMINATOR :i64 = 7000;
pub fn sin(d :i32) -> Fractional {
match d {
0 ..=90 => SINTAB[d as usize],
91 ..=180 => SINTAB[180 - d as usize],
181..=270 => -SINTAB[d as usize - 180],
271..=359 => -SINTAB[360 - d as usize],
_ => sin(d % 360),
}
}
pub fn cos(d :i32) -> Fractional {
match d {
0 ..=90 => SINTAB[90 - d as usize],
91 ..=180 => -SINTAB[90 - (180 - d as usize)],
181..=270 => -SINTAB[90 - (d as usize - 180)],
271..=359 => SINTAB[90 - (360 - d as usize)],
_ => cos(d % 360),
}
}
pub fn tan(d :i32) -> Fractional {
match d {
0 ..=179 => TANTAB[d as usize],
180..=359 => TANTAB[d as usize - 180],
_ => tan(d % 360),
}
}
// hold sin Fractionals from 0 to 89 ...
// luckily with a bit of index tweeking this can also be used for
// cosine values.
lazy_static::lazy_static! {
static ref SINTAB :Vec<Fractional> =
(0..=90).map(|x| _sin(x)).collect();
}
// This table exists only because the sin(α) / cos(α) method
// yields very large unreducable denominators in a lot of cases.
lazy_static::lazy_static! {
static ref TANTAB :Vec<Fractional> =
(0..180).map(|x| _tan(x)).collect();
}
// fractional sin from f64 sin. (From 0° to 90°)
fn _sin(d: u32) -> Fractional {
match d {
0 => Fractional(0, 1),
90 => Fractional(1, 1),
_ => reduce(d, PRECISION, &f64::sin),
}
}
// fractional tan from f64 tan. (From 0° to 179°)
fn _tan(d: u32) -> Fractional {
match d {
0 => Fractional(0, 1),
45 => Fractional(1, 1),
90 => Fractional(1, 0), // although they are both inf and -inf.
135 => -Fractional(1, 1),
_ => reduce(d, PRECISION, &f64::tan),
}
}
// search for a fraction with a denominator less than 10000 that
// provides the minimal precision criteria.
// !! With f = &f64::tan and d close to the inf boundarys of tan
// we get very large numerators because the numerator becomes a
// multiple of the denominator.
fn reduce(d :u32, p :i64, f :&dyn Fn(f64) -> f64) -> Fractional {
// This is undefined behaviour for very large f64, but our f64
// is always between 0.0 and 10000000.0 which should be fine.
let s = (f(f64::to_radians(d as f64)) * p as f64).round() as i64;
let Fractional(n, dn) = Fractional(s, p).reduce();
match dn.abs().cmp(&MAX_DENOMINATOR) {
Ordering::Less => Fractional(n, dn),
_ => reduce(d, p + 1, f),
}
}