bp-core.c

   1 #include "bp.h"
   2 #include "bp-utils.h"
   3
   4
   5 #define NOTFOUND -2
   6 #define CONTINUE -3
   7 #define END -4
   8 #define FOUND -5
   9
  10 #define SBid(i) ((i)>>logSB)
  11 #define SBfirst(i) ((i) & (~(SB-1)))
  12 #define SBlast(i) ((i) | (SB-1))
  13
  14 #define MBid(i) ((i)>>logMB)
  15 #define MBfirst(i) ((i) & (~(MB-1)))
  16 #define MBlast(i) ((i) | (MB-1))
  17 #define max(a,b) \
  18   ({ __typeof__ (a) _a = (a); \
  19   __typeof__ (b) _b = (b); \
  20   _a > _b ? _a : _b; })
  21
  22
  23 #define min(a,b) \
  24   ({ __typeof__ (a) _a = (a); \
  25   __typeof__ (b) _b = (b); \
  26    _a <= _b ? _a : _b; })
  27
  28
  29 pb getpat_preorder(pb *b)
  30 {
  31   return *b;
  32 }
  33
  34 pb getpat_postorder(pb *b)
  35 {
  36   return ~(*b);
  37 }
  38
  39 pb getpat_leaf(pb *b)
  40 {
  41   pb w1,w2,w;
  42   w1 = b[0];
  43   w2 = (w1 << 1) + (b[1] >> (D-1));
  44   w = w1 & (~w2);
  45   return w;
  46 }
  47
  48 pb getpat_inorder(pb *b)
  49 {
  50   pb w1,w2,w;
  51   w1 = b[0];
  52   w2 = (w1 << 1) + (b[1] >> (D-1));
  53   w = (~w1) & w2;
  54   return w;
  55 }
  56
  57 pb getpat_dfuds_leaf(pb *b)
  58 {
  59   pb w1,w2,w;
  60   w1 = b[0];
  61   w2 = (w1 << 1) + (b[1] >> (D-1));
  62   w = (~w1) & (~w2);
  63   return w;
  64 }
  65
  66
  67
  68 ///////////////////////////////////////////
  69 //  bp_depth(bp *b, int s)
  70 //    returns the depth of s
  71 //  The root node has depth 1
  72 ///////////////////////////////////////////
  73 int bp_depth(bp *b, int s)
  74 {
  75   int d;
  76   if (s < 0) return 0;
  77   d = 2 * bp_darray_rank(b->da,s) - (s+1);
  78   return d;
  79 }
  80
  81 int fast_depth(bp *b, int s)
  82 {
  83   int d;
  84   darray *da;
  85   int r,j;
  86
  87   s++;
  88   if ((s & (RRR-1)) != 0) {
  89     //printf("fast_depth:warning s=%d\n",s);
  90     return bp_depth(b,s);
  91   }
  92   da = b->da;
  93   //d = 2 * (da->rl[s>>logR] + da->rm[s>>logRR] + da->rs[s>>logRRR]) - s;
  94
  95   r = da->rl[s>>logR] + da->rm[s>>logRR];
  96   j = (s>>logRRR) & (RR/RRR-1);
  97   while (j > 0) {
  98     r += da->rs[((s>>logRR)<<(logRR-logRRR))+j-1];
  99     j--;
 100   }
 101   d = 2 * r - s;
 102
 103   return d;
 104 }
 105
 106 int search_SB_r(bp *b, int i, int rel)
 107 {
 108   int j,r,n,il;
 109   pb *p,x,w;
 110
 111   n = b->n;
 112   il = min((SBid(i) + 1) << logSB,n);
 113   p = &b->B[i>>logD];
 114   while (i<il) {
 115     x = *p++;
 116     j = i & (D-1);
 117     x <<= j;
 118     j = D-j;
 119     while (j>0) {
 120       w = (x >> (D-ETW)) & ((1<<ETW)-1);
 121       if (rel >= -ETW && rel <= ETW) {
 122         r = fwdtbl[((rel+ETW)<<ETW)+w];
 123         if (r<ETW && r<j) {
 124           if (i+r >= n) return NOTFOUND;
 125           return i+r;
 126         }
 127       }
 128       r = min(j,ETW);
 129       rel -= 2*popcount(w)-r;
 130       x <<= r;
 131       i += r;
 132       j -= r;
 133     }
 134   }
 135   return CONTINUE;
 136 }
 137
 138 int search_MB_r(bp *b, int i, int td)
 139 {
 140   int il,d;
 141   int m,M,n;
 142   pb *B;
 143
 144   B = b->B;
 145   n = b->n;
 146
 147   il = min((MBid(i) + 1) << logMB, n);
 148   for (  ;  i < il;  i+=SB) {
 149 #if (SB % RRR != 0)
 150     d = bp_depth(b,i-1);
 151 #else
 152     d = fast_depth(b,i-1);
 153 #endif
 154     m = d + b->sm[SBid(i)] - SB;
 155     M = d + b->sM[SBid(i)] - 1;
 156     if (m <= td && td <= M) {
 157       return search_SB_r(b,i,td-d);
 158     }
 159   }
 160   if (i >= n) return NOTFOUND;
 161   return CONTINUE;
 162 }
 163
 164 ///////////////////////////////////////////
 165 //  bp_fwd_excess(bp *b,int s, int rel)
 166 //    find the leftmost value depth(s)+rel to the right of s (exclusive)
 167 ///////////////////////////////////////////
 168 int bp_fwd_excess(bp *b,int s, int rel)
 169 {
 170   int i,n;
 171   int d,td;
 172   int m,M;
 173   int m_ofs;
 174   pb *B;
 175
 176   n = b->n;  B = b->B;
 177
 178   i = s+1;
 179
 180   d = search_SB_r(b,i,rel);
 181   if (d >= NOTFOUND) return d;
 182
 183   i = min((SBid(i) + 1) << logSB, n);
 184   td = bp_depth(b,s) + rel;
 185   d = search_MB_r(b,i,td);
 186   if (d >= NOTFOUND) return d;
 187
 188   m_ofs = b->m_ofs;
 189
 190   i = MBid(s) + m_ofs;
 191
 192   while (i > 0) {
 193     if ((i&1) == 0) {
 194       i++;
 195       m = b->mm[i];
 196       M = b->mM[i];
 197
 198       if (m <= td && td <= M) {
 199               while (i < m_ofs) {
 200                       i <<= 1;
 201                       m = b->mm[i];
 202                       M = b->mM[i];
 203                       if (!(m <= td && td <= M)) i++;
 204               }
 205               i -= m_ofs;
 206               i <<= logMB;
 207
 208               return search_MB_r(b,i,td);
 209       };
 210
 211     }
 212     i >>= 1;
 213   }
 214   return NOTFOUND;
 215 }
 216
 217
 218 int degree_SB_slow(bp *b, int i, int t, int rel, int *ans, int ith)
 219 {
 220   int j,r,n,il;
 221   pb *p,x,w,w2;
 222   int d, deg, v;
 223
 224   n = t;
 225   il = min((SBid(i) + 1) << logSB,n);
 226   d = deg = 0;
 227
 228   while (i < il) {
 229     if (bp_getbit(b->B,i)==OP) {
 230       d++;
 231     } else {
 232       d--;
 233     }
 234     if (d < rel) {  // reached the end
 235       if (ith > 0) {
 236         return NOTFOUND;
 237       } else {
 238         *ans = deg;
 239         return END;
 240       }
 241     }
 242     if (d == rel) {  // found the same depth
 243       deg++;
 244       if (deg == ith) {
 245         *ans = i;
 246         return FOUND;
 247       }
 248     }
 249     i++;
 250   }
 251   *ans = deg;
 252   return CONTINUE;
 253 }
 254
 255 int degree_SB(bp *b, int i, int t, int rel, int *ans, int ith)
 256 {
 257   int j,r,n,il;
 258   pb *p,x,w,w2;
 259   int d, deg, v;
 260   int degtmp,itmp;
 261   int ith2,d2;
 262
 263   n = t;
 264   il = min((SBid(i) + 1) << logSB,n);
 265   d = deg = 0;
 266
 267   p = &b->B[i>>logD];
 268   while (i < il) {
 269     x = *p++;
 270     j = i & (D-1);
 271     x <<= j;
 272     j = min(D-j,il-i);
 273     while (j>0) {
 274       w = (x >> (D-ETW)) & ((1<<ETW)-1);
 275       w2 = 0;
 276       if (j < ETW || il-i < ETW) {
 277         r = max(ETW-j,ETW-(il-i));
 278         w2 = (1<<r)-1;
 279       }
 280       v = minmaxtbl_v[0][w | w2];
 281       if (d + v < rel) {
 282         if (ith > 0) {
 283 #if 0
 284           for (r = 0; r < ETW; r++) {
 285             if (w & 0x80) {
 286               d++;
 287             } else {
 288               d--;
 289               if (d < rel) break;
 290               if (d == rel) {
 291                 ith--;
 292                 if (ith == 0) {
 293                   *ans = i + r;
 294                   return FOUND;
 295                 }
 296               }
 297             }
 298             w <<= 1;
 299           }
 300           return NOTFOUND;
 301 #else
 302           r = childtbl2[rel-d+ETW][ith-1][w];
 303           if (r >= 0) {
 304             *ans = i + r;
 305             return FOUND;
 306           }
 307           return NOTFOUND;
 308 #endif
 309         }
 310         r = ETW-1-minmaxtbl_i[0][w | w2];
 311         w2 = (1<<r)-1;
 312         deg += degtbl2[((rel-d+ETW)<<ETW) + (w & (~w2))];
 313         *ans = deg;
 314         return END;
 315       }
 316       if (d + v == rel) {
 317         r = degtbl[w | w2];
 318         deg += r;
 319         if (ith > 0) {
 320           if (ith <= r) {
 321             *ans = i + childtbl[((ith-1)<<ETW) + (w | w2)];
 322             return FOUND;
 323           }
 324           ith -= r;
 325         }
 326       }
 327
 328       r = min(j,ETW);
 329       d += 2*popcount(w)-r;
 330       x <<= r;
 331       i += r;
 332       j -= r;
 333     }
 334   }
 335
 336   *ans = deg;
 337   return CONTINUE;
 338 }
 339
 340 int degree_MB(bp *b, int i, int t, int td, int *ans,int ith)
 341 {
 342   int il,d;
 343   int m,M,n,r;
 344   pb *B;
 345   int deg,degtmp;
 346
 347   d = 0;
 348   B = b->B;
 349   n = t;
 350
 351   il = min((MBid(i) + 1) << logMB,n);
 352   deg = 0;
 353   for (  ;  i+SB-1 < il;  i+=SB) {
 354 #if (SB % RRR != 0)
 355     d = bp_depth(b,i-1);
 356 #else
 357     d = fast_depth(b,i-1);
 358 #endif
 359     m = d + b->sm[SBid(i)] - SB;
 360     if (m < td) {
 361       r = degree_SB(b,i,n,td-d,&degtmp,ith);
 362       if (ith > 0) {
 363         if (r == NOTFOUND) return NOTFOUND;
 364         *ans = degtmp;
 365         return FOUND;
 366       } else {
 367         *ans = deg + degtmp;
 368         return END;
 369       }
 370     }
 371     if (m == td) {
 372       if (ith > 0) {
 373         if (ith <= b->sd[SBid(i)]) break;
 374         ith -= b->sd[SBid(i)];
 375       }
 376       deg += b->sd[SBid(i)];
 377     }
 378   }
 379   if (i < il) {
 380 #if (SB % RRR != 0)
 381     d = bp_depth(b,i-1);
 382 #else
 383     d = fast_depth(b,i-1);
 384 #endif
 385     r = degree_SB(b,i,n,td-d,&degtmp,ith);
 386     if (ith > 0) {
 387       if (r == NOTFOUND) return NOTFOUND;
 388       if (r == FOUND) {
 389         *ans = degtmp;
 390         return FOUND;
 391       }
 392     } else {
 393       deg += degtmp;
 394     }
 395   }
 396   *ans = deg;
 397   if (i >= n) return END;
 398   return CONTINUE;
 399 }
 400
 401 static void partition_range(int s,int t)
 402 {
 403   int i,j,h;
 404
 405   printf("partition [%d,%d] => ",s,t);
 406   h = 1;
 407   while (s <= t) {
 408     if (s & h) {
 409       if (s+h-1 <= t) {
 410         printf("[%d,%d] ",s,s+h-1);
 411         s += h;
 412       }
 413     } else {
 414       if (s+h > t) break;
 415     }
 416     h <<= 1;
 417   }
 418   while (h > 0) {
 419     if (s+h-1 <= t) {
 420       printf("[%d,%d] ",s,s+h-1);
 421       s += h;
 422     }
 423     h >>= 1;
 424   }
 425   printf("\n");
 426 }
 427
 428
 429
 430
 431 ///////////////////////////////////////////
 432 //  bp_fast_degree(bp *b,int s, int t, int ith)
 433 //    returns the number of children of s, to the left of t
 434 //    returns the position of (ith)-th child of s if (ith > 0)
 435 ///////////////////////////////////////////
 436 int bp_fast_degree(bp *b,int s, int t, int ith)
 437 {
 438   int i,j,n;
 439   int d,td;
 440   int m,M;
 441   int m_ofs;
 442   pb *B;
 443   int deg,degtmp;
 444   int sm,tm,ss,h;
 445
 446   n = t;
 447   B = b->B;
 448
 449   deg = 0;
 450
 451   i = s+1;
 452   if (i != SBfirst(i)) {
 453     d = degree_SB(b,i,n,0,&degtmp,ith);
 454     if (ith > 0) {
 455       if (d == NOTFOUND) return -1;
 456       if (d == FOUND) return degtmp;
 457       ith -= degtmp;
 458     }
 459     if (d == END) return degtmp;
 460     deg += degtmp;
 461   }
 462
 463   td = bp_depth(b,s);
 464
 465   i = SBid(i+SB-1) << logSB;
 466
 467   if (i != MBfirst(i)) {
 468     d = degree_MB(b,i,n,td,&degtmp,ith);
 469     if (ith > 0) {
 470       if (d == NOTFOUND) return -1;
 471       if (d == FOUND) return degtmp;
 472       ith -= degtmp;
 473       deg += degtmp;
 474     } else {
 475       deg += degtmp;
 476       if (d == END) {
 477         return deg;
 478       }
 479     }
 480   }
 481
 482 #if 0
 483   // sequential search
 484
 485   sm = MBid(i+MB-1);
 486   tm = MBid((n-1)+1)-1; // the rightmost MB fully contained in [0,n-1]
 487
 488   m_ofs = b->m_ofs;
 489   sm += m_ofs;  tm += m_ofs;
 490   for (i=sm; i<=tm; i++) {
 491     if (b->mm[i] < td) {
 492       break;
 493     }
 494     if (b->mm[i] == td) {
 495       if (ith > 0) {
 496         if (ith <= b->md[i]) break;
 497         ith -= b->md[i];
 498       }
 499       deg += b->md[i];
 500     }
 501   }
 502   ss = i - m_ofs;
 503 #else
 504   sm = MBid(i+MB-1);
 505   tm = MBid((n-1)+1)-1; // the rightmost MB fully contained in [0,n-1]
 506
 507   m_ofs = b->m_ofs;
 508   sm += m_ofs;  tm += m_ofs;
 509   ss = sm;
 510
 511   //partition_range(sm,tm);
 512
 513   //printf("partition [%d,%d] => ",sm,tm);
 514   h = 1;
 515   while (sm <= tm) {
 516     if (sm & h) {
 517       if (sm+h-1 <= tm) {
 518         //printf("[%d,%d] ",sm,sm+h-1);
 519         j = sm / h;
 520         if (b->mm[j] < td) {
 521           h >>= 1;
 522           break;
 523         }
 524         if (b->mm[j] == td) {
 525           if (ith > 0) {
 526             if (ith <= b->md[j]) {
 527               h >>= 1;
 528               break;
 529             }
 530             ith -= b->md[j];
 531           }
 532           deg += b->md[j];
 533         }
 534         sm += h;
 535       }
 536     } else {
 537       if (sm+h > tm) break;
 538     }
 539     h <<= 1;
 540   }
 541   while (h > 0) {
 542     if (sm+h-1 <= tm) {
 543       //printf("[%d,%d] ",sm,sm+h-1);
 544       j = sm / h;
 545       if (ith > 0) {
 546         if (b->mm[j] >= td) {
 547           if (b->mm[j] == td) {
 548             if (ith > b->md[j]) {
 549               ith -= b->md[j];
 550               sm += h;
 551             } else {
 552               deg += b->md[j];
 553             }
 554           } else {
 555             sm += h;
 556           }
 557         }
 558       } else {
 559         if (b->mm[j] >= td) {
 560           if (b->mm[j] == td) {
 561             deg += b->md[j];
 562           }
 563           sm += h;
 564         }
 565       }
 566     }
 567     h >>= 1;
 568   }
 569   //printf("\n");
 570   ss = sm;
 571
 572   ss -= m_ofs;
 573
 574 #endif
 575
 576   ss <<= logMB;
 577
 578   d = degree_MB(b,ss,n,td,&degtmp,ith);
 579   if (ith > 0) {
 580     if (d == NOTFOUND) return -1;
 581     if (d == FOUND) return degtmp;
 582   }
 583   deg += degtmp;
 584   if (d == END) return deg;
 585   return deg;
 586
 587   // unexpected (bug)
 588   printf("degree: ???\n");
 589   return -99;
 590
 591 }
 592
 593 int search_SB_l(bp *b, int i, int rel)
 594 {
 595   int j,r,il;
 596   pb *p,x,w;
 597
 598   il = SBfirst(i);
 599
 600   p = &b->B[i>>logD];
 601   while (i>=il) {
 602     x = *p--;
 603     j = (i & (D-1))+1;
 604     x >>= D-j;
 605     while (j>0) {
 606       w = x & ((1<<ETW)-1);
 607       if (rel >= -ETW && rel <= ETW) {
 608         r = bwdtbl[((rel+ETW)<<ETW)+w];
 609         if (r<ETW && r<j) {
 610           if (i-r < 0) return NOTFOUND;
 611           return i-r-1;
 612         }
 613       }
 614       r = min(j,ETW);
 615       rel += 2*popcount(w)-r;
 616       x >>= r;
 617       i -= r;
 618       j -= r;
 619     }
 620
 621   }
 622   if (i < 0) return NOTFOUND;
 623   return CONTINUE;
 624 }
 625
 626 int search_MB_l(bp *b, int i, int td)
 627 {
 628   int il,d;
 629   int m,M;
 630   pb *B;
 631
 632 #if 0
 633   if (i % SB != SB-1) {
 634     printf("search_MB_l:error!!! i=%d SB=%d\n",i,i%SB);
 635   }
 636 #endif
 637   B = b->B;
 638
 639   il = MBfirst(i);
 640   for (  ;  i >= il;  i-=SB) {
 641 #if (SB % RRR != 0)
 642     d = bp_depth(b,i-SB);
 643 #else
 644     d = fast_depth(b,i-SB);
 645 #endif
 646     m = d + b->sm[SBid(i)] - SB;
 647     M = d + b->sM[SBid(i)] - 1;
 648     if (m <= td && td <= M) {
 649 #if (SB % RRR != 0)
 650       d = bp_depth(b,i);
 651 #else
 652       d = fast_depth(b,i);
 653 #endif
 654       if (d == td) return i;
 655       return search_SB_l(b,i,td-d);
 656     }
 657   }
 658   return CONTINUE;
 659 }
 660
 661 ///////////////////////////////////////////
 662 //  bwd_excess(bp *b,int s, int rel)
 663 //    find the rightmost value depth(s)+rel to the left of s (exclusive)
 664 ///////////////////////////////////////////
 665 int bp_bwd_excess(bp *b,int s, int rel)
 666 {
 667   int i,n;
 668   int d,td;
 669   int m,M;
 670   int m_ofs;
 671   pb *B;
 672   n = b->n;  B = b->B;
 673
 674   i = s;
 675   d = search_SB_l(b,i,rel);
 676   if (d >= NOTFOUND) return d;
 677
 678   i = SBfirst(i) -1;
 679
 680   td = bp_depth(b,s) + rel;
 681
 682   d = search_MB_l(b,i,td);
 683   if (d >= NOTFOUND) return d;
 684
 685   m_ofs = b->m_ofs;
 686   i = (s>>logMB) + m_ofs;
 687   while (i > 0) {
 688     if ((i&1) == 1) {
 689       i--;
 690       m = b->mm[i];
 691       M = b->mM[i];
 692       if (m <= td && td <= M) break;
 693     }
 694     i >>= 1;
 695   }
 696   if (i == 0) {
 697     if (td == 0) return -1;
 698     else return NOTFOUND;
 699   }
 700
 701   while (i < m_ofs) {
 702     i = i*2 + 1;
 703     m = b->mm[i];
 704     M = b->mM[i];
 705     if (!(m <= td && td <= M)) i--;
 706   }
 707   i -= m_ofs;
 708   i = ((i+1)<<logMB)-1;
 709
 710   d = search_MB_l(b,i,td);
 711   if (d >= NOTFOUND) return d;
 712
 713   // unexpected (bug)
 714   printf("bwd_excess: ???\n");
 715   return -99;
 716
 717 }
 718
 719 int rmq_SB(bp *b, int s, int t, int opt, int *dm)
 720 {
 721   int i,d;
 722   int is,ds;
 723   pb *p,x,w,w2;
 724   int j,v,r;
 725   int lr;
 726   int op;
 727
 728   lr = 0;  if (opt & OPT_RIGHT) lr = 1;
 729   op = opt & (OPT_RIGHT | OPT_MAX);
 730
 731   is = s;  ds = d = 0;
 732   i = s+1;
 733
 734 #if SB >= ETW
 735   p = &b->B[i>>logD];
 736   while (i <= t) {
 737     x = *p++;
 738     j = i & (D-1);
 739     x <<= j;
 740     j = min(D-j,t-i+1);
 741     while (j>0) {
 742       w = (x >> (D-ETW)) & ((1<<ETW)-1);
 743       w2 = 0;
 744       if (j < ETW || t-i < ETW-1) {
 745         r = max(ETW-j,ETW-1-(t-i));
 746         w2 = (1<<r)-1;
 747       }
 748
 749       if (op & OPT_MAX) {
 750         v = minmaxtbl_v[op][w & (~w2)];
 751         if (d + v + lr > ds) {
 752           ds = d + v;
 753           is = i + minmaxtbl_i[op][w & (~w2)];
 754         }
 755       } else {
 756         v = minmaxtbl_v[op][w | w2];
 757         if (d + v < ds +lr) {
 758           ds = d + v;
 759           is = i + minmaxtbl_i[op][w | w2];
 760         }
 761       }
 762
 763       r = min(j,ETW);
 764       d += 2*popcount(w)-r;
 765       x <<= r;
 766       i += r;
 767       j -= r;
 768     }
 769   }
 770 #else
 771   while (i <= t) {
 772     if (bp_getbit(b->B,i)==OP) {
 773       d++;
 774       if (op & OPT_MAX) {
 775         if (d + lr > ds) {
 776           ds = d;  is = i;
 777         }
 778       }
 779     } else {
 780       d--;
 781       if (!(op & OPT_MAX)) {
 782         if (d < ds + lr) {
 783           ds = d;  is = i;
 784         }
 785       }
 786     }
 787     i++;
 788   }
 789 #endif
 790   *dm = ds;
 791   return is;
 792 }
 793
 794 int rmq_MB(bp *b, int s, int t, int opt, int *dm)
 795 {
 796   int i,d,m;
 797   int mi,md;
 798   int lr;
 799
 800   lr = 0;  if (opt & OPT_RIGHT) lr = 1;
 801
 802   md = *dm;  mi = -1;
 803   for (i = s;  i <= t;  i++) {
 804 #if (SB % RRR != 0)
 805     d = bp_depth(b,(i<<logSB)-1);
 806 #else
 807     d = fast_depth(b,(i<<logSB)-1);
 808 #endif
 809     if (opt & OPT_MAX) {
 810       m = d + b->sM[i] - 1;
 811       if (m + lr > md) {
 812         md = m;  mi = i;
 813       }
 814     } else {
 815       m = d + b->sm[i] - SB;
 816       if (m < md + lr) {
 817         md = m;  mi = i;
 818       }
 819     }
 820   }
 821   *dm = md;
 822   return mi;
 823 }
 824
 825
 826
 827
 828 ///////////////////////////////////////////
 829 //  bp_rmq(bp *b, int s, int t, int opt)
 830 //    returns the index of leftmost/rightmost minimum/maximum value
 831 //                 in the range [s,t] (inclusive)
 832 //    returns the leftmost minimum if opt == 0
 833 //            the rightmost one if (opt & OPT_RIGHT)
 834 //            the maximum if (opt & OPT_MAX)
 835 ///////////////////////////////////////////
 836 int bp_rmq(bp *b, int s, int t, int opt)
 837 {
 838   int ss, ts;  // SB index of s and t
 839   int sm, tm;  // MB index of s and t
 840   int ds;   // current min value
 841   int is;   // current min index
 842   int ys;   // type of current min index
 843                // 0: is is the index of min
 844                // 1: is is the SB index
 845                // 2: is is the MB index
 846   int m_ofs;
 847   int i,j,il,d,n;
 848   int dm;
 849   int lr;
 850
 851   lr = 0;  if (opt & OPT_RIGHT) lr = 1;
 852
 853   n = b->n;
 854   if (s < 0 || t >= n || s > t) {
 855     printf("rmq: error s=%d t=%d n=%d\n",s,t,n);
 856     return -1;
 857   }
 858   if (s == t) return s;
 859
 860
 861   ////////////////////////////////////////////////////////////
 862   // search the SB of s
 863   ////////////////////////////////////////////////////////////
 864
 865   il = min(SBlast(s),t);
 866   is = rmq_SB(b,s,il,opt,&dm);
 867   if (il == t) {  // scan reached the end of the range
 868     return is;
 869   }
 870   ds = bp_depth(b,s) + dm;  ys = 0;
 871
 872   ////////////////////////////////////////////////////////////
 873   // search the MB of s
 874   ////////////////////////////////////////////////////////////
 875
 876   ss = SBid(s) + 1;
 877   il = min(SBid(MBlast(s)),SBid(t)-1);
 878   dm = ds;
 879   j = rmq_MB(b,ss,il,opt,&dm);
 880   //if (dm < ds + lr) {
 881   if (j >= 0) {
 882     ds = dm;  is = j;  ys = 1;
 883   }
 884
 885   ////////////////////////////////////////////////////////////
 886   // search the tree
 887   ////////////////////////////////////////////////////////////
 888
 889   sm = MBid(s) + 1;
 890   tm = MBid(t) - 1;
 891
 892 #if 0
 893   // sequential search
 894   m_ofs = b->m_ofs;
 895   sm += m_ofs;  tm += m_ofs;
 896   for (i=sm; i<=tm; i++) {
 897     if (opt & OPT_MAX) {
 898       if (b->mM[i] + lr > ds) {
 899         ds = b->mM[i];  is = i - m_ofs;  ys = 2;
 900       }
 901     } else {
 902       if (b->mm[i] < ds + lr) {
 903         ds = b->mm[i];  is = i - m_ofs;  ys = 2;
 904       }
 905     }
 906   }
 907
 908 #else
 909   if (sm <= tm) {
 910     int h;
 911     h = blog(sm ^ tm);
 912
 913     m_ofs = b->m_ofs;
 914     sm += m_ofs;  tm += m_ofs;
 915
 916     if (opt & OPT_MAX) {
 917       if (b->mM[sm] + lr > ds) {
 918         ds = b->mM[sm];  is = sm;  ys = 2;
 919       }
 920       for (i=0; i<=h-2; i++) {
 921         j = sm>>i;
 922         if ((j&1) == 0) {
 923           if (b->mM[j+1] + lr > ds) {
 924             ds = b->mM[j+1];  is = j+1;  ys = 2;
 925           }
 926         }
 927       }
 928       for (i=h-2; i>=0; i--) {
 929         j = tm>>i;
 930         if ((j&1)==1) {
 931           if (b->mM[j-1] + lr > ds) {
 932             ds = b->mM[j-1];  is = j-1;  ys = 2;
 933           }
 934         }
 935       }
 936       if (b->mM[tm] + lr > ds) {
 937         ds = b->mM[tm];  is = tm;  ys = 2;
 938       }
 939       if (ys == 2) {
 940         while (is < m_ofs) {
 941           is <<= 1;
 942           if (b->mM[is+1] + lr > b->mM[is]) is++;
 943         }
 944         is -= m_ofs;
 945       }
 946     } else { // MIN
 947       if (b->mm[sm] < ds + lr) {
 948         ds = b->mm[sm];  is = sm;  ys = 2;
 949       }
 950       for (i=0; i<=h-2; i++) {
 951         j = sm>>i;
 952         if ((j&1) == 0) {
 953           if (b->mm[j+1] < ds + lr) {
 954             ds = b->mm[j+1];  is = j+1;  ys = 2;
 955           }
 956         }
 957       }
 958       for (i=h-2; i>=0; i--) {
 959         j = tm>>i;
 960         if ((j&1)==1) {
 961           if (b->mm[j-1] < ds + lr) {
 962             ds = b->mm[j-1];  is = j-1;  ys = 2;
 963           }
 964         }
 965       }
 966       if (b->mm[tm] < ds + lr) {
 967         ds = b->mm[tm];  is = tm;  ys = 2;
 968       }
 969       if (ys == 2) {
 970         while (is < m_ofs) {
 971           is <<= 1;
 972           if (b->mm[is+1] < b->mm[is] + lr) is++;
 973         }
 974         is -= m_ofs;
 975       }
 976     }
 977   }
 978
 979 #endif
 980
 981   ////////////////////////////////////////////////////////////
 982   // search the MB of t
 983   ////////////////////////////////////////////////////////////
 984
 985
 986   ts = max(SBid(MBfirst(t)),SBid(s)+1);
 987   il = SBid(SBfirst(t)-1);
 988   dm = ds;
 989   j = rmq_MB(b,ts,il,opt,&dm);
 990   //if (dm < ds + lr) {
 991   if (j >= 0) {
 992     ds = dm;  is = j;  ys = 1;
 993   }
 994
 995   ////////////////////////////////////////////////////////////
 996   // search the SB of t
 997   ////////////////////////////////////////////////////////////
 998
 999   i = SBfirst(t);
1000   j = rmq_SB(b,i,t,opt,&dm);
1001   d = bp_depth(b,i) + dm;
1002   if (opt & OPT_MAX) {
1003     if (d + lr > ds) {
1004       ds = d;  is = j;  ys = 0;
1005     }
1006   } else {
1007     if (d < ds + lr) {
1008       ds = d;  is = j;  ys = 0;
1009     }
1010   }
1011
1012   ////////////////////////////////////////////////////////////
1013   // search the rest
1014   ////////////////////////////////////////////////////////////
1015
1016   if (ys == 2) {
1017     ss = SBid(is << logMB);
1018     il = SBid(MBlast(is << logMB));
1019     if (opt & OPT_MAX) {
1020       dm = -n-1;
1021     } else {
1022       dm = n+1;
1023     }
1024     j = rmq_MB(b,ss,il,opt,&dm);
1025     ds = dm;  is = j;  ys = 1;
1026   }
1027
1028   if (ys == 1) {
1029     ss = is << logSB;
1030     il = SBlast(is << logSB);
1031     is = rmq_SB(b,ss,il,opt,&dm);
1032     //ds = bp_depth(b,ss) + dm;  ys = 0;
1033   }
1034
1035   return is;
1036 }
1037