gsl-design.texi

math library from gnu

@node Persistence, Using Return Values, Exceptions and Error handling, Design
@section Persistence

If you make an object foo which uses blocks of memory (e.g. vector,
matrix, histogram) you can provide functions for reading and writing
those blocks,

@smallexample
int gsl_foo_fread (FILE * stream, gsl_foo * v);
int gsl_foo_fwrite (FILE * stream, const gsl_foo * v);
int gsl_foo_fscanf (FILE * stream, gsl_foo * v);
int gsl_foo_fprintf (FILE * stream, const gsl_foo * v, const char *format);
@end smallexample

@noindent
Only dump out the blocks of memory, not any associated parameters such
as lengths.  The idea is for the user to build higher level input/output
facilities using the functions the library provides.  The fprintf/fscanf
versions should be portable between architectures, while the binary
versions should be the "raw" version of the data. Use the functions

@smallexample
int gsl_block_fread (FILE * stream, gsl_block * b);
int gsl_block_fwrite (FILE * stream, const gsl_block * b);
int gsl_block_fscanf (FILE * stream, gsl_block * b);
int gsl_block_fprintf (FILE * stream, const gsl_block * b, const char *format);
@end smallexample

@noindent
or

@smallexample
int gsl_block_raw_fread (FILE * stream, double * b, size_t n, size_t stride);
int gsl_block_raw_fwrite (FILE * stream, const double * b, size_t n, size_t stri
de);
int gsl_block_raw_fscanf (FILE * stream, double * b, size_t n, size_t stride);
int gsl_block_raw_fprintf (FILE * stream, const double * b, size_t n, size_t str
ide, const char *format);
@end smallexample

@noindent
to do the actual reading and writing.

@node Using Return Values, Variable Names, Persistence, Design
@section Using Return Values

Always assign a return value to a variable before using it.  This allows
easier debugging of the function, and inspection and modification of the
return value.  If the variable is only needed temporarily then enclose
it in a suitable scope.

For example, instead of writing,

@example
a = f(g(h(x,y)))
@end example

@noindent
use temporary variables to store the intermediate values,
@example
@{
  double u = h(x,y);
  double v = g(u);
  a = f(v);
@}
@end example

@noindent
These can then be inspected more easily in the debugger, and breakpoints
can be placed more precisely.  The compiler will eliminate the temporary
variables automatically when the program is compiled with optimization.

@node Variable Names, Datatype widths, Using Return Values, Design
@section Variable Names

Try to follow existing conventions for variable names,

@table @code
@item dim
number of dimensions
@item w
pointer to workspace 
@item state
pointer to state variable (use @code{s} if you need to save characters)
@item result
pointer to result (output variable)
@item abserr
absolute error 
@item relerr
relative error
@item epsabs
absolute tolerance 
@item epsrel
relative tolerance
@item size
the size of an array or vector e.g. double array[size] 
@item stride
the stride of a vector
@item size1
the number of rows in a matrix
@item size2
the number of columns in a matrix
@item n
general integer number, e.g. number of elements of array, in fft, etc
@item r
random number generator (gsl_rng)
@end table

@node Datatype widths, size_t, Variable Names, Design
@section Datatype widths

Be aware that in ANSI C the type @code{int} is only guaranteed to
provide 16-bits. It may provide more, but is not guaranteed to.
Therefore if you require 32 bits you must use @code{long int}, which
will have 32 bits or more.  Of course, on many platforms the type
@code{int} does have 32 bits instead of 16 bits but we have to code to
the ANSI standard rather than a specific platform.

@node size_t, Arrays vs Pointers, Datatype widths, Design
@section size_t

All objects (blocks of memory, etc) should be measured in terms of a
@code{size_t} type.  Therefore any iterations (e.g. @code{for(i=0; i<N;
i++)}) should also use an index of type @code{size_t}.  

Don't mix @code{int} and @code{size_t}.  They are @emph{not}
interchangeable.

If you need to write a descending loop you have to be careful because
@code{size_t} is unsigned, so instead of

@example
for (i = N - 1; i >= 0; i--) @{ ... @} /* DOESN'T WORK */
@end example

@noindent
use something like

@example
for (i = N; i-- > 0;) @{ ... @}
@end example

@noindent
to avoid problems with wrap-around at @code{i=0}.  Note that the
post-decrement ensures that the loop variable is tested before it
reaches zero.  Beware that @code{i} will wraparound on exit from the
loop.  (This could also be written as @code{for (i = N; i--;)} since
the test for @code{i>0} is equivalent to @code{i!=0} for an unsigned
integer)

If you really want to avoid confusion use a separate variable to invert
the loop order,
@example
for (i = 0; i < N; i++) @{ j = N - i; ... @}
@end example

Note (BJG). Originally, I suggested using

@example
for (i = N; i > 0 && i--;) @{ ... @}
@end example

