the quicksort algorithm in Numerical Recipes in C (2nd edition)

[Edward Hua]

Hi,

I'm wondering if anybody has ever copied the quicksort algorithm from
the book Numerical Recipes in C: The Art of Scientific Computing (2nd
ed.), by Press, Teukolsky, Vetterling, and Flannery, in Chapter 8.
Quicsort is supposed to sort an array of, say, doubles, in ascending
order.

I tried to copy this algorithm line by line to test it. It seems that
there's an error in this algorithm as given in this book. Namely, the
first two entry (0th entry) of the array is never moved, and the 2nd entry
(1st entry) is also set to 0 after the sorting was supposed to be
completed. All other entries seemed to be correctly updated.

I'm wondering if anybody on this forum's ever had the same problem when
using this algorithm taken from this book. If so, could you tell me what
the problem is and how to fix it? I'm a little time-pressed to investigate
this myself. Thanks in advance.

-Ed
------------------------------------------------------------------
To give an example of this error, I have an initial array of 20 doubles as
follows:

sample_array[0] = 5.754104
sample_array[1] = 0.973984
sample_array[2] = 14.586500
sample_array[3] = 11.796852
sample_array[4] = 4.032750
sample_array[5] = 0.618612
sample_array[6] = 8.638782
sample_array[7] = 2.309558
sample_array[8] = 4.144175
sample_array[9] = 15.159169
sample_array[10] = 1.541796
sample_array[11] = 18.063583
sample_array[12] = 4.739992
sample_array[13] = 4.905477
sample_array[14] = 8.836610
sample_array[15] = 15.108313
sample_array[16] = 18.243906
sample_array[17] = 14.667340
sample_array[18] = 5.801667
sample_array[19] = 2.857230

And after the sorting, it becomes:

sample_array[0]=5.754104
sample_array[1]=0.000000
sample_array[2]=0.618612
sample_array[3]=0.973984
sample_array[4]=1.541796
sample_array[5]=2.309558
sample_array[6]=2.857230
sample_array[7]=4.032750
sample_array[8]=4.144175
sample_array[9]=4.739992
sample_array[10]=4.905477
sample_array[11]=5.801667
sample_array[12]=8.638782
sample_array[13]=8.836610
sample_array[14]=11.796852
sample_array[15]=14.586500
sample_array[16]=14.667340
sample_array[17]=15.108313
sample_array[18]=15.159169
sample_array[19]=18.063583

[Edward Hua]

Please ignore this question. I figured it out.

-Ed

On Wed, 26 Oct 2005, Edward Hua wrote:

Hi,

I'm wondering if anybody has ever copied the quicksort algorithm from
the book Numerical Recipes in C: The Art of Scientific Computing (2nd
ed.), by Press, Teukolsky, Vetterling, and Flannery, in Chapter 8.
Quicsort is supposed to sort an array of, say, doubles, in ascending
order.

I tried to copy this algorithm line by line to test it. It seems that
there's an error in this algorithm as given in this book. Namely, the
first two entry (0th entry) of the array is never moved, and the 2nd entry
(1st entry) is also set to 0 after the sorting was supposed to be
completed. All other entries seemed to be correctly updated.

I'm wondering if anybody on this forum's ever had the same problem when
using this algorithm taken from this book. If so, could you tell me what
the problem is and how to fix it? I'm a little time-pressed to investigate
this myself. Thanks in advance.

-Ed
------------------------------------------------------------------
To give an example of this error, I have an initial array of 20 doubles as
follows:

sample_array[0] = 5.754104
sample_array[1] = 0.973984
sample_array[2] = 14.586500
sample_array[3] = 11.796852
sample_array[4] = 4.032750
sample_array[5] = 0.618612
sample_array[6] = 8.638782
sample_array[7] = 2.309558
sample_array[8] = 4.144175
sample_array[9] = 15.159169
sample_array[10] = 1.541796
sample_array[11] = 18.063583
sample_array[12] = 4.739992
sample_array[13] = 4.905477
sample_array[14] = 8.836610
sample_array[15] = 15.108313
sample_array[16] = 18.243906
sample_array[17] = 14.667340
sample_array[18] = 5.801667
sample_array[19] = 2.857230

And after the sorting, it becomes:

sample_array[0]=5.754104
sample_array[1]=0.000000
sample_array[2]=0.618612
sample_array[3]=0.973984
sample_array[4]=1.541796
sample_array[5]=2.309558
sample_array[6]=2.857230
sample_array[7]=4.032750
sample_array[8]=4.144175
sample_array[9]=4.739992
sample_array[10]=4.905477
sample_array[11]=5.801667
sample_array[12]=8.638782
sample_array[13]=8.836610
sample_array[14]=11.796852
sample_array[15]=14.586500
sample_array[16]=14.667340
sample_array[17]=15.108313
sample_array[18]=15.159169
sample_array[19]=18.063583

