Add tangens and improve fractional rep

fractional/src/trigonometry.rs

@@ -1,10 +1,12 @@
 //
 // Some trigonometic functions with Fractions results.
-// Currently only sin and cos are implemented. As I was unable
-// to find a really good integral approximation for them I
-// implement them as tables which are predefined using the
-// floating point function f64::sin and then transformed into
-// a fraction of a given PRECISION.
+// 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>
 //
@@ -23,75 +25,94 @@
 // 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
 
-const PRECISION :i64 = 100000;
+// 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;
 
-#[inline]
-pub fn rad(d: u32) -> f64 {
-    use std::f64::consts::PI;
-    d as f64 * PI / 180.0
-}
-
-pub fn sin(d: i32) -> Fractional {
-    // hold sin Fractionals from 0 to 89 ...
-    lazy_static::lazy_static! {
-        static ref SINTAB :Vec<Fractional> =
-            (0..90).map(|x| _sin(x)).collect();
+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),
     }
+}
 
-    // fractional sin from f64 sin. (From 1° to 89°)
-    fn _sin(d: u32) -> Fractional {
-        match d {
-            0 => Fractional(0, 1),
-            _ => {
-                // This is undefined behaviour for very large f64, but our f64
-                // is always between 0.0 and 10000.0 which should be fine.
-                let s = (f64::sin(rad(d)) * PRECISION as f64).round() as i64;
-                Fractional(s, PRECISION).reduce()
-            }
-        }
+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 {
-        90        => Fractional(1, 1),
-        180       => SINTAB[0],
-        270       => -Fractional(1, 1),
-        1..=89    => SINTAB[d as usize],
-        91..=179  => SINTAB[180 - d as usize],
-        181..=269 => -SINTAB[d as usize - 180],
-        271..=359 => -SINTAB[360 - d as usize],
-        _         => sin(d % 360),
+        0  ..=179 => TANTAB[d as usize],
+        180..=359 => TANTAB[d as usize - 180],
+        _         => tan(d % 360),
     }
 }
 
-pub fn cos(d: i32) -> Fractional {
-    lazy_static::lazy_static! {
-        static ref COSTAB :Vec<Fractional> =
-            (0..90).map(|x| _cos(x)).collect();
-    }
+// 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();
+}
 
-    fn _cos(d: u32) -> Fractional {
-        match d {
-            0 => Fractional(1, 1),
-            _ => {
-                let s = (f64::cos(rad(d)) * PRECISION as f64).round() as i64;
-                Fractional(s, PRECISION).reduce()
-            }
-        }
+// 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 {
-        90 | 270  => Fractional(0, 1),
-        180       => -COSTAB[0],
-        1..=89    => COSTAB[d as usize],
-        91..=179  => -COSTAB[180 - d as usize],
-        181..=269 => -COSTAB[d as usize - 180],
-        271..=359 => COSTAB[360 - d as usize],
-        _         => cos(d % 360),
+        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),
     }
 }