which makes the test for @code{i>0} explicit and leaves @code{i=0} on
exit from the loop.  However, it is slower as there is an additional
branch which prevents unrolling.  Thanks to J. Seward for pointing
this out.

Note: As a matter of style, please use post-increment (@code{i++}) and
post-decrement (@code{i--}) operators by default and only use
pre-increment (@code{++i}) and pre-decrement (@code{--i}) operators
where specifically needed.

@node Arrays vs Pointers, Pointers, size_t, Design
@section Arrays vs Pointers

A function can be declared with either pointer arguments or array
arguments.  The C standard considers these to be equivalent. However, it
is useful to distinguish between the case of a pointer, representing a
single object which is being modified, and an array which represents a
set of objects with unit stride (that are modified or not depending on
the presence of @code{const}).  For vectors, where the stride is not
required to be unity, the pointer form is preferred.

@smallexample
/* real value, set on output */
int foo (double * x);                           

/* real vector, modified */
int foo (double * x, size_t stride, size_t n);  

/* constant real vector */
int foo (const double * x, size_t stride, size_t n);  

/* real array, modified */
int bar (double x[], size_t n);                 

/* real array, not modified */
int baz (const double x[], size_t n);           
@end smallexample

@node  Pointers, Constness, Arrays vs Pointers, Design
@section Pointers

Avoid dereferencing pointers on the right-hand side of an expression where
possible.  It's better to introduce a temporary variable.  This is
easier for the compiler to optimise and also more readable since it
avoids confusion between the use of @code{*} for multiplication and
dereferencing.

@example
while (fabs (f) < 0.5)
@{
  *e = *e - 1;
  f  *= 2;
@}
@end example

@noindent
is better written as,

@example
@{ 
  int p = *e;

  while (fabs(f) < 0.5)
    @{
     p--;
     f *= 2;
    @}

  *e = p;
@}
@end example

@node Constness, Pseudo-templates, Pointers, Design
@section Constness

Use @code{const} in function prototypes wherever an object pointed to by
a pointer is constant (obviously).  For variables which are meaningfully
constant within a function/scope use @code{const} also.  This prevents
you from accidentally modifying a variable which should be constant
(e.g. length of an array, etc).  It can also help the compiler do
optimization.  These comments also apply to arguments passed by value
which should be made @code{const} when that is meaningful.

@node Pseudo-templates, Arbitrary Constants, Constness, Design
@section Pseudo-templates

There are some pseudo-template macros available in @file{templates_on.h}
and @file{templates_off.h}.  See a directory link @file{block} for
details on how to use them.  Use sparingly, they are a bit of a
nightmare, but unavoidable in places.

In particular, the convention is: templates are used for operations on
"data" only (vectors, matrices, statistics, sorting).  This is intended
to cover the case where the program must interface with an external
data-source which produces a fixed type. e.g. a big array of char's
produced by an 8-bit counter.

All other functions can use double, for floating point, or the
appropriate integer type for integers (e.g. unsigned long int for random
numbers).  It is not the intention to provide a fully templated version
of the library. 

That would be "putting a quart into a pint pot". To summarize, almost
everything should be in a "natural type" which is appropriate for
typical usage, and templates are there to handle a few cases where it is
unavoidable that other data-types will be encountered.

For floating point work "double" is considered a "natural type".  This
sort of idea is a part of the C language.

@node  Arbitrary Constants, Test suites, Pseudo-templates, Design
@section Arbitrary Constants

Avoid arbitrary constants.

For example, don't hard code "small" values like '1e-30', '1e-100' or
@code{10*GSL_DBL_EPSILON} into the routines.  This is not appropriate
for a general purpose library.

Compute values accurately using IEEE arithmetic.  If errors are
potentially significant then error terms should be estimated reliably
and returned to the user, by analytically deriving an error propagation
formula, not using guesswork.

A careful consideration of the algorithm usually shows that arbitrary
constants are unnecessary, and represent an important parameter which
should be accessible to the user.

For example, consider the following code:

@example
if (residual < 1e-30) @{
   return 0.0;  /* residual is zero within round-off error */
@}
@end example

@noindent
This should be rewritten as,

@example
   return residual;
@end example

@noindent
in order to allow the user to determine whether the residual is 
significant or not.

The only place where it is acceptable to use constants like
@code{GSL_DBL_EPSILON} is in function approximations, (e.g. taylor
series, asymptotic expansions, etc).  In these cases it is not an
arbitrary constant, but an inherent part of the algorithm.

@node Test suites, Compilation, Arbitrary Constants, Design
@section Test suites

The implementor of each module should provide a reasonable test suite
for the routines.

The test suite should be a program that uses the library and checks the
result against known results, or invokes the library several times and
does a statistical analysis on the results (for example in the case of
random number generators).

Ideally the one test program per directory should aim for 100% path
coverage of the code.  Obviously it would be a lot of work to really