r2rs-nmath 0.1.1

// Translation of nmath's chebyshev
//
// Copyright (C) 2020 neek-sss
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <https://www.gnu.org/licenses/>.

use num_traits::Float;

/// "chebyshev_init" determines the number of terms for the
/// double precision orthogonal series "dos" needed to insure
/// the error is no larger than "eta".  Ordinarily eta will be
/// chosen to be one-tenth machine precision.
pub fn chebyshev_init(dos: &[f64], nos: usize, eta: f64) -> usize {
    let mut i = 0; /* just to avoid compiler warnings */
    let mut ii = 0;
    let mut err = 0.0;

    if nos < 1 {
        return 0;
    }

    err = 0.0;
    ii = 1;
    while ii <= nos {
        i = nos - ii;
        err += (dos[i]).abs();
        if err > eta {
            return i;
        }
        ii += 1
    }
    return i;
}

/// "chebyshev_eval" evaluates the n-term Chebyshev series
/// "a" at "x".
///
/// Based on the Fortran routine dcsevl by W. Fullerton.
/// Adapted from R. Broucke, Algorithm 446, CACM., 16, 254 (1973).
pub fn chebyshev_eval(x: f64, a: &[f64], n: usize) -> f64 {
    let mut b0 = 0.0;
    let mut b1 = 0.0;
    let mut b2 = 0.0;
    let mut twox = 0.0;
    let mut i = 0;

    if n < 1 || n > 1000 {
        return f64::nan();
    }

    if x < -1.1 || x > 1.1 {
        return f64::nan();
    }

    twox = x * 2.0;
    b1 = 0.0;
    b2 = b1;
    b0 = 0.0;
    i = 1;
    while i <= n {
        b2 = b1;
        b1 = b0;
        b0 = twox * b1 - b2 + a[n - i];
        i += 1
    }
    return (b0 - b2) * 0.5;
}