How to fix the GSL API

["E. Robert Tisdale" <edwin at netwood dot net>]

/*
I am suggesting changes to the GSL API
as illustrated in the attached ANSI C program.
I compiled the program with my GNU C compiler
and ran it under the Linux operating system.

        $ uname -m -r -s -v
        Linux 2.0.36 #1 Tue Oct 13 22:17:11 EDT 1998 i686
        $ gcc --version
        2.7.2.3
        $ gcc -DGSL_DEBUG -O2 -o gsl gsl.c
        $ gsl
        0 1 2 3 4 5 6 7 8 9
        4 = sizeof(struct gsl_handle_private)
        4 = sizeof(gsl_handle)
        20 = sizeof(struct gsl_vector_private)
        20 = sizeof(gsl_vector)
        $

The actual data representation is hidden in opaque objects.
Zero overhead access to data members is provided
through inline functions or C preprocessor macros.
Instead of allocating vector view objects
from the free store (heap)
they are allocated from automatic (stack) memory
just like subvector view objects.
The gsl_handle handle abstraction
is substituted for a pointer to a gsl_block object.
Negative and zero as well as positive strides are supported.
Empty arrays, vector views and subvector views are OK.
Type definitions for gsl_offset, gsl_extent and gsl_stride
liberate GSL developers to choose integral types
other than size_t for offsets and extents
or ptrdiff_t for strides.
Defining C preprocessor macro GSL_DEBUG
includes optional run time error checking code.
The error checking code helps application programmers
to detect programming errors at run time during testing.
After the application has been thoroughly tested,
error checking can be turned off
so that it doesn't affect performance.
*/
        #include<stdio.h>
        #include<stdlib.h>

        typedef int ptrdiff_t;

        struct gsl_handle_private {
          double*       pD;
          };

        typedef struct {        /* opaque data structure */
          int dummy[sizeof(struct gsl_handle_private)/sizeof(int)];
          } gsl_handle;
        typedef    size_t gsl_offset;
        typedef    size_t gsl_extent;
        typedef ptrdiff_t gsl_stride;
        #ifdef GSL_DEBUG
        typedef int       gsl_keeper;
        #endif GSL_DEBUG

        inline void
        gsl_handle_initialize(gsl_handle* ph, double* pd) {
          struct gsl_handle_private* php = (struct gsl_handle_private*)ph;
          php->pD = pd; }

        inline gsl_handle
        gsl_handle_allocate(size_t n) {
          gsl_handle    h;
          void* p = (0 < n)? malloc(n*sizeof(double)): NULL;
          gsl_handle_initialize(&h, (double*)p);
          return h; }

        inline void
        gsl_handle_free(gsl_handle* ph) {
          struct gsl_handle_private* php = (struct gsl_handle_private*)ph;
          free((void*)(php->pD));
          php->pD = (double*)NULL; }

        inline double
        gsl_handle_get(const gsl_handle* ph, gsl_offset j) {
          struct gsl_handle_private* php = (struct gsl_handle_private*)ph;
          return php->pD[j]; }

        inline void
        gsl_handle_put(const gsl_handle* ph, gsl_offset j, double x) {
          struct gsl_handle_private* php = (struct gsl_handle_private*)ph;
          php->pD[j] = x; }

        struct gsl_vector_private {
          gsl_handle    H;
          gsl_offset    O;
          gsl_extent    N;
          gsl_stride    S;
        #ifdef GSL_DEBUG
          gsl_keeper    K;      /* optional data member */
        #endif GSL_DEBUG
          };

        typedef struct {        /* opaque data structure */
          int   dummy[sizeof(struct gsl_vector_private)/sizeof(int)];
          } gsl_vector;

        inline void
        gsl_vector_initialize(gsl_vector* pv,
          gsl_handle h, gsl_offset o, gsl_extent n, gsl_stride s) {
          struct gsl_vector_private* pvp = (struct gsl_vector_private*)pv;
          pvp->H = h; pvp->O = o; pvp->N = n; pvp->S = s;
        #ifdef GSL_DEBUG
          pvp->K = (gsl_keeper)0;
        #endif GSL_DEBUG
          }

        inline gsl_vector
        gsl_vector_allocate(gsl_extent n) {
          gsl_vector    v;
          gsl_vector_initialize(&v,
          gsl_handle_allocate(n), (gsl_offset)0, n, (gsl_stride)1);
        #ifdef GSL_DEBUG
          { struct gsl_vector_private* pvp = (struct gsl_vector_private*)&v;
            pvp->K = (gsl_keeper)1; }
        #endif GSL_DEBUG
          return v; }

        inline void
        gsl_vector_free(gsl_vector* pv) {
          struct gsl_vector_private* pvp = (struct gsl_vector_private*)pv;
        #ifdef GSL_DEBUG
          if (pvp->K)
            gsl_handle_free(&pvp->H);
          else
            fprintf(stderr, "In function gsl_vector_free(gsl_vector*),\n"
              "attempt to free data array from subvector!\n");
        #else
          gsl_handle_free(&pvp->H);
        #endif GSL_DEBUG
          }

        inline gsl_handle
        gsl_vector_handle(const gsl_vector* pv) {
          struct gsl_vector_private* pvp = (struct gsl_vector_private*)pv;
          return pvp->H; }

        inline gsl_offset
        gsl_vector_offset(const gsl_vector* pv) {
          struct gsl_vector_private* pvp = (struct gsl_vector_private*)pv;
          return pvp->O; }

        inline gsl_extent
        gsl_vector_extent(const gsl_vector* pv) {
          struct gsl_vector_private* pvp = (struct gsl_vector_private*)pv;
          return pvp->N; }

        inline gsl_stride
        gsl_vector_stride(const gsl_vector* pv) {
          struct gsl_vector_private* pvp = (struct gsl_vector_private*)pv;
          return pvp->S; }

        inline gsl_vector
        gsl_subvector(const gsl_vector* pv,
          gsl_offset j, gsl_extent n, gsl_stride s) {
          gsl_vector    u;
          gsl_vector_initialize(&u,
            gsl_vector_handle(pv),
            gsl_vector_offset(pv) + (gsl_stride)j*gsl_vector_stride(pv),
            n,
            s*gsl_vector_stride(pv));
          return u; }

        inline double
        gsl_vector_get(const gsl_vector* pv, gsl_offset j) {
          gsl_handle    h = gsl_vector_handle(pv);
          gsl_offset    o = gsl_vector_offset(pv);
          gsl_stride    s = gsl_vector_stride(pv);
          return gsl_handle_get(&h, o + (gsl_stride)j*s); }

        inline void
        gsl_vector_put(const gsl_vector* pv, gsl_offset j, double x) {
          gsl_handle    h = gsl_vector_handle(pv);
          gsl_offset    o = gsl_vector_offset(pv);
          gsl_stride    s = gsl_vector_stride(pv);
          gsl_handle_put(&h, o + (gsl_stride)j*s, x); }

        inline void
        gsl_vector_fill(gsl_vector* pv, double x) {
          gsl_handle    h = gsl_vector_handle(pv);
          gsl_offset    k = gsl_vector_offset(pv);
          gsl_extent    n = gsl_vector_extent(pv);
          gsl_stride    s = gsl_vector_stride(pv);
          gsl_offset    j = 0;
          for (j = 0; j < n; j++) {
            gsl_handle_put(&h, k, x);
            k += s;
            }
          }

        inline void
        gsl_vector_ramp(gsl_vector* pv) {
          gsl_handle    h = gsl_vector_handle(pv);
          gsl_offset    k = gsl_vector_offset(pv);
          gsl_extent    n = gsl_vector_extent(pv);
          gsl_stride    s = gsl_vector_stride(pv);
          gsl_offset    j = 0;
          for (j = 0; j < n; j++) {
            gsl_handle_put(&h, k, (double)j);
            k += s;
            }
          }

        int
        main() {
          const
          gsl_extent    n = 10;
          gsl_vector    v = gsl_vector_allocate(n);
          gsl_vector_ramp(&v);
          { gsl_offset  j = 0;
            for (j = 0; j < n; j++)
              printf("%g ", gsl_vector_get(&v, j));
            printf("\n");
            }
          printf("%d = sizeof(struct gsl_handle_private)\n",
            sizeof(struct gsl_handle_private));
          printf("%d = sizeof(gsl_handle)\n",
            sizeof(gsl_handle));
          printf("%d = sizeof(struct gsl_vector_private)\n",
            sizeof(struct gsl_vector_private));
          printf("%d = sizeof(gsl_vector)\n",
            sizeof(gsl_vector));
          return 0;
          }