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
153
154
155
156
157
158
159
160
161
162

//
// 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 ffi;
use enums;
use std::f64::consts::PI;
use c_vec::CSlice;

pub struct ChebSeries {
    c: *mut ffi::gsl_cheb_series,
    data: CSlice<f64>
}

impl ChebSeries {
    pub fn new(n: usize) -> Option<ChebSeries> {
        let tmp = unsafe { ffi::gsl_cheb_alloc(n) };

        if tmp.is_null() {
            None
        } else {
            unsafe {
                Some(ChebSeries {
                    c: tmp,
                    data: CSlice::new((*tmp).c, (*tmp).order as usize + 1)
                })
            }
        }
    }

    /// 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.
    pub fn init<T>(&mut self, func: ::function<T>, a: f64, b: f64, param: &mut T) -> enums::Value {
        if a >= b {
            rgsl_error!("null function interval [a,b]", ::Value::Domain);
            ::Value::Failure
        } else {
            unsafe {
                (*self.c).a = a;
                (*self.c).b = b;

                let bma = 0.5 * (b - a);
                let bpa = 0.5 * (b + a);
                let fac = 2.0 / ((*self.c).order as f64 + 1.0);

                let mut tmp_vec = CSlice::new((*self.c).f, (*self.c).order_sp as usize + 1);
                for k in 0..((*self.c).order + 1) {
                    let y = (PI * (k as f64 + 0.5) / ((*self.c).order as f64 + 1f64)).cos();
                    tmp_vec.as_mut()[k as usize] = func(y * bma + bpa, param);
                }

                for j in 0..((*self.c).order + 1) {
                    let mut sum = 0f64;

                    for k in 0..((*self.c).order + 1) {
                        sum += tmp_vec.as_ref()[k as usize] * (PI * j as f64 * (k as f64 + 0.5) / ((*self.c).order as f64 + 1f64)).cos();
                    }
                    self.data.as_mut()[j as usize] = fac * sum;
                }
            }
            ::Value::Success
        }
    }

    /// This function returns the order of Chebyshev series cs.
    pub fn order(&self) -> usize {
        unsafe { ffi::gsl_cheb_order(self.c) }
    }

    /// This function returns the size of the Chebyshev coefficient array c[] for the Chebyshev series cs.
    pub fn size(&self) -> usize {
        unsafe { ffi::gsl_cheb_size(self.c) }
    }

    /// This function returns a pointer to the coefficient array c[] location in memory for the Chebyshev series cs.
    pub fn as_slice<'r>(&'r self) -> &'r [f64] {
        self.data.as_ref()
    }

    /// This function returns a pointer to the coefficient array c[] location in memory for the Chebyshev series cs.
    pub fn as_mut_slice<'r>(&'r mut self) -> &'r mut [f64] {
        self.data.as_mut()
    }

    /// This function evaluates the Chebyshev series cs at a given point x.
    pub fn eval(&self, x: f64) -> f64 {
        unsafe { ffi::gsl_cheb_eval(self.c, 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.
    pub fn eval_err(&self, x: f64, result: &mut f64, abs_err: &mut f64) -> enums::Value {
        unsafe { ffi::gsl_cheb_eval_err(self.c, x, result, abs_err) }
    }

    /// This function evaluates the Chebyshev series cs at a given point x, to (at most) the given order order.
    pub fn eval_n(&self, order: usize, x: f64) -> f64 {
        unsafe { ffi::gsl_cheb_eval_n(self.c, 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.
    pub fn eval_n_err(&self, order: usize, x: f64, result: &mut f64, abs_err: &mut f64) -> enums::Value {
        unsafe { ffi::gsl_cheb_eval_n_err(self.c, order, x, 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.
    pub fn calc_deriv(&self, deriv: &ChebSeries) -> enums::Value {
        unsafe { ffi::gsl_cheb_calc_deriv(deriv.c, self.c) }
    }

    /// 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.
    pub fn calc_integ(&self, integ: &ChebSeries) -> enums::Value {
        unsafe { ffi::gsl_cheb_calc_integ(integ.c, self.c) }
    }
}

impl Drop for ChebSeries {
    fn drop(&mut self) {
        unsafe { ffi::gsl_cheb_free(self.c) };
        self.c = ::std::ptr::null_mut();
    }
}

impl ffi::FFI<ffi::gsl_cheb_series> for ChebSeries {
    fn wrap(c: *mut ffi::gsl_cheb_series) -> ChebSeries {
        unsafe {
            ChebSeries {
                c: c,
                data: CSlice::new((*c).c, (*c).order as usize + 1)
            }
        }
    }

    fn unwrap(c: &ChebSeries) -> *mut ffi::gsl_cheb_series {
        c.c
    }
}