bpcore.c

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