bp.c

   1 #include "bp.h"\r
   2 \r
   3 //#define CHECK\r
   4 #define RANDOM\r
   5 \r
   6 int msize=0;\r
   7 #define mymalloc(p,n,f) {p =(__typeof__(p)) malloc((n)*sizeof(*p)); msize += (f)*(n)*sizeof(*p); /* if (f) printf("malloc %d bytes at line %d total %d\n",(n)*sizeof(*p),__LINE__,msize);  */ if ((p)==NULL) {printf("not enough memory (%d bytes) in line %d\n",msize,__LINE__); exit(1);};}\r
   8 \r
   9 int postorder_select_bsearch(bp *b,int s);\r
  10 \r
  11 int naive_depth(bp *b, int s)\r
  12 {\r
  13   int i,d;\r
  14   if (s < 0) return 0;\r
  15   d = 0;\r
  16   for (i=0; i<=s; i++) {\r
  17     if (getbit(b->B,i)==OP) {\r
  18       d++;\r
  19     } else {\r
  20       d--;\r
  21     }\r
  22   }\r
  23   return d;\r
  24 }\r
  25 \r
  26 void printbp(bp *b, int s, int t)\r
  27 {\r
  28   int i,c,d;\r
  29   d = 0;\r
  30   for (i=s; i<=t; i++) {\r
  31     if (getbit(b->B,i)==OP) {\r
  32       c = '(';\r
  33       d++;\r
  34     } else {\r
  35       c = ')';\r
  36       d--;\r
  37     }\r
  38     putchar(c);\r
  39   }\r
  40 }\r
  41 \r
  42 int *matchtbl,*parenttbl;\r
  43 void make_naivetbl(pb *B,int n)\r
  44 {\r
  45   int i,v;\r
  46   int *stack,s;\r
  47 \r
  48   mymalloc(matchtbl,n,0);\r
  49   mymalloc(parenttbl,n,0);\r
  50   mymalloc(stack,n,0);\r
  51 \r
  52   for (i=0; i<n; i++) matchtbl[i] = parenttbl[i] = -1;\r
  53 \r
  54   s = 0;\r
  55   v = 0;\r
  56   stack[0] = -1;\r
  57   for (i=0; i<n; i++) {\r
  58     if (getbit(B,i)==OP) {\r
  59       v++;\r
  60       if (v > 0) {\r
  61         parenttbl[i] = stack[v-1];\r
  62         stack[v] = i;\r
  63       }\r
  64     } else {\r
  65       if (v > 0) {\r
  66         matchtbl[stack[v]] = i;  // close\r
  67         matchtbl[i] = stack[v];   // open\r
  68       }\r
  69       v--;\r
  70     }\r
  71   }\r
  72 \r
  73   free(stack);\r
  74 }\r
  75 \r
  76 int popCount[1<<ETW];\r
  77 int fwdtbl[(2*ETW+1)*(1<<ETW)];\r
  78 int bwdtbl[(2*ETW+1)*(1<<ETW)];\r
  79 #if 0\r
  80 int mintbl_li[1<<ETW], mintbl_lv[1<<ETW];\r
  81 int mintbl_ri[1<<ETW], mintbl_rv[1<<ETW];\r
  82 int maxtbl_li[1<<ETW], maxtbl_lv[1<<ETW];\r
  83 int maxtbl_ri[1<<ETW], maxtbl_rv[1<<ETW];\r
  84 #endif\r
  85 int minmaxtbl_i[4][1<<ETW], minmaxtbl_v[4][1<<ETW];\r
  86 int degtbl[1<<ETW];\r
  87 int degtbl2[(2*ETW+1)*(1<<ETW)];\r
  88 int childtbl[(ETW)*(1<<ETW)];\r
  89 int childtbl2[2*ETW+1][ETW][(1<<ETW)];\r
  90 int depthtbl[(2*ETW+1)*(1<<ETW)];\r
  91 \r
  92 void make_matchtbl(void)\r
  93 {\r
  94   int i,j,x,r;\r
  95   int m,M;\r
  96   pb buf[1];\r
  97   int deg;\r
  98 \r
  99   for (x = 0; x < (1<<ETW); x++) {\r
 100     setbits(buf,0,ETW,x);\r
 101     for (r=-ETW; r<=ETW; r++) fwdtbl[((r+ETW)<<ETW)+x] = ETW;\r
 102     for (r=-ETW; r<=ETW; r++) bwdtbl[((r+ETW)<<ETW)+x] = ETW;\r
 103     for (r=-ETW; r<=ETW; r++) degtbl2[((r+ETW)<<ETW)+x] = 0;\r
 104     for (r=-ETW; r<=ETW; r++) depthtbl[((r+ETW)<<ETW)+x] = 0;\r
 105 \r
 106     r = 0;\r
 107     for (i=0; i<ETW; i++) {\r
 108       if (getbit(buf,i)==OP) {\r
 109         r++;\r
 110       } else {\r
 111         r--;\r
 112       }\r
 113       if (fwdtbl[((r+ETW)<<ETW)+x] == ETW) fwdtbl[((r+ETW)<<ETW)+x] = i;\r
 114     }\r
 115 \r
 116     r = 0;\r
 117     for (i=ETW-1; i>=0; i--) {\r
 118       if (getbit(buf,i)==OP) {\r
 119         r--;\r
 120       } else {\r
 121         r++;\r
 122       }\r
 123       if (bwdtbl[((r+ETW)<<ETW)+x] == ETW) bwdtbl[((r+ETW)<<ETW)+x] = ETW-1-i;\r
 124     }\r
 125 \r
 126     r = 0;\r
 127     for (i=0; i<ETW; i++) {\r
 128       if (getbit(buf,i)==OP) {\r
 129         r++;\r
 130       } else {\r
 131         r--;\r
 132       }\r
 133       depthtbl[((r+ETW)<<ETW)+x] += (1<<(ETW-1));\r
 134     }\r
 135 \r
 136     r = 0;\r
 137     for (i=0; i<ETW; i++) {\r
 138       if (getbit(buf,i)==OP) r++;\r
 139     }\r
 140     popCount[x] = r;\r
 141 \r
 142     r = 0;  m = 0;  M = 0;\r
 143     m = ETW+1;  M = -ETW-1;\r
 144     //maxtbl_lv[x] = -ETW-1;\r
 145     //mintbl_lv[x] = ETW+1;\r
 146     minmaxtbl_v[OPT_MAX | OPT_LEFT][x] = -ETW-1;\r
 147     minmaxtbl_v[OPT_MIN | OPT_LEFT][x] = ETW+1;\r
 148     deg = 0;\r
 149     for (i=0; i<ETW; i++) {\r
 150       if (getbit(buf,i)==OP) {\r
 151         r++;\r
 152         if (r > M) {\r
 153           M = r;  \r
 154           //maxtbl_li[x] = i;  maxtbl_lv[x] = r;\r
 155           minmaxtbl_i[OPT_MAX | OPT_LEFT][x] = i;\r
 156           minmaxtbl_v[OPT_MAX | OPT_LEFT][x] = r;\r
 157         }\r
 158       } else {\r
 159         r--;\r
 160         if (r == m) {\r
 161           deg++;\r
 162           childtbl[((deg-1)<<ETW) + x] = i;\r
 163         }\r
 164         if (r < m) {\r
 165           m = r;\r
 166           //mintbl_li[x] = i;  mintbl_lv[x] = r;\r
 167           minmaxtbl_i[OPT_MIN | OPT_LEFT][x] = i;\r
 168           minmaxtbl_v[OPT_MIN | OPT_LEFT][x] = r;\r
 169           deg = 1;\r
 170           childtbl[((deg-1)<<ETW) + x] = i;\r
 171         }\r
 172       }\r
 173       if (r <= m) degtbl2[((r+ETW)<<ETW)+x]++;\r
 174     }\r
 175     degtbl[x] = deg;\r
 176 \r
 177     r = 0;  m = 0;  M = 0;\r
 178     //maxtbl_rv[x] = -ETW-1;\r
 179     //mintbl_rv[x] = ETW+1;\r
 180     minmaxtbl_v[OPT_MAX | OPT_RIGHT][x] = -ETW-1;\r
 181     minmaxtbl_v[OPT_MIN | OPT_RIGHT][x] = ETW+1;\r
 182     for (i=0; i<ETW; i++) {\r
 183       if (getbit(buf,i)==OP) {\r
 184         r++;\r
 185         if (r >= M) {\r
 186           M = r;  \r
 187           //maxtbl_ri[x] = i;  maxtbl_rv[x] = r;\r
 188           minmaxtbl_i[OPT_MAX | OPT_RIGHT][x] = i;\r
 189           minmaxtbl_v[OPT_MAX | OPT_RIGHT][x] = r;\r
 190         }\r
 191       } else {\r
 192         r--;\r
 193         if (r <= m) {\r
 194           m = r;\r
 195           //mintbl_ri[x] = i;  mintbl_rv[x] = r;\r
 196           minmaxtbl_i[OPT_MIN | OPT_RIGHT][x] = i;\r
 197           minmaxtbl_v[OPT_MIN | OPT_RIGHT][x] = r;\r
 198         }\r
 199       }\r
 200     }\r
 201 \r
 202     for (i = 0; i < ETW; i++) {\r
 203       for (j = -ETW; j <= ETW; j++) {\r
 204         childtbl2[j+ETW][i][x] = -1;\r
 205       }\r
 206     }\r
 207 \r
 208     for (j=-ETW; j<=ETW; j++) {\r
 209       int ith;\r
 210       ith = 0;\r
 211       r = 0;\r
 212       for (i = 0; i < ETW; i++) {\r
 213         if (getbit(buf,i)==OP) {\r
 214           r++;\r
 215         } else {\r
 216           r--;\r
 217           if (r < j) break;\r
 218           if (r == j) {\r
 219             ith++;\r
 220             childtbl2[j+ETW][ith-1][x] = i;\r
 221           }\r
 222         }\r
 223       }\r
 224     }\r
 225   }\r
 226 \r
 227 }\r
 228 \r
 229 \r
 230 int bp_construct(bp *b,int n, pb *B, int opt)\r
 231 {\r
 232   int i,j,d;\r
 233   int m,M,ds;\r
 234   int ns,nm;\r
 235   byte *sm, *sM;\r
 236   byte *sd;\r
 237   int *mm, *mM;\r
 238   int *md;\r
 239   int m_ofs;\r
 240   int r; // # of minimum values\r
 241 \r
 242 #if 0\r
 243   if (SB % D != 0) {\r
 244     printf("warning: SB=%d should be a multiple of D=%d\n",SB,D);\r
 245     // not necessarily?\r
 246   }\r
 247   if (SB % RRR != 0) {\r
 248     printf("warning: SB=%d should be a multiple of RRR=%d\n",SB,RRR);\r
 249   }\r
 250 #endif\r
 251 \r
 252   b->B = B;\r
 253   b->n = n;\r
 254   b->opt = opt;\r
 255   b->idx_size = 0;\r
 256   b->sm = NULL;\r
 257   b->sM = NULL;\r
 258   b->sd = NULL;\r
 259   b->mm = NULL;\r
 260   b->mM = NULL;\r
 261   b->md = NULL;\r
 262   b->da_leaf = NULL;\r
 263   b->da_inorder = NULL;\r
 264   b->da_postorder = NULL;\r
 265   b->da_dfuds_leaf = NULL;\r
 266   mymalloc(b->da,1,0);\r
 267   darray_construct(b->da,n,B, opt & OPT_FAST_PREORDER_SELECT);\r
 268   b->idx_size += b->da->idx_size;\r
 269   //Kim: comment this and the following, they polute the printing of the xpath library\r
 270   //printf("preorder rank/select table: %d bytes (%1.2f bpc)\n",b->da->idx_size,(double)b->da->idx_size*8/n);\r
 271 \r
 272   make_matchtbl();\r
 273 \r
 274   ns = (n+SB-1)/SB;\r
 275   mymalloc(sm, ns, 0);  b->idx_size += ns * sizeof(*sm);\r
 276   mymalloc(sM, ns, 0);  b->idx_size += ns * sizeof(*sM);\r
 277   b->sm = sm;\r
 278   b->sM = sM;\r
 279   if (opt & OPT_DEGREE) {\r
 280     mymalloc(sd, ns, 0);  b->idx_size += ns * sizeof(*sd);\r
 281     b->sd = sd;\r
 282     //printf("SB degree table: %d bytes (%1.2f bpc)\n",ns * sizeof(*sd), (double)ns * sizeof(*sd) * 8/n);\r
 283   }\r
 284   //printf("SB table: %d bytes (%1.2f bpc)\n",ns * sizeof(*sm) * 2, (double)ns * sizeof(*sm)*2 * 8/n);\r
 285 \r
 286   for (i=0; i<n; i++) {\r
 287     if (i % SB == 0) {\r
 288       ds = depth(b,i);\r
 289       m = M = ds;\r
 290       r = 1;\r
 291     } else {\r
 292       d = depth(b,i);\r
 293       if (d == m) r++;\r
 294       if (d < m) {\r
 295         m = d;\r
 296         r = 1;\r
 297       }\r
 298       if (d > M) M = d;\r
 299     }\r
 300     if (i % SB == SB-1 || i==n-1) {\r
 301       ds = depth(b,(i/SB)*SB-1);\r
 302       if (m - ds + SB < 0 || m - ds + SB > 255) {\r
 303         printf("error m=%d ds=%d\n",m,ds);\r
 304       }\r
 305       if (M - ds + 1 < 0 || M - ds + 1 > 255) {\r
 306         printf("error M=%d ds=%d\n",M,ds);\r
 307       }\r
 308       sm[i/SB] = m - ds + SB;\r
 309       sM[i/SB] = M - ds + 1;\r
 310       if (opt & OPT_DEGREE) sd[i/SB] = r;\r
 311     }\r
 312   }\r
 313 \r
 314 #if 0\r
 315   printf("sd: ");\r
 316   for (i=0;i<n/SB;i++) printf("%d ",sd[i]);\r
 317   printf("\n");\r
 318 #endif\r
 319 \r
 320 \r
 321   nm = (n+MB-1)/MB;\r
 322   m_ofs = 1 << blog(nm-1);\r
 323   b->m_ofs = m_ofs;\r
 324 \r
 325   mymalloc(mm, nm + m_ofs, 0);  b->idx_size += (nm+m_ofs) * sizeof(*mm);\r
 326   mymalloc(mM, nm + m_ofs, 0);  b->idx_size += (nm+m_ofs) * sizeof(*mM);\r
 327   b->mm = mm;\r
 328   b->mM = mM;\r
 329   if (opt & OPT_DEGREE) {\r
 330     mymalloc(md, nm + m_ofs, 0);  b->idx_size += (nm+m_ofs) * sizeof(*md);\r
 331     b->md = md;\r
 332     //printf("MB degree table: %d bytes (%1.2f bpc)\n",(nm+m_ofs) * sizeof(*md), (double)(nm+m_ofs) * sizeof(*md) * 8/n);\r
 333   }\r
 334   //printf("MB table: %d bytes (%1.2f bpc)\n",(nm+m_ofs) * sizeof(*mm) * 2, (double)(nm+m_ofs) * sizeof(*mm)*2 * 8/n);\r
 335 \r
 336   for (i=0; i<n; i++) {\r
 337     d = depth(b,i);\r
 338     if (i % MB == 0) {\r
 339       m = M = d;\r
 340       r = 1;\r
 341     } else {\r
 342       if (d == m) r++;\r
 343       if (d < m) {\r
 344         m = d;\r
 345         r = 1;\r
 346       }\r
 347       if (d > M) M = d;\r
 348     }\r
 349     if (i % MB == MB-1 || i==n-1) {\r
 350       mm[m_ofs+ i/MB] = m;\r
 351       mM[m_ofs+ i/MB] = M;\r
 352       if (opt & OPT_DEGREE) md[m_ofs+ i/MB] = r;\r
 353     }\r
 354   }\r
 355   \r
 356   for (j=m_ofs-1; j > 0; j--) {\r
 357     m = 0;\r
 358     if (j*2 < nm + m_ofs) m = mm[j*2];\r
 359     if (j*2+1 < nm + m_ofs) m = min(m,mm[j*2+1]);\r
 360     M = 0;\r
 361     if (j*2 < nm + m_ofs) M = mM[j*2];\r
 362     if (j*2+1 < nm + m_ofs) M = max(M,mM[j*2+1]);\r
 363     mm[j] = m;  mM[j] = M;\r
 364     if (opt & OPT_DEGREE) {\r
 365       d = 0;\r
 366       if (j*2 < nm + m_ofs) d = md[j*2];\r
 367       if (j*2+1 < nm + m_ofs) {\r
 368         if (mm[j*2] == mm[j*2+1]) d += md[j*2+1];\r
 369         if (mm[j*2] > mm[j*2+1]) d = md[j*2+1];\r
 370       }\r
 371       md[j] = d;\r
 372     }\r
 373   }\r
 374   mm[0] = -1;\r
 375   mM[0] = mM[1];\r
 376   if (opt & OPT_DEGREE) {\r
 377     md[0] = -1;\r
 378   }\r
 379 \r
 380 \r
 381 #if 0\r
 382   printf("md: ");\r
 383   for (i=0;i<m_ofs + n/MB;i++) printf("%d ",md[i]);\r
 384   printf("\n");\r
 385 #endif\r
 386 \r
 387   if (opt & OPT_LEAF) {\r
 388     mymalloc(b->da_leaf,1,0);\r
 389     darray_pat_construct(b->da_leaf, n, B, 2, 0x2, opt & OPT_FAST_LEAF_SELECT);\r
 390     //printf("leaf rank/select table: %d bytes (%1.2f bpc)\n",b->da_leaf->idx_size,(double)b->da_leaf->idx_size*8/n);\r
 391     b->idx_size += b->da_leaf->idx_size;  \r
 392   } else {\r
 393     b->da_leaf = NULL;\r
 394   }\r
 395 \r
 396   if (opt & OPT_INORDER) {\r
 397     mymalloc(b->da_inorder,1,0);\r
 398     darray_pat_construct(b->da_inorder, n, B, 2, 0x1, opt & OPT_FAST_INORDER_SELECT);\r
 399     //printf("inorder rank/select table: %d bytes (%1.2f bpc)\n",b->da_inorder->idx_size,(double)b->da_inorder->idx_size*8/n);\r
 400     b->idx_size += b->da_inorder->idx_size;\r
 401   } else {\r
 402     b->da_inorder = NULL;\r
 403   }\r
 404 \r
 405   if (opt & OPT_FAST_POSTORDER_SELECT) {\r
 406     mymalloc(b->da_postorder,1,0);\r
 407     darray_pat_construct(b->da_postorder, n, B, 1, 0x0, (opt & OPT_FAST_POSTORDER_SELECT) | OPT_NO_RANK);\r
 408     //printf("postorder rank/select table: %d bytes (%1.2f bpc)\n",b->da_postorder->idx_size,(double)b->da_postorder->idx_size*8/n);\r
 409     b->idx_size += b->da_postorder->idx_size;\r
 410   } else {\r
 411     b->da_postorder = NULL;\r
 412   }\r
 413 \r
 414   if (opt & OPT_DFUDS_LEAF) {\r
 415     mymalloc(b->da_dfuds_leaf,1,0);\r
 416     darray_pat_construct(b->da_dfuds_leaf, n, B, 2, 0x0, opt & OPT_FAST_DFUDS_LEAF_SELECT);\r
 417     //printf("dfuds leaf rank/select table: %d bytes (%1.2f bpc)\n",b->da_dfuds_leaf->idx_size,(double)b->da_dfuds_leaf->idx_size*8/n);\r
 418     b->idx_size += b->da_dfuds_leaf->idx_size;\r
 419   } else {\r
 420     b->da_dfuds_leaf = NULL;\r
 421   }\r
 422 \r
 423   return 0;\r
 424 }\r
 425 \r
 426 // destroyTree: frees the memory of tree.\r
 427 void destroyTree(bp *b) {\r
 428    if (!b) return; // nothing to free\r
 429 \r
 430    destroyDarray(b->da); // destroys da data structure\r
 431    if (b->da) free(b->da);     \r
 432    \r
 433    if (b->sm) free(b->sm);\r
 434    if (b->sM) free(b->sM);\r
 435    if (b->sd) free(b->sd);\r
 436    if (b->mm) free(b->mm);\r
 437    if (b->mM) free(b->mM);\r
 438    if (b->md) free(b->md);\r
 439    \r
 440    destroyDarray(b->da_leaf);\r
 441    if (b->da_leaf) free(b->da_leaf);\r
 442    \r
 443    destroyDarray(b->da_inorder);\r
 444    if (b->da_inorder) free(b->da_inorder);\r
 445    \r
 446    destroyDarray(b->da_postorder);\r
 447    if (b->da_postorder) free(b->da_postorder);\r
 448    \r
 449    destroyDarray(b->da_dfuds_leaf);\r
 450    if (b->da_dfuds_leaf) free(b->da_dfuds_leaf);\r
 451 }\r
 452 \r
 453 \r
 454 // saveTree: saves parentheses data structure to file\r
 455 // By Diego Arroyuelo\r
 456 void saveTree(bp *b, FILE *fp) {\r
 457 \r
 458    if (fwrite(&(b->n), sizeof(int), 1, fp) != 1) {\r
 459       printf("Error: cannot save number of parentheses to file\n");\r
 460       exit(1);\r
 461    }\r
 462 \r
 463    if (fwrite(b->B, sizeof(pb), (b->n+D-1)/D, fp) != ((b->n+D-1)/D)) {\r
 464       printf("Error: cannot save parentheses sequence to file\n");\r
 465       exit(1);\r
 466    }\r
 467 \r
 468    if (fwrite(&(b->opt), sizeof(int), 1, fp) != 1) {\r
 469       printf("Error: cannot save opt in parentheses to file\n");\r
 470       exit(1);\r
 471    }\r
 472 }\r
 473 \r
 474 // loadTree: load parentheses data structure from file\r
 475 // By Diego Arroyuelo\r
 476 void loadTree(bp *b, FILE *fp) {\r
 477 \r
 478    pb *B;\r
 479    int n, opt;\r
 480 \r
 481    if (fread(&n, sizeof(int), 1, fp) != 1) {\r
 482       printf("Error: cannot read number of parentheses from file\n");\r
 483       exit(1);\r
 484    }\r
 485 \r
 486    mymalloc(B,(n+D-1)/D,0);\r
 487    \r
 488    if (fread(B, sizeof(pb), (n+D-1)/D, fp) != ((n+D-1)/D)) {\r
 489       printf("Error: cannot read parentheses sequence from file\n");\r
 490       exit(1);\r
 491    }\r
 492 \r
 493    if (fread(&opt, sizeof(int), 1, fp) != 1) {\r
 494       printf("Error: cannot read opt in parentheses from file\n");\r
 495       exit(1);\r
 496    }\r
 497    \r
 498    bp_construct(b, n, B, opt); \r
 499    \r
 500 }\r
 501 \r
 502 \r
 503 \r
 504 int naive_fwd_excess(bp *b,int s, int rel)\r
 505 {\r
 506   int i,v,n;\r
 507   pb *B;\r
 508   n = b->n;  B = b->B;\r
 509   v = 0;\r
 510   for (i=s+1; i<n; i++) {\r
 511     if (getbit(B,i)==OP) {\r
 512       v++;\r
 513     } else {\r
 514       v--;\r
 515     }\r
 516     if (v == rel) return i;\r
 517   }\r
 518   return -1;\r
 519 }\r
 520 \r
 521 int naive_bwd_excess(bp *b,int s, int rel)\r
 522 {\r
 523   int i,v;\r
 524   pb *B;\r
 525   B = b->B;\r
 526   v = 0;\r
 527   for (i=s; i>=0; i--) {\r
 528     if (getbit(B,i)==OP) {\r
 529       v--;\r
 530     } else {\r
 531       v++;\r
 532     }\r
 533     if (v == rel) return i-1;\r
 534   }\r
 535   return -2;\r
 536 }\r
 537 \r
 538 int naive_search_SB_l(bp *b, int i, int rel)\r
 539 {\r
 540   int il,v;\r
 541 \r
 542   il = (i / SB) * SB;\r
 543   for (; i>=il; i--) {\r
 544     if (getbit(b->B,i)==OP) {\r
 545       rel++;\r
 546     } else {\r
 547       rel--;\r
 548     }\r
 549     if (rel == 0) return i-1;\r
 550   }\r
 551   if (i < 0) return -2;\r
 552   return -3;\r
 553 }\r
 554 \r
 555 int naive_rmq(bp *b, int s, int t,int opt)\r
 556 {\r
 557   int d,i,dm,im;\r
 558 \r
 559   if (opt & OPT_RIGHT) {\r
 560     d = dm = depth(b,t);  im = t;\r
 561     i = t-1;\r
 562     while (i >= s) {\r
 563       if (getbit(b->B,i+1)==CP) {\r
 564         d++;\r
 565         if (opt & OPT_MAX) {\r
 566           if (d > dm) {\r
 567             dm = d;  im = i;\r
 568           }\r
 569         }\r
 570       } else {\r
 571         d--;\r
 572         if (!(opt & OPT_MAX)) {\r
 573           if (d < dm) {\r
 574             dm = d;  im = i;\r
 575           }\r
 576         }\r
 577       }\r
 578       i--;\r
 579     }\r
 580   } else {\r
 581     d = dm = depth(b,s);  im = s;\r
 582     i = s+1;\r
 583     while (i <= t) {\r
 584       if (getbit(b->B,i)==OP) {\r
 585         d++;\r
 586         if (opt & OPT_MAX) {\r
 587           if (d > dm) {\r
 588             dm = d;  im = i;\r
 589           }\r
 590         }\r
 591       } else {\r
 592         d--;\r
 593         if (!(opt & OPT_MAX)) {\r
 594           if (d < dm) {\r
 595             dm = d;  im = i;\r
 596           }\r
 597         }\r
 598       }\r
 599       i++;\r
 600     }\r
 601   }\r
 602   return im;\r
 603 }\r
 604 \r
 605 int root_node(bp *b)\r
 606 {\r
 607   return 0;\r
 608 }\r
 609 \r
 610 \r
 611 int rank_open(bp *b, int s)\r
 612 {\r
 613   return darray_rank(b->da,s);\r
 614 }\r
 615 \r
 616 int rank_close(bp *b, int s)\r
 617 {\r
 618   return s+1 - darray_rank(b->da,s);\r
 619 }\r
 620 \r
 621 int select_open(bp *b, int s)\r
 622 {\r
 623   if (b->opt & OPT_FAST_PREORDER_SELECT) {\r
 624     return darray_select(b->da,s,1);\r
 625   } else {\r
 626     return darray_select_bsearch(b->da,s,getpat_preorder);\r
 627   }\r
 628 }\r
 629 \r
 630 int select_close(bp *b, int s)\r
 631 {\r
 632   if (b->opt & OPT_FAST_POSTORDER_SELECT) {\r
 633     return darray_pat_select(b->da_postorder,s,getpat_postorder);\r
 634   } else {\r
 635     return postorder_select_bsearch(b,s);\r
 636   }\r
 637 }\r
 638 \r
 639 ///////////////////////////////////////////\r
 640 //  find_close(bp *b,int s)\r
 641 //    returns the matching close parenthesis of s\r
 642 ///////////////////////////////////////////\r
 643 int find_close(bp *b,int s)\r
 644 {\r
 645   return fwd_excess(b,s,-1);\r
 646 }\r
 647 \r
 648 ///////////////////////////////////////////\r
 649 //  find_open(bp *b,int s)\r
 650 //    returns the matching open parenthesis of s\r
 651 ///////////////////////////////////////////\r
 652 int find_open(bp *b,int s)\r
 653 {\r
 654   int r;\r
 655   r = bwd_excess(b,s,0);\r
 656   if (r >= -1) return r+1;\r
 657   return -1;\r
 658 }\r
 659 \r
 660 ///////////////////////////////////////////\r
 661 //  parent(bp *b,int s)\r
 662 //    returns the parent of s\r
 663 //            -1 if s is the root\r
 664 ///////////////////////////////////////////\r
 665 int parent(bp *b,int s)\r
 666 {\r
 667   int r;\r
 668   r = bwd_excess(b,s,-2);\r
 669   if (r >= -1) return r+1;\r
 670   return -1;\r
 671 }\r
 672 \r
 673 int enclose(bp *b,int s)\r
 674 {\r
 675   return parent(b,s);\r
 676 }\r
 677 \r
 678 ///////////////////////////////////////////\r
 679 //  level_ancestor(bp *b,int s,int d)\r
 680 //    returns the ancestor of s with relative depth d (d < 0)\r
 681 //            -1 if no such node\r
 682 ///////////////////////////////////////////\r
 683 int level_ancestor(bp *b,int s,int d)\r
 684 {\r
 685   int r;\r
 686   r = bwd_excess(b,s,d-1);\r
 687   if (r >= -1) return r+1;\r
 688   return -1;\r
 689 }\r
 690 \r
 691 ///////////////////////////////////////////\r
 692 //  lca(bp *b, int s, int t)\r
 693 //    returns the lowest common ancestor of s and t\r
 694 ///////////////////////////////////////////\r
 695 int lca(bp *b, int s, int t)\r
 696 {\r
 697   return parent(b,rmq(b,s,t,0)+1);\r
 698 }\r
 699 \r
 700 \r
 701 ///////////////////////////////////////////\r
 702 //  preorder_rank(bp *b,int s)\r
 703 //    returns the preorder (>= 1) of node s (s >= 0)\r
 704 ///////////////////////////////////////////\r
 705 int preorder_rank(bp *b,int s)\r
 706 {\r
 707   return darray_rank(b->da,s);\r
 708 }\r
 709 \r
 710 ///////////////////////////////////////////\r
 711 //  preorder_select(bp *b,int s)\r
 712 //    returns the node with preorder s (s >= 1)\r
 713 //            -1 if no such node\r
 714 ///////////////////////////////////////////\r
 715 int preorder_select(bp *b,int s)\r
 716 {\r
 717   //  no error handling\r
 718   if (b->opt & OPT_FAST_PREORDER_SELECT) {\r
 719     return darray_select(b->da,s,1);\r
 720   } else {\r
 721     return darray_select_bsearch(b->da,s,getpat_preorder);\r
 722   }\r
 723 }\r
 724 \r
 725 ///////////////////////////////////////////\r
 726 //  postorder_rank(bp *b,int s)\r
 727 //    returns the postorder (>= 1) of node s (s >= 0)\r
 728 //            -1 if s-th bit is not OP\r
 729 ///////////////////////////////////////////\r
 730 int postorder_rank(bp *b,int s)\r
 731 {\r
 732   int t;\r
 733   if (inspect(b,s) == CP) return -1;\r
 734   t = find_close(b,s);\r
 735   //  return t+1 - darray_rank(b->da,t);\r
 736   return rank_close(b,t);\r
 737 }\r
 738 \r
 739 int postorder_select_bsearch(bp *b,int s)\r
 740 {\r
 741   int l,r,m;\r
 742 \r
 743   if (s == 0) return -1;\r
 744 \r
 745   if (s > b->da->n - b->da->m) {\r
 746     return -1;\r
 747   }\r
 748   l = 0;  r = b->da->n - 1;\r
 749 \r
 750   while (l < r) {\r
 751     m = (l+r)/2;\r
 752     //printf("m=%d rank=%d s=%d\n",m,m+1 - darray_rank(b->da,m),s);\r
 753     if (m+1 - darray_rank(b->da,m) >= s) {\r
 754       r = m;\r
 755     } else {\r
 756       l = m+1;\r
 757     }\r
 758   }\r
 759   return l;\r
 760 }\r
 761 \r
 762 ///////////////////////////////////////////\r
 763 //  postorder_select(bp *b,int s)\r
 764 //    returns the position of CP of the node with postorder s (>= 1)\r
 765 ///////////////////////////////////////////\r
 766 int postorder_select(bp *b,int s)\r
 767 {\r
 768 #if 0\r
 769   if (b->opt & OPT_FAST_POSTORDER_SELECT) {\r
 770     return darray_pat_select(b->da_postorder,s,getpat_postorder);\r
 771   } else {\r
 772     return postorder_select_bsearch(b->da,s);\r
 773   }\r
 774 #else\r
 775   return select_close(b,s);\r
 776 #endif\r
 777 }\r
 778 \r
 779 ///////////////////////////////////////////\r
 780 //  leaf_rank(bp *b,int s)\r
 781 //    returns the number of leaves to the left of s\r
 782 ///////////////////////////////////////////\r
 783 int leaf_rank(bp *b,int s)\r
 784 {\r
 785   if ((b->opt & OPT_LEAF) == 0) {\r
 786     printf("leaf_rank: error!!!   not supported\n");\r
 787     return -1;\r
 788   }\r
 789   if (s >= b->n-1) {\r
 790     s = b->n-2;\r
 791   }\r
 792   return darray_pat_rank(b->da_leaf,s,getpat_leaf);\r
 793 }\r
 794 \r
 795 ///////////////////////////////////////////\r
 796 //  leaf_select(bp *b,int s)\r
 797 //    returns the position of s-th leaf\r
 798 ///////////////////////////////////////////\r
 799 int leaf_select(bp *b,int s)\r
 800 {\r
 801   if ((b->opt & OPT_LEAF) == 0) {\r
 802     printf("leaf_select: error!!!   not supported\n");\r
 803     return -1;\r
 804   }\r
 805   if (s > b->da_leaf->m) return -1;\r
 806   if (b->opt & OPT_FAST_LEAF_SELECT) {\r
 807     return darray_pat_select(b->da_leaf,s,getpat_leaf);\r
 808   } else {\r
 809     return darray_select_bsearch(b->da_leaf,s,getpat_leaf);\r
 810   }\r
 811 }\r
 812 \r
 813 \r
 814 ///////////////////////////////////////////\r
 815 //  inorder_rank(bp *b,int s)\r
 816 //    returns the number of ")("  (s >= 0)\r
 817 ///////////////////////////////////////////\r
 818 int inorder_rank(bp *b,int s)\r
 819 {\r
 820   if ((b->opt & OPT_INORDER) == 0) {\r
 821     printf("inorder_rank: error!!!   not supported\n");\r
 822     return -1;\r
 823   }\r
 824   if (s >= b->n-1) {\r
 825     s = b->n-2;\r
 826   }\r
 827   return darray_pat_rank(b->da_inorder,s,getpat_inorder);\r
 828 }\r
 829 \r
 830 ///////////////////////////////////////////\r
 831 //  inorder_select(bp *b,int s)\r
 832 //    returns the s-th position of ")("  (s >= 1)\r
 833 ///////////////////////////////////////////\r
 834 int inorder_select(bp *b,int s)\r
 835 {\r
 836   if ((b->opt & OPT_INORDER) == 0) {\r
 837     printf("inorder_select: error!!!   not supported\n");\r
 838     return -1;\r
 839   }\r
 840   if (b->opt & OPT_FAST_INORDER_SELECT) {\r
 841     return darray_pat_select(b->da_inorder,s,getpat_inorder);\r
 842   } else {\r
 843     return darray_select_bsearch(b->da_inorder,s,getpat_inorder);\r
 844   }\r
 845 }\r
 846 \r
 847 ///////////////////////////////////////////\r
 848 //  leftmost_leaf(bp *b, int s)\r
 849 ///////////////////////////////////////////\r
 850 int leftmost_leaf(bp *b, int s)\r
 851 {\r
 852   if ((b->opt & OPT_LEAF) == 0) {\r
 853     printf("leftmost_leaf: error!!!   not supported\n");\r
 854     return -1;\r
 855   }\r
 856   return leaf_select(b,leaf_rank(b,s)+1);\r
 857 }\r
 858 \r
 859 ///////////////////////////////////////////\r
 860 //  rightmost_leaf(bp *b, int s)\r
 861 ///////////////////////////////////////////\r
 862 int rightmost_leaf(bp *b, int s)\r
 863 {\r
 864   int t;\r
 865   if ((b->opt & OPT_LEAF) == 0) {\r
 866     printf("leftmost_leaf: error!!!   not supported\n");\r
 867     return -1;\r
 868   }\r
 869   t = find_close(b,s);\r
 870   return leaf_select(b,leaf_rank(b,t));\r
 871 }\r
 872 \r
 873 \r
 874 \r
 875 ///////////////////////////////////////////\r
 876 //  inspect(bp *b, int s)\r
 877 //    returns OP (==1) or CP (==0) at s-th bit (0 <= s < n)\r
 878 ///////////////////////////////////////////\r
 879 int inspect(bp *b, int s)\r
 880 {\r
 881   if (s < 0 || s >= b->n) {\r
 882     printf("inspect: error s=%d is out of [%d,%d]\n",s,0,b->n-1);\r
 883   }\r
 884   return getbit(b->B,s);\r
 885 }\r
 886 \r
 887 int isleaf(bp *b, int s)\r
 888 {\r
 889   if (inspect(b,s) != OP) {\r
 890     printf("isleaf: error!!! B[%d] = OP\n",s);\r
 891   }\r
 892   if (inspect(b,s+1) == CP) return 1;\r
 893   else return 0;\r
 894 }\r
 895 \r
 896 \r
 897 ///////////////////////////////////////////\r
 898 //  subtree_size(bp *b, int s)\r
 899 //    returns the number of nodes in the subtree of s\r
 900 ///////////////////////////////////////////\r
 901 int subtree_size(bp *b, int s)\r
 902 {\r
 903   return (find_close(b,s) - s + 1) / 2;\r
 904 }\r
 905 \r
 906 ///////////////////////////////////////////\r
 907 //  first_child(bp *b, int s)\r
 908 //    returns the first child\r
 909 //            -1 if s is a leaf\r
 910 ///////////////////////////////////////////\r
 911 int first_child(bp *b, int s)\r
 912 {\r
 913   if (inspect(b,s+1) == CP) return -1;\r
 914   return s+1;\r
 915 }\r
 916 \r
 917 ///////////////////////////////////////////\r
 918 //  next_sibling(bp *b,int s)\r
 919 //    returns the next sibling of parent(s)\r
 920 //            -1 if s is the last child\r
 921 //////////////////////////////////////////\r
 922 int next_sibling(bp *b, int s)\r
 923 {\r
 924   int t;\r
 925   t = find_close(b,s)+1;\r
 926   if (t >= b->n) {\r
 927     printf("next_sibling: error s=%d t=%d\n",s,t);\r
 928   }\r
 929   if (inspect(b,t) == CP) return -1;\r
 930   return t;\r
 931 }\r
 932 \r
 933 ///////////////////////////////////////////\r
 934 //  prev_sibling(bp *b,int s)\r
 935 //    returns the previous sibling of parent(s)\r
 936 //            -1 if s is the first child\r
 937 //////////////////////////////////////////\r
 938 int prev_sibling(bp *b, int s)\r
 939 {\r
 940   int t;\r
 941   if (s < 0) {\r
 942     printf("prev_sibling: error s=%d\n",s);\r
 943   }\r
 944   if (s == 0) return -1;\r
 945   if (inspect(b,s-1) == OP) return -1;\r
 946   t = find_open(b,s-1);\r
 947   return t;\r
 948 }\r
 949 \r
 950 ///////////////////////////////////////////\r
 951 //  deepest_node(bp *b,int s)\r
 952 //    returns the first node with the largest depth in the subtree of s\r
 953 ///////////////////////////////////////////\r
 954 int deepest_node(bp *b,int s)\r
 955 {\r
 956   int t,m;\r
 957   t = find_close(b,s);\r
 958   m = rmq(b,s,t, OPT_MAX);\r
 959   return m;\r
 960 }\r
 961 \r
 962 ///////////////////////////////////////////\r
 963 //  subtree_height(bp *b,int s)\r
 964 //    returns the height of the subtree of s\r
 965 //            0 if s is a leaf\r
 966 ///////////////////////////////////////////\r
 967 int subtree_height(bp *b,int s)\r
 968 {\r
 969   int t;\r
 970   t = deepest_node(b,s);\r
 971   return depth(b,t) - depth(b,s);\r
 972 }\r
 973 \r
 974 int naive_degree(bp *b, int s)\r
 975 {\r
 976   int t,d;\r
 977   d = 0;\r
 978   t = first_child(b,s);\r
 979   while (t >= 0) {\r
 980     d++;\r
 981     t = next_sibling(b,t);\r
 982   }\r
 983   return d;\r
 984 }\r
 985 \r
 986 ///////////////////////////////////////////\r
 987 //  degree(bp *b, int s)\r
 988 //    returns the number of children of s\r
 989 //            0 if s is a leaf\r
 990 ///////////////////////////////////////////\r
 991 int degree(bp *b, int s)\r
 992 {\r
 993   if (b->opt & OPT_DEGREE) {\r
 994     return fast_degree(b,s,b->n,0);\r
 995   } else {\r
 996     return naive_degree(b,s);\r
 997   }\r
 998 }\r
 999 \r
