Commit fef116697a5f39740ccdb8b7ad7442e77dd19ad2

Authored by Georg Hopp
1 parent d70c0220

add trigon functions sin and cos as static tables of Fractional

  1 +[package]
  2 +name = "fractional"
  3 +version = "0.1.0"
  4 +authors = ["Georg Hopp <georg@steffers.org>"]
  5 +edition = "2018"
  6 +
  7 +[dependencies]
  8 +lazy_static = "1.4.0"
... ...
  1 +//
  2 +// Some code to support fractional numbers for full precision rational number
  3 +// calculations.
  4 +// TODO
  5 +// - maybe this could be build as a generic for all integral numbers.
  6 +// (Question, how can I assure that it is build from integral numbers?
  7 +//
  8 +// Georg Hopp <georg@steffers.org>
  9 +//
  10 +// Copyright © 2019 Georg Hopp
  11 +//
  12 +// This program is free software: you can redistribute it and/or modify
  13 +// it under the terms of the GNU General Public License as published by
  14 +// the Free Software Foundation, either version 3 of the License, or
  15 +// (at your option) any later version.
  16 +//
  17 +// This program is distributed in the hope that it will be useful,
  18 +// but WITHOUT ANY WARRANTY; without even the implied warranty of
  19 +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  20 +// GNU General Public License for more details.
  21 +//
  22 +// You should have received a copy of the GNU General Public License
  23 +// along with this program. If not, see <http://www.gnu.org/licenses/>.
  24 +//
  25 +use std::cmp::Ordering;
  26 +use std::ops::{Add,Sub,Neg,Mul,Div};
  27 +use std::fmt;
  28 +use std::convert::{TryFrom, TryInto};
  29 +use std::num::TryFromIntError;
  30 +
  31 +#[derive(Debug, Eq, Clone, Copy)]
  32 +pub struct Fractional (pub i64, pub i64);
  33 +
  34 +#[inline]
  35 +fn hcf(x :i64, y :i64) -> i64 {
  36 + match y {
  37 + 0 => x,
  38 + _ => hcf(y, x % y),
  39 + }
  40 +}
  41 +
  42 +impl Fractional {
  43 + #[inline]
  44 + pub fn gcd(self, other: Self) -> i64 {
  45 + let Fractional(_, d1) = self;
  46 + let Fractional(_, d2) = other;
  47 + (d1 * d2) / hcf(d1, d2)
  48 + }
  49 +
  50 + #[inline]
  51 + pub fn reduce(self) -> Self {
  52 + let Fractional(n, d) = self;
  53 + Self(n / hcf(n, d), d / hcf(n, d))
  54 + }
  55 +
  56 + #[inline]
  57 + pub fn numerator(self) -> i64 {
  58 + self.0
  59 + }
  60 +
  61 + #[inline]
  62 + pub fn denominator(self) -> i64 {
  63 + self.1
  64 + }
  65 +}
  66 +
  67 +impl fmt::Display for Fractional {
  68 + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
  69 + write!(f, "({}/{})", self.0, self.1)
  70 + }
  71 +}
  72 +
  73 +impl From<i64> for Fractional {
  74 + fn from(x: i64) -> Self {
  75 + Self(x, 1)
  76 + }
  77 +}
  78 +
  79 +impl TryFrom<usize> for Fractional {
  80 + type Error = &'static str;
  81 +
  82 + fn try_from(x: usize) -> Result<Self, Self::Error> {
  83 + let v = i64::try_from(x);
  84 + match v {
  85 + Err(_) => Err("Conversion from usize to i32 failed"),
  86 + Ok(_v) => Ok(Self(_v, 1)),
  87 + }
  88 + }
  89 +}
  90 +
  91 +pub fn from_vector(xs: &Vec<i64>) -> Vec<Fractional> {
  92 + xs.iter().map(|x| Fractional(*x, 1)).collect()
  93 +}
  94 +
  95 +impl TryInto<f64> for Fractional {
  96 + type Error = TryFromIntError;
  97 +
  98 + fn try_into(self) -> Result<f64, Self::Error> {
  99 + let n :i32 = self.0.try_into()?;
  100 + let d :i32 = self.1.try_into()?;
  101 + Ok(f64::from(n) / f64::from(d))
  102 + }
  103 +}
  104 +
  105 +impl TryInto<(i32, i32)> for Fractional {
  106 + type Error = TryFromIntError;
  107 +
  108 + fn try_into(self) -> Result<(i32, i32), Self::Error> {
  109 + let a :i32 = (self.0 / self.1).try_into()?;
  110 + let b :i32 = (self.0 % self.1).try_into()?;
  111 + Ok((a, b))
  112 + }
  113 +}
  114 +
  115 +impl PartialEq for Fractional {
  116 + fn eq(&self, other: &Self) -> bool {
  117 + let Fractional(n1, d1) = self;
  118 + let Fractional(n2, d2) = other;
  119 + n1 * (self.gcd(*other) / d1) == n2 * (self.gcd(*other) / d2)
  120 + }
  121 +}
  122 +
  123 +impl PartialOrd for Fractional {
  124 + fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
  125 + Some(self.cmp(other))
  126 + }
  127 +}
  128 +
  129 +impl Ord for Fractional {
  130 + fn cmp(&self, other: &Self) -> Ordering {
  131 + let Fractional(n1, d1) = self;
  132 + let Fractional(n2, d2) = other;
  133 + let x = n1 * (self.gcd(*other) / d1);
  134 + let y = n2 * (self.gcd(*other) / d2);
  135 + x.cmp(&y)
  136 + }
  137 +}
  138 +
  139 +impl Add for Fractional {
  140 + type Output = Self;
  141 +
  142 + fn add(self, other: Self) -> Self {
  143 + let Fractional(n1, d1) = self;
  144 + let Fractional(n2, d2) = other;
  145 + let n = n1 * (self.gcd(other) / d1) + n2 * (self.gcd(other) / d2);
  146 + Self(n, self.gcd(other)).reduce()
  147 + }
  148 +}
  149 +
  150 +impl Sub for Fractional {
  151 + type Output = Self;
  152 +
  153 + fn sub(self, other: Self) -> Self {
  154 + let Fractional(n1, d1) = self;
  155 + let Fractional(n2, d2) = other;
  156 + let n = n1 * (self.gcd(other) / d1) - n2 * (self.gcd(other) / d2);
  157 + Self(n, self.gcd(other)).reduce()
  158 + }
  159 +}
  160 +
  161 +impl Neg for Fractional {
  162 + type Output = Self;
  163 +
  164 + fn neg(self) -> Self {
  165 + let Fractional(n, d) = self;
  166 + Self(-n, d).reduce()
  167 + }
  168 +}
  169 +
  170 +impl Mul for Fractional {
  171 + type Output = Self;
  172 +
  173 + fn mul(self, other :Self) -> Self {
  174 + let Fractional(n1, d1) = self;
  175 + let Fractional(n2, d2) = other;
  176 + Self(n1 * n2, d1 * d2).reduce()
  177 + }
  178 +}
  179 +
  180 +impl Div for Fractional {
  181 + type Output = Self;
  182 +
  183 + fn div(self, other: Self) -> Self {
  184 + let Fractional(n, d) = other;
  185 + self * Fractional(d, n)
  186 + }
  187 +}
  188 +
  189 + /* some stuff that could be tested...
  190 + let x = Fractional(1, 3);
  191 + let y = Fractional(1, 6);
  192 +
  193 + println!(
  194 + "Greatest common denominator of {} and {}: {}", x, y, x.gcd(y));
  195 + println!("Numerator of {}: {}", x, x.numerator());
  196 + println!("Denominator of {}: {}", x, x.denominator());
  197 + assert_eq!(Fractional(1, 3), Fractional(2, 6));
  198 + assert_eq!(Fractional(1, 3), Fractional(1, 3));
  199 + assert_eq!(y < x, true);
  200 + assert_eq!(y > x, false);
  201 + assert_eq!(x == y, false);
  202 + assert_eq!(x == x, true);
  203 + assert_eq!(x + y, Fractional(1, 2));
  204 + println!("{} + {} = {}", x, y, x + y);
  205 + assert_eq!(x - y, Fractional(1, 6));
  206 + println!("{} - {} = {}", x, y, x - y);
  207 + assert_eq!(y - x, Fractional(-1, 6));
  208 + println!("{} - {} = {}", y, x, y - x);
  209 + assert_eq!(-x, Fractional(-1, 3));
  210 + println!("-{} = {}", x, -x);
  211 + assert_eq!(x * y, Fractional(1, 18));
  212 + println!("{} * {} = {}", x, y, x * y);
  213 + assert_eq!(x / y, Fractional(2, 1));
  214 + println!("{} / {} = {}", x, y, x / y);
  215 + assert_eq!(y / x, Fractional(1, 2));
  216 + println!("{} / {} = {}", y, x, y / x);
  217 +
  218 + println!("Fractional from 3: {}", Fractional::from(3));
  219 + let z :f64 = Fractional::into(x);
  220 + println!("Floating point of {}: {}", x, z);
  221 + let (d, r) = Fractional::into(x);
  222 + println!("(div, rest) of {}: ({}, {})", x, d, r);
  223 + */
  224 +
