* Bug #8779 fixed - gsort() did not preserve order of equal elements, in lexicographi...
[scilab.git] / scilab / modules / elementary_functions / src / c / qsort.c
1 /*
2 * See Copyright below
3 * Copyright (c) 1992, 1993
4 * The Regents of the University of California.  All rights reserved.
5 */
6
7 #include <stdlib.h>
8 #include <string.h>
9 /*--------------------------------------------------------------------------*/
10 #include "qsort.h"
11 #include "qsort-int.h"
12 #include "qsort-short.h"
13 #include "qsort-char.h"
14 #include "qsort-double.h"
15 #include "qsort-string.h"
16 #include "core_math.h"
17 /*--------------------------------------------------------------------------*/
18
19 /*--------------------------------------------------------------------------*/
20
21 /*--------------------------------------------------------------------------*/
22 /*      $NetBSD: qsort.c,v 1.5 1995/12/28 08:52:36 thorpej Exp $        */
23 /*-
24 * Copyright (c) 1992, 1993
25 *       The Regents of the University of California.  All rights reserved.
26 *
27 * Redistribution and use in source and binary forms, with or without
28 * modification, are permitted provided that the following conditions
29 * are met:
30 * 1. Redistributions of source code must retain the above copyright
31 *    notice, this list of conditions and the following disclaimer.
32 * 2. Redistributions in binary form must reproduce the above copyright
33 *    notice, this list of conditions and the following disclaimer in the
34 *    documentation and/or other materials provided with the distribution.
35 * 3. All advertising materials mentioning features or use of this software
36 *    must display the following acknowledgement:
37 *       This product includes software developed by the University of
38 *       California, Berkeley and its contributors.
39 * 4. Neither the name of the University nor the names of its contributors
40 *    may be used to endorse or promote products derived from this software
41 *    without specific prior written permission.
42 *
43 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
44 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
45 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
46 * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
47 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
48 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
49 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
50 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
51 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
52 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
53 * SUCH DAMAGE.
54 *
55 *  Modified for Scilab Jean-Philippe Chancelier  to keep a permutation index
56 *  Modified for Scilab by Serge Steer to make it stable when permutation index is computed.
57 */
58 /*--------------------------------------------------------------------------*/
59 /*
60 * Qsort routine from Bentley & McIlroy's "Engineering a Sort Function".
61 * Software---Practice and Experience, 23(11):1249-1265
62 */
63 /*--------------------------------------------------------------------------*/
64 void sciqsort(char *a, char *tab, int flag, int n, int es, int es1, int (*cmp)(), int (*swapcode)(), int (*lswapcodeind)())
65 {
66     char *pa, *pb, *pc, *pd, *pl, *pm, *pn;
67     char *taba, *tabb, *tabc, *tabd, *tabl, *tabm, *tabn;
68     int d, dind, r, r1,  swap_cnt;
69
70 loop:
71     swap_cnt = 0;
72     if (n < 7)   /* Insertion sort on smallest arrays */
73     {
74         for (pm = a + es, tabm = tab + es1 ; pm < (char *) a + n * es; pm += es, tabm += es1 )
75         {
76             for (pl = pm, tabl = tabm ; pl > (char *) a && cmp(pl - es, pl, tabl - es1, tabl, flag) > 0;  pl -= es, tabl -= es1)
77             {
78                 swapind(tabl, tabl - es1);
79                 swap(pl, pl - es);
80             }
81         }
82         return;
83     }
84     /*Determine the pivot */
85     pm = a + (n / 2) * es;/* Small arrays, middle element */
86     tabm = tab + (n / 2) * es1 ;
87
88     pn = a + (n - 1) * es;
89     tabn = tab + (n - 1) * es1;
90     if (n > 7)
91     {
92         pl = a;
93         tabl = tab;
94         if (n > 40)  /* Big arrays, pseudomedian of 9 */
95         {
96             dind = (n / 8) * es1;
97             d =   (n / 8) * es;
98             med3(pl, tabl, pl, pl + d, pl + 2 * d, tabl, tabl + dind, tabl + 2 * dind, cmp);
99             med3(pm, tabm, pm - d, pm, pm + d, tabm - dind, tabm, tabm + dind, cmp);
100             med3(pn, tabn, pn - 2 * d, pn - d, pn, tabn - 2 * dind, tabn - dind, tabn, cmp);
101         }
102         med3(pm, tabm, pl, pm, pn, tabl, tabm, tabn, cmp);
103     }
104     /* Put it at the first position */
105     /* Partionning */
106     if (cmp(pn, a, tabn, tab, flag))
107     {
108         swapind(tab, tabm);
109         swap(a, pm);
110     }
111
112     /* pointers on data array */
113     pa = pb = a + es;/* pa and pb start from the beginning of the array */
114     pc = pd = a + (n - 1) * es;/* pc and pd start from the end of the array */
115
116     /* similar pointers for index array */
117     taba = tabb = tab + es1;
118     tabc = tabd = tab + (n - 1) * es1;
119
120     /* here we have
121     |a  |pa                            | pc|
122     |a  |pb                            | pd|
123     |*a |                ?             | ? |
124     */
125     for (;;)
126     {
127         /* increase the pointer pb while it points on values lesser  than the pivot (pointer a) */
128         while (pb <= pc && (r = cmp(pb, a, tabb, tab, flag)) <= 0)
129         {
130             if (r == 0)  /*The pivot and  value pointed to by pb are equal */
131             {
132                 /* store the equal value at the location pa and increase pa */
133                 swap_cnt = 1;
134                 swapind(taba, tabb);
135                 taba += es1;
136                 swap(pa, pb);
137                 pa += es;
138             }
139             pb += es;/* next number */
140             tabb += es1;
141         }
142
143         /* here pb points on a value greater than the pivot */
144         /* decrease the pointer pc while it points on a value greater than the pivot (pointer a) */
145         while (pb <= pc && (r = cmp(pc, a, tabc, tab, flag)) >= 0)
146         {
147             if (r == 0)  /*The pivot and  value pointed to by pc are equal */
148             {
149                 /* store the equal value at the location pd and decrease pd */
150                 swap_cnt = 1;
151                 swapind(tabc, tabd);
152                 tabd -= es1;
153                 swap(pc, pd);
154                 pd -= es;
155             }
156             pc -= es;
157             tabc -= es1;
158         }
159         /* here pc points on a value lesser than the pivot */
160         if (pb > pc)
161         {
162             /* here we have
163             |a      |pa      |pc|pb     pd|      $|
164             |   =*a |    <*a |       >*a  | =*a  $|
165             */
166             /* partition is done */
167             break;
168         }
169         /*here
170         pc  points on a value lesser than the pivot
171         and
172         pb points on a value greater than the pivot
173         swap the values
174         */
175         swapind(tabb, tabc);
176         tabb += es1;
177         tabc -= es1;
178         swap(pb, pc);
179         swap_cnt = 1;
180         /* increase pb and decrease pc */
181         pb += es;
182         pc -= es;
183         /* here we have
184         |a      |pa      |pb       pc|       pd|      $|
185         |   =*a |    <*a |     ?     |    >*a  | =*a  $|
186         */
187     }
188
189     if (swap_cnt == 0)    /* Switch to insertion sort */
190     {
191         for (pm = a + es, tabm = tab + es1 ; pm < (char *) a + n * es; pm += es, tabm += es1)
192         {
193             for (pl = pm, tabl = tabm ; pl > (char *) a && cmp(pl - es, pl, tabl - es1, tabl, flag) > 0;  pl -= es, tabl -= es1)
194             {
195                 swapind(tabl, tabl - es1);
196                 swap(pl, pl - es);
197             }
198         }
199         return;
200     }
201     /* put the equal values in the middle */
202     pn = a + n * es;
203     r = (int)Min(pa - (char *)a, pb - pa);
204     vecswap(a, pb - r, r);
205
206     tabn = tab + n * es1 ;
207     r1 = (int)Min(taba - (char *) tab, tabb - taba);
208     vecswapind(tab, tabb - r1, r1);
209
210     r = (int)Min(pd - pc, pn - pd - es);
211     vecswap(pb, pn - r, r);
212
213     r1 = (int)Min(tabd - tabc, tabn - tabd - es1 );
214     vecswapind(tabb, tabn - r1, r1);
215
216     if ((r = (int)(pb - pa)) > es )
217         /* recall  sciqsort for the lower part */
218     {
219         sciqsort(a, tab, flag, r / es, es, es1, cmp, swapcode, lswapcodeind);
220     }
221     if ((r = (int)(pd - pc)) > es)
222     {
223         /* Iterate rather than recurse to save stack space */
224         a = pn - r;
225         tab = tabn - (tabd - tabc);
226         n = r / es;
227         goto loop;
228     }
229 }
230 /*--------------------------------------------------------------------------*/
231 int swapcodeint(char * parmi, char * parmj, int n, int incr)
232 {
233     int i = n;
234     register int *pi = (int *) (parmi);
235     register int *pj = (int *) (parmj);
236     register int inc1 = incr / sizeof(int);
237     do
238     {
239         register int t = *pi;
240         *pi = *pj;
241         *pj = t;
242         pi += inc1;
243         pj += inc1;
244     }
245     while (--i > 0);
246     return(0);
247 }
248 /*--------------------------------------------------------------------------*/