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