[Ben Pfaff]

Edward Hua <ey**@ece.corne ll.edu> writes:

I'm wondering if anybody has ever copied the quicksort algorithm from
the book Numerical Recipes in C: The Art of Scientific Computing (2nd
ed.), by Press, Teukolsky, Vetterling, and Flannery, in Chapter 8.

Copying code from that book is a mistake. The explanatory text
can be insightful, but the code is not well written.
--
"When I have to rely on inadequacy, I prefer it to be my own."
--Richard Heathfield

[Rich Gibbs]

Edward Hua said the following, on 10/26/05 10:15:

Hi,

I'm wondering if anybody has ever copied the quicksort algorithm from
the book Numerical Recipes in C: The Art of Scientific Computing (2nd
ed.), by Press, Teukolsky, Vetterling, and Flannery, in Chapter 8.
Quicsort is supposed to sort an array of, say, doubles, in ascending
order.

I tried to copy this algorithm line by line to test it. It seems that
there's an error in this algorithm as given in this book. Namely, the
first two entry (0th entry) of the array is never moved, and the 2nd entry
(1st entry) is also set to 0 after the sorting was supposed to be
completed. All other entries seemed to be correctly updated.

I'm wondering if anybody on this forum's ever had the same problem when
using this algorithm taken from this book. If so, could you tell me what
the problem is and how to fix it? I'm a little time-pressed to investigate
this myself. Thanks in advance.

This really is not a C question as you have posed it, but I have a
suspicion that you have a C-related problem.

The original _Numerical Recipes_ book used FORTRAN, I think. I don't
have my copy of it, or of _Numerical Recipes in C_, handy, but I seem to
recall that in the C version the authors used some funky macro hacks so
that they could use 1-origin arrays (a la FORTRAN), rather than C's
0-origin arrays.

If you have an "off-by-one" problem related to this, it might account
for both of the symptoms you are seeing: the zeroth element is never
looked at; and the element 'array[20]' (which of course doesn't exist as
far as C is concerned) _is_ used, and might be zero by chance.

--
Rich Gibbs
ri*****@gmail.c om
"You can observe a lot by watching." -- Yogi Berra

[Edward Hua]

Ok. Is there a more authoritative book/on-line refrence that has an
archive of C functions (like those in NRiC)? Thanks.

-Ed

On Wed, 26 Oct 2005, Ben Pfaff wrote:

Edward Hua <ey**@ece.corne ll.edu> writes:

I'm wondering if anybody has ever copied the quicksort algorithm from
the book Numerical Recipes in C: The Art of Scientific Computing (2nd
ed.), by Press, Teukolsky, Vetterling, and Flannery, in Chapter 8.

Copying code from that book is a mistake. The explanatory text
can be insightful, but the code is not well written.
--
"When I have to rely on inadequacy, I prefer it to be my own."
--Richard Heathfield

[pete]

Edward Hua wrote:

how to fix it?

#define BYTE_SWAP(A, B) \
{ \
p1 = (A); \
p2 = (B); \
end = p2 + size; \
do { \
swap = *p1; \
*p1++ = *p2; \
*p2++ = swap; \
} while (p2 != end); \
}

void q2sort(void *base, size_t nmemb, size_t size,
int (*compar)(const void *, const void *))
{
unsigned char *left;
size_t middle, last, right;
struct {
size_t bytes;
void *base;
} stack[CHAR_BIT * sizeof nmemb], *stack_ptr;
unsigned char *p1, *p2, *end, swap;

stack -> bytes = nmemb * size;
stack -> base = base;
stack_ptr = stack + 1;
do {
--stack_ptr;
if (stack_ptr -> bytes > size) {
left = stack_ptr -> base;
right = last = stack_ptr -> bytes - size;
middle = size;
while (compar(left, left + middle) > 0 && middle != right) {
middle += size;
}
while (compar(left + last, left) > 0) {
last -= size;
}
while (last > middle) {
BYTE_SWAP(left + middle, left + last);
do {
middle += size;
} while (compar(left, left + middle) > 0);
do {
last -= size;
} while (compar(left + last, left) > 0);
}
BYTE_SWAP(left, left + last);
if (right - last > last) {
stack_ptr -> base = left + last + size;
stack_ptr++ -> bytes = right - last;
stack_ptr -> base = left;
stack_ptr++ -> bytes = last;
} else {
stack_ptr++ -> bytes = last;
stack_ptr -> base = left + last + size;
stack_ptr++ -> bytes = right - last;
}
}
} while (stack_ptr != stack);
}

--
pete