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

//
// A rust binding for the GSL library by Guillaume Gomez (guillaume1.gomez@gmail.com)
//

/*!
# Chebyshev Approximations

This chapter describes routines for computing Chebyshev approximations to univariate functions. A
Chebyshev approximation is a truncation of the series f(x) = \sum c_n T_n(x), where the Chebyshev
polynomials T_n(x) = \cos(n \arccos x) provide an orthogonal basis of polynomials on the interval
[-1,1] with the weight function 1 / \sqrt{1-x^2}. The first few Chebyshev polynomials are,
T_0(x) = 1, T_1(x) = x, T_2(x) = 2 x^2 - 1.

For further information see Abramowitz & Stegun, Chapter 22.

## Definitions

The approximation is made over the range [a,b] using order+1 terms, including the coefficient
`c[0]`. The series is computed using the following convention,

f(x) = (c_0 / 2) + \sum_{n=1} c_n T_n(x)

which is needed when accessing the coefficients directly.

## References and Further Reading

The following paper describes the use of Chebyshev series,

R. Broucke, `Ten Subroutines for the Manipulation of Chebyshev Series [C1] (Algorithm 446)`.
Communications of the ACM 16(4), 254–256 (1973)
!*/

use crate::Value;
use ffi::FFI;

ffi_wrapper!(ChebSeries, *mut sys::gsl_cheb_series, gsl_cheb_free);

impl ChebSeries {
    #[doc(alias = "gsl_cheb_alloc")]
    pub fn new(n: usize) -> Option<Self> {
        let tmp = unsafe { sys::gsl_cheb_alloc(n) };

        if tmp.is_null() {
            None
        } else {
            Some(Self::wrap(tmp))
        }
    }

    /// This function computes the Chebyshev approximation cs for the function f over the range
    /// (a,b) to the previously specified order. The computation of the Chebyshev approximation is
    /// an O(n^2) process, and requires n function evaluations.
    #[doc(alias = "gsl_cheb_init")]
    pub fn init<F: Fn(f64) -> f64>(&mut self, f: F, a: f64, b: f64) -> Result<(), Value> {
        let function = wrap_callback!(f, F);

        let ret = unsafe { sys::gsl_cheb_init(self.unwrap_unique(), &function, a, b) };
        result_handler!(ret, ())
    }

    /// This function returns the order of Chebyshev series cs.
    #[doc(alias = "gsl_cheb_order")]
    pub fn order(&self) -> usize {
        unsafe { sys::gsl_cheb_order(self.unwrap_shared()) }
    }

    /// This function returns the size of the Chebyshev coefficient array c[] for the Chebyshev
    /// series cs.
    #[doc(alias = "gsl_cheb_size")]
    pub fn size(&self) -> usize {
        unsafe { sys::gsl_cheb_size(self.unwrap_shared()) }
    }

    /// This function evaluates the Chebyshev series cs at a given point x.
    #[doc(alias = "gsl_cheb_eval")]
    pub fn eval(&self, x: f64) -> f64 {
        unsafe { sys::gsl_cheb_eval(self.unwrap_shared(), x) }
    }

    /// This function computes the Chebyshev series cs at a given point x, estimating both the
    /// series result and its absolute error abserr. The error estimate is made from the first
    /// neglected term in the series.
    ///
    /// Returns `(result, abs_err)`.
    #[doc(alias = "gsl_cheb_eval_err")]
    pub fn eval_err(&self, x: f64) -> Result<(f64, f64), Value> {
        let mut result = 0.;
        let mut abs_err = 0.;

        let ret =
            unsafe { sys::gsl_cheb_eval_err(self.unwrap_shared(), x, &mut result, &mut abs_err) };
        result_handler!(ret, (result, abs_err))
    }

    /// This function evaluates the Chebyshev series cs at a given point x, to (at most) the given
    /// order order.
    #[doc(alias = "gsl_cheb_eval_n")]
    pub fn eval_n(&self, order: usize, x: f64) -> f64 {
        unsafe { sys::gsl_cheb_eval_n(self.unwrap_shared(), order, x) }
    }

    /// This function evaluates a Chebyshev series cs at a given point x, estimating both the series
    /// result and its absolute error abserr, to (at most) the given order order. The error estimate
    /// is made from the first neglected term in the series.
    ///
    /// Returns `(result, abs_err)`.
    #[doc(alias = "gsl_cheb_eval_n_err")]
    pub fn eval_n_err(&self, order: usize, x: f64) -> Result<(f64, f64), Value> {
        let mut result = 0.;
        let mut abs_err = 0.;

        let ret = unsafe {
            sys::gsl_cheb_eval_n_err(self.unwrap_shared(), order, x, &mut result, &mut abs_err)
        };
        result_handler!(ret, (result, abs_err))
    }

    /// This function computes the derivative of the series cs, storing the derivative coefficients
    /// in the previously allocated deriv. The two series cs and deriv must have been allocated with
    /// the same order.
    #[doc(alias = "gsl_cheb_calc_deriv")]
    pub fn calc_deriv(&self, deriv: &mut ChebSeries) -> Result<(), Value> {
        let ret = unsafe { sys::gsl_cheb_calc_deriv(deriv.unwrap_unique(), self.unwrap_shared()) };
        result_handler!(ret, ())
    }

    /// This function computes the integral of the series cs, storing the integral coefficients in
    /// the previously allocated integ. The two series cs and integ must have been allocated with
    /// the same order. The lower limit of the integration is taken to be the left hand end of the
    /// range a.
    #[doc(alias = "gsl_cheb_calc_integ")]
    pub fn calc_integ(&self, integ: &mut ChebSeries) -> Result<(), Value> {
        let ret = unsafe { sys::gsl_cheb_calc_integ(integ.unwrap_unique(), self.unwrap_shared()) };
        result_handler!(ret, ())
    }
}