... ...
  1 +//
  2 +// Lib for fractional calculations.
  3 +//
  4 +// Georg Hopp <georg@steffers.org>
  5 +//
  6 +// Copyright © 2019 Georg Hopp
  7 +//
  8 +// This program is free software: you can redistribute it and/or modify
  9 +// it under the terms of the GNU General Public License as published by
  10 +// the Free Software Foundation, either version 3 of the License, or
  11 +// (at your option) any later version.
  12 +//
  13 +// This program is distributed in the hope that it will be useful,
  14 +// but WITHOUT ANY WARRANTY; without even the implied warranty of
  15 +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  16 +// GNU General Public License for more details.
  17 +//
  18 +// You should have received a copy of the GNU General Public License
  19 +// along with this program. If not, see <http://www.gnu.org/licenses/>.
  20 +//
  21 +extern crate lazy_static;
  22 +
  23 +pub mod fractional;
  24 +pub mod trigonometry;
  25 +
  26 +use fractional::{Fractional};
... ...
  1 +//
  2 +// Test our fractional crate / module...
  3 +//
  4 +// Georg Hopp <georg@steffers.org>
  5 +//
  6 +// Copyright © 2019 Georg Hopp
  7 +//
  8 +// This program is free software: you can redistribute it and/or modify
  9 +// it under the terms of the GNU General Public License as published by
  10 +// the Free Software Foundation, either version 3 of the License, or
  11 +// (at your option) any later version.
  12 +//
  13 +// This program is distributed in the hope that it will be useful,
  14 +// but WITHOUT ANY WARRANTY; without even the implied warranty of
  15 +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  16 +// GNU General Public License for more details.
  17 +//
  18 +// You should have received a copy of the GNU General Public License
  19 +// along with this program. If not, see <http://www.gnu.org/licenses/>.
  20 +//
  21 +use std::convert::{TryFrom, TryInto};
  22 +use std::num::TryFromIntError;
  23 +use std::f64::consts::PI as FPI;
  24 +
  25 +use fractional::fractional::{Fractional, from_vector};
  26 +use fractional::trigonometry::{sin, cos, PI};
  27 +
  28 +struct Vector {
  29 + x: Fractional,
  30 + y: Fractional,
  31 + z: Fractional,
  32 +}
  33 +
  34 +fn mean(v: &Vec<i64>) -> Result<Fractional, TryFromIntError> {
  35 + let r = v.iter().fold(0, |acc, x| acc + x);
  36 + let l = i64::try_from(v.len())?;
  37 + Ok(Fractional(r, l))
  38 +}
  39 +
  40 +fn main() {
  41 + let a = vec![3, 6, 1, 9];
  42 + let b = from_vector(&a);
  43 + let c = mean(&a).unwrap(); // This might fail if the len of the
  44 + // vector (usize) does not fit into i32.
  45 + let d :f64 = c.try_into().unwrap();
  46 + let e :f64 = Fractional::try_into(c).unwrap();
  47 +
  48 + println!(" [i32] : {:?}" , a);
  49 + println!(" [Fractional] : {:?}" , b);
  50 + println!(" mean of [i32] : {}" , c);
  51 + println!(" as f64 : {}" , d);
  52 + println!(" and as f64 : {}" , e);
  53 + println!(" again as f64 : {}" , TryInto::<f64>::try_into(c).unwrap());
  54 + println!(" Rust π : {}" , FPI);
  55 + println!(" π : {} {}" , TryInto::<f64>::try_into(PI).unwrap(), PI);
  56 + println!(" π as tuple : {:?}" , TryInto::<(i32, i32)>::try_into(PI).unwrap());
  57 + println!(" Rust π² : {}" , FPI * FPI);
  58 + println!(" π² : {} {}" , TryInto::<f64>::try_into(PI * PI).unwrap(), PI * PI);
  59 + println!(" sin 9 : {}" , sin(9));
  60 + println!(" sin 9 : {}" , TryInto::<f64>::try_into(sin(9)).unwrap());
  61 + println!(" Rust sin 9 : {}" , f64::sin(9.0 * FPI / 180.0));
  62 + println!(" sin 17 : {}" , sin(17));
  63 + println!(" sin 17 : {}" , TryInto::<f64>::try_into(sin(17)).unwrap());
  64 + println!(" Rust sin 17 : {}" , f64::sin(17.0 * FPI / 180.0));
  65 + println!(" sin 31 : {}" , sin(31));
  66 + println!(" sin 31 : {}" , TryInto::<f64>::try_into(sin(31)).unwrap());
  67 + println!(" Rust sin 31 : {}" , f64::sin(31.0 * FPI / 180.0));
  68 + println!(" sin 45 : {}" , sin(45));
  69 + println!(" sin 45 : {}" , TryInto::<f64>::try_into(sin(45)).unwrap());
  70 + println!(" Rust sin 45 : {}" , f64::sin(45.0 * FPI / 180.0));
  71 + println!(" sin 73 : {}" , sin(73));
  72 + println!(" sin 73 : {}" , TryInto::<f64>::try_into(sin(73)).unwrap());
  73 + println!(" Rust sin 73 : {}" , f64::sin(73.0 * FPI / 180.0));
  74 + println!(" sin 123 : {}" , sin(123));
  75 + println!(" sin 123 : {}" , TryInto::<f64>::try_into(sin(123)).unwrap());
  76 + println!(" Rust sin 123 : {}" , f64::sin(123.0 * FPI / 180.0));
  77 + println!(" sin 213 : {}" , sin(213));
  78 + println!(" sin 213 : {}" , TryInto::<f64>::try_into(sin(213)).unwrap());
  79 + println!(" Rust sin 213 : {}" , f64::sin(213.0 * FPI / 180.0));
  80 + println!(" sin 312 : {}" , sin(312));
  81 + println!(" sin 312 : {}" , TryInto::<f64>::try_into(sin(312)).unwrap());
  82 + println!(" Rust sin 312 : {}" , f64::sin(312.0 * FPI / 180.0));
  83 + println!(" sin 876 : {}" , sin(876));
  84 + println!(" sin 876 : {}" , TryInto::<f64>::try_into(sin(876)).unwrap());
  85 + println!(" Rust sin 876 : {}" , f64::sin(876.0 * FPI / 180.0));
  86 + println!(" cos 9 : {}" , cos(9));
  87 + println!(" cos 9 : {}" , TryInto::<f64>::try_into(cos(9)).unwrap());
  88 + println!(" Rust cos 9 : {}" , f64::cos(9.0 * FPI / 180.0));
  89 + println!(" cos 17 : {}" , cos(17));
  90 + println!(" cos 17 : {}" , TryInto::<f64>::try_into(cos(17)).unwrap());
  91 + println!(" Rust cos 17 : {}" , f64::cos(17.0 * FPI / 180.0));
  92 + println!(" cos 31 : {}" , cos(31));
  93 + println!(" cos 31 : {}" , TryInto::<f64>::try_into(cos(31)).unwrap());
  94 + println!(" Rust cos 31 : {}" , f64::cos(31.0 * FPI / 180.0));
  95 + println!(" cos 45 : {}" , cos(45));
  96 + println!(" cos 45 : {}" , TryInto::<f64>::try_into(cos(45)).unwrap());
  97 + println!(" Rust cos 45 : {}" , f64::cos(45.0 * FPI / 180.0));
  98 + println!(" cos 73 : {}" , cos(73));
  99 + println!(" cos 73 : {}" , TryInto::<f64>::try_into(cos(73)).unwrap());
  100 + println!(" Rust cos 73 : {}" , f64::cos(73.0 * FPI / 180.0));
  101 + println!(" cos 123 : {}" , cos(123));
  102 + println!(" cos 123 : {}" , TryInto::<f64>::try_into(cos(123)).unwrap());
  103 + println!(" Rust cos 123 : {}" , f64::cos(123.0 * FPI / 180.0));
  104 + println!(" cos 213 : {}" , cos(213));
  105 + println!(" cos 213 : {}" , TryInto::<f64>::try_into(cos(213)).unwrap());
  106 + println!(" Rust cos 213 : {}" , f64::cos(213.0 * FPI / 180.0));
  107 + println!(" cos 312 : {}" , cos(312));
  108 + println!(" cos 312 : {}" , TryInto::<f64>::try_into(cos(312)).unwrap());
  109 + println!(" Rust cos 312 : {}" , f64::cos(312.0 * FPI / 180.0));
  110 + println!(" cos 876 : {}" , cos(876));
  111 + println!(" cos 876 : {}" , TryInto::<f64>::try_into(cos(876)).unwrap());
  112 + println!(" Rust cos 876 : {}" , f64::cos(876.0 * FPI / 180.0));
  113 +}
