text: 'use strict'
  
  module.exports = boxIntersectIter
  
  var pool = require('typedarray-pool')
  var bits = require('bit-twiddle')
  var bruteForce = require('./brute')
  var bruteForcePartial = bruteForce.partial
  var bruteForceFull = bruteForce.full
  var sweep = require('./sweep')
  var findMedian = require('./median')
  var genPartition = require('./partition')
  
  //Twiddle parameters
  var BRUTE_FORCE_CUTOFF    = 128       //Cut off for brute force search
  var SCAN_CUTOFF           = (1<<22)   //Cut off for two way scan
  var SCAN_COMPLETE_CUTOFF  = (1<<22)  
  
  //Partition functions
  var partitionInteriorContainsInterval = genPartition(
    '!(lo>=p0)&&!(p1>=hi)', 
    ['p0', 'p1'])
  
  var partitionStartEqual = genPartition(
    'lo===p0',
    ['p0'])
  
  var partitionStartLessThan = genPartition(
    'lo<p0',
    ['p0'])
  
  var partitionEndLessThanEqual = genPartition(
    'hi<=p0',
    ['p0'])
  
  var partitionContainsPoint = genPartition(
    'lo<=p0&&p0<=hi',
    ['p0'])
  
  var partitionContainsPointProper = genPartition(
    'lo<p0&&p0<=hi',
    ['p0'])
  
  //Frame size for iterative loop
  var IFRAME_SIZE = 6
  var DFRAME_SIZE = 2
  
  //Data for box statck
  var INIT_CAPACITY = 1024
  var BOX_ISTACK  = pool.mallocInt32(INIT_CAPACITY)
  var BOX_DSTACK  = pool.mallocDouble(INIT_CAPACITY)
  
  //Initialize iterative loop queue
  function iterInit(d, count) {
    var levels = (8 * bits.log2(count+1) * (d+1))|0
    var maxInts = bits.nextPow2(IFRAME_SIZE*levels)
    if(BOX_ISTACK.length < maxInts) {
      pool.free(BOX_ISTACK)
      BOX_ISTACK = pool.mallocInt32(maxInts)
    }
    var maxDoubles = bits.nextPow2(DFRAME_SIZE*levels)
    if(BOX_DSTACK < maxDoubles) {
      pool.free(BOX_DSTACK)
      BOX_DSTACK = pool.mallocDouble(maxDoubles)
    }
  }
  
  //Append item to queue
  function iterPush(ptr,
    axis, 
    redStart, redEnd, 
    blueStart, blueEnd, 
    state, 
    lo, hi) {
  
    var iptr = IFRAME_SIZE * ptr
    BOX_ISTACK[iptr]   = axis
    BOX_ISTACK[iptr+1] = redStart
    BOX_ISTACK[iptr+2] = redEnd
    BOX_ISTACK[iptr+3] = blueStart
    BOX_ISTACK[iptr+4] = blueEnd
    BOX_ISTACK[iptr+5] = state
  
    var dptr = DFRAME_SIZE * ptr
    BOX_DSTACK[dptr]   = lo
    BOX_DSTACK[dptr+1] = hi
  }
  
  //Special case:  Intersect single point with list of intervals
  function onePointPartial(
    d, axis, visit, flip,
    redStart, redEnd, red, redIndex,
    blueOffset, blue, blueId) {
  
    var elemSize = 2 * d
    var bluePtr  = blueOffset * elemSize
    var blueX    = blue[bluePtr + axis]
  
  red_loop:
    for(var i=redStart, redPtr=redStart*elemSize; i<redEnd; ++i, redPtr+=elemSize) {
      var r0 = red[redPtr+axis]
      var r1 = red[redPtr+axis+d]
      if(blueX < r0 || r1 < blueX) {
        continue
      }
      if(flip && blueX === r0) {
        continue
      }
      var redId = redIndex[i]
      for(var j=axis+1; j<d; ++j) {
        var r0 = red[redPtr+j]
        var r1 = red[redPtr+j+d]
        var b0 = blue[bluePtr+j]
        var b1 = blue[bluePtr+j+d]
        if(r1 < b0 || b1 < r0) {
          continue red_loop
        }
      }
      var retval
      if(flip) {
        retval = visit(blueId, redId)
      } else {
        retval = visit(redId, blueId)
      }
      if(retval !== void 0) {
        return retval
      }
    }
  }
  
  //Special case:  Intersect one point with list of intervals
  function onePointFull(
    d, axis, visit,
    redStart, redEnd, red, redIndex,
    blueOffset, blue, blueId) {
  
    var elemSize = 2 * d
    var bluePtr  = blueOffset * elemSize
    var blueX    = blue[bluePtr + axis]
  
  red_loop:
    for(var i=redStart, redPtr=redStart*elemSize; i<redEnd; ++i, redPtr+=elemSize) {
      var redId = redIndex[i]
      if(redId === blueId) {
        continue
      }
      var r0 = red[redPtr+axis]
      var r1 = red[redPtr+axis+d]
      if(blueX < r0 || r1 < blueX) {
        continue
      }
      for(var j=axis+1; j<d; ++j) {
        var r0 = red[redPtr+j]
        var r1 = red[redPtr+j+d]
        var b0 = blue[bluePtr+j]
        var b1 = blue[bluePtr+j+d]
        if(r1 < b0 || b1 < r0) {
          continue red_loop
        }
      }
      var retval = visit(redId, blueId)
      if(retval !== void 0) {
        return retval
      }
    }
  }
  
  //The main box intersection routine
  function boxIntersectIter(
    d, visit, initFull,
    xSize, xBoxes, xIndex,
    ySize, yBoxes, yIndex) {
  
    //Reserve memory for stack
    iterInit(d, xSize + ySize)
  
    var top  = 0
    var elemSize = 2 * d
    var retval
  
    iterPush(top++,
        0,
        0, xSize,
        0, ySize,
        initFull ? 16 : 0, 
        -Infinity, Infinity)
    if(!initFull) {
      iterPush(top++,
        0,
        0, ySize,
        0, xSize,
        1, 
        -Infinity, Infinity)
    }
  
    while(top > 0) {
      top  -= 1
  
      var iptr = top * IFRAME_SIZE
      var axis      = BOX_ISTACK[iptr]
      var redStart  = BOX_ISTACK[iptr+1]
      var redEnd    = BOX_ISTACK[iptr+2]
      var blueStart = BOX_ISTACK[iptr+3]
      var blueEnd   = BOX_ISTACK[iptr+4]
      var state     = BOX_ISTACK[iptr+5]
  
      var dptr = top * DFRAME_SIZE
      var lo        = BOX_DSTACK[dptr]
      var hi        = BOX_DSTACK[dptr+1]
  
      //Unpack state info
      var flip      = (state & 1)
      var full      = !!(state & 16)
  
      //Unpack indices
      var red       = xBoxes
      var redIndex  = xIndex
      var blue      = yBoxes
      var blueIndex = yIndex
      if(flip) {
        red         = yBoxes
        redIndex    = yIndex
        blue        = xBoxes
        blueIndex   = xIndex
      }
  
      if(state & 2) {
        redEnd = partitionStartLessThan(
          d, axis,
          redStart, redEnd, red, redIndex,
          hi)
        if(redStart >= redEnd) {
          continue
        }
      }
      if(state & 4) {
        redStart = partitionEndLessThanEqual(
          d, axis,
          redStart, redEnd, red, redIndex,
          lo)
        if(redStart >= redEnd) {
          continue
        }
      }
      
      var redCount  = redEnd  - redStart
      var blueCount = blueEnd - blueStart
  
      if(full) {
        if(d * redCount * (redCount + blueCount) < SCAN_COMPLETE_CUTOFF) {
          retval = sweep.scanComplete(
            d, axis, visit, 
            redStart, redEnd, red, redIndex,
            blueStart, blueEnd, blue, blueIndex)
          if(retval !== void 0) {
            return retval
          }
          continue
        }
      } else {
        if(d * Math.min(redCount, blueCount) < BRUTE_FORCE_CUTOFF) {
          //If input small, then use brute force
          retval = bruteForcePartial(
              d, axis, visit, flip,
              redStart,  redEnd,  red,  redIndex,
              blueStart, blueEnd, blue, blueIndex)
          if(retval !== void 0) {
            return retval
          }
          continue
        } else if(d * redCount * blueCount < SCAN_CUTOFF) {
          //If input medium sized, then use sweep and prune
          retval = sweep.scanBipartite(
            d, axis, visit, flip, 
            redStart, redEnd, red, redIndex,
            blueStart, blueEnd, blue, blueIndex)
          if(retval !== void 0) {
            return retval
          }
          continue
        }
      }
      
      //First, find all red intervals whose interior contains (lo,hi)
      var red0 = partitionInteriorContainsInterval(
        d, axis, 
        redStart, redEnd, red, redIndex,
        lo, hi)
  
      //Lower dimensional case
      if(redStart < red0) {
  
        if(d * (red0 - redStart) < BRUTE_FORCE_CUTOFF) {
          //Special case for small inputs: use brute force
          retval = bruteForceFull(
            d, axis+1, visit,
            redStart, red0, red, redIndex,
            blueStart, blueEnd, blue, blueIndex)
          if(retval !== void 0) {
            return retval
          }
        } else if(axis === d-2) {
          if(flip) {
            retval = sweep.sweepBipartite(
              d, visit,
              blueStart, blueEnd, blue, blueIndex,
              redStart, red0, red, redIndex)
          } else {
            retval = sweep.sweepBipartite(
              d, visit,
              redStart, red0, red, redIndex,
              blueStart, blueEnd, blue, blueIndex)
          }
          if(retval !== void 0) {
            return retval
          }
        } else {
          iterPush(top++,
            axis+1,
            redStart, red0,
            blueStart, blueEnd,
            flip,
            -Infinity, Infinity)
          iterPush(top++,
            axis+1,
            blueStart, blueEnd,
            redStart, red0,
            flip^1,
            -Infinity, Infinity)
        }
      }
  
      //Divide and conquer phase
      if(red0 < redEnd) {
  
        //Cut blue into 3 parts:
        //
        //  Points < mid point
        //  Points = mid point
        //  Points > mid point
        //
        var blue0 = findMedian(
          d, axis, 
          blueStart, blueEnd, blue, blueIndex)
        var mid = blue[elemSize * blue0 + axis]
        var blue1 = partitionStartEqual(
          d, axis,
          blue0, blueEnd, blue, blueIndex,
          mid)
  
        //Right case
        if(blue1 < blueEnd) {
          iterPush(top++,
            axis,
            red0, redEnd,
            blue1, blueEnd,
            (flip|4) + (full ? 16 : 0),
            mid, hi)
        }
  
        //Left case
        if(blueStart < blue0) {
          iterPush(top++,
            axis,
            red0, redEnd,
            blueStart, blue0,
            (flip|2) + (full ? 16 : 0),
            lo, mid)
        }
  
        //Center case (the hard part)
        if(blue0 + 1 === blue1) {
          //Optimization: Range with exactly 1 point, use a brute force scan
          if(full) {
            retval = onePointFull(
              d, axis, visit,
              red0, redEnd, red, redIndex,
              blue0, blue, blueIndex[blue0])
          } else {
            retval = onePointPartial(
              d, axis, visit, flip,
              red0, redEnd, red, redIndex,
              blue0, blue, blueIndex[blue0])
          }
          if(retval !== void 0) {
            return retval
          }
        } else if(blue0 < blue1) {
          var red1
          if(full) {
            //If full intersection, need to handle special case
            red1 = partitionContainsPoint(
              d, axis,
              red0, redEnd, red, redIndex,
              mid)
            if(red0 < red1) {
              var redX = partitionStartEqual(
                d, axis,
                red0, red1, red, redIndex,
                mid)
              if(axis === d-2) {
                //Degenerate sweep intersection:
                //  [red0, redX] with [blue0, blue1]
                if(red0 < redX) {
                  retval = sweep.sweepComplete(
                    d, visit,
                    red0, redX, red, redIndex,
                    blue0, blue1, blue, blueIndex)
                  if(retval !== void 0) {
                    return retval
                  }
                }
  
                //Normal sweep intersection:
                //  [redX, red1] with [blue0, blue1]
                if(redX < red1) {
                  retval = sweep.sweepBipartite(
                    d, visit,
                    redX, red1, red, redIndex,
                    blue0, blue1, blue, blueIndex)
                  if(retval !== void 0) {
                    return retval
                  }
                }
              } else {
                if(red0 < redX) {
                  iterPush(top++,
                    axis+1,
                    red0, redX,
                    blue0, blue1,
                    16,
                    -Infinity, Infinity)
                }
                if(redX < red1) {
                  iterPush(top++,
                    axis+1,
                    redX, red1,
                    blue0, blue1,
                    0,
                    -Infinity, Infinity)
                  iterPush(top++,
                    axis+1,
                    blue0, blue1,
                    redX, red1,
                    1,
                    -Infinity, Infinity)
                }
              }
            }
          } else {
            if(flip) {
              red1 = partitionContainsPointProper(
                d, axis,
                red0, redEnd, red, redIndex,
                mid)
            } else {
              red1 = partitionContainsPoint(
                d, axis,
                red0, redEnd, red, redIndex,
                mid)
            }
            if(red0 < red1) {
              if(axis === d-2) {
                if(flip) {
                  retval = sweep.sweepBipartite(
                    d, visit,
                    blue0, blue1, blue, blueIndex,
                    red0, red1, red, redIndex)
                } else {
                  retval = sweep.sweepBipartite(
                    d, visit,
                    red0, red1, red, redIndex,
                    blue0, blue1, blue, blueIndex)
                }
              } else {
                iterPush(top++,
                  axis+1,
                  red0, red1,
                  blue0, blue1,
                  flip,
                  -Infinity, Infinity)
                iterPush(top++,
                  axis+1,
                  blue0, blue1,
                  red0, red1,
                  flip^1,
                  -Infinity, Infinity)
              }
            }
          }
        }
      }
    }
  }