sys/src/cmd/disk/sacfs/ssort6.c

   1 #include <u.h>
   2 #include <libc.h>
   3 #include "ssort.h"
   4
   5 #define pred(i, h) ((t=(i)-(h))<0?  t+n: t)
   6 #define succ(i, h) ((t=(i)+(h))>=n? t-n: t)
   7
   8 enum
   9 {
  10         BUCK = ~(~0u>>1),       /* high bit */
  11         MAXI = ~0u>>1,          /* biggest int */
  12 };
  13
  14 typedef int     Elem;
  15 static  void    qsort2(Elem*, Elem*, int n);
  16 static  int     ssortit(int a[], int p[], int s[], int q[], int n, int h, int *pe, int nbuck);
  17 static  void    lift(int si, int q[], int i);
  18 int     sharedlen(int i, int j, int s[], int q[]);
  19
  20 int
  21 ssort(int a[], int s[], int n)
  22 {
  23         int i, l;
  24         int c, cc, ncc, lab, cum, nbuck;
  25         int k;
  26         int *p = 0;
  27         int result = 0;
  28         int *q = 0;
  29         int *al;
  30         int *pl;
  31
  32 #       define finish(r)  if(1){result=r; goto out;}else
  33
  34         for(k=0,i=0; i<n; i++)
  35                 if(a[i] > k)
  36                         k = a[i];       /* max element */
  37         k++;
  38         if(k>n)
  39                 finish(2);
  40
  41         nbuck = 0;
  42         p = malloc(n*sizeof(int));
  43         if(p == 0)
  44                 finish(1);
  45
  46         if(s) {                                 /* shared lengths */
  47                 q = malloc(((n+1)>>1)*sizeof(int));
  48                 if(q == 0)
  49                         finish(1);
  50                 for(i=0; i<n; i++)
  51                         s[i] = q[i>>1] = MAXI;
  52                 q[i>>1] = MAXI;
  53         }
  54
  55         pl = p + n - k;
  56         al = a;
  57         memset(pl, -1, k*sizeof(int));
  58
  59         for(i=0; i<n; i++) {            /* (1) link */
  60                 l = a[i];
  61                 al[i] = pl[l];
  62                 pl[l] = i;
  63         }
  64
  65         for(i=0; i<k; i++)              /* check input - no holes */
  66                 if(pl[i]<0)
  67                         finish(2);
  68
  69
  70         lab = 0;                        /* (2) create p and label a */
  71         cum = 0;
  72         i = 0;
  73         for(c = 0; c < k; c++){
  74                 for(cc = pl[c]; cc != -1; cc = ncc){
  75                         ncc = al[cc];
  76                         al[cc] = lab;
  77                         cum++;
  78                         p[i++] = cc;
  79                 }
  80                 if(lab + 1 == cum) {
  81                         i--;
  82                 } else {
  83                         p[i-1] |= BUCK;
  84                         nbuck++;
  85                 }
  86                 if(s) {
  87                         s[lab] = 0;
  88                         lift(0, q, lab);
  89                 }
  90                 lab = cum;
  91         }
  92
  93         ssortit(a, p, s, q, n, 1, p+i, nbuck);
  94         memcpy(a, p, n*sizeof(int));
  95
  96 out:
  97         free(p);
  98         free(q);
  99         return result;
 100 }
 101
 102 /*
 103  * calculate the suffix array for the bytes in buf,
 104  * terminated by a unique end marker less than any character in buf
 105  * returns the index of the identity permutation,
 106  * and -1 if there was an error.
 107  */
 108 int
 109 ssortbyte(uchar buf[], int p[], int s[], int n)
 110 {
 111         int *a, *q, buckets[256*256];
 112         int i, last, lastc, cum, c, cc, ncc, lab, id, nbuck;
 113
 114         a = malloc((n+1)*sizeof(int));
 115         if(a == nil)
 116                 return -1;
 117
 118         q = nil;
 119         if(s) {                                 /* shared lengths */
 120                 q = malloc(((n+2)>>1)*sizeof(int));
 121                 if(q == nil){
 122                         free(a);
 123                         return -1;
 124                 }
 125                 for(i=0; i<n+1; i++)
 126                         s[i] = q[i>>1] = MAXI;
 127                 q[i>>1] = MAXI;
 128         }
 129
 130         memset(buckets, -1, sizeof(buckets));
 131         c = buf[n-1] << 8;
 132         last = c;
 133         for(i = n - 2; i >= 0; i--){
 134                 c = (buf[i] << 8) | (c >> 8);
 135                 a[i] = buckets[c];
 136                 buckets[c] = i;
 137         }
 138
 139         /*
 140          * end of string comes before anything else
 141          */
 142         a[n] = 0;
 143         if(s) {
 144                 s[0] = 0;
 145                 lift(0, q, 0);
 146         }
 147
 148         lab = 1;
 149         cum = 1;
 150         i = 0;
 151         lastc = 1;              /* won't match c & 0xff00 for any c */
 152         nbuck = 0;
 153         for(c = 0; c < 256*256; c++) {
 154                 /*
 155                  * last character is followed by unique end of string
 156                  */
 157                 if(c == last) {
 158                         a[n-1] = lab;
 159                         if(s) {
 160                                 s[lab] = 0;
 161                                 lift(0, q, lab);
 162                         }
 163                         cum++;
 164                         lab++;
 165                         lastc = c & 0xff00;
 166                 }
 167
 168                 for(cc = buckets[c]; cc != -1; cc = ncc) {
 169                         ncc = a[cc];
 170                         a[cc] = lab;
 171                         cum++;
 172                         p[i++] = cc;
 173                 }
 174                 if(lab == cum)
 175                         continue;
 176                 if(lab + 1 == cum)
 177                         i--;
 178                 else {
 179                         p[i - 1] |= BUCK;
 180                         nbuck++;
 181                 }
 182                 if(s) {
 183                         cc = (c & 0xff00) == lastc;
 184                         s[lab] = cc;
 185                         lift(cc, q, lab);
 186                 }
 187                 lastc = c & 0xff00;
 188                 lab = cum;
 189         }
 190
 191         id = ssortit(a, p, s, q, n+1, 2, p+i, nbuck);
 192         free(a);
 193         free(q);
 194         return id;
 195 }
 196
 197 /*
 198  * can get an interval for the shared lengths from [h,3h) by recording h
 199  * rather than h + sharedlen(..) when relabelling.  if so, no calls to lift are needed.
 200  */
 201 static int
 202 ssortit(int a[], int p[], int shared[], int q[], int n, int h, int *pe, int nbuck)
 203 {
 204         int *s, *ss, *packing, *sorting;
 205         int v, sv, vv, packed, lab, t, i;
 206
 207         for(; h < n && p < pe; h=2*h) {
 208                 packing = p;
 209                 nbuck = 0;
 210
 211                 for(sorting = p; sorting < pe; sorting = s){
 212                         /*
 213                          * find length of stuff to sort
 214                          */
 215                         lab = a[*sorting];
 216                         for(s = sorting; ; s++) {
 217                                 sv = *s;
 218                                 v = a[succ(sv & ~BUCK, h)];
 219                                 if(v & BUCK)
 220                                         v = lab;
 221                                 a[sv & ~BUCK] = v | BUCK;
 222                                 if(sv & BUCK)
 223                                         break;
 224                         }
 225                         *s++ &= ~BUCK;
 226                         nbuck++;
 227
 228                         qsort2(sorting, a, s - sorting);
 229
 230                         v = a[*sorting];
 231                         a[*sorting] = lab;
 232                         packed = 0;
 233                         for(ss = sorting + 1; ss < s; ss++) {
 234                                 sv = *ss;
 235                                 vv = a[sv];
 236                                 if(vv == v) {
 237                                         *packing++ = ss[-1];
 238                                         packed++;
 239                                 } else {
 240                                         if(packed) {
 241                                                 *packing++ = ss[-1] | BUCK;
 242                                         }
 243                                         lab += packed + 1;
 244                                         if(shared) {
 245                                                 v = h + sharedlen(v&~BUCK, vv&~BUCK, shared, q);
 246                                                 shared[lab] = v;
 247                                                 lift(v, q, lab);
 248                                         }
 249                                         packed = 0;
 250                                         v = vv;
 251                                 }
 252                                 a[sv] = lab;
 253                         }
 254                         if(packed) {
 255                                 *packing++ = ss[-1] | BUCK;
 256                         }
 257                 }
 258                 pe = packing;
 259         }
 260
 261         /*
 262          * reconstuct the permutation matrix
 263          * return index of the entire string
 264          */
 265         v = a[0];
 266         for(i = 0; i < n; i++)
 267                 p[a[i]] = i;
 268
 269         return v;
 270 }
 271
 272 /* Propagate a new entry s[i], with value si, into q[]. */
 273 static void
 274 lift(int si, int q[], int i)
 275 {
 276         for(i >>= 1; q[i] > si; i &= ~-i)               /* squash the low 1-bit */
 277                 q[i] = si;
 278 }
 279
 280 /*
 281  * Find in s[] the minimum value over the interval i<=k<=j, using
 282  * tree q[] to do logarithmic, rather than linear search
 283  */
 284 int
 285 sharedlen(int i, int j, int s[], int q[])
 286 {
 287         int k, v;
 288         int min = MAXI;
 289         int max = 0;
 290
 291         if(i > j) {             /* swap i & j */
 292                 i ^= j;
 293                 j ^= i;
 294                 i ^= j;
 295         }
 296         i++;
 297         while(i <= j && min > max) {
 298                 k = i & -i;
 299                 if(i & 1)
 300                         v = s[i];
 301                 else
 302                         v = q[i>>1];
 303                 if(i+k > j+1) {
 304                         if(v > max)
 305                                 max = v;
 306                         if(s[i] < min)
 307                                 min = s[i];
 308                         i++;
 309                 } else {
 310                         if(v < min)
 311                                 min = v;
 312                         i += k;
 313                 }
 314         }
 315         return min;
 316 }
 317
 318 /*
 319  * qsort specialized for sorting permutations based on successors
 320  */
 321 static void
 322 vecswap2(Elem *a, Elem *b, int n)
 323 {
 324         while (n-- > 0) {
 325                 Elem t = *a;
 326                 *a++ = *b;
 327                 *b++ = t;
 328         }
 329 }
 330
 331 #define swap2(a, b) { t = *(a); *(a) = *(b); *(b) = t; }
 332 #define ptr2char(i, asucc) (asucc[*(i)])
 333
 334 static Elem*
 335 med3(Elem *a, Elem *b, Elem *c, Elem *asucc)
 336 {
 337         Elem va, vb, vc;
 338
 339         if ((va=ptr2char(a, asucc)) == (vb=ptr2char(b, asucc)))
 340                 return a;
 341         if ((vc=ptr2char(c, asucc)) == va || vc == vb)
 342                 return c;
 343         return va < vb ?
 344                   (vb < vc ? b : (va < vc ? c : a))
 345                 : (vb > vc ? b : (va < vc ? a : c));
 346 }
 347
 348 static void
 349 inssort(Elem *a, Elem *asucc, int n)
 350 {
 351         Elem *pi, *pj, t;
 352
 353         for (pi = a + 1; --n > 0; pi++)
 354                 for (pj = pi; pj > a; pj--) {
 355                         if(ptr2char(pj-1, asucc) <= ptr2char(pj, asucc))
 356                                 break;
 357                         swap2(pj, pj-1);
 358                 }
 359 }
 360
 361 static void
 362 qsort2(Elem *a, Elem *asucc, int n)
 363 {
 364         Elem d, r, partval;
 365         Elem *pa, *pb, *pc, *pd, *pl, *pm, *pn, t;
 366
 367         if (n < 15) {
 368                 inssort(a, asucc, n);
 369                 return;
 370         }
 371         pl = a;
 372         pm = a + (n >> 1);
 373         pn = a + (n-1);
 374         if (n > 30) { /* On big arrays, pseudomedian of 9 */
 375                 d = (n >> 3);
 376                 pl = med3(pl, pl+d, pl+2*d, asucc);
 377                 pm = med3(pm-d, pm, pm+d, asucc);
 378                 pn = med3(pn-2*d, pn-d, pn, asucc);
 379         }
 380         pm = med3(pl, pm, pn, asucc);
 381         swap2(a, pm);
 382         partval = ptr2char(a, asucc);
 383         pa = pb = a + 1;
 384         pc = pd = a + n-1;
 385         for (;;) {
 386                 while (pb <= pc && (r = ptr2char(pb, asucc)-partval) <= 0) {
 387                         if (r == 0) {
 388                                 swap2(pa, pb);
 389                                 pa++;
 390                         }
 391                         pb++;
 392                 }
 393                 while (pb <= pc && (r = ptr2char(pc, asucc)-partval) >= 0) {
 394                         if (r == 0) {
 395                                 swap2(pc, pd);
 396                                 pd--;
 397                         }
 398                         pc--;
 399                 }
 400                 if (pb > pc)
 401                         break;
 402                 swap2(pb, pc);
 403                 pb++;
 404                 pc--;
 405         }
 406         pn = a + n;
 407         r = pa-a;
 408         if(pb-pa < r)
 409                 r = pb-pa;
 410         vecswap2(a, pb-r, r);
 411         r = pn-pd-1;
 412         if(pd-pc < r)
 413                 r = pd-pc;
 414         vecswap2(pb, pn-r, r);
 415         if ((r = pb-pa) > 1)
 416                 qsort2(a, asucc, r);
 417         if ((r = pd-pc) > 1)
 418                 qsort2(a + n-r, asucc, r);
 419 }