1 /*
2 * Mathlib : A C Library of Special Functions
3 * Copyright (C) 1998 Ross Ihaka
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, a copy is available at
17 * http://www.r-project.org/Licenses/
18 *
19 * SYNOPSIS
20 *
21 * int chebyshev_init(double *dos, int nos, double eta)
22 * double chebyshev_eval(double x, double *a, int n)
23 *
24 * DESCRIPTION
25 *
26 * "chebyshev_init" determines the number of terms for the
27 * double precision orthogonal series "dos" needed to insure
28 * the error is no larger than "eta". Ordinarily eta will be
29 * chosen to be one-tenth machine precision.
30 *
31 * "chebyshev_eval" evaluates the n-term Chebyshev series
32 * "a" at "x".
33 *
34 * NOTES
35 *
36 * These routines are translations into C of Fortran routines
37 * by W. Fullerton of Los Alamos Scientific Laboratory.
38 *
39 * Based on the Fortran routine dcsevl by W. Fullerton.
40 * Adapted from R. Broucke, Algorithm 446, CACM., 16, 254 (1973).
41 */
42
43 #include "nmath.h"
44
45 /* NaNs propagated correctly */
46
47
chebyshev_init(double * dos,int nos,double eta)48 int attribute_hidden chebyshev_init(double *dos, int nos, double eta)
49 {
50 int i, ii;
51 double err;
52
53 if (nos < 1)
54 return 0;
55
56 err = 0.0;
57 i = 0; /* just to avoid compiler warnings */
58 for (ii=1; ii<=nos; ii++) {
59 i = nos - ii;
60 err += fabs(dos[i]);
61 if (err > eta) {
62 return i;
63 }
64 }
65 return i;
66 }
67
68
chebyshev_eval(double x,const double * a,const int n)69 double attribute_hidden chebyshev_eval(double x, const double *a, const int n)
70 {
71 double b0, b1, b2, twox;
72 int i;
73
74 if (n < 1 || n > 1000) ML_ERR_return_NAN;
75
76 if (x < -1.1 || x > 1.1) ML_ERR_return_NAN;
77
78 twox = x * 2;
79 b2 = b1 = 0;
80 b0 = 0;
81 for (i = 1; i <= n; i++) {
82 b2 = b1;
83 b1 = b0;
84 b0 = twox * b1 - b2 + a[n - i];
85 }
86 return (b0 - b2) * 0.5;
87 }
88