a gsl_interp_poly type candidate submited [code maturity is low]

a gsl_interp_poly type candidate submited [code maturity is low]

[Dan, Ho-Jin]

[-- Attachment #1: Type: text/plain, Size: 997 bytes --]

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.




[-- Attachment #2: interp_poly.c --]
[-- Type: text/plain, Size: 3041 bytes --]

#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;
}
   

Re: a gsl_interp_poly type candidate submited [code maturity is low]

[Brian Gough]

Dan, Ho-Jin writes:
 > Hi there, I couldn't seek the global polynomial interpolation for
 > distinct n points similar to POLINT in the slatec library.

Thanks, that looks good. Can you add the derivative and integral
methods too to complete the interface?

Brian

Re: a gsl_interp_poly type candidate submited [code maturity is low]

[Dan, Ho-Jin]

Okay, I'll do the work soon.
have a nice day...

ho-jin

Brian Gough wrote:

>Dan, Ho-Jin writes:
> > Hi there, I couldn't seek the global polynomial interpolation for
> > distinct n points similar to POLINT in the slatec library.
>
>Thanks, that looks good. Can you add the derivative and integral
>methods too to complete the interface?
>
>Brian
>
>

a gsl_interp_poly type candidate submited [code maturity is low]

[Dan, Ho-Jin]

[-- Attachment #1: Type: text/plain, Size: 997 bytes --]

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.




[-- Attachment #2: interp_poly.c --]
[-- Type: text/plain, Size: 3041 bytes --]

#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;
}
   

Re: a gsl_interp_poly type candidate submited [code maturity is low]

[Dan, Ho-Jin]

Okay, I'll do the work soon.
have a nice day...

ho-jin

Brian Gough wrote:

>Dan, Ho-Jin writes:
> > Hi there, I couldn't seek the global polynomial interpolation for
> > distinct n points similar to POLINT in the slatec library.
>
>Thanks, that looks good. Can you add the derivative and integral
>methods too to complete the interface?
>
>Brian
>
>

Re: a gsl_interp_poly type candidate submited [code maturity is low]

[Brian Gough]

Dan, Ho-Jin writes:
 > Hi there, I couldn't seek the global polynomial interpolation for
 > distinct n points similar to POLINT in the slatec library.

Thanks, that looks good. Can you add the derivative and integral
methods too to complete the interface?

Brian

a gsl_interp_poly type candidate submited [code maturity is low]

[Dan, Ho-Jin]

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.

Re: a gsl_interp_poly type candidate submited [code maturity is low]

[Brian Gough]

Dan, Ho-Jin writes:
 > Hi there, I couldn't seek the global polynomial interpolation for
 > distinct n points similar to POLINT in the slatec library.

Thanks, that looks good. Can you add the derivative and integral
methods too to complete the interface?

Brian