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 /*--------------------------------------------------------------------------*/