bp.c

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