gsl-design.texi

math library from gnu

best reference books in the field, rather than lesser known text-books
or introductory books which happen to be available (e.g. from
undergraduate studies).  For example, references concerning algorithms
should be to Knuth, references concerning statistics should be to
Kendall & Stuart, references concerning special functions should be to
Abramowitz & Stegun (Handbook of Mathematical Functions AMS-55), etc.
Whereever possible refer to Abramowitz & Stegun rather than other
reference books because it is a public domain work, so it is
inexpensive and freely redistributable.

The standard references have a better chance of being available in an
accessible library for the user.  If they are not available and the user
decides to buy a copy in order to look up the reference then this also
gives them the best quality book which should also cover the largest
number of other references in the GSL Manual.  If many different books
were to be referenced this would be an expensive and inefficient use of
resources for a user who needs to look up the details of the algorithms.
Reference books also stay in print much longer than text books, which
are often out-of-print after a few years.

Similarly, cite original papers wherever possible.  Be sure to keep
copies of these for your own reference (e.g. when dealing with bug
reports) or to pass on to future maintainers.

If you need help in tracking down references, ask on the
@code{gsl-discuss} mailing list.  There is a group of volunteers with
access to good libraries who have offered to help GSL developers get
copies of papers.

@c [JT section: written by James Theiler

@c And we furthermore promise to try as hard as possible to document
@c the software: this will ideally involve discussion of why you might want
@c to use it, what precisely it does, how precisely to invoke it, 
@c how more-or-less it works, and where we learned about the algorithm,
@c and (unless we wrote it from scratch) where we got the code.
@c We do not plan to write this entire package from scratch, but to cannibalize
@c existing mathematical freeware, just as we expect our own software to
@c be cannibalized.]

To write mathematics in the texinfo file you can use the @code{@@math}
command with @emph{simple} TeX commands. These are automatically
surrounded by @code{$...$} for math mode. For example,

@example
to calculate the coefficient @@math@{\alpha@} use the function...
@end example

@noindent
will be correctly formatted in both online and TeX versions of the
documentation.

Note that you cannot use the special characters @{ and @}
inside the @code{@@math} command because these conflict between TeX
and Texinfo.  This is a problem if you want to write something like
@code{\sqrt@{x+y@}}.  

To work around it you can preceed the math command with a special
macro @code{@@c} which contains the explicit TeX commands you want to
use (no restrictions), and put an ASCII approximation into the
@code{@@math} command (you can write @code{@@@{} and
@code{@@@}} there for the left and right braces).  The explicit TeX
commands are used in the TeX ouput and the argument of @code{@@math}
in the plain info output.  

Note that the @code{@@c@{@}} macro must go at the end of the
preceeding line, because everything else after it is ignored---as far
as texinfo is concerned it's actually a 'comment'. The comment
command @@c has been modified to capture a TeX expression which is
output by the next @@math command. For ordinary comments use the @@comment
command.

For example,

@example
this is a test @@c@{$\sqrt@{x+y@}$@}
@@math@{\sqrt@@@{x+y@@@}@}
@end example

@noindent
is equivalent to @code{this is a test $\sqrt@{x+y@}$} in plain TeX
and @code{this is a test @@math@{\sqrt@@@{x+y@@@}@}} in Info.  

It looks nicer if some of the more cryptic TeX commands are given
a C-style ascii version, e.g.

@example
for @@c@{$x \ge y$@}
@@math@{x >= y@}
@end example

@noindent
will be appropriately displayed in both TeX and Info.


@node Namespace, Header files, Documentation, Design
@section Namespace

Use @code{gsl_} as a prefix for all exported functions and variables.

Use @code{GSL_} as a prefix for all exported macros.

All exported header files should have a filename with the prefix @code{gsl_}.

All installed libraries should have a name like libgslhistogram.a

Any installed executables (utility programs etc) should have the prefix
@code{gsl-} (with a hyphen, not an underscore).

All function names, variables, etc should be in lower case.  Macros and
preprocessor variables should be in upper case.

Some common conventions in variable and function names:

@table @code
@item p1
plus 1, e.g. function @code{log1p(x)} or a variable like @code{kp1}, @math{=k+1}.

@item m1
minus 1, e.g. function @code{expm1(x)} or a variable like @code{km1}, @math{=k-1}.
@end table

@node Header files, Target system, Namespace, Design
@section Header files

Installed header files should be idempotent, i.e. surround them by the
preprocessor conditionals like the following,

@example
#ifndef __GSL_HISTOGRAM_H__
#define __GSL_HISTOGRAM_H__
...
#endif /* __GSL_HISTOGRAM_H__ */
@end example

@node Target system, Function Names, Header files, Design
@section Target system

The target system is ANSI C, with a full Standard C Library, and IEEE
arithmetic.

@node Function Names, Object-orientation, Target system, Design
@section Function Names

Each module has a name, which prefixes any function names in that
module, e.g. the module gsl_fft has function names like
gsl_fft_init. The modules correspond to subdirectories of the library
source tree.

@node Object-orientation, Comments, Function Names, Design
@section Object-orientation

The algorithms should be object oriented, but only to the extent that is
easy in portable ANSI C.  The use of casting or other tricks to simulate
inheritance is not desirable, and the user should not have to be aware
of anything like that.  This means many types of patterns are ruled
out.  However, this is not considered a problem -- they are too
complicated for the library.  

Note: it is possible to define an abstract base class easily in C, using
function pointers.  See the rng directory for an example. 

When reimplementing public domain fortran code, please try to introduce
the appropriate object concepts as structs, rather than translating the
code literally in terms of arrays.  The structs can be useful just
within the file, you don't need to export them to the user.

For example, if a fortran program repeatedly uses a subroutine like,

