/*
integration.c
Ruby/GSL: Ruby extension library for GSL (GNU Scientific Library)
(C) Copyright 2001-2006 by Yoshiki Tsunesada
Ruby/GSL is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License.
This library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY.
*/
#include "rb_gsl_config.h"
#include "rb_gsl_array.h"
#include "rb_gsl_function.h"
#include "rb_gsl_integration.h"
#include "rb_gsl_common.h"
#ifndef CHECK_WORKSPACE
#define CHECK_WORKSPACE(x) if(CLASS_OF(x)!=cgsl_integration_workspace)\
rb_raise(rb_eTypeError,\
"wrong argument type %s (Integration::Workspace expected)",\
rb_class2name(CLASS_OF(x)));
#endif
#define EPSABS_DEFAULT 0.0
#define EPSREL_DEFAULT 1e-10
#define LIMIT_DEFAULT 1000
#define KEY_DEFAULT GSL_INTEG_GAUSS61
static VALUE cgsl_integration_qaws_table, cgsl_integration_qawo_table;
static VALUE cgsl_integration_workspace;
static int get_a_b(int argc, VALUE *argv, int argstart, double *a, double *b);
static int get_epsabs_epsrel(int argc, VALUE *argv, int argstart,
double *epsabs, double *epsrel);
static int get_a_b_epsabs_epsrel(int argc, VALUE *argv, int argstart,
double *a, double *b, double *epsabs,
double *epsrel);
static int get_limit_key_workspace(int argc, VALUE *argv, int argstart,
size_t *limit, int *key,
gsl_integration_workspace **w);
static int get_limit_workspace(int argc, VALUE *argv, int argstart,
size_t *limit,
gsl_integration_workspace **w);
static int get_epsabs_epsrel_limit_workspace(int argc, VALUE *argv, int argstart,
double *epsabs, double *epsrel,
size_t *limit,
gsl_integration_workspace **w);
static int get_a_b(int argc, VALUE *argv, int argstart, double *a, double *b)
{
int itmp;
VALUE aa, bb;
if (argstart >= argc) return argstart;
if (TYPE(argv[argstart]) == T_ARRAY) {
aa = rb_ary_entry(argv[argstart], 0);
bb = rb_ary_entry(argv[argstart], 1);
Need_Float(aa); Need_Float(bb);
*a = RFLOAT(aa)->value;
*b = RFLOAT(bb)->value;
itmp = argstart + 1;
} else {
Need_Float(argv[argstart]); Need_Float(argv[argstart+1]);
*a = NUM2DBL(argv[argstart]);
*b = NUM2DBL(argv[argstart+1]);
itmp = argstart + 2;
}
return itmp;
}
static int get_epsabs_epsrel(int argc, VALUE *argv, int argstart,
double *epsabs, double *epsrel)
{
int itmp;
VALUE aa, bb;
*epsabs = EPSABS_DEFAULT;
*epsrel = EPSREL_DEFAULT;
if (argstart >= argc) return argstart;
if (TYPE(argv[argstart]) == T_ARRAY) {
aa = rb_ary_entry(argv[argstart], 0);
bb = rb_ary_entry(argv[argstart], 1);
Need_Float(aa); Need_Float(bb);
*epsabs = NUM2DBL(aa);
*epsrel = NUM2DBL(bb);
itmp = 1;
} else {
Need_Float(argv[argstart]); Need_Float(argv[argstart+1]);
*epsabs = NUM2DBL(argv[argstart]);
*epsrel = NUM2DBL(argv[argstart+1]);
itmp = 2;
}
return argstart + itmp;
}
static int get_a_b_epsabs_epsrel(int argc, VALUE *argv, int argstart,
double *a, double *b, double *epsabs,
double *epsrel)
{
int itmp;
*epsabs = EPSABS_DEFAULT;
*epsrel = EPSREL_DEFAULT;
itmp = get_a_b(argc, argv, argstart, a, b);
itmp = get_epsabs_epsrel(argc, argv, itmp, epsabs, epsrel);
return itmp;
}
static int get_limit_key_workspace(int argc, VALUE *argv, int argstart,
size_t *limit, int *key,
gsl_integration_workspace **w)
{
int flag = 0;
switch (argc-argstart) {
case 3:
CHECK_FIXNUM(argv[argstart]);
CHECK_FIXNUM(argv[argstart+1]);
CHECK_WORKSPACE(argv[argstart+2]);
*limit = FIX2INT(argv[argstart]);
*key = FIX2INT(argv[argstart+1]);
Data_Get_Struct(argv[argstart+2], gsl_integration_workspace, *w);
flag = 0;
break;
case 1:
CHECK_FIXNUM(argv[argstart]);
*key = FIX2INT(argv[argstart]);
*limit = LIMIT_DEFAULT;
*w = gsl_integration_workspace_alloc(*limit);
flag = 1;
break;
case 2:
if (TYPE(argv[argc-1]) == T_FIXNUM) {
CHECK_FIXNUM(argv[argc-2]);
*limit = FIX2INT(argv[argc-2]);
*key = FIX2INT(argv[argc-1]);
*w = gsl_integration_workspace_alloc(*limit);
flag = 1;
} else {
CHECK_FIXNUM(argv[argc-2]);
CHECK_WORKSPACE(argv[argc-1]);
*key = FIX2INT(argv[argc-2]);
Data_Get_Struct(argv[argc-1], gsl_integration_workspace, *w);
*limit = (*w)->limit;
flag = 0;
}
break;
case 0:
*key = KEY_DEFAULT;
*limit = LIMIT_DEFAULT;
*w = gsl_integration_workspace_alloc(*limit);
flag = 1;
break;
default:
rb_raise(rb_eArgError, "wrong number of arguments");
break;
}
if (*w == NULL) rb_raise(rb_eRuntimeError, "something wrong with workspace");
return flag;
}
static int get_limit_workspace(int argc, VALUE *argv, int argstart,
size_t *limit,
gsl_integration_workspace **w)
{
int flag = 0;
switch (argc-argstart) {
case 2:
CHECK_FIXNUM(argv[argstart]);
*limit = FIX2INT(argv[argstart]);
CHECK_WORKSPACE(argv[argstart+1]);
Data_Get_Struct(argv[argstart+1], gsl_integration_workspace, *w);
flag = 0;
break;
case 0:
*limit = LIMIT_DEFAULT;
*w = gsl_integration_workspace_alloc(*limit);
flag = 1;
break;
case 1:
switch (TYPE(argv[argstart])) {
case T_FIXNUM:
case T_BIGNUM:
CHECK_FIXNUM(argv[argstart]);
*limit = FIX2INT(argv[argstart]);
*w = gsl_integration_workspace_alloc(*limit);
flag = 1;
break;
default:
CHECK_WORKSPACE(argv[argc-1]);
Data_Get_Struct(argv[argc-1], gsl_integration_workspace, *w);
*limit = (*w)->limit;
flag = 0;
break;
}
break;
default:
rb_raise(rb_eArgError, "wrong number of arguments");
break;
}
if (*w == NULL) rb_raise(rb_eRuntimeError, "something wrong with workspace");
return flag;
}
static int get_epsabs_epsrel_limit_workspace(int argc, VALUE *argv, int argstart,
double *epsabs, double *epsrel,
size_t *limit,
gsl_integration_workspace **w)
{
int flag = 0, itmp;
itmp = argstart;
*epsabs = EPSABS_DEFAULT;
*epsrel = EPSREL_DEFAULT;
*limit = LIMIT_DEFAULT;
switch (argc-itmp) {
case 0:
*w = gsl_integration_workspace_alloc(*limit);
flag = 1;
break;
case 1:
if (TYPE(argv[itmp]) == T_ARRAY) {
get_epsabs_epsrel(argc, argv, itmp, epsabs, epsrel);
*w = gsl_integration_workspace_alloc(*limit);
flag = 1;
} else {
flag = get_limit_workspace(argc, argv, itmp, limit, w);
}
break;
case 2:
case 3:
switch (TYPE(argv[itmp])) {
case T_ARRAY:
itmp = get_epsabs_epsrel(argc, argv, itmp, epsabs, epsrel);
flag = get_limit_workspace(argc, argv, itmp, limit, w);
break;
case T_FLOAT:
get_epsabs_epsrel(argc, argv, itmp, epsabs, epsrel);
*w = gsl_integration_workspace_alloc(*limit);
flag = 1;
break;
default:
flag = get_limit_workspace(argc, argv, itmp, limit, w);
break;
}
break;
case 4:
itmp = get_epsabs_epsrel(argc, argv, itmp, epsabs, epsrel);
flag = get_limit_workspace(argc, argv, itmp, limit, w);
break;
default:
rb_raise(rb_eArgError, "wrong number of arguments");
break;
}
if (*w == NULL) rb_raise(rb_eRuntimeError, "something wrong with workspace");
return flag;
}
static VALUE rb_gsl_integration_qng(int argc, VALUE *argv, VALUE obj)
{
double a, b, epsabs = EPSABS_DEFAULT, epsrel = EPSREL_DEFAULT;
double result, abserr;
size_t neval;
gsl_function *F = NULL;
int status, itmp;
if (argc < 1) rb_raise(rb_eArgError,
"wrong number of arguments (%d for >= 1)", argc);
switch (TYPE(obj)) {
case T_MODULE: case T_CLASS: case T_OBJECT:
CHECK_FUNCTION(argv[0]);
Data_Get_Struct(argv[0], gsl_function, F);
itmp = get_a_b_epsabs_epsrel(argc, argv, 1, &a, &b, &epsabs, &epsrel);
break;
default:
itmp = get_a_b_epsabs_epsrel(argc, argv, 0, &a, &b, &epsabs, &epsrel);
Data_Get_Struct(obj, gsl_function, F);
break;
}
status = gsl_integration_qng(F, a, b, epsabs, epsrel,
&result, &abserr, &neval);
return rb_ary_new3(4, rb_float_new(result), rb_float_new(abserr),
INT2FIX(neval), INT2FIX(status));
}
static VALUE rb_gsl_integration_qag(int argc, VALUE *argv, VALUE obj)
{
double a, b, epsabs = EPSABS_DEFAULT, epsrel = EPSREL_DEFAULT;
double result, abserr;
size_t limit = LIMIT_DEFAULT;
gsl_function *F = NULL;
gsl_integration_workspace *w = NULL;
int key = KEY_DEFAULT, status, intervals, itmp, flag = 0;
if (argc < 1) rb_raise(rb_eArgError,
"wrong number of arguments (%d for >= 1)", argc);
switch (TYPE(obj)) {
case T_MODULE: case T_CLASS: case T_OBJECT:
CHECK_FUNCTION(argv[0]);
Data_Get_Struct(argv[0], gsl_function, F);
if (argc == 3) {
CHECK_FIXNUM(argv[2]);
get_a_b(argc, argv, 1, &a, &b);
key = FIX2INT(argv[2]);
w = gsl_integration_workspace_alloc(limit);
flag = 1;
} else if (argc == 4) {
CHECK_FIXNUM(argv[3]);
get_a_b(argc, argv, 1, &a, &b);
key = FIX2INT(argv[3]);
w = gsl_integration_workspace_alloc(limit);
flag = 1;
} else {
itmp = get_a_b_epsabs_epsrel(argc, argv, 1, &a, &b, &epsabs, &epsrel);
flag = get_limit_key_workspace(argc, argv, itmp, &limit, &key, &w);
}
break;
default:
if (argc == 2) {
if (FIXNUM_P(argv[1])) {
key = FIX2INT(argv[1]);
w = gsl_integration_workspace_alloc(limit);
flag = 1;
} else if (rb_obj_is_kind_of(argv[1], cgsl_integration_workspace)) {
Data_Get_Struct(argv[1], gsl_integration_workspace, w);
flag = 0;
} else {
rb_raise(rb_eTypeError, "Key or Workspace expected");
}
itmp = get_a_b(argc, argv, 0, &a, &b);
} else if (argc == 3) {
if (FIXNUM_P(argv[2])) {
key = FIX2INT(argv[2]);
w = gsl_integration_workspace_alloc(limit);
flag = 1;
} else if (rb_obj_is_kind_of(argv[2], cgsl_integration_workspace)) {
Data_Get_Struct(argv[2], gsl_integration_workspace, w);
flag = 0;
} else {
rb_raise(rb_eTypeError, "Key or Workspace expected");
}
itmp = get_a_b(argc, argv, 0, &a, &b);
} else {
itmp = get_a_b_epsabs_epsrel(argc, argv, 0, &a, &b, &epsabs, &epsrel);
flag = get_limit_key_workspace(argc, argv, itmp, &limit, &key, &w);
}
Data_Get_Struct(obj, gsl_function, F);
break;
}
status = gsl_integration_qag(F, a, b, epsabs, epsrel, limit, key, w,
&result, &abserr);
intervals = w->size;
if (flag == 1) gsl_integration_workspace_free(w);
return rb_ary_new3(4, rb_float_new(result), rb_float_new(abserr),
INT2FIX(intervals), INT2FIX(status));
}
static VALUE rb_gsl_integration_qags(int argc, VALUE *argv, VALUE obj)
{
double a, b, epsabs = EPSABS_DEFAULT, epsrel = EPSREL_DEFAULT;
double result, abserr;
size_t limit = LIMIT_DEFAULT;
gsl_function *F = NULL;
gsl_integration_workspace *w = NULL;
int status, intervals, flag = 0, itmp;
switch (TYPE(obj)) {
case T_MODULE: case T_CLASS: case T_OBJECT:
CHECK_FUNCTION(argv[0]);
Data_Get_Struct(argv[0], gsl_function, F);
itmp = get_a_b(argc, argv, 1, &a, &b);
break;
default:
Data_Get_Struct(obj, gsl_function, F);
itmp = get_a_b(argc, argv, 0, &a, &b);
break;
}
flag = get_epsabs_epsrel_limit_workspace(argc, argv, itmp, &epsabs, &epsrel,
&limit, &w);
status = gsl_integration_qags(F, a, b, epsabs, epsrel, limit, w,
&result, &abserr);
intervals = w->size;
if (flag == 1) gsl_integration_workspace_free(w);
return rb_ary_new3(4, rb_float_new(result), rb_float_new(abserr),
INT2FIX(intervals), INT2FIX(status));
}
static VALUE rb_gsl_integration_qagp(int argc, VALUE *argv, VALUE obj)
{
double epsabs, epsrel;
double result, abserr;
size_t limit;
gsl_function *F = NULL;
gsl_vector *v = NULL;
gsl_integration_workspace *w = NULL;
int status, intervals, flag = 0, flag2 = 0, itmp;
switch (TYPE(obj)) {
case T_MODULE: case T_CLASS: case T_OBJECT:
CHECK_FUNCTION(argv[0]);
Data_Get_Struct(argv[0], gsl_function, F);
itmp = 1;
break;
default:
Data_Get_Struct(obj, gsl_function, F);
itmp = 0;
break;
}
if (TYPE(argv[itmp]) == T_ARRAY) {
v = make_cvector_from_rarray(argv[itmp]);
flag2 = 1;
} else {
Data_Get_Vector(argv[itmp], v);
flag2 = 0;
}
itmp += 1;
flag = get_epsabs_epsrel_limit_workspace(argc, argv, itmp, &epsabs, &epsrel,
&limit, &w);
status = gsl_integration_qagp(F, v->data, v->size, epsabs, epsrel, limit, w,
&result, &abserr);
intervals = w->size;
if (flag == 1) gsl_integration_workspace_free(w);
if (flag2 == 1) gsl_vector_free(v);
return rb_ary_new3(4, rb_float_new(result), rb_float_new(abserr),
INT2FIX(intervals), INT2FIX(status));
}
/* (-infty --- +infty) */
static VALUE rb_gsl_integration_qagi(int argc, VALUE *argv, VALUE obj)
{
double epsabs, epsrel;
double result, abserr;
size_t limit;
gsl_function *F = NULL;
gsl_integration_workspace *w = NULL;
int status, intervals, flag = 0, itmp;
switch (TYPE(obj)) {
case T_MODULE: case T_CLASS: case T_OBJECT:
CHECK_FUNCTION(argv[0]);
Data_Get_Struct(argv[0], gsl_function, F);
itmp = 1;
break;
default:
Data_Get_Struct(obj, gsl_function, F);
itmp = 0;
break;
}
flag = get_epsabs_epsrel_limit_workspace(argc, argv, itmp, &epsabs, &epsrel,
&limit, &w);
status = gsl_integration_qagi(F, epsabs, epsrel, limit, w,
&result, &abserr);
intervals = w->size;
if (flag == 1) gsl_integration_workspace_free(w);
return rb_ary_new3(4, rb_float_new(result), rb_float_new(abserr),
INT2FIX(intervals), INT2FIX(status));
}
/* (a --- +infty) */
static VALUE rb_gsl_integration_qagiu(int argc, VALUE *argv, VALUE obj)
{
double a, epsabs, epsrel;
double result, abserr;
size_t limit;
gsl_function *F = NULL;
gsl_integration_workspace *w = NULL;
int status, intervals, flag = 0, itmp;
switch (TYPE(obj)) {
case T_MODULE: case T_CLASS: case T_OBJECT:
CHECK_FUNCTION(argv[0]);
Data_Get_Struct(argv[0], gsl_function, F);
itmp = 1;
break;
default:
Data_Get_Struct(obj, gsl_function, F);
itmp = 0;
break;
}
Need_Float(argv[itmp]);
a = NUM2DBL(argv[itmp]);
itmp += 1;
flag = get_epsabs_epsrel_limit_workspace(argc, argv, itmp, &epsabs, &epsrel,
&limit, &w);
status = gsl_integration_qagiu(F, a, epsabs, epsrel, limit, w,
&result, &abserr);
intervals = w->size;
if (flag == 1) gsl_integration_workspace_free(w);
return rb_ary_new3(4, rb_float_new(result), rb_float_new(abserr),
INT2FIX(intervals), INT2FIX(status));
}
/* (-infty --- b) */
static VALUE rb_gsl_integration_qagil(int argc, VALUE *argv, VALUE obj)
{
double b, epsabs, epsrel;
double result, abserr;
size_t limit;
gsl_function *F = NULL;
gsl_integration_workspace *w = NULL;
int status, intervals, flag = 0, itmp;
switch (TYPE(obj)) {
case T_MODULE: case T_CLASS: case T_OBJECT:
CHECK_FUNCTION(argv[0]);
Data_Get_Struct(argv[0], gsl_function, F);
itmp = 1;
break;
default:
Data_Get_Struct(obj, gsl_function, F);
itmp = 0;
break;
}
Need_Float(argv[itmp]);
b = NUM2DBL(argv[itmp]);
flag = get_epsabs_epsrel_limit_workspace(argc, argv, itmp+1, &epsabs, &epsrel,
&limit, &w);
Data_Get_Struct(obj, gsl_function, F);
status = gsl_integration_qagil(F, b, epsabs, epsrel, limit, w,
&result, &abserr);
intervals = w->size;
if (flag == 1) gsl_integration_workspace_free(w);
return rb_ary_new3(4, rb_float_new(result), rb_float_new(abserr),
INT2FIX(intervals), INT2FIX(status));
}
static VALUE rb_gsl_integration_qawc(int argc, VALUE *argv, VALUE obj)
{
double a, b, c, epsabs, epsrel;
double result, abserr;
size_t limit;
gsl_function *F = NULL;
gsl_integration_workspace *w = NULL;
int status, intervals, itmp, flag = 0;
switch (TYPE(obj)) {
case T_MODULE: case T_CLASS: case T_OBJECT:
CHECK_FUNCTION(argv[0]);
Data_Get_Struct(argv[0], gsl_function, F);
itmp = 1;
break;
default:
Data_Get_Struct(obj, gsl_function, F);
itmp = 0;
break;
}
itmp = get_a_b(argc, argv, itmp, &a, &b);
if (argc-itmp <= 0) rb_raise(rb_eArgError, "The pole is not given");
Need_Float(argv[itmp]);
c = NUM2DBL(argv[itmp]);
flag = get_epsabs_epsrel_limit_workspace(argc, argv, itmp+1, &epsabs, &epsrel,
&limit, &w);
status = gsl_integration_qawc(F, a, b, c, epsabs, epsrel, limit, w, &result, &abserr);
intervals = w->size;
if (flag == 1) gsl_integration_workspace_free(w);
return rb_ary_new3(4, rb_float_new(result), rb_float_new(abserr), INT2FIX(intervals),
INT2FIX(status));
}
VALUE rb_gsl_integration_qaws_table_alloc(int argc, VALUE *argv, VALUE klass)
{
gsl_integration_qaws_table *t = NULL;
VALUE alpha, beta, mu, nu;
if (TYPE(argv[0]) == T_ARRAY) {
alpha = rb_ary_entry(argv[0], 0);
beta = rb_ary_entry(argv[0], 1);
mu = rb_ary_entry(argv[0], 2);
nu = rb_ary_entry(argv[0], 3);
} else {
Need_Float(argv[0]); Need_Float(argv[1]);
CHECK_FIXNUM(argv[2]); CHECK_FIXNUM(argv[3]);
alpha = argv[0];
beta = argv[1];
mu = argv[2];
nu = argv[3];
}
t = gsl_integration_qaws_table_alloc(NUM2DBL(alpha), NUM2DBL(beta),
FIX2INT(mu), FIX2INT(nu));
return Data_Wrap_Struct(klass, 0, gsl_integration_qaws_table_free, t);
}
static VALUE rb_gsl_integration_qaws_table_set(int argc, VALUE *argv, VALUE obj)
{
gsl_integration_qaws_table *t = NULL;
double alpha, beta;
int mu, nu, type;
if (argc != 1 && argc != 4)
rb_raise(rb_eArgError, "wrong number of argument (%d for 1 or 3)", argc);
type = TYPE(argv[0]);
Data_Get_Struct(obj, gsl_integration_qaws_table, t);
if (type == T_FIXNUM || type == T_BIGNUM || type == T_FLOAT) {
alpha = NUM2DBL(argv[0]);
beta = NUM2DBL(argv[1]);
mu = FIX2INT(argv[2]);
nu = FIX2INT(argv[3]);
} else if (type == T_ARRAY) {
alpha = NUM2DBL(rb_ary_entry(argv[0], 0));
beta = NUM2DBL(rb_ary_entry(argv[0], 1));
mu = FIX2INT(rb_ary_entry(argv[0], 2));
nu = FIX2INT(rb_ary_entry(argv[0], 3));
} else {
rb_raise(rb_eTypeError, "wrong argument type %s", rb_class2name(CLASS_OF(argv[0])));
}
gsl_integration_qaws_table_set(t, alpha, beta, mu, nu);
return obj;
}
static VALUE rb_gsl_integration_qaws_table_to_a(VALUE obj)
{
gsl_integration_qaws_table *t = NULL;
VALUE ary;
Data_Get_Struct(obj, gsl_integration_qaws_table, t);
ary = rb_ary_new2(4);
rb_ary_store(ary, 0, rb_float_new(t->alpha));
rb_ary_store(ary, 1, rb_float_new(t->beta));
rb_ary_store(ary, 2, INT2FIX(t->mu));
rb_ary_store(ary, 3, INT2FIX(t->nu));
return ary;
}
static gsl_integration_qaws_table* make_qaws_table(VALUE ary);
static VALUE rb_gsl_ary_to_integration_qaws_table(VALUE ary)
{
gsl_integration_qaws_table *t = NULL;
t = make_qaws_table(ary);
return Data_Wrap_Struct(cgsl_integration_qaws_table,
0, gsl_integration_qaws_table_free, t);
}
static gsl_integration_qaws_table* make_qaws_table(VALUE ary)
{
double alpha, beta;
int mu, nu;
alpha = NUM2DBL(rb_ary_entry(ary, 0));
beta = NUM2DBL(rb_ary_entry(ary, 1));
mu = FIX2INT(rb_ary_entry(ary, 2));
nu = FIX2INT(rb_ary_entry(ary, 3));
return gsl_integration_qaws_table_alloc(alpha, beta, mu, nu);
}
static VALUE rb_gsl_integration_qaws(int argc, VALUE *argv, VALUE obj)
{
double a, b, epsabs, epsrel;
double result, abserr;
size_t limit;
gsl_function *F = NULL;
gsl_integration_workspace *w = NULL;
gsl_integration_qaws_table *t = NULL;
int status, intervals, itmp, flag = 0, flagt = 0;
switch (TYPE(obj)) {
case T_MODULE: case T_CLASS: case T_OBJECT:
if (argc < 2) rb_raise(rb_eArgError, "too few arguments");
CHECK_FUNCTION(argv[0]);
Data_Get_Struct(argv[0], gsl_function, F);
itmp = 1;
break;
default:
if (argc < 1) rb_raise(rb_eArgError, "too few arguments");
Data_Get_Struct(obj, gsl_function, F);
itmp = 0;
break;
}
itmp = get_a_b(argc, argv, itmp, &a, &b);
if (TYPE(argv[itmp]) == T_ARRAY) {
flagt = 1;
t = make_qaws_table(argv[itmp]);
} else {
flagt = 0;
if (!rb_obj_is_kind_of(argv[itmp], cgsl_integration_qaws_table))
rb_raise(rb_eTypeError, "Integration::QAWS_Table expected");
Data_Get_Struct(argv[itmp], gsl_integration_qaws_table, t);
}
flag = get_epsabs_epsrel_limit_workspace(argc, argv, itmp+1, &epsabs, &epsrel,
&limit, &w);
status = gsl_integration_qaws(F, a, b, t, epsabs, epsrel, limit, w, &result, &abserr);
intervals = w->size;
if (flag == 1) gsl_integration_workspace_free(w);
if (flagt == 1) gsl_integration_qaws_table_free(t);
return rb_ary_new3(4, rb_float_new(result), rb_float_new(abserr), INT2FIX(intervals),
INT2FIX(status));
}
static gsl_integration_qawo_table* make_qawo_table(VALUE ary);
static VALUE rb_gsl_integration_qawo_table_alloc(int argc, VALUE *argv,
VALUE klass)
{
gsl_integration_qawo_table *t = NULL;
double omega, L;
enum gsl_integration_qawo_enum sine;
size_t n;
if (argc != 1 && argc != 4)
rb_raise(rb_eArgError, "wrong nubmer of arguments (%d for 1 or 4)", argc);
if (TYPE(argv[0]) == T_ARRAY) {
omega = NUM2DBL(rb_ary_entry(argv[0], 0));
L = NUM2DBL(rb_ary_entry(argv[0], 1));
sine = FIX2INT(rb_ary_entry(argv[0], 2));
n = FIX2INT(rb_ary_entry(argv[0], 3));
} else {
omega = NUM2DBL(argv[0]);
L = NUM2DBL(argv[1]);
sine = FIX2INT(argv[2]);
n = FIX2INT(argv[3]);
}
t = gsl_integration_qawo_table_alloc(omega, L, sine, n);
return Data_Wrap_Struct(klass, 0, gsl_integration_qawo_table_free, t);
}
static VALUE rb_gsl_integration_qawo_table_to_a(VALUE obj)
{
gsl_integration_qawo_table *t = NULL;
VALUE ary;
Data_Get_Struct(obj, gsl_integration_qawo_table, t);
ary = rb_ary_new2(4);
rb_ary_store(ary, 0, rb_float_new(t->omega));
rb_ary_store(ary, 1, rb_float_new(t->L));
rb_ary_store(ary, 2, INT2FIX(t->sine));
rb_ary_store(ary, 3, INT2FIX(t->n));
return ary;
}
static VALUE rb_gsl_ary_to_integration_qawo_table(VALUE ary)
{
gsl_integration_qawo_table *t = NULL;
t = make_qawo_table(ary);
return Data_Wrap_Struct(cgsl_integration_qawo_table,
0, gsl_integration_qawo_table_free, t);
}
static gsl_integration_qawo_table* make_qawo_table(VALUE ary)
{
double omega, L;
enum gsl_integration_qawo_enum sine;
size_t n;
omega = NUM2DBL(rb_ary_entry(ary, 0));
L = NUM2DBL(rb_ary_entry(ary, 1));
sine = FIX2INT(rb_ary_entry(ary, 2));
n = FIX2INT(rb_ary_entry(ary, 3));
return gsl_integration_qawo_table_alloc(omega, L, sine, n);
}
static VALUE rb_gsl_integration_qawo_table_set(int argc, VALUE *argv, VALUE obj)
{
gsl_integration_qawo_table *t = NULL;
double omega, L;
enum gsl_integration_qawo_enum sine;
int type;
if (argc != 1 && argc != 3)
rb_raise(rb_eArgError, "wrong number of argument (%d for 1 or 3)", argc);
type = TYPE(argv[0]);
Data_Get_Struct(obj, gsl_integration_qawo_table, t);
if (type == T_FIXNUM || type == T_BIGNUM || type == T_FLOAT) {
omega = NUM2DBL(argv[0]);
L = NUM2DBL(argv[1]);
sine = FIX2INT(argv[2]);
} else if (type == T_ARRAY) {
omega = NUM2DBL(rb_ary_entry(argv[0], 0));
L = NUM2DBL(rb_ary_entry(argv[0], 1));
sine = FIX2INT(rb_ary_entry(argv[0], 2));
} else {
rb_raise(rb_eTypeError, "wrong argument type %s", rb_class2name(CLASS_OF(argv[0])));
}
gsl_integration_qawo_table_set(t, omega, L, sine);
return obj;
}
static VALUE rb_gsl_integration_qawo_table_set_length(VALUE obj, VALUE L)
{
gsl_integration_qawo_table *t = NULL;
Need_Float(L);
Data_Get_Struct(obj, gsl_integration_qawo_table, t);
gsl_integration_qawo_table_set_length(t, NUM2DBL(L));
return obj;
}
static int get_qawo_table(VALUE tt, gsl_integration_qawo_table **t);
static VALUE rb_gsl_integration_qawo(int argc, VALUE *argv, VALUE obj)
{
double a, epsabs, epsrel;
double result, abserr;
size_t limit;
gsl_function *F = NULL;
gsl_integration_workspace *w = NULL;
gsl_integration_qawo_table *t = NULL;
int status, intervals, itmp, flag = 0, flagt = 0;
switch (TYPE(obj)) {
case T_MODULE: case T_CLASS: case T_OBJECT:
if (argc < 2) rb_raise(rb_eArgError, "too few arguments");
CHECK_FUNCTION(argv[0]);
Data_Get_Struct(argv[0], gsl_function, F);
itmp = 1;
break;
default:
if (argc < 1) rb_raise(rb_eArgError, "too few arguments");
Data_Get_Struct(obj, gsl_function, F);
itmp = 0;
break;
}
Need_Float(argv[itmp]);
a = NUM2DBL(argv[itmp]);
flagt = get_qawo_table(argv[argc-1], &t);
flag = get_epsabs_epsrel_limit_workspace(argc-1, argv, itmp+1, &epsabs, &epsrel,
&limit, &w);
status = gsl_integration_qawo(F, a, epsabs, epsrel, limit, w, t, &result, &abserr);
intervals = w->size;
if (flag == 1) gsl_integration_workspace_free(w);
if (flagt == 1) gsl_integration_qawo_table_free(t);
return rb_ary_new3(4, rb_float_new(result), rb_float_new(abserr), INT2FIX(intervals),
INT2FIX(status));
}
static int get_qawo_table(VALUE tt,
gsl_integration_qawo_table **t)
{
int flagt;
if (TYPE(tt) == T_ARRAY) {
flagt = 1;
*t = make_qawo_table(tt);
} else {
flagt = 0;
if (!rb_obj_is_kind_of(tt, cgsl_integration_qawo_table))
rb_raise(rb_eTypeError, "Integration::QAWO_Table expected");
Data_Get_Struct(tt, gsl_integration_qawo_table, *t);
}
return flagt;
}
static VALUE rb_gsl_integration_qawf(int argc, VALUE *argv, VALUE obj)
{
double a, epsabs = EPSREL_DEFAULT;
double result, abserr;
size_t limit = LIMIT_DEFAULT;
gsl_function *F = NULL;
gsl_integration_workspace *w = NULL, *cw = NULL;
gsl_integration_qawo_table *t = NULL;
int status, intervals, flag = 0, flagt = 0, itmp;
VALUE *vtmp;
switch (TYPE(obj)) {
case T_MODULE: case T_CLASS: case T_OBJECT:
if (argc < 2) rb_raise(rb_eArgError, "too few arguments");
CHECK_FUNCTION(argv[0]);
Data_Get_Struct(argv[0], gsl_function, F);
itmp = 1;
break;
default:
if (argc < 1) rb_raise(rb_eArgError, "too few arguments");
Data_Get_Struct(obj, gsl_function, F);
itmp = 0;
break;
}
Need_Float(argv[itmp]);
a = NUM2DBL(argv[itmp]);
itmp += 1;
if (TYPE(argv[itmp]) == T_FLOAT) {
epsabs = NUM2DBL(argv[itmp]);
itmp += 1;
}
vtmp = argv + itmp;
flagt = get_qawo_table(argv[argc-1], &t);
switch (argc - 1 - itmp) {
case 0:
w = gsl_integration_workspace_alloc(limit);
cw = gsl_integration_workspace_alloc(limit);
flag = 1;
break;
case 1:
CHECK_FIXNUM(vtmp[0]);
limit = FIX2INT(vtmp[0]);
w = gsl_integration_workspace_alloc(limit);
cw = gsl_integration_workspace_alloc(limit);
flag = 1;
break;
case 2:
CHECK_WORKSPACE(vtmp[0]); CHECK_WORKSPACE(vtmp[1]);
Data_Get_Struct(vtmp[0], gsl_integration_workspace, w);
Data_Get_Struct(vtmp[1], gsl_integration_workspace, cw);
flag = 0;
break;
case 3:
CHECK_FIXNUM(vtmp[0]);
CHECK_WORKSPACE(vtmp[1]); CHECK_WORKSPACE(vtmp[2]);
limit = FIX2INT(vtmp[0]);
Data_Get_Struct(vtmp[1], gsl_integration_workspace, w);
Data_Get_Struct(vtmp[2], gsl_integration_workspace, cw);
flag = 0;
break;
default:
rb_raise(rb_eArgError, "wrong number of arguments");
break;
}
status = gsl_integration_qawf(F, a, epsabs, limit, w, cw, t, &result, &abserr);
intervals = w->size;
if (flag == 1) {
gsl_integration_workspace_free(w);
gsl_integration_workspace_free(cw);
}
if (flagt == 1) gsl_integration_qawo_table_free(t);
return rb_ary_new3(4, rb_float_new(result), rb_float_new(abserr),
INT2FIX(intervals), INT2FIX(status));
}
static void rb_gsl_integration_define_symbols(VALUE module)
{
rb_define_const(module, "GAUSS15", INT2FIX(GSL_INTEG_GAUSS15));
rb_define_const(module, "GAUSS21", INT2FIX(GSL_INTEG_GAUSS21));
rb_define_const(module, "GAUSS31", INT2FIX(GSL_INTEG_GAUSS31));
rb_define_const(module, "GAUSS41", INT2FIX(GSL_INTEG_GAUSS41));
rb_define_const(module, "GAUSS51", INT2FIX(GSL_INTEG_GAUSS51));
rb_define_const(module, "GAUSS61", INT2FIX(GSL_INTEG_GAUSS61));
rb_define_const(module, "COSINE", INT2FIX(GSL_INTEG_COSINE));
rb_define_const(module, "SINE", INT2FIX(GSL_INTEG_SINE));
}
static VALUE rb_gsl_integration_workspace_alloc(int argc, VALUE *argv,
VALUE klass)
{
size_t limit;
if (argc == 1) limit = FIX2INT(argv[0]);
else limit = LIMIT_DEFAULT;
return Data_Wrap_Struct(klass, 0,
gsl_integration_workspace_free,
gsl_integration_workspace_alloc(limit));
}
static VALUE rb_gsl_integration_workspace_limit(VALUE obj)
{
gsl_integration_workspace *w = NULL;
Data_Get_Struct(obj, gsl_integration_workspace, w);
return INT2FIX(w->limit);
}
static VALUE rb_gsl_integration_workspace_size(VALUE obj)
{
gsl_integration_workspace *w = NULL;
Data_Get_Struct(obj, gsl_integration_workspace, w);
return INT2FIX(w->size);
}
static VALUE rb_gsl_integration_workspace_nrmax(VALUE obj)
{
gsl_integration_workspace *w = NULL;
Data_Get_Struct(obj, gsl_integration_workspace, w);
return INT2FIX(w->nrmax);
}
static VALUE rb_gsl_integration_workspace_i(VALUE obj)
{
gsl_integration_workspace *w = NULL;
Data_Get_Struct(obj, gsl_integration_workspace, w);
return INT2FIX(w->i);
}
static VALUE rb_gsl_integration_workspace_maximum_level(VALUE obj)
{
gsl_integration_workspace *w = NULL;
Data_Get_Struct(obj, gsl_integration_workspace, w);
return INT2FIX(w->maximum_level);
}
static VALUE rb_gsl_integration_workspace_to_a(VALUE obj)
{
gsl_integration_workspace *w = NULL;
Data_Get_Struct(obj, gsl_integration_workspace, w);
return rb_ary_new3(5, INT2FIX(w->limit), INT2FIX(w->size), INT2FIX(w->nrmax),
INT2FIX(w->i), INT2FIX(w->maximum_level));
}
static VALUE rb_gsl_integration_workspace_alist(VALUE obj)
{
gsl_integration_workspace *w = NULL;
gsl_vector_view *v = NULL;
Data_Get_Struct(obj, gsl_integration_workspace, w);
v = rb_gsl_make_vector_view(w->alist, w->limit, 1);
return Data_Wrap_Struct(cgsl_vector_view_ro, 0, free, v);
}
static VALUE rb_gsl_integration_workspace_blist(VALUE obj)
{
gsl_integration_workspace *w = NULL;
gsl_vector_view *v = NULL;
Data_Get_Struct(obj, gsl_integration_workspace, w);
v = rb_gsl_make_vector_view(w->blist, w->limit, 1);
return Data_Wrap_Struct(cgsl_vector_view_ro, 0, free, v);
}
static VALUE rb_gsl_integration_workspace_rlist(VALUE obj)
{
gsl_integration_workspace *w = NULL;
gsl_vector_view *v = NULL;
Data_Get_Struct(obj, gsl_integration_workspace, w);
v = rb_gsl_make_vector_view(w->rlist, w->limit, 1);
return Data_Wrap_Struct(cgsl_vector_view_ro, 0, free, v);
}
static VALUE rb_gsl_integration_workspace_elist(VALUE obj)
{
gsl_integration_workspace *w = NULL;
gsl_vector_view *v = NULL;
Data_Get_Struct(obj, gsl_integration_workspace, w);
v = rb_gsl_make_vector_view(w->elist, w->limit, 1);
return Data_Wrap_Struct(cgsl_vector_view_ro, 0, free, v);
}
void Init_gsl_integration(VALUE module)
{
VALUE mgsl_integ;
mgsl_integ = rb_define_module_under(module, "Integration");
rb_gsl_integration_define_symbols(mgsl_integ);
rb_define_method(cgsl_function, "integration_qng", rb_gsl_integration_qng, -1);
rb_define_method(cgsl_function, "integration_qag", rb_gsl_integration_qag, -1);
rb_define_method(cgsl_function, "integration_qags", rb_gsl_integration_qags, -1);
rb_define_method(cgsl_function, "integration_qagp", rb_gsl_integration_qagp, -1);
rb_define_method(cgsl_function, "integration_qagi", rb_gsl_integration_qagi, -1);
rb_define_method(cgsl_function, "integration_qagiu", rb_gsl_integration_qagiu, -1);
rb_define_method(cgsl_function, "integration_qagil", rb_gsl_integration_qagil, -1);
rb_define_method(cgsl_function, "integration_qawc", rb_gsl_integration_qawc, -1);
rb_define_alias(cgsl_function, "qng", "integration_qng");
rb_define_alias(cgsl_function, "qag", "integration_qag");
rb_define_alias(cgsl_function, "qags", "integration_qags");
rb_define_alias(cgsl_function, "qagp", "integration_qagp");
rb_define_alias(cgsl_function, "qagi", "integration_qagi");
rb_define_alias(cgsl_function, "qagiu", "integration_qagiu");
rb_define_alias(cgsl_function, "qagil", "integration_qagil");
rb_define_alias(cgsl_function, "qawc", "integration_qawc");
cgsl_integration_qaws_table = rb_define_class_under(mgsl_integ, "QAWS_Table",
cGSL_Object);
rb_define_singleton_method(cgsl_integration_qaws_table, "alloc",
rb_gsl_integration_qaws_table_alloc, -1);
/* rb_define_singleton_method(cgsl_integration_qaws_table, "new",
rb_gsl_integration_qaws_table_alloc, -1);*/
rb_define_method(cgsl_integration_qaws_table, "to_a",
rb_gsl_integration_qaws_table_to_a, 0);
rb_define_method(cgsl_integration_qaws_table, "set",
rb_gsl_integration_qaws_table_set, -1);
rb_define_method(rb_cArray, "to_gsl_integration_qaws_table",
rb_gsl_ary_to_integration_qaws_table, 0);
rb_define_alias(rb_cArray, "to_qaws_table", "to_gsl_integration_qaws_table");
rb_define_method(cgsl_function, "integration_qaws", rb_gsl_integration_qaws, -1);
rb_define_alias(cgsl_function, "qaws", "integration_qaws");
cgsl_integration_qawo_table = rb_define_class_under(mgsl_integ, "QAWO_Table",
cGSL_Object);
rb_define_singleton_method(cgsl_integration_qawo_table, "alloc",
rb_gsl_integration_qawo_table_alloc, -1);
/* rb_define_singleton_method(cgsl_integration_qawo_table, "new",
rb_gsl_integration_qawo_table_alloc, -1);*/
rb_define_method(cgsl_integration_qawo_table, "to_a",
rb_gsl_integration_qawo_table_to_a, 0);
rb_define_method(rb_cArray, "to_gsl_integration_qawo_table",
rb_gsl_ary_to_integration_qawo_table, 0);
rb_define_method(cgsl_integration_qawo_table, "set",
rb_gsl_integration_qawo_table_set, -1);
rb_define_method(cgsl_integration_qawo_table, "set_length",
rb_gsl_integration_qawo_table_set_length, 1);
rb_define_method(cgsl_function, "integration_qawo", rb_gsl_integration_qawo, -1);
rb_define_method(cgsl_function, "integration_qawf", rb_gsl_integration_qawf, -1);
rb_define_alias(cgsl_function, "qawo", "integration_qawo");
rb_define_alias(cgsl_function, "qawf", "integration_qawf");
cgsl_integration_workspace = rb_define_class_under(mgsl_integ,
"Workspace", cGSL_Object);
/* rb_define_singleton_method(cgsl_integration_workspace, "new",
rb_gsl_integration_workspace_alloc, -1);*/
rb_define_singleton_method(cgsl_integration_workspace, "alloc",
rb_gsl_integration_workspace_alloc, -1);
rb_define_method(cgsl_integration_workspace, "limit",
rb_gsl_integration_workspace_limit, 0);
rb_define_method(cgsl_integration_workspace, "size",
rb_gsl_integration_workspace_size, 0);
rb_define_method(cgsl_integration_workspace, "nrmax",
rb_gsl_integration_workspace_nrmax, 0);
rb_define_method(cgsl_integration_workspace, "i",
rb_gsl_integration_workspace_i, 0);
rb_define_method(cgsl_integration_workspace, "maximum_level",
rb_gsl_integration_workspace_maximum_level, 0);
rb_define_method(cgsl_integration_workspace, "to_a",
rb_gsl_integration_workspace_to_a, 0);
rb_define_method(cgsl_integration_workspace, "alist",
rb_gsl_integration_workspace_alist, 0);
rb_define_method(cgsl_integration_workspace, "blist",
rb_gsl_integration_workspace_blist, 0);
rb_define_method(cgsl_integration_workspace, "rlist",
rb_gsl_integration_workspace_rlist, 0);
rb_define_method(cgsl_integration_workspace, "elist",
rb_gsl_integration_workspace_elist, 0);
/*****/
rb_define_singleton_method(mgsl_integ, "qng", rb_gsl_integration_qng, -1);
rb_define_singleton_method(mgsl_integ, "qag", rb_gsl_integration_qag, -1);
rb_define_singleton_method(mgsl_integ, "qags", rb_gsl_integration_qags, -1);
rb_define_singleton_method(mgsl_integ, "qagp", rb_gsl_integration_qagp, -1);
rb_define_singleton_method(mgsl_integ, "qagi", rb_gsl_integration_qagi, -1);
rb_define_singleton_method(mgsl_integ, "qagiu", rb_gsl_integration_qagiu, -1);
rb_define_singleton_method(mgsl_integ, "qagil", rb_gsl_integration_qagil, -1);
rb_define_singleton_method(mgsl_integ, "qawc", rb_gsl_integration_qawc, -1);
rb_define_singleton_method(mgsl_integ, "qaws", rb_gsl_integration_qaws, -1);
rb_define_singleton_method(mgsl_integ, "qawo", rb_gsl_integration_qawo, -1);
rb_define_singleton_method(mgsl_integ, "qawf", rb_gsl_integration_qawf, -1);
}
#undef EPSABS_DEFAULT
#undef EPSREL_DEFAULT
#undef LIMIT_DEFAULT
#undef KEY_DEFAULT
#ifdef CHECK_WORKSPACE
#undef CHECK_WORKSPACE
#endif
syntax highlighted by Code2HTML, v. 0.9.1