... ...
  1 +//
  2 +// Test our fractional crate / module...
  3 +//
  4 +// Georg Hopp <georg@steffers.org>
  5 +//
  6 +// Copyright © 2019 Georg Hopp
  7 +//
  8 +// This program is free software: you can redistribute it and/or modify
  9 +// it under the terms of the GNU General Public License as published by
  10 +// the Free Software Foundation, either version 3 of the License, or
  11 +// (at your option) any later version.
  12 +//
  13 +// This program is distributed in the hope that it will be useful,
  14 +// but WITHOUT ANY WARRANTY; without even the implied warranty of
  15 +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  16 +// GNU General Public License for more details.
  17 +//
  18 +// You should have received a copy of the GNU General Public License
  19 +// along with this program. If not, see <http://www.gnu.org/licenses/>.
  20 +//
  21 +use crate::{Fractional};
  22 +
  23 +pub const PI :Fractional = Fractional(355, 113); // This is a really close
  24 + // fractional approximation
  25 + // for pi
  26 +
  27 +const PRECISION :i64 = 100000;
  28 +
  29 +#[inline]
  30 +pub fn rad(d: u32) -> f64 {
  31 + use std::f64::consts::PI;
  32 + d as f64 * PI / 180.0
  33 +}
  34 +
  35 +pub fn sin(d: i32) -> Fractional {
  36 + // hold sin Fractionals from 0 to 89 ...
  37 + lazy_static::lazy_static! {
  38 + static ref SINTAB :Vec<Fractional> =
  39 + (0..90).map(|x| _sin(x)).collect();
  40 + }
  41 +
  42 + // fractional sin from f64 sin. (From 1° to 89°)
  43 + fn _sin(d: u32) -> Fractional {
  44 + match d {
  45 + 0 => Fractional(0, 1),
  46 + _ => {
  47 + // This is undefined behaviour for very large f64, but our f64
  48 + // is always between 0.0 and 10000.0 which should be fine.
  49 + let s = (f64::sin(rad(d)) * PRECISION as f64).round() as i64;
  50 + Fractional(s, PRECISION).reduce()
  51 + }
  52 + }
  53 + }
  54 +
  55 + match d {
  56 + 90 => Fractional(1, 1),
  57 + 180 => SINTAB[0],
  58 + 270 => -Fractional(1, 1),
  59 + 1..=89 => SINTAB[d as usize],
  60 + 91..=179 => SINTAB[180 - d as usize],
  61 + 181..=269 => -SINTAB[d as usize - 180],
  62 + 271..=359 => -SINTAB[360 - d as usize],
  63 + _ => sin(d % 360),
  64 + }
  65 +}
  66 +
  67 +pub fn cos(d: i32) -> Fractional {
  68 + lazy_static::lazy_static! {
  69 + static ref COSTAB :Vec<Fractional> =
  70 + (0..90).map(|x| _cos(x)).collect();
  71 + }
  72 +
  73 + fn _cos(d: u32) -> Fractional {
  74 + match d {
  75 + 0 => Fractional(1, 1),
  76 + _ => {
  77 + let s = (f64::cos(rad(d)) * PRECISION as f64).round() as i64;
  78 + Fractional(s, PRECISION).reduce()
  79 + }
  80 + }
  81 + }
  82 +
  83 + match d {
  84 + 90 | 270 => Fractional(0, 1),
  85 + 180 => -COSTAB[0],
  86 + 1..=89 => COSTAB[d as usize],
  87 + 91..=179 => -COSTAB[180 - d as usize],
  88 + 181..=269 => -COSTAB[d as usize - 180],
  89 + 271..=359 => COSTAB[360 - d as usize],
  90 + _ => cos(d % 360),
  91 + }
  92 +}
... ...
Please register or login to post a comment