Interpolation Search

Interpolation search is a value-aware search algorithm within computation and algorithms for locating a numerical key x inside a sorted array A[0 .. n − 1] by estimating the target's array position proportionally from the boundary key values A[low] and A[high]. First described by W. Wesley Peterson in 1957 in the context of open-addressed hash tables and sequential index files, and rigorously analyzed by Andrew Chi-Chih Yao and F. Frances Yao (1976) and by Yehoshua Perl, Alon Itai, and Haim Avni (1978), interpolation search models the intuitive strategy used when looking up a name in a telephone directory: rather than opening the book to the midpoint as in divide-and-conquer binary search, one probes near the front for a surname beginning with "B" and near the back for "W".

When the n sorted keys are drawn independently from a uniform distribution (or any smooth distribution transformed via its cumulative distribution function), each linear interpolation probe squares the relative width of the remaining uncertainty interval, reducing the expected number of array probes from ⌊log₂ n⌋ + 1 down to log₂ log₂ n + O(1)—an asymptotic average-case complexity of O(log log n). Conversely, on adversarially skewed or geometrically growing keys where linear interpolation repeatedly steps by a single index, unmodified interpolation search degrades to O(n) worst-case time unless safeguarded by hybrid alternation with binary or exponential search or by Dan E. Willard's O(√(log n)) Santoro–Sidney interpolation-binary scheme.

Proportional Position Estimation

Let A[0 .. n − 1] be a sorted array of distinct numerical keys and let [low, high] denote the current active search interval satisfying A[low] ≤ x ≤ A[high]. Under the hypothesis that key values grow approximately linearly with array index across [low, high], the fraction of the index interval (pos − low) / (high − low) should match the fraction of the key interval (x − A[low]) / (A[high] − A[low]). Solving for the integer probe coordinate pos yields the linear interpolation probe formula:

At each iteration, the algorithm checks that low ≤ high and A[low] ≤ x ≤ A[high] ( guarding against division by zero when A[low] = A[high]):

1. Evaluate pos using 64-bit or 128-bit integer arithmetic (or floating-point division) so the product (x − A[low]) · (high − low) cannot overflow.
2. Compare the target x against A[pos]:
   • If A[pos] = x, return pos.
   • If A[pos] < x, narrow the search interval to the upper subrange by setting low ← pos + 1.
   • If A[pos] > x, narrow the search interval to the lower subrange by setting high ← pos − 1.

Order Statistics of Uniform Keys and Conditional Variance

Why does linear interpolation shrink the active interval from size m = high − low + 1 to size O(√m) per probe on average? Suppose m keys are drawn independently and uniformly from a continuous interval [a, b] and sorted in ascending order. By the classical theory of Beta order statistics, conditioned on the endpoint values A[low] = a and A[high] = b, the number of keys K falling below the target value x follows a binomial distribution Bin(m − 2, p) where p = (x − a) / (b − a) is the normalized interpolation ratio.

A binomial random variable K ~ Bin(N, p) has mean μ = Np (which is precisely the index offset selected by the interpolation probe pos − low) and standard deviation:

By Chebyshev's inequality and Chernoff tail bounds, the true rank of x inside [low, high] deviates from the interpolated probe index pos by only O(√m) indices in expectation. Consequently, when the next probe (or a secondary bracket check) isolates the subinterval containing x, the effective number of candidate elements drops from m to approximately √m.

Doubly Logarithmic Convergence and the Yao–Yao Lower Bound

Because each interpolation step maps an uncertainty scale of m_k elements to an expected scale of m_k+1 ≈ √m_k = m_k^1/2, after k successive probes an initial array of size n contracts to expected size:

Reducing the active interval down to a constant number of elements (m_k ≤ 2) requires n^2−k ≤ 2, or 2^−k log₂ n ≤ 1. Taking base-2 logarithms a second time gives:

Perl, Itai, and Avni (1978) proved via Martingale potential functions on the conditional variance that the expected number of probes for interpolation search on n independent uniform keys is at most log₂ log₂ n + 3. In a companion theoretical landmark, Andrew Yao and Frances Yao (1976) proved that any search algorithm using arbitrary arithmetic functions of the probed keys must perform at least log₂ log₂ n + O(1) probes on average over uniform keys, establishing that simple linear interpolation is asymptotically optimal. For an array of n = 2³² ≈ 4.29 × 10⁹ uniform 32-bit integers, log₂ log₂(2³²) = 5, so interpolation search inspects roughly 5 to 6 array slots compared with 32 probes for binary search.

Adversarial Skew and Worst-Case O(n) Degradation