1000 int naive_child(bp *b, int s, int d)\r
1001 {\r
1002   int t,i;\r
1003   t = first_child(b,s);\r
1004   for (i = 1; i < d; i++) {\r
1005     if (t == -1) break;\r
1006     t = next_sibling(b,t);\r
1007   }\r
1008   return t;\r
1009 }\r
1010 \r
1011 ///////////////////////////////////////////\r
1012 //  child(bp *b, int s, int d)\r
1013 //    returns the d-th child of s  (1 <= d <= degree(s))\r
1014 //            -1 if no such node\r
1015 ///////////////////////////////////////////\r
1016 int child(bp *b, int s, int d)\r
1017 {\r
1018   int r;\r
1019   if (b->opt & OPT_DEGREE) {\r
1020     //return find_open(b,fast_degree(b,s,b->n,d));\r
1021     if (d==1) return first_child(b,s);\r
1022     r = fast_degree(b,s,b->n,d-1)+1;\r
1023     if (inspect(b,r) == CP) return -1;\r
1024     return r;\r
1025   } else {\r
1026     return naive_child(b,s,d);\r
1027   }\r
1028     \r
1029 }\r
1030 \r
1031 \r
1032 int naive_child_rank(bp *b, int t)\r
1033 {\r
1034   int v,d;\r
1035   d = 0;\r
1036   while (t != -1) {\r
1037     d++;\r
1038     t = prev_sibling(b,t);\r
1039   }\r
1040   return d;\r
1041 }\r
1042 \r
1043 ///////////////////////////////////////////\r
1044 //  child_rank(bp *b, int t)\r
1045 //    returns d if t is the d-th child of the parent of t (d >= 1)\r
1046 //            1 if t is the root\r
1047 ///////////////////////////////////////////\r
1048 int child_rank(bp *b, int t)\r
1049 {\r
1050   int r;\r
1051   if (t == root_node(b)) return 1;\r
1052   if (b->opt & OPT_DEGREE) {\r
1053     r = parent(b,t);\r
1054     return fast_degree(b,r,t,0)+1;\r
1055   } else {\r
1056     return naive_child_rank(b,t);\r
1057   }\r
1058 }\r
1059 \r
1060 \r
1061 \r
1062 ///////////////////////////////////////////\r
1063 //  is_ancestor(bp *b, int s, int t)\r
1064 //     returns 1 if s is an ancestor of t\r
1065 //             0 otherwise\r
1066 ///////////////////////////////////////////\r
1067 int is_ancestor(bp *b, int s, int t)\r
1068 {\r
1069   int v;\r
1070   v = find_close(b,s);\r
1071   if (s <= t && t <= v) return 1;\r
1072   return 0;\r
1073 }\r
1074 \r
1075 ///////////////////////////////////////////\r
1076 //  distance(bp *b, int s, int t)\r
1077 //    returns the length of the shortest path from s to t in the tree\r
1078 ///////////////////////////////////////////\r
1079 int distance(bp *b, int s, int t)\r
1080 {\r
1081   int v,d;\r
1082   v = lca(b,s,t);\r
1083   d = depth(b,v);\r
1084   return (depth(b,s) - d) + (depth(b,t) - d);\r
1085 }\r
1086 \r
1087 ///////////////////////////////////////////\r
1088 //  level_next(bp *b, int d)\r
1089 ///////////////////////////////////////////\r
1090 int level_next(bp *b,int s)\r
1091 {\r
1092   int t;\r
1093   t = fwd_excess(b,s,0);\r
1094   return t;\r
1095 }\r
1096 \r
1097 ///////////////////////////////////////////\r
1098 //  level_prev(bp *b, int d)\r
1099 ///////////////////////////////////////////\r
1100 int level_prev(bp *b,int s)\r
1101 {\r
1102   int t;\r
1103   t = bwd_excess(b,s,0);\r
1104   return t;\r
1105 }\r
1106 \r
1107 ///////////////////////////////////////////\r
1108 //  level_leftmost(bp *b, int d)\r
1109 ///////////////////////////////////////////\r
1110 int level_leftmost(bp *b, int d)\r
1111 {\r
1112   int t;\r
1113   if (d < 1) return -1;\r
1114   if (d == 1) return 0;\r
1115   t = fwd_excess(b,0,d);\r
1116   return t;\r
1117 }\r
1118 \r
1119 ///////////////////////////////////////////\r
1120 //  level_rigthmost(bp *b, int d)\r
1121 ///////////////////////////////////////////\r
1122 int level_rigthmost(bp *b, int d)\r
1123 {\r
1124   int t;\r
1125   if (d < 1) return -1;\r
1126   if (d == 1) return 0;\r
1127   t = bwd_excess(b,0,d-1);\r
1128   return find_open(b,t);\r
1129 }\r
1130 \r
1131 ///////////////////////////////////////////\r
1132 //  leaf_size(bp *b, int s)\r
1133 ///////////////////////////////////////////\r
1134 int leaf_size(bp *b, int s)\r
1135 {\r
1136   int t;\r
1137   if ((b->opt & OPT_LEAF) == 0) {\r
1138     printf("leaf_size: error!!!   not supported\n");\r
1139     return -1;\r
1140   }\r
1141   t = find_close(b,s);\r
1142   return leaf_rank(b,t) - leaf_rank(b,s);\r
1143 }\r