"Numerical Recipes in C" is nonportable code

[Rob Carriere]

In article <1673 at dataio.Data-IO.COM> bright at dataio.Data-IO.COM (Walter Bright)
writes:
>In article <531 at accelerator.eng.ohio-state.edu> rob at kaa.eng.ohio-state.edu 
>(Rob Carriere) writes:
><In article <13258 at mimsy.UUCP< chris at mimsy.UUCP (Chris Torek) writes:
><< [ still on the b = malloc( foo );  bb = b - 1; code in NumRecipes ]
>< [ claiming astonishment in accordance with the Least of.. Principle ]
>On a segmented architecture, like 8086's, malloc can and does return
>a value that is a pointer to the beginning of a segment. 
Yes, and as I've written people who e-mailed me this argument, that
spells out ``broken compiler'' if it gives problems ('cause you can't
malloc more than 64K that way).  Of course on a machine with large
segments, the counterargument doesn't quite hold water either, so...
>[ unItelligent CPU explanation deleted ]
>	array = malloc(MAX * sizeof(array[0]));
>	for (p = &array[MAX-1]; p >= &array[0]; p--)

No!  That wasn't the problem!  (Wish it was, that'd be easy to
avoid!).

The problem is that the authors of Numerical Recipes (NR) observe,
correctly, that many numerical problems are naturally non-zero based.
This gives you the choice between carrying around boatloads of index
arithmatic (inefficient and error-prone), or making non-zero based
arrays.  They opt for the latter, in the following way:

float *my_vec;  /* this is going to be a vector */
int nl, nh;
...
my_vec = vector( nl, nh );  /* allocates a vector with lowest valid
                               index nl, and highest valid index nh
                            */
...
my_vec[3] = foo(bar);
...

Where we have:

float *vector( nl, nh )
     int nl;
     int nh;
{
    float *v;

    v = (float *)malloc( ( nh-nl +1 )* sizeof(float) );
    if( v == 0 ) nrerror( "Allocation error in vector()" );

    return v - nl;
}

This is quite a bit more disciplined than the example above; it is
also quite bit more fundamental.  Fortunately, as far as I've checked
at least, NR only uses vectors and matrices with either 0 or unit
offset, so on broken architectures you could always do 

                  malloc( (nh   + 1 )* sizeof(float) );

    return v;

This would waste a float per vector, and a pointer-to-float plus n
floats for an n-by-something matrix.  Ugly, but it works.  (and we
*are* the throw-away culture after all :-)

Rob Carriere