ext/poly_source.c in gsl-1.14.7 vs ext/poly_source.c in gsl-1.15.3
- old
+ new
@@ -464,13 +464,20 @@
default:
rb_raise(rb_eArgError,
"wrong number of arguments (3 numbers or 1 array or 1 vector)");
break;
}
- r = gsl_vector_alloc(2);
- gsl_vector_set(r, 0, x0);
- gsl_vector_set(r, 1, x1);
+ // If n == 0, we want to return an empty gsl_vector, but gsl_vectors can'y be
+ // empty, so just return an empty Array.
+ if(n == 0) {
+ return rb_ary_new();
+ }
+ r = gsl_vector_alloc(n);
+ switch(n) {
+ case 2: gsl_vector_set(r, 1, x1); /* fall through */
+ case 1: gsl_vector_set(r, 0, x0);
+ }
return Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, r);
}
static VALUE FUNCTION(rb_gsl_poly,complex_solve_quadratic)(int argc, VALUE *argv, VALUE obj)
{
@@ -505,13 +512,20 @@
default:
rb_raise(rb_eArgError,
"wrong number of arguments (3 numbers or 1 array or 1 vector)");
break;
}
- r = gsl_vector_complex_alloc(2);
- gsl_vector_complex_set(r, 0, z0);
- gsl_vector_complex_set(r, 1, z1);
+ // If n == 0, we want to return an empty gsl_vector, but gsl_vectors can'y be
+ // empty, so just return an empty Array.
+ if(n == 0) {
+ return rb_ary_new();
+ }
+ r = gsl_vector_complex_alloc(n);
+ switch(n) {
+ case 2: gsl_vector_complex_set(r, 1, z1); /* fall through */
+ case 1: gsl_vector_complex_set(r, 0, z0);
+ }
return Data_Wrap_Struct(cgsl_vector_complex, 0, gsl_vector_complex_free, r);
}
static VALUE FUNCTION(rb_gsl_poly,solve_cubic)(int argc, VALUE *argv, VALUE obj)
{
@@ -545,14 +559,16 @@
default:
rb_raise(rb_eArgError,
"wrong number of arguments (3 numbers or 1 array or 1 vector)");
break;
}
- r = gsl_vector_alloc(3);
- gsl_vector_set(r, 0, x0);
- gsl_vector_set(r, 1, x1);
- gsl_vector_set(r, 2, x2);
+ r = gsl_vector_alloc(n);
+ switch(n) {
+ case 3: gsl_vector_set(r, 2, x2); /* fall through */
+ case 2: gsl_vector_set(r, 1, x1); /* fall through */
+ case 1: gsl_vector_set(r, 0, x0);
+ }
return Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, r);
}
static VALUE FUNCTION(rb_gsl_poly,complex_solve_cubic)(int argc, VALUE *argv, VALUE obj)
{
@@ -586,14 +602,16 @@
default:
rb_raise(rb_eArgError,
"wrong number of arguments (3 numbers or 1 array or 1 vector)");
break;
}
- r = gsl_vector_complex_alloc(3);
- gsl_vector_complex_set(r, 0, z0);
- gsl_vector_complex_set(r, 1, z1);
- gsl_vector_complex_set(r, 2, z2);
+ r = gsl_vector_complex_alloc(n);
+ switch(n) {
+ case 3: gsl_vector_complex_set(r, 2, z2); /* fall through */
+ case 2: gsl_vector_complex_set(r, 1, z1); /* fall through */
+ case 1: gsl_vector_complex_set(r, 0, z0);
+ }
return Data_Wrap_Struct(cgsl_vector_complex, 0, gsl_vector_complex_free, r);
}
#ifdef HAVE_POLY_SOLVE_QUARTIC
static VALUE FUNCTION(rb_gsl_poly,solve_quartic)(int argc, VALUE *argv, VALUE obj)
@@ -631,15 +649,17 @@
default:
rb_raise(rb_eArgError,
"wrong number of arguments (3 numbers or 1 array or 1 vector)");
break;
}
- r = gsl_vector_alloc(4);
- gsl_vector_set(r, 0, x0);
- gsl_vector_set(r, 1, x1);
- gsl_vector_set(r, 2, x2);
- gsl_vector_set(r, 3, x3);
+ r = gsl_vector_alloc(n);
+ switch(n) {
+ case 4: gsl_vector_set(r, 3, x3); /* fall through */
+ case 3: gsl_vector_set(r, 2, x2); /* fall through */
+ case 2: gsl_vector_set(r, 1, x1); /* fall through */
+ case 1: gsl_vector_set(r, 0, x0);
+ }
return Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, r);
}
static VALUE FUNCTION(rb_gsl_poly,complex_solve_quartic)(int argc, VALUE *argv, VALUE obj)
{
@@ -676,15 +696,17 @@
default:
rb_raise(rb_eArgError,
"wrong number of arguments (3 numbers or 1 array or 1 vector)");
break;
}
- r = gsl_vector_complex_alloc(4);
- gsl_vector_complex_set(r, 0, z0);
- gsl_vector_complex_set(r, 1, z1);
- gsl_vector_complex_set(r, 2, z2);
- gsl_vector_complex_set(r, 3, z3);
+ r = gsl_vector_complex_alloc(n);
+ switch(n) {
+ case 4: gsl_vector_complex_set(r, 0, z0); /* fall through */
+ case 3: gsl_vector_complex_set(r, 1, z1); /* fall through */
+ case 2: gsl_vector_complex_set(r, 2, z2); /* fall through */
+ case 1: gsl_vector_complex_set(r, 3, z3);
+ }
return Data_Wrap_Struct(cgsl_vector_complex, 0, gsl_vector_complex_free, r);
}
#endif
/* singleton method */
@@ -694,11 +716,13 @@
gsl_poly_complex_workspace *w = NULL;
GSL_TYPE(gsl_vector) *v = NULL;
gsl_vector *a = NULL, *z = NULL;
gsl_complex c;
gsl_vector_complex *r = NULL;
- int status, i, flag = 0;
+ int i, flag = 0;
+ // local variable "status" declared and set, but never used
+ //int status;
switch (argc) {
case 1:
break;
case 2:
@@ -745,11 +769,11 @@
flag = 0;
} else {
w = gsl_poly_complex_workspace_alloc(size);
flag = 1;
}
- status = gsl_poly_complex_solve(a->data, size, w, z->data);
+ /*status =*/ gsl_poly_complex_solve(a->data, size, w, z->data);
if (flag == 1) gsl_poly_complex_workspace_free(w);
gsl_vector_free(a);
r = gsl_vector_complex_alloc(size - 1);
for (i = 0; i < size - 1; i++) {
@@ -781,19 +805,23 @@
a1 = FUNCTION(gsl_vector,get)(v, 1);
a0 = FUNCTION(gsl_vector,get)(v, 0);
d = a1*a1 - 4.0*a2*a0; /* determinant */
if (d >= 0.0) {
n = gsl_poly_solve_quadratic(a2, a1, a0, &x0, &x1);
- r = gsl_vector_alloc(2);
- gsl_vector_set(r, 0, x0);
- gsl_vector_set(r, 1, x1);
+ r = gsl_vector_alloc(n);
+ switch(n) {
+ case 2: gsl_vector_set(r, 1, x1); /* fall through */
+ case 1: gsl_vector_set(r, 0, x0);
+ }
return Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, r);
} else {
n = gsl_poly_complex_solve_quadratic(a2, a1, a0, &z0, &z1);
- r2 = gsl_vector_complex_alloc(2);
- gsl_vector_complex_set(r2, 0, z0);
- gsl_vector_complex_set(r2, 1, z1);
+ r2 = gsl_vector_complex_alloc(n);
+ switch(n) {
+ case 2: gsl_vector_complex_set(r2, 1, z1); /* fall through */
+ case 1: gsl_vector_complex_set(r2, 0, z0);
+ }
return Data_Wrap_Struct(cgsl_vector_complex, 0, gsl_vector_complex_free, r2);
}
}
@@ -801,25 +829,25 @@
static VALUE FUNCTION(rb_gsl_poly,complex_solve_quadratic2)(VALUE obj)
{
GSL_TYPE(gsl_vector) *v = NULL;
gsl_vector_complex *r = NULL;
double a2, a1, a0;
- double d; /* determinant */
int n;
gsl_complex z0, z1;
Data_Get_Struct(obj, GSL_TYPE(gsl_vector), v);
if (v->size < 3) {
rb_raise(rb_eArgError, "the order of the object is less than 3.");
}
a2 = FUNCTION(gsl_vector,get)(v, 2); /* coefficients */
a1 = FUNCTION(gsl_vector,get)(v, 1);
a0 = FUNCTION(gsl_vector,get)(v, 0);
- d = a1*a1 - 4.0*a2*a0; /* determinant */
n = gsl_poly_complex_solve_quadratic(a2, a1, a0, &z0, &z1);
- r = gsl_vector_complex_alloc(2);
- gsl_vector_complex_set(r, 0, z0);
- gsl_vector_complex_set(r, 1, z1);
+ r = gsl_vector_complex_alloc(n);
+ switch(n) {
+ case 2: gsl_vector_complex_set(r, 1, z1); /* fall through */
+ case 1: gsl_vector_complex_set(r, 0, z0);
+ }
return Data_Wrap_Struct(cgsl_vector_complex, 0, gsl_vector_complex_free, r);
}
/* x**3 + a2 x**2 + a1 x + a0 = 0 */
static VALUE FUNCTION(rb_gsl_poly,solve_cubic2)(VALUE obj)
@@ -837,14 +865,16 @@
a3 = FUNCTION(gsl_vector,get)(v, 3);
a2 = FUNCTION(gsl_vector,get)(v, 2)/a3; /* coefficients */
a1 = FUNCTION(gsl_vector,get)(v, 1)/a3;
a0 = FUNCTION(gsl_vector,get)(v, 0)/a3;
n = gsl_poly_solve_cubic(a2, a1, a0, &x0, &x1, &x2);
- r = gsl_vector_alloc(3);
- gsl_vector_set(r, 0, x0);
- gsl_vector_set(r, 1, x1);
- gsl_vector_set(r, 2, x2);
+ r = gsl_vector_alloc(n);
+ switch(n) {
+ case 3: gsl_vector_set(r, 2, x2); /* fall through */
+ case 2: gsl_vector_set(r, 1, x1); /* fall through */
+ case 1: gsl_vector_set(r, 0, x0);
+ }
return Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, r);
}
static VALUE FUNCTION(rb_gsl_poly,complex_solve_cubic2)(VALUE obj)
{
@@ -860,14 +890,16 @@
a3 = FUNCTION(gsl_vector,get)(v, 3);
a2 = FUNCTION(gsl_vector,get)(v, 2)/a3; /* coefficients */
a1 = FUNCTION(gsl_vector,get)(v, 1)/a3;
a0 = FUNCTION(gsl_vector,get)(v, 0)/a3;
n = gsl_poly_complex_solve_cubic(a2, a1, a0, &z0, &z1, &z2);
- r = gsl_vector_complex_alloc(3);
- gsl_vector_complex_set(r, 0, z0);
- gsl_vector_complex_set(r, 1, z1);
- gsl_vector_complex_set(r, 2, z2);
+ r = gsl_vector_complex_alloc(n);
+ switch(n) {
+ case 3: gsl_vector_complex_set(r, 2, z2); /* fall through */
+ case 2: gsl_vector_complex_set(r, 1, z1); /* fall through */
+ case 1: gsl_vector_complex_set(r, 0, z0);
+ }
return Data_Wrap_Struct(cgsl_vector_complex, 0, gsl_vector_complex_free, r);
}
#ifdef HAVE_POLY_SOLVE_QUARTIC
/* a4 x**4 + a3 x**3 + a2 x**2 + a1 x + a0 = 0 */
@@ -932,11 +964,13 @@
gsl_vector *a = NULL, *z = NULL;
size_t size, size2, i;
gsl_poly_complex_workspace *w = NULL;
gsl_complex c;
gsl_vector_complex *r = NULL;
- int n, flag = 0;
+ int flag = 0;
+ // local variable "status" declared and set, but never used
+ //int status;
Data_Get_Struct(obj, GSL_TYPE(gsl_vector), v);
size = v->size;
size2 = 2*(size - 1);
z = gsl_vector_alloc(size2);
@@ -951,11 +985,11 @@
} else {
w = gsl_poly_complex_workspace_alloc(size);
flag = 1;
}
- n = gsl_poly_complex_solve(a->data, size, w, z->data);
+ /*status =*/ gsl_poly_complex_solve(a->data, size, w, z->data);
r = gsl_vector_complex_alloc(size - 1);
for (i = 0; i < size - 1; i++) {
c.dat[0] = gsl_vector_get(z, i*2);
c.dat[1] = gsl_vector_get(z, i*2+1);
@@ -1582,20 +1616,16 @@
static VALUE rb_gsl_poly_wfit(int argc, VALUE *argv, VALUE obj)
{
gsl_multifit_linear_workspace *space = NULL;
gsl_matrix *X = NULL, *cov = NULL;
gsl_vector *x, *y = NULL, *w, *c = NULL;
- gsl_vector_view xx, yy, ww;
size_t order, i, j;
double chisq, val;
int status, flag = 0;
VALUE vc, vcov;
if (argc != 4 && argc != 5)
rb_raise(rb_eArgError, "wrong number of arguments (%d for 4 or 5)", argc);
- x = &xx.vector;
- w = &ww.vector;
- y = &yy.vector;
Data_Get_Vector(argv[0], x);
Data_Get_Vector(argv[1], w);
Data_Get_Vector(argv[2], y);
order = NUM2INT(argv[3]);
if (argc == 5) {
@@ -1613,10 +1643,10 @@
for (j = 1; j <= order; j++) {
val *= gsl_vector_get(x, i);
gsl_matrix_set(X, i, j, val);
}
}
- status = gsl_multifit_linear(X, y, c, cov, &chisq, space);
+ status = gsl_multifit_wlinear(X, w, y, c, cov, &chisq, space);
if (flag == 1) gsl_multifit_linear_free(space);
vc = Data_Wrap_Struct(cgsl_poly, 0, gsl_vector_free, c);
vcov = Data_Wrap_Struct(cgsl_matrix, 0, gsl_matrix_free, cov);
gsl_matrix_free(X);
return rb_ary_new3(4, vc, vcov, rb_float_new(chisq), INT2FIX(status));