@example
SUBROUTINE  RESIZE (X, K, ND, K1)
@end example

@noindent
where X(K,D) represents a grid to be resized to X(K1,D) you can make
this more readable by introducing a struct,

@smallexample
struct grid @{
    int nd;  /* number of dimensions */
    int k;   /* number of bins */
    double * x;   /* partition of axes, array of size x[k][nd] */
@}

void
resize_grid (struct grid * g, int k_new)
@{
...
@}
@end smallexample

@noindent
Similarly, if you have a frequently recurring code fragment within a
single file you can define a static or static inline function for it.
This is typesafe and saves writing out everything in full.


@node Comments, Minimal structs, Object-orientation, Design
@section Comments

Follow the GNU Coding Standards.  A relevant quote is,

``Please write complete sentences and capitalize the first word.  If a
lower-case identifier comes at the beginning of a sentence, don't
capitalize it!  Changing the spelling makes it a different identifier.
If you don't like starting a sentence with a lower case letter, write
the sentence differently (e.g., "The identifier lower-case is ...").''

@node Minimal structs, Algorithm decomposition, Comments, Design
@section Minimal structs

We prefer to make structs which are @dfn{minimal}.  For example, if a
certain type of problem can be solved by several classes of algorithm
(e.g. with and without derivative information) it is better to make
separate types of struct to handle those cases.  i.e. run time type
identification is not desirable.

@node Algorithm decomposition, Memory allocation and ownership, Minimal structs, Design
@section Algorithm decomposition

Iterative algorithms should be decomposed into an INITIALIZE, ITERATE,
TEST form, so that the user can control the progress of the iteration
and print out intermediate results.  This is better than using
call-backs or using flags to control whether the function prints out
intermediate results.  In fact, call-backs should not be used -- if they
seem necessary then it's a sign that the algorithm should be broken down
further into individual components so that the user has complete control
over them.  

For example, when solving a differential equation the user may need to
be able to advance the solution by individual steps, while tracking a
realtime process.  This is only possible if the algorithm is broken down
into step-level components.  Higher level decompositions would not give
sufficient flexibility.

@node Memory allocation and ownership, Memory layout, Algorithm decomposition, Design
@section Memory allocation and ownership

Functions which allocate memory on the heap should end in _alloc
(e.g. gsl_foo_alloc) and be deallocated by a corresponding _free function
(gsl_foo_free).

Be sure to free any memory allocated by your function if you have to
return an error in a partially initialized object.

Don't allocate memory 'temporarily' inside a function and then free it
before the function returns.  This prevents the user from controlling
memory allocation.  All memory should be allocated and freed through
separate functions and passed around as a "workspace" argument.  This
allows memory allocation to be factored out of tight loops.

To avoid confusion over ownership, workspaces should not own each
other or contain other workspaces.  For clarity and ease of use in
different contexts, they should be allocated from integer arguments
rather than derived from other structs.

@node Memory layout, Linear Algebra Levels, Memory allocation and ownership, Design
@section Memory layout

We use flat blocks of memory to store matrices and vectors, not C-style
pointer-to-pointer arrays.  The matrices are stored in row-major order
-- i.e. the column index (second index) moves continuously through memory. 

@node Linear Algebra Levels, Error estimates, Memory layout, Design
@section Linear Algebra Levels

Functions using linear algebra are divided into two levels:

For purely "1d" functions we use the C-style arguments (double *,
stride, size) so that it is simpler to use the functions in a normal C
program, without needing to invoke all the gsl_vector machinery.

The philosophy here is to minimize the learning curve. If someone only
needs to use one function, like an fft, they can do so without having
to learn about gsl_vector.  

This leads to the question of why we don't do the same for matrices.
In that case the argument list gets too long and confusing, with
(size1, size2, tda) for each matrix and potential ambiguities over row
vs column ordering. In this case, it makes sense to use gsl_vector and
gsl_matrix, which take care of this for the user.

So really the library has two levels -- a lower level based on C types
for 1d operations, and a higher level based on gsl_matrix and
gsl_vector for general linear algebra.

Of course, it would be possible to define a vector version of the
lower level functions too. So far we have not done that because it was
not essential -- it could be done but it is easy enough to get by
using the C arguments, by typing v->data, v->stride, v->size instead.
A gsl_vector version of low-level functions would mainly be a
convenience.

Please use BLAS routines internally within the library whenever possible
for efficiency.

@node Error estimates, Exceptions and Error handling, Linear Algebra Levels, Design
@section Error estimates

In the special functions error bounds are given as twice the expected
``gaussian'' error.  i.e. 2-sigma, so the result is inside the error
98% of the time.  People expect the true value to be within +/- the
quoted error (this wouldn't be the case 32% of the time for 1 sigma).
Obviously the errors are not gaussian but a factor of two works well
in practice.

@node Exceptions and Error handling, Persistence, Error estimates, Design
@section Exceptions and Error handling

The basic error handling procedure is the return code (see gsl_errno.h
for a list of allowed values).  Use the GSL_ERROR macro to mark an
error.  The current definition of this macro is not ideal but it can be
changed at compile time. 

You should always use the GSL_ERROR macro to indicate an error, rather
than just returning an error code. The macro allows the user to trap
errors using the debugger (by setting a breakpoint on the function
gsl_error).  

The only circumstances where GSL_ERROR should not be used are where the
return value is "indicative" rather than an error -- for example, the
iterative routines use the return code to indicate the success or
failure of an iteration. By the nature of an iterative algorithm
"failure" (a return code of GSL_CONTINUE) is a normal occurrence and
there is no need to use GSL_ERROR there.

Be sure to free any memory allocated by your function if you return an
error (in particular for errors in partially initialized objects).