DEFUN_DLD(GSL_OCTAVE_NAME, args, nargout, "\
  -*- texinfo -*-\n\
@deftypefn {Loadable Function} {@var{y} =} GSL_OCTAVE_NAME (@var{n}, @var{x})\n\
@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} GSL_OCTAVE_NAME (@dots{})\n\
\n\
GSL_FUNC_DOCSTRING
\n\
@var{err} contains an estimate of the absolute error in the value @var{y}.\n\
\n\
This function is from the GNU Scientific Library,\n\
see @url{http://www.gnu.org/software/gsl/} for documentation.\n\
@end deftypefn\n\
")
{
    int i;
    dim_vector dv;
    
    gsl_set_error_handler (octave_gsl_errorhandler);
    
    if(args.length() != 2) {
	print_usage ();
	return octave_value();
    }
    if(!args(0).is_real_type() || !args(1).is_real_type()) {
        error("The arguments must be real.");
	print_usage ();	    
	return octave_value();
    }

    // Nice combinatorial explosion here
    NDArray n = args(0).array_value();
    NDArray x = args(1).array_value();    
    if(n.length() == x.length()) {
	dv = x.dims();
    	NDArray y(dv);
	int len = x.length();
	if(nargout < 2) {
	    for(i = 0; i < len; i++) {
		y.xelem(i) = GSL_FUNC_NAME (static_cast<int>(n.xelem(i)),
					    x.xelem(i));
	    }
	    return octave_value(y);	    
	} else {
	    NDArray err(dv);
	    gsl_sf_result result;
	    octave_value_list retval;
	    for(i = 0; i < len; i++) {
		GSL_FUNC_NAME_e (static_cast<int>(n.xelem(i)),
				 x.xelem(i), &result);
		y.xelem(i) = result.val;
		err.xelem(i) = result.err;
	    }
	    retval(1) = octave_value(err);
	    retval(0) = octave_value(y);
	    return retval;
	}
    } else if(n.length() == 1) {
	dv = x.dims();
    	NDArray y(dv);
	int len = x.length();
	int nint = static_cast<int>(n.xelem(0));
	if(nargout < 2) {
	    for(i = 0; i < len; i++) {
		y.xelem(i) = GSL_FUNC_NAME (nint,
					    x.xelem(i));
	    }
	    return octave_value(y);	    
	} else {
	    NDArray err(dv);
	    gsl_sf_result result;
	    octave_value_list retval;
	    for(i = 0; i < len; i++) {
		GSL_FUNC_NAME_e (nint, x.xelem(i), &result);
		y.xelem(i) = result.val;
		err.xelem(i) = result.err;
	    }
	    retval(1) = octave_value(err);
	    retval(0) = octave_value(y);
	    return retval;
	}
    } else if(x.length() == 1) {
	dv = n.dims();
    	NDArray y(dv);
	int len = n.length();
	double xdouble = x.xelem(0);
	if(nargout < 2) {
	    for(i = 0; i < len; i++) {
		y.xelem(i) = GSL_FUNC_NAME (static_cast<int>(n.xelem(i)),
					    xdouble);
	    }
	    return octave_value(y);	    
	} else {
	    NDArray err(dv);
	    gsl_sf_result result;
	    octave_value_list retval;
	    for(i = 0; i < len; i++) {
		GSL_FUNC_NAME_e (static_cast<int>(n.xelem(i)),
				 xdouble, &result);
		y.xelem(i) = result.val;
		err.xelem(i) = result.err;
	    }
	    retval(1) = octave_value(err);
	    retval(0) = octave_value(y);
	    return retval;
	}
    } else {
	error("First and second argument must either have the same size, or one of them must be scalar.");
	print_usage ();	
    }

    return octave_value();

}


/*
;;; Local Variables: ***
;;; mode: C++ ***
;;; End: ***
*/