diff --git a/include/divsufsort_private.h b/include/divsufsort_private.h
index 7e261c1..d65485d 100644
--- a/include/divsufsort_private.h
+++ b/include/divsufsort_private.h
@@ -197,7 +197,7 @@ sssort(const sauchar_t *Td, const saidx_t *PA,
        saidx_t depth, saidx_t n, saint_t lastsuffix);
 /* trsort.c */
 void
-trsort(saidx_t *ISA, saidx_t *SA, saidx_t n, saidx_t depth);
+trsort(saidx_t *ISA, saidx_t *SA, saidx_t n, saidx_t depth, saidx_t *buf, saidx_t bufsize);
 
 
 #ifdef __cplusplus
diff --git a/lib/divsufsort.c b/lib/divsufsort.c
index 9f64b4f..de729ce 100644
--- a/lib/divsufsort.c
+++ b/lib/divsufsort.c
@@ -156,9 +156,10 @@ sort_typeBstar(const sauchar_t *T, saidx_t *SA,
     }
 
     /* Construct the inverse suffix array of type B* suffixes using trsort. */
-    trsort(ISAb, SA, m, 1);
+    buf = ISAb + m, bufsize = n - (2 * m);
+    trsort(ISAb, SA, m, 1, buf, bufsize);
 
-    /* Set the sorted order of tyoe B* suffixes. */
+    /* Set the sorted order of type B* suffixes. */
     for(i = n - 1, j = m, c0 = T[n - 1]; 0 <= i;) {
       for(--i, c1 = c0; (0 <= i) && ((c0 = T[i]) >= c1); --i, c1 = c0) { }
       if(0 <= i) {
diff --git a/lib/trsort.c b/lib/trsort.c
index 6fe3e67..997f8e0 100644
--- a/lib/trsort.c
+++ b/lib/trsort.c
@@ -328,13 +328,13 @@ tr_introsort(saidx_t *ISA, const saidx_t *ISAd,
              saidx_t *SA, saidx_t *first, saidx_t *last,
              trbudget_t *budget) {
 #define STACK_SIZE TR_STACKSIZE
-  struct { const saidx_t *a; saidx_t *b, *c; saint_t d, e; }stack[STACK_SIZE];
+  struct { const saidx_t *a; saidx_t *b, *c; saint_t d; }stack[STACK_SIZE];
   saidx_t *a, *b, *c;
   saidx_t t;
   saidx_t v, x = 0;
   saidx_t incr = ISAd - ISA;
   saint_t limit, next;
-  saint_t ssize, trlink = -1;
+  saint_t ssize, partial_copy_flag = 0;
 
   for(ssize = 0, limit = tr_ilg(last - first);;) {
 
@@ -353,39 +353,39 @@ tr_introsort(saidx_t *ISA, const saidx_t *ISAd,
 
         /* push */
         if(1 < (b - a)) {
-          STACK_PUSH5(NULL, a, b, 0, 0);
-          STACK_PUSH5(ISAd - incr, first, last, -2, trlink);
-          trlink = ssize - 2;
+          STACK_PUSH(NULL, a, b, partial_copy_flag);
+          STACK_PUSH(ISAd - incr, first, last, -2);
+          partial_copy_flag = 0;
         }
         if((a - first) <= (last - b)) {
           if(1 < (a - first)) {
-            STACK_PUSH5(ISAd, b, last, tr_ilg(last - b), trlink);
+            STACK_PUSH(ISAd, b, last, tr_ilg(last - b));
             last = a, limit = tr_ilg(a - first);
           } else if(1 < (last - b)) {
             first = b, limit = tr_ilg(last - b);
           } else {
-            STACK_POP5(ISAd, first, last, limit, trlink);
+            STACK_POP(ISAd, first, last, limit);
           }
         } else {
           if(1 < (last - b)) {
-            STACK_PUSH5(ISAd, first, a, tr_ilg(a - first), trlink);
+            STACK_PUSH(ISAd, first, a, tr_ilg(a - first));
             first = b, limit = tr_ilg(last - b);
           } else if(1 < (a - first)) {
             last = a, limit = tr_ilg(a - first);
           } else {
-            STACK_POP5(ISAd, first, last, limit, trlink);
+            STACK_POP(ISAd, first, last, limit);
           }
         }
       } else if(limit == -2) {
         /* tandem repeat copy */
         a = stack[--ssize].b, b = stack[ssize].c;
-        if(stack[ssize].d == 0) {
+        if(partial_copy_flag == 0) {
+          partial_copy_flag = stack[ssize].d;
           tr_copy(ISA, SA, first, a, b, last, ISAd - ISA);
         } else {
-          if(0 <= trlink) { stack[trlink].d = -1; }
           tr_partialcopy(ISA, SA, first, a, b, last, ISAd - ISA);
         }
-        STACK_POP5(ISAd, first, last, limit, trlink);
+        STACK_POP(ISAd, first, last, limit);
       } else {
         /* sorted partition */
         if(0 <= *first) {
@@ -401,26 +401,26 @@ tr_introsort(saidx_t *ISA, const saidx_t *ISAd,
           /* push */
           if(trbudget_check(budget, a - first)) {
             if((a - first) <= (last - a)) {
-              STACK_PUSH5(ISAd, a, last, -3, trlink);
+              STACK_PUSH(ISAd, a, last, -3);
               ISAd += incr, last = a, limit = next;
             } else {
               if(1 < (last - a)) {
-                STACK_PUSH5(ISAd + incr, first, a, next, trlink);
+                STACK_PUSH(ISAd + incr, first, a, next);
                 first = a, limit = -3;
               } else {
                 ISAd += incr, last = a, limit = next;
               }
             }
           } else {
-            if(0 <= trlink) { stack[trlink].d = -1; }
+            partial_copy_flag = 1;
             if(1 < (last - a)) {
               first = a, limit = -3;
             } else {
-              STACK_POP5(ISAd, first, last, limit, trlink);
+              STACK_POP(ISAd, first, last, limit);
             }
           }
         } else {
-          STACK_POP5(ISAd, first, last, limit, trlink);
+          STACK_POP(ISAd, first, last, limit);
         }
       }
       continue;
@@ -460,75 +460,75 @@ tr_introsort(saidx_t *ISA, const saidx_t *ISAd,
         if((a - first) <= (last - b)) {
           if((last - b) <= (b - a)) {
             if(1 < (a - first)) {
-              STACK_PUSH5(ISAd + incr, a, b, next, trlink);
-              STACK_PUSH5(ISAd, b, last, limit, trlink);
+              STACK_PUSH(ISAd + incr, a, b, next);
+              STACK_PUSH(ISAd, b, last, limit);
               last = a;
             } else if(1 < (last - b)) {
-              STACK_PUSH5(ISAd + incr, a, b, next, trlink);
+              STACK_PUSH(ISAd + incr, a, b, next);
               first = b;
             } else {
               ISAd += incr, first = a, last = b, limit = next;
             }
           } else if((a - first) <= (b - a)) {
             if(1 < (a - first)) {
-              STACK_PUSH5(ISAd, b, last, limit, trlink);
-              STACK_PUSH5(ISAd + incr, a, b, next, trlink);
+              STACK_PUSH(ISAd, b, last, limit);
+              STACK_PUSH(ISAd + incr, a, b, next);
               last = a;
             } else {
-              STACK_PUSH5(ISAd, b, last, limit, trlink);
+              STACK_PUSH(ISAd, b, last, limit);
               ISAd += incr, first = a, last = b, limit = next;
             }
           } else {
-            STACK_PUSH5(ISAd, b, last, limit, trlink);
-            STACK_PUSH5(ISAd, first, a, limit, trlink);
+            STACK_PUSH(ISAd, b, last, limit);
+            STACK_PUSH(ISAd, first, a, limit);
             ISAd += incr, first = a, last = b, limit = next;
           }
         } else {
           if((a - first) <= (b - a)) {
             if(1 < (last - b)) {
-              STACK_PUSH5(ISAd + incr, a, b, next, trlink);
-              STACK_PUSH5(ISAd, first, a, limit, trlink);
+              STACK_PUSH(ISAd + incr, a, b, next);
+              STACK_PUSH(ISAd, first, a, limit);
               first = b;
             } else if(1 < (a - first)) {
-              STACK_PUSH5(ISAd + incr, a, b, next, trlink);
+              STACK_PUSH(ISAd + incr, a, b, next);
               last = a;
             } else {
               ISAd += incr, first = a, last = b, limit = next;
             }
           } else if((last - b) <= (b - a)) {
             if(1 < (last - b)) {
-              STACK_PUSH5(ISAd, first, a, limit, trlink);
-              STACK_PUSH5(ISAd + incr, a, b, next, trlink);
+              STACK_PUSH(ISAd, first, a, limit);
+              STACK_PUSH(ISAd + incr, a, b, next);
               first = b;
             } else {
-              STACK_PUSH5(ISAd, first, a, limit, trlink);
+              STACK_PUSH(ISAd, first, a, limit);
               ISAd += incr, first = a, last = b, limit = next;
             }
           } else {
-            STACK_PUSH5(ISAd, first, a, limit, trlink);
-            STACK_PUSH5(ISAd, b, last, limit, trlink);
+            STACK_PUSH(ISAd, first, a, limit);
+            STACK_PUSH(ISAd, b, last, limit);
             ISAd += incr, first = a, last = b, limit = next;
           }
         }
       } else {
-        if((1 < (b - a)) && (0 <= trlink)) { stack[trlink].d = -1; }
+        if(1 < (b - a)) { partial_copy_flag = 1; }
         if((a - first) <= (last - b)) {
           if(1 < (a - first)) {
-            STACK_PUSH5(ISAd, b, last, limit, trlink);
+            STACK_PUSH(ISAd, b, last, limit);
             last = a;
           } else if(1 < (last - b)) {
             first = b;
           } else {
-            STACK_POP5(ISAd, first, last, limit, trlink);
+            STACK_POP(ISAd, first, last, limit);
           }
         } else {
           if(1 < (last - b)) {
-            STACK_PUSH5(ISAd, first, a, limit, trlink);
+            STACK_PUSH(ISAd, first, a, limit);
             first = b;
           } else if(1 < (a - first)) {
             last = a;
           } else {
-            STACK_POP5(ISAd, first, last, limit, trlink);
+            STACK_POP(ISAd, first, last, limit);
           }
         }
       }
@@ -536,14 +536,507 @@ tr_introsort(saidx_t *ISA, const saidx_t *ISAd,
       if(trbudget_check(budget, last - first)) {
         limit = tr_ilg(last - first), ISAd += incr;
       } else {
-        if(0 <= trlink) { stack[trlink].d = -1; }
-        STACK_POP5(ISAd, first, last, limit, trlink);
+        partial_copy_flag = 1;
+        STACK_POP(ISAd, first, last, limit);
       }
     }
   }
 #undef STACK_SIZE
 }
 
+static
+void
+tr_first_stage_introsort(saidx_t *ISA, const saidx_t *ISAd,
+                         saidx_t *SA, saidx_t *LA,
+                         saidx_t *first, saidx_t *last,
+                         trbudget_t *budget) {
+#define STACK_SIZE TR_STACKSIZE
+  struct { const saidx_t *a; saidx_t *b, *c; saint_t d; }stack[STACK_SIZE];
+  saidx_t *a, *b;
+  saidx_t t;
+  saidx_t v, x = 0;
+  saidx_t incr = ISAd - ISA;
+  saint_t limit, next;
+  saint_t ssize;
+
+  for(ssize = 0, limit = tr_ilg(last - first);;) {
+
+    if(limit < 0) {
+
+        /* sorted partition */
+        while(*first >= 0) {
+          if(++first >= last) break;
+        }
+        if(first < last) {
+          a = first; do { *a = ~*a; } while(*++a < 0);
+          b = a++;
+
+          /* push */
+          if((ISA[*b] != ISAd[*b]) && trbudget_check(budget, a - first)) {
+            next = tr_ilg(a - first);
+            if((a - first) <= (last - a)) {
+              STACK_PUSH(ISAd, a, last, -3);
+              ISAd += incr, last = a, limit = next;
+            } else {
+              if(1 < (last - a)) {
+                STACK_PUSH(ISAd, first, a, next);
+                first = a, limit = -3;
+              } else {
+                ISAd += incr, last = a, limit = next;
+              }
+            }
+          } else {
+            if((ISA[*b] == ISAd[*b]) && trbudget_check(budget, a - first)) {
+              LA[a - SA - 1] = ISAd - ISA;
+            } else {
+              LA[a - SA - 1] = 0;
+            }
+            LA[a - SA - 2] = first - SA;
+            a[-1] = ~a[-1];
+            if(1 < (last - a)) {
+              first = a, limit = -3;
+            } else {
+              STACK_POP(ISAd, first, last, limit);
+            }
+          }
+        } else {
+          STACK_POP(ISAd, first, last, limit);
+        }
+      continue;
+    }
+
+    if((last - first) <= TR_INSERTIONSORT_THRESHOLD) {
+      tr_insertionsort(ISAd, first, last);
+      limit = -3;
+      continue;
+    }
+
+    if(limit-- == 0) {
+      tr_heapsort(ISAd, first, last - first);
+      for(a = last - 1; first < a; a = b) {
+        for(x = ISAd[*a], b = a - 1; (first <= b) && (ISAd[*b] == x); --b) { *b = ~*b; }
+      }
+      limit = -3;
+      continue;
+    }
+
+    /* choose pivot */
+    a = tr_pivot(ISAd, first, last);
+    SWAP(*first, *a);
+    v = ISAd[*first];
+
+    /* partition */
+    tr_partition(ISAd, first, first + 1, last, &a, &b, v);
+
+      /* push */
+      if((1 < (b - a)) && (ISA[*a] != v) && (trbudget_check(budget, b - a))) {
+        next = tr_ilg(b - a);
+        if((a - first) <= (last - b)) {
+          if((last - b) <= (b - a)) {
+            if(1 < (a - first)) {
+              STACK_PUSH(ISAd + incr, a, b, next);
+              STACK_PUSH(ISAd, b, last, limit);
+              last = a;
+            } else if(1 < (last - b)) {
+              STACK_PUSH(ISAd + incr, a, b, next);
+              first = b;
+            } else {
+              ISAd += incr, first = a, last = b, limit = next;
+            }
+          } else if((a - first) <= (b - a)) {
+            if(1 < (a - first)) {
+              STACK_PUSH(ISAd, b, last, limit);
+              STACK_PUSH(ISAd + incr, a, b, next);
+              last = a;
+            } else {
+              STACK_PUSH(ISAd, b, last, limit);
+              ISAd += incr, first = a, last = b, limit = next;
+            }
+          } else {
+            STACK_PUSH(ISAd, b, last, limit);
+            STACK_PUSH(ISAd, first, a, limit);
+            ISAd += incr, first = a, last = b, limit = next;
+          }
+        } else {
+          if((a - first) <= (b - a)) {
+            if(1 < (last - b)) {
+              STACK_PUSH(ISAd + incr, a, b, next);
+              STACK_PUSH(ISAd, first, a, limit);
+              first = b;
+            } else if(1 < (a - first)) {
+              STACK_PUSH(ISAd + incr, a, b, next);
+              last = a;
+            } else {
+              ISAd += incr, first = a, last = b, limit = next;
+            }
+          } else if((last - b) <= (b - a)) {
+            if(1 < (last - b)) {
+              STACK_PUSH(ISAd, first, a, limit);
+              STACK_PUSH(ISAd + incr, a, b, next);
+              first = b;
+            } else {
+              STACK_PUSH(ISAd, first, a, limit);
+              ISAd += incr, first = a, last = b, limit = next;
+            }
+          } else {
+            STACK_PUSH(ISAd, first, a, limit);
+            STACK_PUSH(ISAd, b, last, limit);
+            ISAd += incr, first = a, last = b, limit = next;
+          }
+        }
+      } else {
+        if(1 < (b - a)) {
+          if((ISA[*a] == v) && trbudget_check(budget, b - a)) {
+            LA[b - SA - 1] = ISAd - ISA;
+          } else {
+            LA[b - SA - 1] = 0;
+          }
+          LA[b - SA - 2] = a - SA;
+          b[-1] = ~b[-1];
+        }
+        if((a - first) <= (last - b)) {
+          if(1 < (a - first)) {
+            STACK_PUSH(ISAd, b, last, limit);
+            last = a;
+          } else if(1 < (last - b)) {
+            first = b;
+          } else {
+            STACK_POP(ISAd, first, last, limit);
+          }
+        } else {
+          if(1 < (last - b)) {
+            STACK_PUSH(ISAd, first, a, limit);
+            first = b;
+          } else if(1 < (a - first)) {
+            last = a;
+          } else {
+            STACK_POP(ISAd, first, last, limit);
+          }
+        }
+      }
+  }
+#undef STACK_SIZE
+}
+
+static
+void
+tr_second_stage_introsort(saidx_t *ISA, const saidx_t *ISAd,
+                          saidx_t *SA, saidx_t *LA,
+                          saidx_t *first, saidx_t *last,
+                          trbudget_t *budget) {
+#define STACK_SIZE TR_STACKSIZE
+  struct { const saidx_t *a; saidx_t *b, *c; saint_t d; }stack[STACK_SIZE];
+  saidx_t *a, *b, *c;
+  saidx_t t;
+  saidx_t v, x = 0;
+  saidx_t incr;
+  saint_t limit, next;
+  saint_t ssize, partial_copy_flag = 0;
+
+  /* update ranks */
+  a = last - 1;
+  b = last;
+  if(*a < 0) {
+    v = a - SA;
+    *a = ~*a;
+    a = SA + LA[v - 1];
+    if(LA[v] > 0) {
+      LA[v - 1] = b - SA;
+      b = a;
+    }
+  }
+  while(--a >= first) {
+    if(0 <= *a) {
+      do { ISA[*a] = a - SA; } while((--a >= first) && (0 <= *a));
+      if(a <= first) { break; }
+    }
+    v = a - SA;
+    c = SA + LA[v - 1];
+    ISA[*a = ~*a] = v;
+    do { ISA[*--a] = v; } while(a > c);
+    if(LA[v] > 0) {
+      LA[v - 1] = b - SA;
+      b = a;
+    }
+  }
+  if(b >= last) return;
+
+  for(ssize = 0, first = b, limit = -5;;) {
+
+    if(limit < 0) {
+      if(limit == -3) {
+        /* sorted partition */
+        while(*first >= 0) {
+          if(++first >= last) break;
+        }
+        if(first < last) {
+          a = first; do { *a = ~*a; } while(*++a < 0);
+          b = a++;
+
+          /* push */
+          if((ISA[*b] != ISAd[*b]) && (LA[ISAd[*b]] > 0) && trbudget_check(budget, a - first)) {
+            incr = LA[ISAd[*b]];
+            next = tr_ilg(a - first);
+            if((a - first) <= (last - a)) {
+              STACK_PUSH(ISAd, a, last, -3);
+              ISAd += incr, last = a, limit = next;
+            } else {
+              if(1 < (last - a)) {
+                STACK_PUSH(ISAd + incr, first, a, next);
+                first = a, limit = -3;
+              } else {
+                ISAd += incr, last = a, limit = next;
+              }
+            }
+          } else {
+            if((ISA[*b] == ISAd[*b]) && trbudget_check(budget, a - first)) {
+              LA[a - SA - 1] = ISAd - ISA;
+            } else {
+              partial_copy_flag = 1;
+              LA[a - SA - 1] = 0;
+            }
+            *(a - 1) = ~*(a - 1);
+            LA[a - SA - 2] = first - SA;
+            if(1 < (last - a)) {
+              first = a, limit = -3;
+            } else {
+              STACK_POP(ISAd, first, last, limit);
+            }
+          }
+        } else {
+          STACK_POP(ISAd, first, last, limit);
+        }
+      } else if(limit == -2) {
+        /* tandem repeat copy */
+        a = stack[--ssize].b, b = stack[ssize].c;
+        if(partial_copy_flag == 0) {
+          partial_copy_flag = stack[ssize].d;
+          tr_copy(ISA, SA, first, a, b, last, ISAd - ISA);
+        } else {
+          tr_partialcopy(ISA, SA, first, a, b, last, ISAd - ISA);
+        }
+        STACK_POP(ISAd, first, last, limit);
+      } else if(limit == -4) {
+        /* update ranks */
+        a = last - 1, b = last;
+        if(*a < 0) {
+          v = a - SA;
+          *a = ~*a;
+          a = SA + LA[v - 1];
+          if(LA[v] > 0) {
+            LA[v - 1] = b - SA;
+            b = a;
+          }
+        }
+        while(--a >= first) {
+          if(0 <= *a) {
+            do { ISA[*a] = a - SA; } while((--a >= first) && (0 <= *a));
+            if(a <= first) { break; }
+          }
+          v = a - SA;
+          c = SA + LA[v - 1];
+          ISA[*a = ~*a] = v;
+          do { ISA[*--a] = v; } while(a > c);
+          if(LA[v] > 0) {
+            LA[v - 1] = b - SA;
+            b = a;
+          }
+        }
+        if(b < last) {
+          first = b;
+          limit = -5;
+        } else {
+          STACK_POP(ISAd, first, last, limit);
+        }
+      } else /* if(limit == -5) */ {
+          v = ISA[*first];
+          b = SA + LA[v - 1];
+          a = SA + v + 1;
+          if(b < last) {
+            if((last - b) < (a - first)) {
+              STACK_PUSH(ISAd, first, a, -5);
+              first = b;
+              continue;
+            }
+            STACK_PUSH(ISAd, b, last, -5);
+          }
+          ISAd = ISA + LA[v];
+          last = a;
+
+          /* tandem repeat partition */
+          tr_partition(ISAd, first, first, last, &a, &b, last - SA - 1);
+  
+          /* update ranks */
+          if(a < last) {
+            for(c = first, v = a - SA - 1; c < a; ++c) { ISA[*c] = v; }
+            if(1 < (a - first)) { LA[a - SA - 1] = LA[last - SA - 1]; }
+          }
+          if(b < last) {
+            for(c = a, v = b - SA - 1; c < b; ++c) { ISA[*c] = v; }
+            if(1 < (b - a)) { LA[b - SA - 1] = LA[last - SA - 1]; }
+          }
+  
+          /* push */
+          if(1 < (b - a)) {
+            STACK_PUSH(NULL, a, b, partial_copy_flag);
+            STACK_PUSH(ISAd, first, last, -2);
+            partial_copy_flag = 0;
+          }
+          if((a - first) <= (last - b)) {
+            if(1 < (a - first)) {
+              STACK_PUSH(ISAd, b, last, -4);
+              STACK_PUSH(ISAd, b, last, tr_ilg(last - b));
+              STACK_PUSH(ISAd, first, a, -4);
+              last = a, limit = tr_ilg(a - first);
+            } else if(1 < (last - b)) {
+              STACK_PUSH(ISAd, b, last, -4);
+              first = b, limit = tr_ilg(last - b);
+            } else {
+              STACK_POP(ISAd, first, last, limit);
+            }
+          } else {
+            if(1 < (last - b)) {
+              STACK_PUSH(ISAd, first, a, -4);
+              STACK_PUSH(ISAd, first, a, tr_ilg(a - first));
+              STACK_PUSH(ISAd, b, last, -4);
+              first = b, limit = tr_ilg(last - b);
+            } else if(1 < (a - first)) {
+              STACK_PUSH(ISAd, first, a, -4);
+              last = a, limit = tr_ilg(a - first);
+            } else {
+              STACK_POP(ISAd, first, last, limit);
+            }
+          }
+      }
+      continue;
+    }
+
+    if((last - first) <= TR_INSERTIONSORT_THRESHOLD) {
+      tr_insertionsort(ISAd, first, last);
+      limit = -3;
+      continue;
+    }
+
+    if(limit-- == 0) {
+      tr_heapsort(ISAd, first, last - first);
+      for(a = last - 1; first < a; a = b) {
+        for(x = ISAd[*a], b = a - 1; (first <= b) && (ISAd[*b] == x); --b) { *b = ~*b; }
+      }
+      limit = -3;
+      continue;
+    }
+
+    /* choose pivot */
+    a = tr_pivot(ISAd, first, last);
+    SWAP(*first, *a);
+    v = ISAd[*first];
+
+    /* partition */
+    tr_partition(ISAd, first, first + 1, last, &a, &b, v);
+
+      /* push */
+      if((1 < (b - a)) && (ISA[*a] != v) && (LA[v] > 0) && (trbudget_check(budget, b - a))) {
+        incr = LA[v];
+        next = tr_ilg(b - a);
+        if((a - first) <= (last - b)) {
+          if((last - b) <= (b - a)) {
+            if(1 < (a - first)) {
+              STACK_PUSH(ISAd + incr, a, b, next);
+              STACK_PUSH(ISAd, b, last, limit);
+              last = a;
+            } else if(1 < (last - b)) {
+              STACK_PUSH(ISAd + incr, a, b, next);
+              first = b;
+            } else {
+              ISAd += incr, first = a, last = b, limit = next;
+            }
+          } else if((a - first) <= (b - a)) {
+            if(1 < (a - first)) {
+              STACK_PUSH(ISAd, b, last, limit);
+              STACK_PUSH(ISAd + incr, a, b, next);
+              last = a;
+            } else {
+              STACK_PUSH(ISAd, b, last, limit);
+              ISAd += incr, first = a, last = b, limit = next;
+            }
+          } else {
+            STACK_PUSH(ISAd, b, last, limit);
+            STACK_PUSH(ISAd, first, a, limit);
+            ISAd += incr, first = a, last = b, limit = next;
+          }
+        } else {
+          if((a - first) <= (b - a)) {
+            if(1 < (last - b)) {
+              STACK_PUSH(ISAd + incr, a, b, next);
+              STACK_PUSH(ISAd, first, a, limit);
+              first = b;
+            } else if(1 < (a - first)) {
+              STACK_PUSH(ISAd + incr, a, b, next);
+              last = a;
+            } else {
+              ISAd += incr, first = a, last = b, limit = next;
+            }
+          } else if((last - b) <= (b - a)) {
+            if(1 < (last - b)) {
+              STACK_PUSH(ISAd, first, a, limit);
+              STACK_PUSH(ISAd + incr, a, b, next);
+              first = b;
+            } else {
+              STACK_PUSH(ISAd, first, a, limit);
+              ISAd += incr, first = a, last = b, limit = next;
+            }
+          } else {
+            STACK_PUSH(ISAd, first, a, limit);
+            STACK_PUSH(ISAd, b, last, limit);
+            ISAd += incr, first = a, last = b, limit = next;
+          }
+        }
+      } else {
+        if(1 < (b - a)) {
+          if((ISA[*a] == v) && trbudget_check(budget, b - a)) {
+            LA[b - SA - 1] = ISAd - ISA;
+          } else {
+            LA[b - SA - 1] = 0;
+            partial_copy_flag = 1;
+          }
+          *(b - 1) = ~*(b - 1);
+          LA[b - SA - 2] = a - SA;
+        }
+        if((a - first) <= (last - b)) {
+          if(1 < (a - first)) {
+            STACK_PUSH(ISAd, b, last, limit);
+            last = a;
+          } else if(1 < (last - b)) {
+            first = b;
+          } else {
+            STACK_POP(ISAd, first, last, limit);
+          }
+        } else {
+          if(1 < (last - b)) {
+            STACK_PUSH(ISAd, first, a, limit);
+            first = b;
+          } else if(1 < (a - first)) {
+            last = a;
+          } else {
+            STACK_POP(ISAd, first, last, limit);
+          }
+        }
+      }
+  }
+#undef STACK_SIZE
+}
+
+static INLINE
+void
+tr_two_stage_introsort(saidx_t *ISA, const saidx_t *ISAd,
+                       saidx_t *SA, saidx_t *LA,
+                       saidx_t *first, saidx_t *last,
+                       trbudget_t *budget) {
+
+  tr_first_stage_introsort(ISA, ISAd, SA, LA, first, last, budget);
+  tr_second_stage_introsort(ISA, ISAd, SA, LA, first, last, budget);
+}
+
 
 
 /*---------------------------------------------------------------------------*/
@@ -552,7 +1045,7 @@ tr_introsort(saidx_t *ISA, const saidx_t *ISAd,
 
 /* Tandem repeat sort */
 void
-trsort(saidx_t *ISA, saidx_t *SA, saidx_t n, saidx_t depth) {
+trsort(saidx_t *ISA, saidx_t *SA, saidx_t n, saidx_t depth, saidx_t *buf, saidx_t bufsize) {
   saidx_t *ISAd;
   saidx_t *first, *last;
   trbudget_t budget;
@@ -571,8 +1064,12 @@ trsort(saidx_t *ISA, saidx_t *SA, saidx_t n, saidx_t depth) {
         last = SA + ISA[t] + 1;
         if(1 < (last - first)) {
           budget.count = 0;
-          tr_introsort(ISA, ISAd, SA, first, last, &budget);
-          if(budget.count != 0) { unsorted += budget.count; }
+          if(bufsize < (last - first)) {
+            tr_introsort(ISA, ISAd, SA, first, last, &budget);
+          } else {
+            tr_two_stage_introsort(ISA, ISAd, SA, buf - (first - SA), first, last, &budget);
+          }
+          if(budget.count != 0) { unsorted += budget.count; bufsize = 0; }
           else { skip = first - last; }
         } else if((last - first) == 1) {
           skip = -1;