bpcore.c

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