This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
a gsl_interp_poly type candidate submited [code maturity is low]
- From: "Dan, Ho-Jin" <hjdan at sys713 dot kaist dot ac dot kr>
- To: gsl-discuss at sources dot redhat dot com
- Date: Mon, 26 Nov 2001 04:06:17 +0900
- Subject: a gsl_interp_poly type candidate submited [code maturity is low]
Hi there,
I couldn't seek the global polynomial interpolation for distinct n
points similar
to POLINT in the slatec library. The global polynomial interpolation
generally
is not good for interpolating data since its error behavior is not so
good. But
it is frequently used as a basic building block in other routines. (I
need also
another for complex data set). I have to submit the low mature code in
order to
confirm my work.
1. I am a novice in developing a free software, so I confuse how to
submit codes,
patches and manuals(I can handle a document of the latex format).
2. in the vim editor, how to set the gnu-recommended c-indenting.
3. If a submitted code doesn't seem to be useful, please suggest another
proxy
4. Is the name of "polynomial" type good or too general?
5. One of Higher-level interfaces will generate coefs. of a polynomial for
gsl_poly_eval(), gsl_poly_solve() etc.
6. complex counterparts will be submited.
note. GNU GPL lisence will be attached.
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_interp.h>
/* Note 1: gsl_interp_polynomial type added
* is polynomial name(analogous to slatec) right?
* another candidate is a divdif (divided difference)
* Note 2: This code is under development
* Note 3: ref: K. E. Atkinson, An introduction to Numerical Analysis 2nd ed.,
* John Wiley & Sons, Inc., Singapore, pp. 138--147
* Dan, Ho-Jin ( mailto:hjdan@sys713.kaist.ac.kr )
* Nov. 26, 2001
*/
/* for ctags ref.
* gsl_interp_alloc
* gsl_spline_alloc
* gsl_interp
* cspline_state_t
* gsl_interp_cspline
*/
typedef struct {
double* d;
} polynomial_state_t;
void* polynomial_alloc(size_t size)
{
polynomial_state_t* state = (polynomial_state_t*)malloc(
sizeof(polynomial_state_t));
if (state == 0) {
GSL_ERROR_NULL("failed to allocate space for polynomial_state_t",
GSL_ENOMEM);
}
state->d = (double*)malloc(sizeof(double)*size);
if (state->d == 0) {
free(state);
GSL_ERROR_NULL("failed to allocate space for polynomial_state_t.d",
GSL_ENOMEM);
}
return state;
}
int polynomial_init(void *vstate,
const double xa[], const double ya[], size_t size)
{
polynomial_state_t* state = (polynomial_state_t*) vstate;
int i, j;
/* newton's divide difference */
state->d[0] = ya[0];
for (j=size-1; j>=1; j--) {
state->d[j] = (ya[j] - ya[j-1])/(xa[j] - xa[j-1]);
}
for (i=2; i<size; i++) {
for (j = size-1; j>=i; j--) {
state->d[j] = (state->d[j] - state->d[j-1])/(xa[j] - xa[j-i]);
}
}
/*
for (i=0; i<size; i++)
printf("[%d]%e\n", i, state->d[i]);
*/
return GSL_SUCCESS;
}
int polynomial_eval(const void* vstate,
const double xa[], const double ya[], size_t size, double x,
gsl_interp_accel *acc, double * y)
{
polynomial_state_t* state = (polynomial_state_t*) vstate;
int i;
*y = state->d[size - 1];
for (i=size-2; i>=0; i--) {
*y = state->d[i] + (x - xa[i])*(*y);
}
return GSL_SUCCESS;
}
void polynomial_free(void* vstate)
{
polynomial_state_t* state;
free(state->d);
free(state);
}
static const gsl_interp_type polynomial_type =
{
"polynomial",
2,
&polynomial_alloc, /* alloc, not applicable */
&polynomial_init,
&polynomial_eval,
NULL, /* currenly unavailable */
NULL, /* currenly unavailable */
NULL, /* currenly unavailable */
&polynomial_free, /* free, not applicable */
};
const gsl_interp_type * gsl_interp_polynomial = &polynomial_type;
int main()
{
int i;
double xi, yi;
double x[5] = {2., 2.1, 2.2, 2.3, 2.4};
double y[5] = {1.414214, 1.449138, 1.483240, 1.516575, 1.549193};
printf("#m=0,S=2\n");
for (i=0; i<5; i++)
{
printf("%g %g\n", x[i], y[i]);
}
printf("#m=1,S=0\n");
{
gsl_interp* divdif = gsl_interp_alloc(gsl_interp_polynomial, 5);
gsl_interp_init(divdif, x, y, 5);
for (xi = x[0]; xi < x[4]; xi+=0.01)
{
yi = gsl_interp_eval(divdif, x, y, xi, 0);
printf("%g %g\n", xi, yi);
}
gsl_interp_free(divdif);
}
return 0;
}