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