When the sorted array violates uniformity—for example, when keys follow a geometric progression A[j] = 10^j (1, 10, 100, …, 10ⁿ⁻¹) and the target is x = 10ⁿ⁻²—the denominator A[high] − A[low] = 10ⁿ⁻¹ − 1 is dominated by the final element. For every target smaller than 10ⁿ⁻², the ratio (x − A[low]) / (A[high] − A[low]) is less than 0.01, so for any array of size n ≤ 100 the formula evaluates to pos = low + 0 = low. Setting low ← pos + 1 advances the left boundary by only a single element per step, degrading interpolation search to n − 1 probes—equivalent to linear search.

Quadratic Binary Search and Interpolation-Exponential Hybrids

To eliminate the O(n) worst case while retaining O(log log n) average performance on uniform keys, two classical modifications are used:

1. Quadratic Binary Search (Gonnet, Rogers, and George, 1980) and Interpolation-Exponential Search: After computing the interpolated probe pos and observing, say, A[pos] < x, instead of re-interpolating across the entire remaining suffix [pos + 1, high], the algorithm walks outward from pos in multiples of √m (or in powers of two pos + 2⁰√m, pos + 2¹√m, …) to bracket x within a narrow window before recursing.
2. Santoro–Sidney and Willard's θ-Search (1985): Nicola Santoro and Jeffrey B. Sidney (1985) and Dan E. Willard (1985) interleaved each interpolation probe with a binary probe (or a guarded fraction probe inside [low + θm, high − θm]), proving that the worst-case number of probes is bounded by O(log n) (or O(√(log n)) under Willard's adaptive fraction schedule) on arbitrary key distributions while still finishing in O(log log n) expected probes on uniform data.

Division Latency vs. Memory-Hierarchy Latency and Learned Indexes

On modern superscalar CPUs, a 64-bit integer division (idiv) requires 15 to 40 clock cycles, whereas binary search uses a single 1-cycle right shift (shr). For tiny arrays that fit inside the L1 cache (n ≤ 256), binary search often finishes faster in wall-clock time despite doing more probes. Once n exceeds the L3 cache (n ≥ 10⁷) or resides on SSD/NVMe storage alongside structures like the B-tree, van Emde Boas tree, or y-fast trie, each random probe incurs a 60–100 ns DRAM miss or microsecond I/O stall, and eliminating 80% of memory probes via interpolation (such as Three-Point Interpolation Search, TIP, or piecewise linear learned indexes) delivers substantial throughput gains.

How does interpolation search avoid integer overflow when computing the probe index?

Because the product (x − A[low]) · (high − low) multiplies a key difference by an array length, both of which may be 32-bit or 64-bit integers, production implementations either promote both operands to a wider integer type (such as 128-bit integer arithmetic for 64-bit keys) before dividing by (A[high] − A[low]), or compute the normalized ratio as an IEEE 754 double-precision float in [0.0, 1.0] and multiply by (high − low) before clamping the result to [low, high].

Can interpolation search be applied to non-uniform distributions or lexicographic strings?

Yes, provided an inexpensive monotonic transformation maps the keys to an approximately uniform space. If keys follow a known cumulative distribution function F(x), applying the probability integral transform u = F(x) converts the keys into uniform [0, 1] variates before interpolating. For fixed-prefix ASCII or UTF-8 strings, interpreting the first 8 bytes of each string as a big-endian 64-bit integer enables numerical interpolation across lexicographically sorted string tables, though common prefixes require stripping shared leading bytes at each narrowing step.

Why must low be updated to pos + 1 rather than pos when A[pos] < x?

When x − A[low] is small relative to A[high] − A[low], the floor operation in the interpolation formula frequently rounds the offset down to 0, yielding pos = low. If the algorithm set low ← pos instead of low ← pos + 1 after confirming A[pos] < x, the next iteration would recompute the identical probe pos = low and enter an infinite loop. Advancing low ← pos + 1 and high ← pos − 1 guarantees that every unsuccessful probe eliminates at least one element and that the search terminates in finite steps.

1. Peterson, W. W. (1957). "Addressing for random-access storage." IBM Journal of Research and Development, 1(2), 130–146. https://doi.org/10.1147/rd.12.0130
2. Yao, A. C.-C., & Yao, F. F. (1976). "The complexity of searching an ordered random table." In 17th Annual Symposium on Foundations of Computer Science (SFCS 1976), 173–177. https://doi.org/10.1109/SFCS.1976.32
3. Perl, Y., Itai, A., & Avni, H. (1978). "Interpolation search—a log log N search." Communications of the ACM, 21(7), 550–553. https://doi.org/10.1145/359545.359557
4. Gonnet, G. H., Rogers, L. D., & George, J. A. (1980). "An algorithmic and complexity analysis of interpolation search." Acta Informatica, 13(1), 39–52. https://doi.org/10.1007/BF00288534
5. Willard, D. E. (1985). "Searching unindexed and nonuniformly generated files in log log N time." SIAM Journal on Computing, 14(4), 1013–1029. https://doi.org/10.1137/0214071