Science the GNU Way, Part I
In my past several articles, I've looked at various packages to do all kinds of science. Sometimes, however, there just isn't a tool to solve a particular problem. That's the great thing about science. There is always something new to discover and study. But, this means it's up to you to develop the software tools you need to do your analysis. This article takes a look at the GNU Scientific Library, or GSL. This library is the Swiss Army library of routines that you will find useful in your work.
First, you need to to get a copy of GSL and install it on your system. Because it is part of the GNU Project, it is hosted at http://www.gnu.org/s/gsl. You always can download and build from the source code, but all major distributions should have packages available. For example, on Debianbased systems, you need to install the package libgsl0dev to develop your code and gslbin to run that code. GSL is meant for C and C++, so you also need a compiler. Most of you probably already are familiar with GCC, so I stick with that here.
The next step is actually the last step. I'm looking at compiling now so that I can focus on all the tools available in GSL. All the header files for GSL are stored in a subdirectory named gsl. So, for example, if you wanted to include the math header file, you would use:
#include <gsl/gsl_math.h>
All the functions are stored in a single library file called libgsl.a or libgsl.so. You also need a library to handle basic linear algebra. GSL provides one (in the file libgslcblas.so), but you can use your own. Additionally, you need to link in the math library. So, the final compile and link command should look like this: Garrick, one line below.
gcc o hello_world hello_world.c lgsl lgslcblas lm
There are optional inline versions for some of the performancecritical
functions in GSL. To use these, you need to include
DHAVE_INLINE
with
your compile command. To try to help with portability issues, GSL offers
some functions that exist only on certain platforms (and not others). As an
example, the BSD math library has a function called
hypot
. GSL offers
its own version, called gsl_hypot
, that you can use on nonBSD
platforms. Some functions have both a general algorithm as well
as optimized versions for specific platforms. This way, if you are running
on a SPARC, for example, you can select a version optimized for SPARC if it
exists.
One of the first things you likely will want to do is check
whether you are getting correct results from your code or if there were
errors. GSL has a number of functions and data structures available in the
header file named gsl_errno.h. Functions return a value of zero if everything
is fine. If there were any problems in trying to complete the requested
action, a nonzero value is returned. This could be an actual error
condition, like a wrong data type or memory error, or it could be a
condition like not being able to converge to within the requested accuracy
in the function call. This is why you always need to check the return value
for all GSL function calls. The actual values returned in an error
condition are error codes, defined in the file gsl_errno.h. They are
defined as macros that start with GSL_
. Examples
include the following:

GSL_EDOM
— domain error, used by functions when an argument doesn't fall into the domain over which the function is defined. 
GSL_ERANGE
— range error, either an overflow or underflow. 
GSL_ENOMEM
— no memory available.
The library will use values only up to 1024. Values above this are available
for use in your own code. There also are string versions of these error
codes available. You can translate the error code to its text value with
the function gsl_errno()
.
Now that you know how to compile your program and what to do with errors,
let's start looking at what kind of work you can do with GSL. Basic
mathematical functions are defined in the file gsl_math.h. The set of
mathematical constants from the BSD math library are provided by this part
of GSL. All of the constants start with M_
. Here are
a few of them:

M_PI
— pi. 
M_SQRT2
— the square root of 2. 
M_EULER
— Euler's constant.
There also are capabilities for dealing with infinities and nonnumbers. Three macros define the values themselves:

GSL_POSINF
— positive infinity. 
GSL_NEGINF
— negative infinity. 
GSL_NAN
— not a number.
There also are functions to test variables:

gsl_isnan
— is it not a number? 
gsl_isinf
— is it infinite? 
gsl_finite
— is it finite?
There is a macro to find the sign of a number.
GSL_SIGN(x)
returns the
sign of x: 1 if it is positive and –1 if it is negative. If you are
interested in seeing whether a number is even or odd, two macros
are defined: GSL_IS_ODD(x)
and
GSL_IS_EVEN(x)
. These return 1 if the
condition is true and 0 if it is not.
A series of elementary functions are part of the BSD math library. GSL provides versions of these for platforms that don't have native versions, including items like:

gsl_hypot
— calculate hypotenuse. 
gsl_asinh
,gsl_acosh
,gsl_atanh
— the arc hyperbolic trig functions.
If you are calculating the power of a number, you would use
gsl_pow_int(x,n)
, which gives you x to the power of n. There are
specific versions for powers less than 10. So if you wanted to find the
cube of a number, you would use gsl_pow_3
. These are very efficient and
highly optimized. You even can inline these specialized functions when
HAVE_INLINE
is defined.
Several macros are defined to help you find
the maximum or minimum of numbers, based on data type. The basic
GSL_MAX(a,b)
and GSL_MIN(a,b)
simply return either the maximum or
minimum of the two numbers a and b. GSL_MAX_DBL
and
GSL_MIN_DBL
find
the maximum and minimum of two doubles using an inline function.
GSL_MAX_INT
and GSL_MIN_INT
do the same for integer arguments.
When you
do any kind of numerical calculation on a computer, errors
always are introduced by roundoff and truncation. This is because you can't exactly
reproduce numbers on a finite binary system. But, what if you want to
compare two numbers and see whether they are approximately the same? GSL
provides the function gsl_fcmp(x,y,epsilon)
. This function compares the
two doubles x and y, and checks to see if they are within epsilon of each
other. If they are within this range, the function returns 0. If x < y, it
returns –1, and it returns 1 if x > y.
Complex numbers are used in many scientific fields. Within GSL, complex data
types are defined in the header file gsl_complex.h, and relevant functions
are defined in gsl_complex_math.h. To store complex numbers, the data
type gsl_complex
is defined. This is a struct that stores the two
portions. You can set the values with the functions
gsl_complex_rect(x,y)
or gsl_complex_polar(x,y)
. In the first, this
represents x+iy; whereas
in the second, x is the radius, and y is the angle in a polar representation.
You can pull out the real and imaginary parts of a complex number with the
macros GSL_REAL
and GSL_IMAG
. There is a function available to find the
absolute value of a complex number,
gsl_complex_abs(x)
, where x is of
type gsl_complex. Because complex numbers actually are built up of two
parts, even basic arithmetic is not simple. To do basic math, you can use
the following:

gsl_complex_add(a,b)

gsl_complex_sub(a,b)

gsl_complex_mul(a,b)

gsl_complex_div(a,b)
You can calculate the conjugate with
gsl_complex_conjugate(a)
and the
inverse with gsl_complex_inverse(a)
.
There are functions for basic
mathematical functions. To calculate the square root, you would use
gsl_complex_sqrt(x)
. To calculate the logarithm, you would use
gsl_complex_log(x)
. Several others are available too.
Trigonometric functions are provided, like
gsl_complex_sin(x)
. There also
are functions for hyperbolic trigonometric functions, along with the relevant
inverse functions.
Now that you have the basics down, my next article will explore all the actual scientific calculations you can do. I'll look at statistics, linear algebra, random numbers and many other topics.
Science image via Shutterstock.com.