bpcore.c

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