The R Project SVN R

Rev

Rev 42302 | Rev 77685 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 42302 Rev 42303
Line -... Line 1...
- 
 
 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 
 
- 
 
 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 
 
- 
 
 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 
}