bpcore.c

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