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
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152

use crate::{Elementary, Error, Func, Function, Round};

/// types that implement the [Integrate](crate::Integrate) trait can safely be integrated within
/// the domain ℝ.
pub trait Integrate {
    /// Method will return an instance of [Integral](crate::Integral) that can be taylored to
    /// specific usecases.
    ///
    /// Example:
    /// ```rust
    /// let function = Function::from("sin(x)");
    ///
    /// let mut integral = function.integrate();

    /// // specify bounds and precision for the integral
    /// integral
    ///     .set_lower_bound(0.)
    ///     .set_upper_bound(PI / 2.)
    ///     .set_precision(20000);
    ///
    /// // evaluate the integral
    /// let value = integral.evaluate().unwrap();
    /// // note that the value of the evaluated integral must be unwrapped if using the `integrate()`
    /// // method because the method cannot guarantee that bounds have been set at the point of
    /// // evaluating. The evaluate_integral() method which is implemented for any instance with the
    /// // Integrate trait is safer and is guaranteed to yield a valid result.
    ///
    /// // round the value to 10 decimal places
    /// value.round_to(10);
    ///
    /// assert_eq!(value, 1.);
    /// ```
    /// Also note that if the precision is not specified, it will default to 1000.
    fn integrate(&self) -> Integral;

    /// Method will return the value of the definite integral of the function evaluated from the
    /// provided lower and upper bound.
    ///
    /// Example:
    /// ```
    /// let function = Function::from("cos(x)");
    ///
    /// // the evaluate_integral() method will automatically round the result to five decimal points.
    /// // This is because higher precision cannot be guaranteed with using the standard precision set
    /// // for the method. Provided that the function is defined for all values between the lower and
    /// // upper bounds, the method will always return a valid result.
    /// let value = function.evaluate_integral(0., PI);
    ///
    /// assert_eq!(value, 0.);
    /// ```
    fn evaluate_integral(&self, lower_bound: f64, upper_bound: f64) -> f64;
}

const STANDARD_PRECISION: usize = 1000;

/// See [Integrate documentation](crate::Integrate) for usage and examples
pub struct Integral {
    function: Func,
    lower_bound: Option<f64>,
    upper_bound: Option<f64>,
    precision: usize,
}

impl Integral {
    pub fn vacant(function: Func) -> Self {
        Self {
            function,
            lower_bound: None,
            upper_bound: None,
            precision: STANDARD_PRECISION,
        }
    }

    pub fn set_lower_bound(&mut self, lower_bound: f64) -> &mut Self {
        self.lower_bound = Some(lower_bound);
        self
    }

    pub fn set_upper_bound(&mut self, upper_bound: f64) -> &mut Self {
        self.upper_bound = Some(upper_bound);
        self
    }

    pub fn set_precision(&mut self, precision: usize) -> &mut Self {
        self.precision = precision;
        self
    }

    pub fn evaluate(&self) -> Result<f64, Error> {
        if let (Some(lower_bound), Some(upper_bound)) = (self.lower_bound, self.upper_bound) {
            Ok(simpsons_rule(
                &self.function,
                lower_bound,
                upper_bound,
                self.precision,
            ))
        } else {
            Err(Error::InternalError(String::from(
                "Bounds of integration must be set in order to evaluate the integral",
            )))
        }
    }
}

/// See [Integrate](crate::Integrate) for usage and examples.
impl Integrate for Elementary {
    fn integrate(&self) -> Integral {
        Integral::vacant(self.clone().call())
    }
    /// Evaluating the integral gives a value of the integral with eight decimal places
    fn evaluate_integral(&self, lower_bound: f64, upper_bound: f64) -> f64 {
        unsafe {
            let mut value = self
                .integrate()
                .set_lower_bound(lower_bound)
                .set_upper_bound(upper_bound)
                .evaluate()
                .unwrap_unchecked(); // this unwrap will never fail because the upper and lower bounds
                                     // will always be set
            value.with_significant_figures(5)
        }
    }
}

/// See [Integrate](crate::Integrate) for usage and examples.
impl Integrate for Function {
    fn integrate(&self) -> Integral {
        self.elementary().integrate()
    }
    fn evaluate_integral(&self, lower_bound: f64, upper_bound: f64) -> f64 {
        self.elementary()
            .evaluate_integral(lower_bound, upper_bound)
    }
}

fn simpsons_rule(funciton: &Func, lower_bound: f64, upper_bound: f64, precision: usize) -> f64 {
    // note that n must be an even number for Simpson's rule to work
    let n = precision * 2;
    let dx = (upper_bound - lower_bound) / n as f64;
    let mut sum: f64 = (1..n)
        .map(|x| {
            if x % 2 == 0 {
                2. * funciton(lower_bound + x as f64 * dx)
            } else {
                4. * funciton(lower_bound + x as f64 * dx)
            }
        })
        .sum();
    sum += funciton(lower_bound) + funciton(upper_bound);

    sum * dx / 3.
}