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