bp.c

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