/src/r-source/src/nmath/chebyshev.c
Line | Count | Source |
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 | | * https://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 | | |
48 | | attribute_hidden int chebyshev_init(double *dos, int nos, double eta) |
49 | 0 | { |
50 | 0 | int i, ii; |
51 | 0 | double err; |
52 | |
|
53 | 0 | if (nos < 1) |
54 | 0 | return 0; |
55 | | |
56 | 0 | err = 0.0; |
57 | 0 | i = 0; /* just to avoid compiler warnings */ |
58 | 0 | for (ii=1; ii<=nos; ii++) { |
59 | 0 | i = nos - ii; |
60 | 0 | err += fabs(dos[i]); |
61 | 0 | if (err > eta) { |
62 | 0 | return i; |
63 | 0 | } |
64 | 0 | } |
65 | 0 | return i; |
66 | 0 | } |
67 | | |
68 | | |
69 | | attribute_hidden double chebyshev_eval(double x, const double *a, const int n) |
70 | 0 | { |
71 | 0 | double b0, b1, b2, twox; |
72 | 0 | int i; |
73 | |
|
74 | 0 | if (n < 1 || n > 1000) ML_WARN_return_NAN; |
75 | |
|
76 | 0 | if (x < -1.1 || x > 1.1) ML_WARN_return_NAN; |
77 | |
|
78 | 0 | twox = x * 2; |
79 | 0 | b2 = b1 = 0; |
80 | 0 | b0 = 0; |
81 | 0 | for (i = 1; i <= n; i++) { |
82 | 0 | b2 = b1; |
83 | 0 | b1 = b0; |
84 | 0 | b0 = twox * b1 - b2 + a[n - i]; |
85 | 0 | } |
86 | 0 | return (b0 - b2) * 0.5; |
87 | 0 | } |