Strengthening Invariants for Efficient Computation ⋆

Yanhong A. Liu^{a,1}, Scott D. Stoller^{a}, Tim Teitelbaum^{b}

aComputer Science Department, State University of New York at Stony Brook, Stony Brook, NY 11794, USA

bDepartment of Computer Science, Cornell University, Ithaca, NY 14853, USA

Abstract

This paper presents program analyses and transformations for strengthening invariants for the purpose of efficient computation. Finding the stronger invariants corresponds to discovering a general class of auxiliary information for any incremental computation problem. Combining the techniques with previous techniques for caching intermediate results, we obtain a systematic approach that transforms non-incremental programs into efficient incremental programs that use and maintain useful auxiliary information as well as useful intermediate results. The use of auxiliary information allows us to achieve a greater degree of incrementality than otherwise possible. Applications of the approach include strength reduction in optimizing compilers and finite differencing in transformational programming.

1 Introduction

Efficient computation via incremental computation. In essence, every program computes by fixed-point iteration, expressed as recursive functions or loops. This is why loop optimizations are so important. A loop body can be regarded as a program f parameterized by an induction variable x that is incremented on each iteration by a change operation ⊕. Efficient iterative

⋆ This work was supported in part by ONR under grants N00014-92-J-1973, N00014-99-1-0132, and N00014-99-1-0358 and by NSF under grants CCR-9503319, CCR-9711253, and CCR-9876058. This article is a revised and extended version of a paper that appeared in Conference Record of the 23rd Annual ACM Symposium on Principles of Programming Languages, St. Petersburg Beach, Florida, January 1996. 1 Corresponding author. E-mail: liu@cs.sunysb.edu

Preprint submitted to Elsevier Preprint 4 April 2001

computation relies on effective use of state, i.e., computing the result of each iteration incrementally using stored results of previous iterations. This is why strength reduction [3], finite differencing [64], and related techniques are crucial for performance.

Given a program f and an input change operation ⊕, a program f′ that computes f(x ⊕ y) efficiently by using the result of the previous computation of f(x) is called an incremental version of f under ⊕. Often, information other than the result of f(x) needs to be maintained and used for efficient incremental computation of f(x ⊕ y). We call a program that computes such information an extended version of f. Thus, the goal of computing loops efficiently corresponds to constructing an extended version of a program f and deriving an incremental version of the extended version under an input change operation ⊕.

In general, incremental computation aims to solve a problem on a sequence of inputs that differ only slightly from one another, making use of the previously computed output in computing a new output, instead of computing the new output from scratch. Incremental computation is a fundamental issue relevant throughout computer software, e.g., optimizing compilers [2,3,17,23,78], transformational program development [8,20,62,65,77], and interactive systems [5,6,10,22,33,41,71,72]. Numerous techniques for incremental computation have been developed, e.g., [3,4,25,34–36,56,64,68,70,73,76,79,86].

Strengthening invariants for incrementalization. We are engaged in an ambitious effort to derive incremental extended programs automatically (or semi-automatically) from non-incremental programs written in standard programming languages.This approach contrasts with many other approaches that aim to evaluate non-incremental programs incrementally. We call this approach incrementalization. We have partitioned the core of the problem into three subproblems:

• P1. Exploiting the result, i.e., the return value, of f(x).

• P2. Caching, maintaining, and exploiting intermediate results of the computation f(x).

• P3. Discovering, computing, maintaining, and exploiting auxiliary information about x, i.e., values not computed by f(x).

Our current approaches to problems P1 and P2 are described in [56] and [54], respectively. In this paper, we address problem P3 and contribute:

• A novel proposal for finding auxiliary information.

• A comprehensive methodology for deriving incremental programs that addresses all three subproblems.

Since auxiliary information is not computed by the original body of computation f(x), adding it strengthens the invariants that hold over the iterative computation that uses f(x).

Some approaches to efficient computation have exploited specific kinds of auxiliary information for strengthening invariants, e.g., auxiliary arithmetic associated with some classical strength-reduction rules [3], auxiliary maps maintained by finite differencing rules for aggregate primitives in SETL [64] and INC [86], and auxiliary data structures for problems with certain properties like stable decomposition [70]. However, systematic discovery of auxiliary information for arbitrary programs has been a subject completely open for study.

Auxiliary information is, by definition, useful information about x that is not computed by f(x). Where, then, can one find it? The key insight of our proposal is:

A. Consider, as candidate auxiliary information for f, all intermediate results of an incremental version of f that depend only on x; such an incremental version can be obtained using some techniques we developed for solving P1 and P2. P2 is included here so that candidate auxiliary information useful for efficiently maintaining intermediate results is also included.

How can one discover which pieces of candidate auxiliary information are useful and how they can be used? We propose:

B. Extend f with all candidate auxiliary information, and then apply some techniques used in our methods for P1 and P2 to obtain an extended version and an incremental extended version that together compute, exploit, and maintain only useful intermediate results and useful auxiliary information.

Thus, on the one hand, one can regard the method for P3 in this paper as an extension to methods for P1 and P2. On the other hand, one can regard methods for P1 and P2 (suitably revised for their different applications here) as aids for solving P3. The modular components complement one another to form a comprehensive principled approach for incrementalization and therefore also for achieving efficient iterative computation generally. Although the entire approach seems complex, each module or step is simple.

We summarize here the essence of our methods:

P1. In [56], we gave a systematic transformational approach for deriving an incremental version f′ of a program f under an input change ⊕. The basic idea is to identify in the computation of f(x ⊕ y) those subcomputations that are also performed in the computation of f(x) and whose values can be retrieved from the cached result r of f(x). The computation of f(x ⊕ y) is symbolically transformed to avoid re-performing these subcomputations by replacing them with corresponding retrievals. This efficient way of computing f(x ⊕ y) is captured in the definition of f′(x, y, r).

P2. In [54], we gave a method, called cache-and-prune, for statically transforming programs to cache all intermediate results useful for incremental computation. The basic idea is to (I) extend the program f to a program f̄ that returns all intermediate results, (II) incrementalize the program f̄ under ⊕ to obtain an incremental version f̄′ of f̄ using our method for P1, and (III) analyze the dependencies in f̄′, then prune the extended program f̄to a program f̂ that returns only the useful intermediate results, and prune the program f̄′ to obtain a program f̂′ that incrementally maintains only the useful intermediate results.

P3. This paper presents a two-phase method that discovers a general class of auxiliary information for any incremental computation problem. The two phases correspond to A and B above. For Phase A, we have developed an embedding analysis that helps avoid including redundant information in an extended version, and we have exploited a forward dependence analysis that helps identify candidate auxiliary information. All the program analyses and transformations used in this method are combined with considerations for caching intermediate results, so we obtain incremental extended programs that exploit and maintain intermediate results as well as auxiliary information.

We illustrate our approach by applying it to problems in list processing, VLSI design, graph algorithms, and other application areas.

The rest of this paper is organized as follows. Section 2 formulates the problem. Section 3 discusses discovering candidate auxiliary information. Section 4 describes how candidate auxiliary information is used. A number of examples are given in Section 6. Finally, we discuss related work and conclude in Section 7.

2 Formulating the problem

We use a simple first-order, call-by-value functional programming language. The expressions of the language are given by the following grammar:

e ::= v  variable
| c(e_{1}, ..., e_{n})  constructor application
| p(e_{1}, ..., e_{n})  primitive function application
| f(e_{1}, ..., e_{n})  function application
| if e_{1} then e_{2} else e_{3} conditional expression
| let v = e_{1} in e_{2}  binding expression

A program is a set of mutually recursive function definitions of the form

f(v_{1}, ..., v_{n}) ≜ e (1)

and a function f_{0} that is to be evaluated with some input x = ⟨x_{1}, ..., x_{n}⟩. Figure 1 gives some example definitions.

cmp compares sum of odd and product of even positions of list x. cmp(x) ≜ sum(odd(x)) ≤ prod(even(x))

odd(x) ≜ if null(x) then nil else cons(car(x), even(cdr(x)))
even(x) ≜ if null(x) then nil else odd(cdr(x))

sum(x) ≜ if null(x) then 0 else car(x) + sum(cdr(x))
prod(x) ≜ if null(x) then 1 else car(x) ∗ prod(cdr(x))

Fig. 1. Example function definitions.

An input change operation ⊕ to a function f_{0} combines an old input x = ⟨x_{1}, ..., x_{n}⟩ and a change y = ⟨y_{1}, ..., y_{m}⟩ to form a new input x′ = ⟨x′_{1}, ..., x′_{n}⟩ = x ⊕ y, where each x′_{i} is some function of x_{j}’s and y_{k}’s. For example, an input change operation to function cmp in Figure 1 may be defined by x′ = x ⊕ y = cons(y, x).

We use an asymptotic cost model for measuring time complexity and write t(f(v_{1}, ..., v_{n})) to denote the asymptotic time of computing f(v_{1}, ..., v_{n}). Thus, assuming all primitive functions take constant time, it suffices to consider only the values of function applications as candidate information to cache. Of course, maintaining extra information takes extra space. Our primary goal is to improve the asymptotic running time of the incremental computation. We save space by maintaining only information useful for achieving this.

Given a program f_{0} and an input change operation ⊕, we use the approach in [56] to derive an incremental version f′_{0} of f_{0} under ⊕, such that, if f_{0}(x) = r, then whenever f_{0}(x ⊕ y) returns a value, f′_{0}(x, y, r) returns the same value and is asymptotically at least as fast. While f_{0}(x) abbreviates f_{0}(x_{1}, ..., x_{n}), and f_{0}(x ⊕ y) abbreviates f_{0}(⟨x_{1}, ..., x_{n}⟩ ⊕ ⟨y_{1}, ..., y_{m}⟩), f′_{0}(x, y, r) abbreviates f′_{0}(x_{1}, ..., x_{n}, y_{1}, ..., y_{m}, r). Note that some of the parameters of f′_{0} may be dead and eliminated [56]. For example, for function sum of Figure 1 and input change operation x ⊕ y = cons(y, x), function sum′ in Figure 2 is derived.

In order to use also intermediate results of f_{0}(x) to compute f_{0}(x⊕y) possibly faster, we use the approach in [54] to cache useful intermediate results of f_{0} and obtain a program that incrementally computes the return value and maintains these intermediate results. For example, for function cmp of Figure 1 and input change operation x ⊕ ⟨y_{1}, y_{2}⟩ = cons(y_{1}, cons(y_{2}, x)), the intermediate results sum(odd(x)) and prod(even(x)) are cached, and functions ĉmp and ĉmp′ in Figure 2 are obtained.

Sometimes, auxiliary information other than the intermediate results of f_{0}(x) is needed to compute f_{0}(x ⊕ y) quickly. For example, for function cmp of Figure 1 and input change operation x⊕y = cons(y, x), the values of sum(even(x)) and prod(odd(x)), in addition to the intermediate values sum(odd(x)) and prod(even(x)), are crucial for computing cmp(cons(y, x)) incrementally but are not computed in cmp(x). Using the method in this paper, we can derive functions c̃mp and c̃mp′ in Figure 2 that compute these pieces of auxiliary information, use them in computing cmp(cons(y, x)), and maintain them as well. Function c̃mp′ computes incrementally using only O(1) time. We use this example as a running example.

If sum(x) = r, then sum′(y, r) = sum(cons(y, x)). For x of length n, sum′(y, r) takes time O(1); sum(cons(y, x)) takes time O(n). sum′(y, r) ≜ y + r

cmp(x) = 1st(_{cmp}(x)).

For x of length n, _{cmp}(x) takes time O(n); cmp(x) takes time O(n).

If _{cmp}(x) = r̂, then _{cmp}′(y_{1}, y_{2}, r̂) = _{cmp}(cons(y_{1}, cons(y_{2}, x))).

For x of length n, _{cmp}′(y_{1}, y_{2}, r̂) takes time O(1); _{cmp}(cons(y_{1}, cons(y_{2}, x))) takes time O(n).

_{cmp}(x) ≜ let v_{1} = sum(odd(x)) in let v_{2} = prod(even(x)) in < v_{1} ≤ v_{2}, v_{1}, v_{2} >

_{cmp}′(y_{1}, y_{2}, r̂) ≜ let v_{1} = y_{1}+2nd(r̂) in let v_{2} = y_{2}∗3rd(r̂) in < v_{1} ≤ v_{2}, v_{1}, v_{2} >

cmp(x) = 1st(_{cmp}(x)).

For x of length n, _{cmp}(x) takes time O(n); cmp(x) takes time O(n).

If _{cmp}(x) = _{r}, then _{cmp}′(y, _{r}) = _{cmp}(cons(y, x)).

For x of length n, _{cmp}′(y, _{r}) takes time O(1); _{cmp}(cons(y, x)) takes time O(n).

_{cmp}(x) ≜ let v_{1} = odd(x) in let u_{1} = sum(v_{1}) in let v_{2} = even(x) in let u_{2} = prod(v_{2}) in <u_{1} ≤ u_{2}, u_{1}, u_{2}, sum(v_{2}), prod(v_{1})>

_{cmp}′(y, _{r}) ≜ <\begin{matrix}y+4th(_{r}) ≤ 5th(_{r}), \\ y+4th(_{r}), 5th(_{r}), 2nd(_{r}), y∗3rd(_{r})> \end{matrix}

Fig. 2. Resulting function definitions.

Notation. We use <> to construct tuples that bundle intermediate results and auxiliary information with the original return value of a function. We use selectors 1st, 2nd, 3rd, ... to select the first, second, third, ... elements of such a tuple.

We use x to denote the previous input to f_{0}; r, the cached result of f_{0}(x); y, the input change parameter; x′, the new input x ⊕ y; and f′_{0}, an incremental version of f_{0} under ⊕. We let f̄_{0} return all intermediate results of f_{0}, and let f̌_{0} return candidate auxiliary information for f_{0} under ⊕. We use ¯ ; ˇf_{0} to denote a function that returns all intermediate results and candidate auxiliary information;¯ ; ř, the cached result of ¯ ; ˇf_{0}(x); and ¯ ; ˇf_{0}′, an incremental version of ¯ ; ˇf_{0} under ⊕. Finally, we use f̃_{0} to denote a function that returns only the useful intermediate results and auxiliary information; r̃, the cached result of f̃_{0}(x); and f̃_{0}′, a function that incrementally maintains only the useful intermediate results and auxiliary information. Note that (useful) intermediate results include the original return value. Table 1 summarizes the notation.

function  return value  denoted as incremental function
f_{0}  original value  r  f′_{0}
f̄_{0}  all i.r.  r̄
f̌_{0}  candidate a.i.  ř
\begin{matrix} & ¯ \\ ˇf_{0} \end{matrix}  all i.r. & candidate a.i.  \begin{matrix} & ¯ \\ ř \end{matrix}  \begin{matrix} & ¯ \\ ˇf_{0} \end{matrix}′
_{f_{0}}  useful i.r. & useful a.i.  _{r}  _{f_{0}}′

Table 1 Notation.

For convenience, we assume that bound variables have distinct names.

3 Phase A: Discovering candidate auxiliary information

Auxiliary information is, by definition, useful information not computed by the original program f_{0}, so it cannot be obtained directly from f_{0}. However, auxiliary information is information depending only on x that can speed up the computation of f_{0}(x ⊕ y). Seeking to obtain such information systematically, we come to the idea that when computing f_{0}(x ⊕ y), for example in the manner of incremental version f′_{0}(x, y, r), there are often subcomputations that depend only on x and r, but not on y, and whose values cannot be retrieved from the return value or intermediate results of f_{0}(x). If the values of these subcomputations were available, then we might make f′_{0} faster.

To obtain such candidate auxiliary information, the basic idea is to transform f_{0}(x ⊕ y) similarly as for deriving f′_{0} and to collect subcomputations in the transformed f_{0}(x ⊕ y) that depend only on x and whose values cannot be retrieved from the return value or intermediate results of f_{0}(x). Note that computing intermediate results of f_{0}(x) incrementally, with their corresponding auxiliary information, is often crucial for efficient incremental computation. Thus, we modify the basic idea just described so that it starts with f̄_{0}(x ⊕ y) instead of f_{0}(x ⊕ y).

Phase A has three steps. Step 1 extends f_{0} to a function f̄_{0} that caches all intermediate results. Step 2 transforms f̄_{0}(x ⊕ y) into a function f̄_{0}^{‵}, similar to f̄_{0}′, that exposes candidate auxiliary information. Step 3 constructs a function f̌_{0} that computes only the candidate auxiliary information in f̄_{0}^{‵}.

3.1 Step A.1: Caching all intermediate results

Extending f_{0} to cache all intermediate results uses the transformations in Stage I of [54]. It first performs a straightforward extension transformation to embed all intermediate results in the final return value and then performs administrative simplifications.

Certain improvements, suggested in [55] but not given in [55] or [54], can be made to the extension transformation. In particular, we can avoid caching redundant intermediate results, i.e., values of function applications that are already embedded in the values of their enclosing computations, since these omitted values can be retrieved from the results of the enclosing applications. These improvements are more important for discovering auxiliary information, since the resulting program should be much simpler and therefore easier to treat in subsequent analyses and transformations. These improvements also benefit the modified version of this extension transformation used in Step A.3.

We first briefly describe the extension transformation in [54]; then, we describe an embedding analysis that leads to the desired improvements to the extension transformation.

Extension transformation. For each function definition f(v_{1}, ..., v_{n}) ≜ e, we construct a function definition

f̄(v_{1}, ..., v_{n}) ≜ Ext[[e]] (2)

where Ext[[e]] extends an expression e to return the values of all function calls made in computing e, i.e., it considers subexpressions of e in applicative and left-to-right order, introduces bindings that name the results of function calls, builds up tuples of these values together with the values of the original subexpressions, and passes these values from subcomputations to enclosing computations.

The definition of Ext is given in Figure 3. We assume that each introduced binding uses a fresh variable name. For a constructed tuple <>, while we use 1st to return the first element, which is the original return value, we use rst to return a tuple of the remaining elements, which are the corresponding intermediate results. We use an infix operation @ to concatenate two tuples. For transforming a conditional expression, the transformation Pad[[e]] generates a tuple of ’s of length equal to the number of function applications in e, where

is a dummy constant that just occupies a spot. The length of the tuple generated by Pad[[e]] can easily be determined by static inspection of e. The use of Pad ensures that each possible intermediate result appears in a fixed position independent of values of the conditions in conditional expressions.

Ext[[v]]  = < v >
Ext[[g(e_{1}, ..., e_{n})]] where g is c or p

= let v_{1} =Ext[[e_{1}]] in ... let v_{n} =Ext[[e_{n}]] in

< g(1st(v_{1}), ..., 1st(v_{n})) > @ rst(v_{1}) @ ... @ rst(v_{n})

Ext[[f(e_{1}, ..., e_{n})]] = let v_{1} =Ext[[e_{1}]] in ... let v_{n} =Ext[[e_{n}]] in

let v =f̄(1st(v_{1}), ..., 1st(v_{n})) in

< 1st(v) > @ rst(v_{1}) @ ... @ rst(v_{n}) @ < v >

Ext[[if e_{1} then e_{2} else e_{3}]]

= let v_{1} =Ext[[e_{1}]] in

if 1st(v_{1}) \begin{matrix}then & \begin{matrix}let v_{2} =Ext[[e_{2}]] in \\ < 1st(v_{2}) > @ rst(v_{1}) @ rst(v_{2}) @ Pad[[e_{3}]] \end{matrix} \\ else & \begin{matrix}let v_{3} =Ext[[e_{3}]] in \\ < 1st(v_{3}) > @ rst(v_{1}) @ Pad[[e_{2}]] @ rst(v_{3}) \end{matrix} \end{matrix}

Ext[[let v =e_{1} in e_{2}]]

= let v_{1} =Ext[[e_{1}]] in

let v =1st(v_{1}) in \begin{matrix}let v_{2} =Ext[[e_{2}]] in \\ < 1st(v_{2}) > @ rst(v_{1}) @ rst(v_{2}) \end{matrix}

Fig. 3. Definition of Ext.

It is easy to show by induction that each rule Ext[[e]] = e′ guarantees that 1st(e′) = e. Thus, 1st(f̄(v_{1}, ..., v_{n})) = f(v_{1}, ..., v_{n}). Essentially, f̄ performs the same computation as f except that f̄ builds a tree of intermediate results computed in computing f. The original value returned by any subcomputation is always a leftmost component. Each leftmost component is consumed as in the original computation to produce the original value. Values of function calls are also kept in the other components, as first established by the last < v > in the rule for function applications and propagated by @rst(v_{i}) in the other rules.

Administrative simplifications are performed on the resulting functions to simplify tuple operations for passing intermediate results, unwind binding expressions that become unnecessary as a result of simplifying their subexpressions, and lift bindings out of enclosing expressions whenever possible to enhance readability. For example, for functions sum and prod in Figure 1, we obtain functions \overline{sum} and \overline{prod}, respectively, in (7), where intermediate values of recursive calls to sum and prod, respectively, are returned as well as the original return values.

The following improvements, not given in [55] or [54], can be made to the above brute-force caching of all intermediate results. First, before applying the extension transformation, common subcomputations in both branches of a conditional expression are lifted out of the conditional. This simplifies programs in general. For caching all intermediate results, this lifting saves the extension transformation from caching values of common subcomputations at different positions in different branches, which makes it easier to reason about using these values for incremental computation. The same effect can be achieved by explicitly allocating, for values of common subcomputations in different branches, the same slot in each corresponding branch.

Next, we concentrate on major improvements. These improvements are based on an embedding analysis.

Embedding analysis. First, we define and compute embedding predicates Mf and Me, where “is embedded in” means “can be retrieved from”. We use Mf(f, i) to indicate whether the value of v_{i} is embedded in the value of f(v_{1}, ..., v_{n}), and we use Me(e, v) to indicate whether the value of variable v is embedded in the value of expression e. These predicates must satisfy the following safety requirements:

\begin{matrix}if Mf(f, i) = true, then there exists a function f_{i}^{−1} \\ such that, if u = f(v_{1}, ..., v_{n}), then v_{i} = f_{i}^{−1}(u) \\ if Me(e, v) = true, then there exists a function e_{v}^{−1} \\ such that, if u = e, then v = e_{v}^{−1}(u) \end{matrix} (3)

For each function definition f(v_{1}, ..., v_{n}) ≜ e_{f}, we define Mf(f, i) = Me(e_{f}, v_{i}), and we define Me recursively as in Figure 4. For a primitive function p, ∃p_{i}^{−1} denotes true if p has an inverse for the ith argument, and false otherwise. For a conditional expression, if_{e_{2}e_{3}}^{e_{1}} denotes true if the value of e_{1} can be determined statically or inferred from the value of if e_{1} then e_{2} else e_{3}, and false otherwise. For example, if_{e_{2}e_{3}}^{e_{1}} is true if e_{1} is T (for true) or F (for false), or if the two branches of the conditional expression return applications of different constructors. For a Boolean expression e_{1}, e_{1} ⊢Me(e, v) means that whenever the value of e_{1} is T, the value of v is embedded in the value of e. In order that the embedding analysis does not obviate useful caching, it considers a value to be embedded only if the value can be retrieved from the value of its immediately enclosing computation in constant time; in particular, this constraint applies to the retrievals when ∃p_{i}^{−1} or if_{e_{2}e_{3}}^{e_{1}} is true.

Me(u, v)  = \{ \begin{matrix}true & if v = u \\ false & otherwise \end{matrix}
Me(c(e_{1}, ..., e_{n}), v)  = Me(e_{1}, v) ∨ ... ∨ Me(e_{n}, v)
Me(p(e_{1}, ..., e_{n}), v)  = (∃p_{1}^{−1} ∧ Me(e_{1}, v)) ∨ ... ∨ (∃p_{n}^{−1} ∧ Me(e_{n}, v))
Me(f(e_{1}, ..., e_{n}), v)  = (Mf(f, 1) ∧ Me(e_{1}, v)) ∨ ... ∨ (Mf(f, n) ∧ Me(e_{n}, v))
Me(if e_{1} then e_{2} else e_{3}, v) = if_{e_{2}e_{3}}^{e_{1}} ∧ (e_{1} ⊢Me(e_{2}, v)) ∧ (¬e_{1} ⊢Me(e_{3}, v))
Me(let u = e_{1} in e_{2}, v)  = Me(e_{2}, v) ∨ (Me(e_{1}, v) ∧ Me(e_{2}, u))

Fig. 4. Definition of Me.

We can easily show by induction that the safety requirements (3) are satisfied. To compute Mf, we start with Mf(f, i) = true for every f and i and iterate using the above definitions to compute the greatest fixed point in the pointwise extension of the Boolean domain with false ⊑ true. The iteration always terminates since these definitions are monotonic and the domain is finite.

Next, we compute embedding tags. For each function definition f(v_{1}, ..., v_{n}) ≜ e_{f}, we associate an embedding tag Mtag(e) with each occurrence of each subexpression e of e_{f}, indicating whether the value of e is embedded in the value of e_{f}. Mtag can be defined in a similar fashion to Me. We define Mtag(e_{f}) = true, and define the true values of Mtag for subexpressions e of e_{f} as in Figure 5; the tags of other subexpressions of e_{f} are defined to be false. These tags can be computed directly once the above embedding

if Mtag(c(e_{1}, ..., e_{n})) = true  then Mtag(e_{i}) = true, for i = 1..n
if Mtag(p(e_{1}, ..., e_{n})) = true  then Mtag(e_{i}) = true if ∃p_{i}^{−1}, for i = 1..n
if Mtag(f(e_{1}, ..., e_{n})) = true  then Mtag(e_{i}) = true if Mf(f, i), for i = 1..n
if Mtag(if e_{1} then e_{2} else e_{3}) = true then Mtag(e_{i}) = true if if_{e_{2}e_{3}}^{e_{1}}, for i = 1, 2, 3
if Mtag(let v =e_{1} in e_{2}) = true  then Mtag(e_{2}) = true; Mtag(e_{1}) = true if Me(e_{2}, v)

Fig. 5. Definition of Mtag.

predicates are computed.

Finally, we use the embedding tags to compute, for each function f, an embedding-all property Mall(f) indicating whether all intermediate results of f are embedded in the value of f. We define, for each function f(v_{1}, ..., v_{n}) ≜ e_{f},

Mall(f) = ⋀ ; all function applications g(e_{1}, ..., e_{n}) occurring in e_{f} Mtag(g(e_{1}, ..., e_{m})) ∧ Mall(g) (4) where Mtag is with respect to e_{f}. To compute Mall, we start with Mall(f) = true for all f and iterate using the definition in (4) until the greatest fixed point is reached. This fixed point exists for similar reasons as for Mf.

Embedding analysis is related to the transmission analysis for compile-time garbage collection [39], in that both analyze values that appear in the result of the enclosing computation. While embedding analysis determines whether a value appears in the result, transmission analysis determines what components of a value appear in certain components of the result. The domain formed with true and false for embedding analysis is finite, so it allows an efficient linear-time analysis, which is also precise for our usage here. The components in transmission analysis are described as nested patterns that form an infinite domain, causing the analysis to be imprecise or inefficient in general. Embedding analysis is also related to the propagation analysis for update optimization [81], but that analysis only determines whether a value is the same as the result of the enclosing computation.

Improved extension transformation. The above embedding analysis is used to improve the extension transformation as follows.

First, if Mall(f) = true, i.e., if all intermediate results of f are embedded in the value of f, then we do not construct an extended function for f. This makes the transformation for caching all intermediate results idempotent.

If there is a function not all of whose intermediate results are embedded in its return value, then an extended function for it needs to be defined as in (2). We modify the definition of Ext[[f(e_{1}, ..., e_{n})]] as follows. If Mall(f) = true, which includes the case where f does not contain function applications, then, due to the first improvement, f is not extended, so we reference the value of f directly:

Ext[[f(e_{1}, ..., e_{n})]] = let v_{1} =Ext[[e_{1}]] in ... let v_{n} =Ext[[e_{n}]] in let v =f(1st(v_{1}), ..., 1st(v_{n})) in < v > @ rst(v_{1}) @ ... @ rst(v_{n}) @ < v > (5)

Furthermore, if Mall(f) = true, and Mtag(f(e_{1}, ..., e_{n})) = true, i.e, the value of f(e_{1}, ..., e_{n}) is embedded in the value of its enclosing application, then we avoid caching the value of f separately:

Ext[[f(e_{1}, ..., e_{n})]] = let v_{1} =Ext[[e_{1}]] in ... let v_{n} =Ext[[e_{n}]] in <f(1st(v_{1}), ..., 1st(v_{n}))> @ rst(v_{1}) @ ... @ rst(v_{n})(6)

To summarize, the transformation Ext remains the same as in Figure 3 except that the rule for a function application f(e_{1}, ..., e_{n}) is replaced with the following: if Mall(f) = true and Mtag(f(e_{1}, ..., e_{n})) = true, then define Ext[[f(e_{1}, ..., e_{n})]] as in (6); else if Mall(f) = true but Mtag(f(e_{1}, ..., e_{n})) = false, then define Ext[[f(e_{1}, ..., e_{n})]] as in (5); otherwise define Ext[[f(e_{1}, ..., e_{n})]] as in Figure 3. Note that function applications f(e_{1}, ..., e_{n}) such that Mall(f) = true and Mtag(f(e_{1}, ..., e_{n})) = true should not be counted by Pad. The lengths of tuples generated by Pad can still be statically determined.

For function cmp in Figure 1, this improved extension transformation yields the following functions:

\begin{matrix}\overline{cmp}(x) ≜ let v_{1} = odd(x) in let u_{1} = \overline{sum}(v_{1}) in let v_{2} = even(x) in let u_{2} = \overline{prod}(v_{2}) in < 1st(u_{1}) ≤ 1st(u_{2}), v_{1}, u_{1}, v_{2}, u_{2} > \\ \overline{sum}(x) ≜ if null(x) then < 0, > else \begin{matrix}let v_{1} = \overline{sum}(cdr(x)) in \\ < car(x) + 1st(v_{1}), v_{1} > \end{matrix} \\ \overline{prod}(x) ≜ if null(x) then < 1, > else \begin{matrix}let v_{1} = \overline{prod}(cdr(x)) in \\ < car(x) ∗ 1st(v_{1}), v_{1} > \end{matrix} \end{matrix} (7)

where \overline{cmp} is extended to return values of odd and even, as well as recursive calls to sum and prod. Functions odd and even are not extended, since all their intermediate results are embedded in their return values.

3.2 Step A.2: Exposing auxiliary information by incrementalization

This step transforms f̄_{0}(x ⊕ y) to expose subcomputations depending only on x and whose values cannot be retrieved from the cached result of f̄_{0}(x). It uses analyses and transformations similar to those in [56] that derive an incremental program f̄_{0}′(x, y, r̄), by expanding subcomputations of f̄_{0}(x ⊕ y) depending on both x and y and replacing those depending only on x by retrievals from r̄ when possible.

Our goal here is not to quickly retrieve values from r̄, but to find potentially useful auxiliary information, i.e., subcomputations depending on x (and r̄) but not y whose values cannot be retrieved from r̄. Thus, time considerations in [56] are dropped here but are picked up after Step A.3, as discussed in Section 5.

In particular, in [56], a (recursive) application of a function f is replaced by an application of an incremental version f′ only if a fast retrieval from some cached result of the previous computation can be used as the argument for the parameter of f′ that corresponds to a cached result. For example, if an incremental version f′(x, y, r) is introduced to compute f(x⊕y) incrementally for r = f(x), then in [56], a function application f(g(x) ⊕ h(y)) is replaced by an application of f′ only if some fast retrieval p(r) for the value of f(g(x)) can be used as the argument for the parameter r of f′(x, y, r), in which case the application is replaced by f′(g(x), h(y), p(r)). In Step A.2 here, an application of f is replaced by an application of f′ also when a retrieval cannot be found; in this case, the value needed for the cache parameter is computed directly, so for this example, the application f(g(x) ⊕ h(y)) is replaced by f′(g(x), h(y), f(g(x))). It is easy to see that, in this case, f(g(x)) becomes a piece of candidate auxiliary information.

Since the functions obtained from this step may be different from the incremental functions f′ obtained in [56], we denote them by f^{‵}.

For function \overline{cmp} in (7) and input change operation x ⊕ y = cons(y, x), we transform the computation of \overline{cmp}(cons(y, x)), with \overline{cmp}(x) = r̄:

1. unfold \overline{cmp}(cons(y, x))

= let v_{1} = odd(cons(y, x)) in let u_{1} = \overline{sum}(v_{1}) in let v_{2} = even(cons(y, x)) in let u_{2} = \overline{prod}(v_{2}) in < 1st(u_{1}) ≤ 1st(u_{2}), v_{1}, u_{1}, v_{2}, u_{2} >

2. unfold odd, \overline{sum}, even and simplify

= let v′_{1} = even(x) in let u′_{1} = \overline{sum}(v′_{1}) in let v_{2} = odd(x) in let u_{2} = \overline{prod}(v_{2}) in < y+1st(u′_{1}) ≤ 1st(u_{2}), cons(y, v′_{1}), <y+1st(u′_{1}), u′_{1} >, v_{2}, u_{2} >

3. replace applications of even and odd by retrievals

= let v′_{1} = 4th(r̄) in let u′_{1} = \overline{sum}(v′_{1}) in let v_{2} = 2nd(r̄) in let u_{2} = \overline{prod}(v_{2}) in < y+1st(u′_{1}) ≤ 1st(u_{2}), cons(y, v′_{1}), <y+1st(u′_{1}), u′_{1} >, v_{2}, u_{2} >

Simplification, i.e., unwinding the bindings for v′_{1} and v_{2} here, yields the following function \overline{cmp}^{‵} such that, if \overline{cmp}(x) = r̄, then \overline{cmp}^{‵}(y, r̄) = \overline{cmp}(cons(y, x)):

\overline{cmp}^{‵}(y, r̄) ≜ let u′_{1} = \overline{sum}(4th(r̄)) in let u_{2} = \overline{prod}(2nd(r̄)) in <y+1st(u′_{1})≤1st(u_{2}), cons(y, 4th(r̄)), < y+1st(u′_{1}), u′_{1} >, 2nd(r̄), u_{2} > (8)

3.3 Step A.3: Collecting candidate auxiliary information

This step collects candidate auxiliary information, i.e., intermediate results of f̄_{0}^{‵}(x, y, r̄) that depend only on x and r̄, and yields function f̌_{0}(x, r̄). It is similar to Step A.1 in that both collect intermediate results; they differ in that Step A.1 collects all intermediate results, while this step collects only those that depend only on x and r̄.

Forward dependence analysis. First, we use a forward dependence analysis to identify subcomputations of f̄_{0}^{‵}(x, y, r̄) that depend only on x and r̄. The analysis is in the same spirit as binding-time analysis [38,45] for partial evaluation, if we regard the arguments corresponding to x and r̄ as static and the rest as dynamic. We compute the following sets, called forward dependency sets, directly.

For each function f(v_{1}, ..., v_{n}) ≜ e_{f}, we compute a set Sf(f) that contains the indices of the arguments of f such that, in all uses of f, the values of these arguments depend only on x and r̄, and, for each occurrence of each subexpression e of e_{f}, we compute a set Se(e) that contains the free variables in e that depend only on x and r̄. The recursive definitions of these sets are given in Figure 6, where F V (e) denotes the set of free variables in e and is

For each function f(v_{1}, ..., v_{n}) ≜ e_{f}, define Se(e_{f}) = {v_{i} | i ∈Sf(f)} and, for each subexpression e of e_{f},
if e is c(e_{1}, ..., e_{n}) or p(e_{1}, ..., e_{n}) then Se(e_{1}) = ... = Se(e_{n}) = Se(e)
if e is f_{1}(e_{1}, ..., e_{n})  then Se(e_{1}) = ... = Se(e_{n}) = Se(e), Sf(f_{1})={i | F V (e_{i}) ⊆Se(e)} ∩ Sf(f_{1})
if e is if e_{1} then e_{2} else e_{3}  then Se(e_{1}) = Se(e_{2}) = Se(e_{3}) = Se(e)
if e is let v = e_{1} in e_{2}  then Se(e_{1}) = Se(e), Se(e_{2}) = \begin{cases}Se(e) ∪ \{v\} & if F V (e_{1}) ⊆ Se(e) \\ Se(e) \ \{v\} & otherwise \end{cases}

Fig. 6. Definition of Sf and Se.

defined as follows:

F V (v)  = {v}
F V (g(e_{1}, ..., e_{n})), where g is c, p, or f = F V (e_{1}) ∪ ... ∪ F V (e_{n})
F V (if e_{1} then e_{2} else e_{3})  = F V (e_{1}) ∪ F V (e_{2}) ∪ F V (e_{3})
F V (let v = e_{1} in e_{2})  = F V (e_{1}) ∪ (F V (e_{2}) \ {v})

To compute these sets, we start with Sf(f̄_{0}^{‵}) containing the indices of the arguments of f̄_{0}^{‵} corresponding to x and r̄, and, for all other functions f, Sf(f) containing the indices of all arguments of f, and iterate until a greatest fixed point with respect to the point-wise extension of subset ordering is reached. This iteration always terminates since, for each function f, f has a fixed arity, Sf(f) decreases, and a lower bound ∅ exists.

For the running example, we obtain Sf(\overline{cmp}^{‵}) = \{2\} and Sf(\overline{sum}) = Sf(\overline{prod}) = \{1\}. For every subexpression e in the definition of \overline{cmp}^{‵}(y, r̄), r̄ ∈ Se(e). For every subexpression e in the definitions of \overline{sum}(x) and \overline{prod}(x), Se(e) = \{x\}.

Collection transformation. Next, we use a collection transformation to collect the candidate auxiliary information. The main difference between this collection transformation and the extension transformation in Step A.1 is that, in the former, the value originally computed by a subexpression is returned only if it depends only on x and r̄, while in the latter, the value originally computed by a subexpression is always returned.

For each function f(v_{1}, ..., v_{n}) ≜ e called in the program for f̄_{0}^{‵} and such that Sf(f) ≠ ∅, we construct a function definition

f̌(v_{i_{1}}, ..., v_{i_{k}}) ≜ Col[[e]] (9)

where Sf(f) = {i_{1}, ..., i_{k}} and 1 ≤ i_{1} < ... < i_{k} ≤ n. Col[[e]] collects the results of intermediate function applications in e that have been statically determined to depend only on x and r̄. Note, however, that an improvement similar to that in Step A.1 is made, namely, we avoid constructing such a collected version for f if Sf(f) = {1, ..., n} and Mall(f) = true.

The transformation Col always first examines whether its argument expression e has been determined to depend only on x and r̄, i.e., F V (e) ⊆ Se(e). If so, Col[[e]] = Ext[[e]], where Ext is the improved extension transformation defined in Step A.1. Otherwise, Col[[e]] is defined as in Figure 7, where Pǎd[[e]] generates a tuple of ’s of length equal to the number of the function applications in e, except that function applications f(e_{1}, ..., e_{n}) such that Sf(f) = ∅, or Sf(f) = {1, ..., n} but Mall(f) = true and Mtag(f(e_{1}, ..., e_{n})) = true are not counted. Note that if e has been determined to depend only on x and r̄, then 1st(Col[[e]]) is the original value of e, and rst(Col[[e]]) contains values of intermediate function applications that depend only on x and r̄; otherwise, Col[[e]] contains only values of intermediate function applications that depend only on x and r̄.

Col[[v]]  = <>
Col[[g(e_{1}, ..., e_{n})]] where g is c or p = Col[[e_{1}]] @ ... @ Col[[e_{n}]]
Col[[f(e_{1}, ..., e_{n})]] = let v_{1} =Col[[e_{1}]] in ... let v_{n} =Col[[e_{n}]] in e′_{1} @ ... @e′_{n} @ e′

where e′_{i} = \begin{cases}rst(v_{i}) & if i ∈ Sf(f) \\ v_{i} & otherwise \end{cases}

e′ = \begin{cases}<> & if Sf(f)= ∅ \\ <f̌(1st(v_{i_{1}}), ...,1st(v_{i_{k}}))> & otherwise \end{cases}

where Sf(f)= {i_{1}, ..., i_{k}} and 1≤ i_{1} <...<i_{k} ≤ n

Col[[if e_{1} then e_{2} else e_{3}]]

= let v_{1} =Col[[e_{1}]] in

if 1st(v_{1}) \begin{matrix}then & \begin{matrix}let v_{2} =Col[[e_{2}]] in \\ rst(v_{1}) @ v_{2} @ Pǎd[[e_{3}]] \end{matrix} \\ else & \begin{matrix}let v_{3} =Col[[e_{3}]] in \\ rst(v_{1}) @ Pǎd[[e_{2}]] @ v_{3} \end{matrix} \end{matrix}

if F V (e_{1}) ⊆ Se(e_{1})

= let v_{1} =Col[[e_{1}]] in let v_{2} =Col[[e_{2}]] in let v_{3} =Col[[e_{3}]] in v_{1} @ v_{2} @ v_{3}  otherwise

Col[[let v =e_{1} in e_{2}]]

= let v_{1} =Col[[e_{1}]] in

let v =1st(v_{1}) in \begin{matrix}let v_{2} =Col[[e_{2}]] in \\ rst(v_{1}) @ v_{2} \end{matrix}

if F V (e_{1}) ⊆ Se(e_{1})

= let v_{1} =Col[[e_{1}]] in \begin{matrix}let v_{2} =Col[[e_{2}]] in \\ v_{1} @ v_{2} \end{matrix}  otherwise

Fig. 7. Definition of Col when F V (e) ⊈ Se(e).

Although this forward dependence analysis is equivalent to binding time analysis in partial evaluation [37], the application here is different. In partial evaluation, the goal is to obtain a residual program that is specialized on a given set of static arguments and takes only the dynamic arguments, while here, we construct a program that computes only on the “static” arguments. In a sense, what is computed by our resulting program corresponds to what is computed by the first of the two stages from the staging transformation [40]. The resulting program obtained here is similar to the slice obtained from forward slicing [84]. However, our forward dependence analysis finds parts of a program that depend only on certain information, while forward slicing finds parts of a program that depend possibly on certain information. Furthermore, our resulting program also returns all intermediate results on the arguments of interest.

For function \overline{cmp}^{‵} in (8), collecting all intermediate results that depend only on its second parameter yields

čmp(r̄) ≜ < \overline{sum}(4th(r̄)), \overline{prod}(2nd(r̄)) > (10)

We can see that computing čmp(r̄) is no slower than computing cmp(x). We will see that this guarantees that incremental computation using the program obtained at the end is at least as fast as computing cmp from scratch.

4 Phase B: Using auxiliary information

Phase B determines which pieces of the collected candidate auxiliary information are useful for incremental computation of f_{0}(x ⊕ y) and exactly how they can be used. The basic idea is to merge the candidate auxiliary information with the original computation of f_{0}(x), derive an incremental version for the merged program, and determine the least information useful for computing the value of f_{0}(x ⊕ y) in that incremental version.

However, we want the incremental computation of f_{0}(x ⊕ y) to have access to the auxiliary information in addition to the intermediate results of f_{0}(x). Thus, we merge the candidate auxiliary information in f̌_{0}(x, r̄) with f̄_{0}(x) instead of f_{0}(x). After deriving an incremental version for the resulting program, we prune out the useless auxiliary information and the useless intermediate results.

Phase B has three steps. Step 1 merges f̌_{0} with f̄_{0} to form a function ¯ ; ˇf_{0} that returns candidate auxiliary information as well as all intermediate results. It also determines a projection Π_{0} that projects the return value of f_{0} out of ¯ ; ˇf_{0}. Step 2 incrementalizes ¯ ; ˇf_{0} under ⊕ to obtain an incremental version ¯ ; ˇf_{0}′. Step 3 prunes out of ¯ ; ˇf_{0} and ¯ ; ˇf_{0}′ the intermediate results and auxiliary information that are not useful.

4.1 Step B.1: Combining intermediate results and auxiliary information

To merge the candidate auxiliary information with f̄_{0}, we could simply attach it onto f̄_{0} by defining ¯ ; ˇf_{0} to be the pair of f̄_{0} and f̌_{0}:

¯ ; ˇf_{0}(x) ≜ let r̄ = f̄_{0}(x) in let ř = f̌_{0}(x, r̄) in < r̄, ř >

and use the projection Π_{0}(¯ ; ř) = 1st(1st(¯ ; ř)) to project out the original return value of f_{0}. However, we can do better by using a transformation to integrate the computation of f̌_{0} more tightly into the computation of f̄_{0}, as opposed to carrying out two disjoint computations. The integrated computation is usually more efficient; so is its incremental version.

We do not describe the integration in detail. Basically, it uses traditional transformation techniques [14] like those used in tupling tactic [24,66,15]. We require only that Π_{0}(¯ ; ˇf_{0}(x)) always project out 1st(f̄_{0}(x)), which is the value of f_{0}(x), and that the values of all other components of f̄_{0}(x) and f̌_{0}(x, r̄) are embedded in the value of ¯ ; ˇf_{0}(x). This allows re-arranging the components in the return value.

For functions \overline{cmp} in (7) and čmp in (10), we first define a function

\begin{matrix}¯ & ¯ \\ ˇ cmp \end{matrix}(x) ≜ let r̄ = \overline{cmp}(x) in let ř = čmp(r̄) in < r̄, ř >

and a projection Π_{0}(¯ ; ř) = 1st(1st(¯ ; ř)). Next, we transform \begin{matrix}¯ & ¯ \\ ˇ cmp \end{matrix}(x) to integrate the computations of \overline{cmp} and čmp,

1. unfold \begin{matrix}¯ & ¯ \\ ˇ cmp \end{matrix}, then \overline{cmp} and čmp

= let r̄ = let v_{1} = odd(x) in let u_{1} = \overline{sum}(v_{1}) in let v_{2} = even(x) in let u_{2} = \overline{prod}(v_{2}) in < 1st(u_{1}) ≤ 1st(u_{2}), v_{1}, u_{1}, v_{2}, u_{2} > in let ř = < \overline{sum}(4th(r̄)), \overline{prod}(2nd(r̄)) > in < r̄, ř >

2. lift bindings for v_{1}, u_{1}, v_{2}, u_{2}, and simplify

= let v_{1} = odd(x) in let u_{1} = \overline{sum}(v_{1}) in let v_{2} = even(x) in let u_{2} = \overline{prod}(v_{2}) in let r̄ = <1st(u_{1}) ≤ 1st(u_{2}), v_{1}, u_{1}, v_{2}, u_{2} > in let ř = <\overline{sum}(v_{2}), \overline{prod}(v_{1}) > in < r̄, ř >

3. unfold bindings for r̄ and ř

= let v_{1} = odd(x) in let u_{1} = \overline{sum}(v_{1}) in let v_{2} = even(x) in let u_{2} = \overline{prod}(v_{2}) in < <1st(u_{1}) ≤ 1st(u_{2}), v_{1}, u_{1}, v_{2}, u_{2} >, < \overline{sum}(v_{2}), \overline{prod}(v_{1}) >>

Simplifying the return value and Π_{0}, we obtain the function

\begin{matrix}¯ & ¯ \\ ˇ cmp \end{matrix}(x) ≜ let v_{1} = odd(x) in let u_{1} = \overline{sum}(v_{1}) in let v_{2} = even(x) in let u_{2} = \overline{prod}(v_{2}) in < 1st(u_{1}) ≤ 1st(u_{2}), v_{1}, u_{1}, v_{2}, u_{2}, \overline{sum}(v_{2}), \overline{prod}(v_{1}) > (11)

and the projection Π_{0}(¯ ; ř) = 1st(¯ ; ř).

4.2 Step B.2: Incrementalization

To derive an incremental version ¯ ; ˇf_{0}′ of ¯ ; ˇf_{0} under ⊕, we can use the method in [56], as sketched in Section 1. Depending on the power expected from the derivation, the method can be made semi-automatic or fully automatic.

For function \begin{matrix}¯ & ¯ \\ ˇ cmp \end{matrix} in (11) and input change operation x ⊕ y = cons(y, x), we derive an incremental version of \begin{matrix}¯ & ¯ \\ ˇ cmp \end{matrix} under ⊕:

1. unfold \begin{matrix}¯ & ¯ \\ ˇ cmp \end{matrix}(cons(y, x))

= let v_{1} = odd(cons(y, x)) in let u_{1} = \overline{sum}(v_{1}) in let v_{2} = even(cons(y, x)) in let u_{2} = \overline{prod}(v_{2}) in < 1st(u_{1}) ≤ 1st(u_{2}), v_{1}, u_{1}, v_{2}, u_{2}, \overline{sum}(v_{2}), \overline{prod}(v_{1}) >

2. unfold odd, \overline{sum}, even, \overline{prod} and simplify

= let v′_{1} = even(x) in let u′_{1} = \overline{sum}(v′_{1}) in let v_{2} = odd(x) in let u_{2} = \overline{prod}(v_{2}) in let u′_{4} = \overline{prod}(v′_{1}) in <\begin{matrix}y + 1st(u′_{1}) ≤ 1st(u_{2}), cons(y, v′_{1}), <y+1st(u′_{1}), u′_{1} >, v_{2}, u_{2}, \\ \overline{sum}(v_{2}), <y∗1st(u′_{4}), u′_{4} >> \end{matrix}

3. replace all applications by retrievals

= let v′_{1} = 4th(\begin{matrix} & ¯ \\ ř \end{matrix}) in let u′_{1} = 6th(\begin{matrix} & ¯ \\ ř \end{matrix}) in let v_{2} = 2nd(\begin{matrix} & ¯ \\ ř \end{matrix}) in let u_{2} = 7th(\begin{matrix} & ¯ \\ ř \end{matrix}) in let u′_{4} = 5th(\begin{matrix} & ¯ \\ ř \end{matrix}) in <\begin{matrix}y + 1st(u′_{1}) ≤ 1st(u_{2}), cons(y, v′_{1}), <y+1st(u′_{1}), u′_{1} >, v_{2}, u_{2}, \\ 3rd(\begin{matrix} & ¯ \\ ř \end{matrix}), <y∗1st(u′_{4}), u′_{4} >> \end{matrix}

Simplification, i.e., unwinding all the bindings, yields the following incremental version \begin{matrix}¯ & ¯ \\ ˇ cmp \end{matrix}′ such that, if \begin{matrix}¯ & ¯ \\ ˇ cmp \end{matrix}(x) =¯ ; ř, then \begin{matrix}¯ & ¯ \\ ˇ cmp \end{matrix}′(y,¯ ; ř) = \begin{matrix}¯ & ¯ \\ ˇ cmp \end{matrix}(cons(y, x)):

\begin{matrix}¯ & ¯ \\ ˇ cmp \end{matrix}′(y,\begin{matrix} & ¯ \\ ř \end{matrix}) ≜ <\begin{matrix}y+1st(6th(\begin{matrix} & ¯ \\ ř \end{matrix})) ≤ 1st(7th(\begin{matrix} & ¯ \\ ř \end{matrix})), \\ cons(y, 4th(\begin{matrix} & ¯ \\ ř \end{matrix})), <y+1st(6th(\begin{matrix} & ¯ \\ ř \end{matrix})), 6th(\begin{matrix} & ¯ \\ ř \end{matrix})>, 2nd(\begin{matrix} & ¯ \\ ř \end{matrix}), 7th(\begin{matrix} & ¯ \\ ř \end{matrix}), \\ 3rd(\begin{matrix} & ¯ \\ ř \end{matrix}), <y∗1st(5th(\begin{matrix} & ¯ \\ ř \end{matrix})), 5th(\begin{matrix} & ¯ \\ ř \end{matrix})>> \end{matrix} (12)

Clearly, \begin{matrix}¯ & ¯ \\ ˇ cmp \end{matrix}′(y,¯ ; ř) computes \begin{matrix}¯ & ¯ \\ ˇ cmp \end{matrix}(cons(y, x)) in only O(1) time.

4.3 Step B.3: Pruning

To prune ¯ ; ˇf_{0} and ¯ ; ˇf_{0}′, we use the analyses and transformations in Stage III of [54]. A backward dependence analysis determines the components of ¯ ; řand subcomputations of ¯ ; ˇf_{0}′ whose values are useful in computing Π_{0}(¯ ; ˇf_{0}′(x, y,¯ ; ř)), which is the value of f_{0}. A pruning transformation replaces useless computations with . Finally, the resulting functions are optimized by eliminating the

components, adjusting the selectors, etc. Improved methods for pruning are described in [48,52].

For functions \begin{matrix}¯ & ¯ \\ ˇ cmp \end{matrix} in (11) and \begin{matrix}¯ & ¯ \\ ˇ cmp \end{matrix}′ in (12), 1st(\begin{matrix}¯ & ¯ \\ ˇ cmp \end{matrix}′(y,¯ ; ř)) depends on 1st(6th(¯ ; ř)) and 1st(7th(¯ ; ř)), which depend on 1st(3rd(¯ ; ř)) and 1st(5th(¯ ; ř)), respectively. All other components are not needed. We obtain

_{cmp}(x) ≜ let v_{1} = odd(x) in let u_{1} = \overline{sum}(v_{1}) in let v_{2} = even(x) in let u_{2} = \overline{prod}(v_{2}) in < \begin{matrix}1st(u_{1}) ≤ 1st(u_{2}), , < 1st(u_{1}), >, , < 1st(u_{2}), >, \\ < 1st(\overline{sum}(v_{2})), >, < 1st(\overline{prod}(v_{1})), >> \end{matrix}

_{cmp}′(y, _{r}) ≜ <\begin{matrix}y+1st(6th(_{r})) ≤ 1st(7th(_{r})), , < y+1st(6th(_{r})), >, , < 1st(7th(_{r})), >, \\ < 1st(3rd(_{r})), >, < y∗1st(5th(_{r})), >> \end{matrix}

Optimizing these functions, i.e., eliminating unneeded components, adjusting indexing, and simplifying tuple constructions and selections, yields the final definitions of c̃mp and c̃mp′, which appear in Figure 2.

5 Discussion

Correctness. Auxiliary information is maintained incrementally, so at the step of discovering it, we should not be concerned with the time complexity of computing it from scratch; this is why time considerations were dropped in Step A.2. However, to make the overall approach effective, we must consider the cost of computing and maintaining the auxiliary information. Here, we simply require that the candidate auxiliary information be computed at least as fast as the original program, i.e., t(f̌_{0}(x, r̄)) ≤ t(f_{0}(x)) for r̄ = f̄_{0}(x). This can be checked after Step A.3. We guarantee this condition by simply dropping pieces of candidate auxiliary information for which it cannot be confirmed. Standard constructions for mechanical time analysis [75,83,49] can be used, although further study is needed, and it is being carried out for both time analysis [29] and space analysis [80]. The trade-off between time and space is a problem open for study.

Suppose Step B.1 projects out the original value using 1st. With the above condition, in a similar way to [54], we can show that, if f_{0}(x) = r, then

1st(f̃_{0}(x)) = r and t(f̃_{0}(x)) ≤ t(f_{0}(x)) (13)

and if f̃_{0}(x) = r̃ and f_{0}(x ⊕ y) = r′, then

\begin{matrix}1st(f̃_{0}′(x, y, r̃)) = r′, f̃_{0}′(x, y, r̃) = f̃_{0}(x ⊕ y), \\ and t(f̃_{0}′(x, y, r̃)) ≤ t(f_{0}(x ⊕ y)). \end{matrix} (14)

i.e., the programs f̃_{0} and f̃_{0}′ preserve the semantics and compute asymptotically at least as fast as f_{0}. Note that f̃_{0}(x) may terminate more often than f_{0}(x), and f̃_{0}′(x, y, r̃) may terminate more often than f_{0}(x ⊕ y), due to the transformations used in Steps B.2 and B.3.

Multi-pass discovery of auxiliary information. The program f̃_{0} can sometimes be computed even faster by maintaining auxiliary information useful for incremental computation of the auxiliary information already in f̃_{0}. We can obtain such auxiliary information of auxiliary information by iterating the above approach. Whether to continue iteration using the approach depends on whether the desired performance improvement is achieved. Therefore, again, analysis of performance of the transformed program is important.

Other auxiliary information. There are cases where the auxiliary information discovered using the above approach is not sufficient for efficient incremental computation. For example, we do not yet know how to discover a heap data structure systematically for computing a minimum element in a set under element deletions. In these cases, classes of special parameterized data structures are often used. Ideally, we could collect them as auxiliary information parameterized with certain classes of data types. Then, we could systematically extend a program to compute such auxiliary information and maintain it incrementally. How to do this precisely is a problem open for study. In the worst case, we could code manually discovered auxiliary information to obtain an extended program f̃_{0}, and then use our systematic approach to derive an incremental version f̃_{0}′ that incrementally computes the new output using the auxiliary information and also maintains the auxiliary information.

Incrementalization for program efficiency improvement. To use our approach for program efficiency improvement, we need to identify expensive computations f and appropriate change operations ⊕.

Expensive computations are generally easy to determine, e.g., function calls, especially calls of non-linear functions, in a language that contains functions, operations on sets and maps in a very-high-level language, or multiplications and exponentiations in a hardware description language. Change operations in an iterative program are simply updates to the parameters of f by the loop body, but change operations in a recursive program, especially for nonlinear recursive functions, are non-trivial to determine. Nevertheless, Liu and Stoller have studied a general method for identifying change operations [51] for recursive functions. The basic idea is to use a minimal input change that is in the opposite direction of change compared to arguments of recursive calls. Using the opposite direction of change yields an increment; using a minimal change allows maximum reuse, i.e., maximum incrementality.

Section 6 gives examples where efficient incremental programs are derived using the method in this paper based on identified expensive computations f and change operations ⊕.

6 Examples

The running example on list processing illustrates the application of our approach to solving explicit incremental problems for, e.g., interactive systems and reactive systems. Other applications include optimizing compilers and transformational programming. This section first presents an example for each of these two applications, based on problems in VLSI design and graph algorithms, respectively. Then, we describe four other examples, taken from problems in games, string processing, combinatorial optimization, and image processing. They help show that this systematic method for discovering auxiliary information is powerful, that the underlying principle is general and applies to other language features, and that it has many application areas.

6.1 Strength reduction in optimizing compilers: binary integer square root

This example is from formal hardware design [60], where a specification of a non-restoring binary integer square root algorithm is transformed into a VLSI circuit design and implementation. There, a strength-reduced program was manually discovered and then proved correct using Nuprl [18]. Here, we show how our method can automatically derive the strength reductions. This is of particular interest in light of the Pentium chip flaw in 1994 [28].

The initial specification of the algorithm is given in (15). Given a binary integer n of l bits, where n > 0 (and l is usually 8, 16, ...), it computes the binary integer square root m of n using the non-restoring method [26,60], which is exact for perfect squares and off by at most 1 for other integers. The method is to first set bit l − 1 of the answer m to 1 and then adjust bits l − 2 to 0 in order: increase m by 2^{i} if m is smaller than \sqrt{n}, and decrease m by 2^{i} if m is larger than \sqrt{n}.

read(n, l); m := 2^{l−1}; for i := l − 2 downto 0 do p := n − m^{2}; if p > 0 then m := m + 2^{i}; else if p < 0 then m := m − 2^{i}; write(m); (15)

In hardware, multiplications and exponentials are much more expensive than additions and shifts (doublings or halvings), so the goal is to replace the former by the latter.

To simplify the presentation, we jump to the heart of the problem, namely, computing n − m^{2} and 2^{i} incrementally in each iteration under the change m′ = m ± 2^{i} and i′ = i − 1. Let f_{0} be

f_{0}(n, m, i) ≜ pair(n − m^{2}, 2^{i})

where pair is a constructor with selectors fst(a, b) = a and snd(a, b) = b, and let input change operation ⊕ be

⟨n′, m′, i′⟩ = ⟨n, m, i⟩ ⊕ ⟨ ⟩ = ⟨n, m ± 2^{i}, i − 1⟩

Step A.1. We cache all intermediate results of f_{0}, obtaining

f̄_{0}(n, m, i) ≜ let v = m^{2} in < pair(n − v, 2^{i}), v >

Step A.2. We transform f̄_{0} under ⊕, obtaining

f̄_{0}^{‵}(n, m, i, r̄) ≜ \begin{matrix}let v = 2nd(r̄) ± 2∗m∗snd(1st(r̄)) + (snd(1st(r̄)))^{2} in \\ < pair(n − v, snd(1st(r̄))/2), v > \end{matrix}

Step A.3. We collect candidate auxiliary information, obtaining

f̌_{0}(n, m, i, r̄) ≜ < 2∗m∗snd(1st(r̄)), (snd(1st(r̄)))^{2} > (16)

Step B.1. We merge the collected candidate auxiliary information with f̄_{0}, obtaining Π_{0}(¯ ; ř) = 1st(¯ ; ř) and

\begin{matrix} & ¯ \\ ˇf_{0} \end{matrix}(n, m, i) ≜ let v = m^{2} in let u = 2^{i} in <pair(n − v, u), v, 2∗m∗u, u^{2} >

Step B.2. We derive an incremental version of ¯ ; ˇf_{0} under ⊕, obtaining

\begin{matrix} & ¯ \\ ˇf_{0} \end{matrix}′(n, m, i,\begin{matrix} & ¯ \\ ř \end{matrix}) ≜ let v = 2nd(\begin{matrix} & ¯ \\ ř \end{matrix}) ± 3rd(\begin{matrix} & ¯ \\ ř \end{matrix}) + 4th(\begin{matrix} & ¯ \\ ř \end{matrix}) in let u = snd(1st(\begin{matrix} & ¯ \\ ř \end{matrix}))/2 in < \begin{matrix}pair(fst(1st(\begin{matrix} & ¯ \\ ř \end{matrix})) ∓ 3rd(\begin{matrix} & ¯ \\ ř \end{matrix}) − 4th(\begin{matrix} & ¯ \\ ř \end{matrix}), u), \\ v, 3rd(\begin{matrix} & ¯ \\ ř \end{matrix})/2 ± 4th(\begin{matrix} & ¯ \\ ř \end{matrix}), 4th(\begin{matrix} & ¯ \\ ř \end{matrix})/4 > \end{matrix}

Step B.3. We prune functions ¯ ; ˇf_{0} and ¯ ; ˇf_{0}′, eliminating their second components and obtaining

_{f_{0}}(n, m, i) ≜ let u = 2^{i} in < pair(n − m^{2}, u), 2∗m∗u, u^{2} > (17)

_{f_{0}}′(n, m, i, _{r}) ≜ < \begin{matrix}pair(fst(1st(_{r})) ∓ 2nd(_{r}) − 3rd(_{r}), snd(1st(_{r}))/2), \\ 2nd(_{r})/2 ± 3rd(_{r}), 3rd(_{r})/4 > \end{matrix} (18)

Thus, we can initialize using (17) and perform the corresponding incremental update in the loop body using (18). The expensive multiplications and exponentials in the loop body have been completely replaced with additions, subtractions, and shifts. We even discover that an unnecessary shift is done in [60]. Thus, a systematic approach such as ours is desirable not only for automating designs and guaranteeing correctness, but also for reducing costs. Explicit incrementalization also enables further optimizations on the overall loop structure, as discussed in [47].

6.2 Promotion and accumulation in transformational programming: path sequence problem

This example was used by Bird to illustrate important program transformation strategies called promotion and accumulation [8]. Given a directed acyclic graph, and a string whose elements are vertices in the graph, the problem is to compute the length of the longest subsequence in the string that forms a path in the graph. We focus on the second half of the example, where an exponential-time recursive solution is improved (incorrectly in [8], corrected in [9]).

Function llp defined below computes the desired length. The input string is given explicitly as the argument to llp. The input graph is represented by a predicate arc such that arc(a, b) is true iff there is an edge from vertex a to vertex b in the graph. The primitive function max returns the maximum of its two arguments.

\begin{matrix}llp(l) ≜ if null(l) then 0 else max(llp(cdr(l)), 1+f(car(l), cdr(l))) \\ f(n, l) ≜ if null(l) then 0 else \begin{matrix}if arc(n, car(l)) then \\ max(f(n, cdr(l)), 1+f(car(l), cdr(l))) \end{matrix} else f(n, cdr(l)) \end{matrix} (19)

The problem is to compute llp incrementally under the input change operation l′ = l ⊕ i = cons(i, l). Using the method described in this paper, we obtain

\begin{matrix}_{llp}(l) ≜ if null(l) then <0> else \begin{matrix}let v = _{f}(car(l), cdr(l)) in \\ <max(llp(cdr(l)), 1+1st(v)), v > \end{matrix} \\ _{f}(n, l) ≜ if null(l) then <0> else let u = _{f}(car(l), cdr(l)) in if arc(n, car(l)) then <max(f(n, cdr(l)), 1+1st(u), u> else <f(n, cdr(l)), u> \end{matrix} (20)

and

\begin{matrix}_{llp}′(i, l, _{r}) ≜ if null(l) then <1, <0>> else \begin{matrix}let v = _{f′}(i, l, 2nd(_{r})) in \\ <max(1st(_{r}), 1+1st(v)), v > \end{matrix} \\ _{f′}(i, l, _{r_{1}}) ≜ if null(cdr(l)) then if arc(i, car(l)) then <1, <0>> else <0, <0>> else let u = _{f′}(i, cdr(l), 2nd(_{r_{1}})) in if arc(i, car(l)) then <max(1st(u), 1+1st(_{r_{1}})), _{r_{1}} > else <1st(u), _{r_{1}} > \end{matrix} (21)

Computing llp(cons(i, l)) from scratch takes exponential time, but computing l̃lp′(i, l, r̃) takes only O(n) time, where n is the length of l, since l̃lp′(i, l, r̃) calls f̃′, which goes through the list l once.

Finally, we use these derived functions to compute the original function llp. Note that llp(l) = 1st(l̃lp(l)) and, if l̃lp(l) = r̃, then l̃lp(cons(i, l)) = l̃lp′(i, l, r̃). Using the definition of l̃lp′ in (21) in this last equation, we obtain

_{llp}(cons(i, l)) = if null(l) then <1, <0>> else let _{r} = _{llp}(l) in let v = _{f′}(i, l, 2nd(_{r})) in <max(1st(_{r}), 1+1st(v)), v >

Using this equation and the base case l̃lp(nil) = < 0 >, we obtain a new definition of l̃lp:

_{llp}(l) ≜ if null(l) then <0> else if null(cdr(l)) then <1, <0>> else let _{r} = _{llp}(cdr(l)) in let v = _{f′}(car(l), cdr(l), 2nd(_{r})) in <max(1st(_{r}), 1+1st(v)), v > (22)

where f̃′ is defined in (21). This new l̃lp recursively considers tails of the input sequence and, for each element at the head of a subsequence, calls f̃′ to use it to extend recursively computed sequence lengths. It takes only O(n^{2}) time, since it calls f̃′ only O(n) times.

6.3 Towers of Hanoi

This problem computes the sequence of steps needed to move a stack of n disks, one at a time, from one peg a, initially with no larger disks on top of smaller ones, to a second peg b via a third peg c, without ever putting a larger disk on top of a smaller one [1,67]. This can be solved using a simple recursion that has one move operation, two recursive calls, and two concatenation operations, but computing the recursive function for n disks takes O(2^{n}) time, due to repeated recursive calls on the same arguments.

Pettorossi and Proietti use tupling [66] to give, to our knowledge, the most recent derivation of a linear-time program for this problem [67]. The transformation involves unfolding functions 9 times and eventually identifying that 6 function calls can be computed in a tuple using a recursion with step size 2. The resulting program consists of 9 moves and 18 concatenations in a recursive equation and 3 moves and 6 concatenations together in five other non-recursive equations.

Our method identifies a change operation that corresponds to step size 1, uses the method for discovering auxiliary information twice, and yields an incremental function that has 3 moves and 6 concatenations. Using this function directly to form a new recursion, we obtain a linear-time program that is one quarter the size of that in [67].

Function hanoi defined below computes the desired length, where skip and move are constructors, and :: is a primitive function for concatenation.

hanoi(n, a, b, c) ≜ if n ≤ 0 then skip else hanoi(n − 1, a, c, b) :: move(a, b) :: hanoi(n − 1, c, b, a) (23)

The problem can be formulated as computing hanoi incrementally under the input change operation ⟨n′, a′, b′, c′⟩ = ⟨n, a, b, c⟩⊕⟨⟩ = ⟨n+1, a, c, b⟩, identified using the method in [51] as the minimum input increment operation. After discovering auxiliary information twice using the method in this paper, we obtain

{ \atop hanoi}(n, a, b, c) ≜ < hanoi(n, a, b, c), hanoi(n, b, c, a), hanoi(n, c, a, b) > (24)

and

{ \atop hanoi}′(n, a, b, c, _{r}) ≜ if n + 1 ≤ 0 then < skip, skip, skip > else < \begin{matrix}1st(_{r}) :: move(a, c) :: 2nd(_{r}), \\ 3rd(_{r}) :: move(c, b) :: 1st(_{r}), \\ 2nd(_{r}) :: move(b, a) :: 3rd(_{r}) > \end{matrix} (25)

Clearly, computing { \atop hanoi}(n + 1, a, c, b) from scratch takes exponential time, but computing { \atop hanoi}′(n, a, b, c, r̃) takes only O(1) time. To compute the origi-

nal hanoi, note that hanoi(n, a, b, c) = 1st({ \atop hanoi}(n, a, b, c)) and, if { \atop hanoi}(n, a, b, c) =

r̃, then { \atop hanoi}(n+1, a, c, b) = { \atop hanoi}′(n, a, b, c, r̃). Using the definition of { \atop hanoi}′

in this last equation, we obtain

{ \atop hanoi}(n + 1, a, c, b) = if n + 1 ≤ 0 then < skip, skip, skip > else \begin{matrix}let _{r} = { \atop hanoi}(n, a, b, c) in \\ < \begin{matrix}1st(_{r}) :: move(a, c) :: 2nd(_{r}), \\ 3rd(_{r}) :: move(c, b) :: 1st(_{r}), \\ 2nd(_{r}) :: move(b, a) :: 3rd(_{r}) > \end{matrix} \end{matrix}

Replacing n + 1 with n, and exchanging b and c, we obtain a new definition

of { \atop hanoi}:

{ \atop hanoi}(n, a, b, c) ≜ if n ≤ 0 then < skip, skip, skip > else \begin{matrix}let _{r} = { \atop hanoi}(n − 1, a, c, b) in \\ < \begin{matrix}1st(_{r}) :: move(a, b) :: 2nd(_{r}), \\ 3rd(_{r}) :: move(b, c) :: 1st(_{r}), \\ 2nd(_{r}) :: move(c, a) :: 3rd(_{r}) > \end{matrix} \end{matrix} (26)

Clearly, we can compute hanoi using the new { \atop hanoi}, which takes O(n) time.

6.4 Longest common subsequence

Suppose we want to compute the length of the longest common subsequence of two given sequences x and y of lengths n and m, respectively [19]. This can be programmed as a straightforward recursive function that takes O(2^{n+m}) time to compute. Using the method in this paper, in a similar fashion to the path sequence problem discussed earlier, we can derive an optimized program that runs in O(n ∗ m) time.

Function lcs below computes the desired length, where x(n) is the nth element of sequence x, and y(m) is the mth element of sequence y.

lcs(n, m) ≜ if n ≤ 0 or m ≤ 0 then 0 else if x(n) = y(m) then lcs(n − 1, m − 1) + 1 else max(lcs(n, m − 1), lcs(n − 1, m)) (27)

We incrementalize lcs under the input change operation ⟨n′, m′⟩ = ⟨n, m⟩ ⊕ ⟨⟩ = ⟨n + 1, m⟩, identified using the method in [51] as a minimum input increment operation, and obtain

_{lcs}(n, m) ≜ if n ≤ 0 or m ≤ 0 then < 0 > else let v = _{lcs}(n, m − 1) in if x(n) = y(m) then < lcs(n − 1, m − 1) + 1, v > else < max(1st(v), lcs(n − 1, m)), v > (28)

and

_{lcs}′(n, m, _{r}) ≜ if n + 1 ≤ 0 or m ≤ 0 then < 0 > else if n ≤ 0 then let v = _{lcs}′(n + 1, m − 1, < 0 >) in if x(n + 1) = y(m) then < 1, v > else < 1st(v), v > else let v = _{lcs}′(n + 1, m − 1, 2nd(_{r})) in if x(n + 1) = y(m) then < 1st(2nd(_{r})) + 1, v > else < max(1st(v), 1st(_{r})), v > (29)

Computing l̃cs(n + 1, m) from scratch takes exponential time, but computing l̃cs′(n, m, r̃) takes only O(m) time. Note that lcs(n, m) = 1st(l̃cs(n, m)) and, if l̃cs(n, m) = r̃, then l̃cs(n + 1, m) = l̃cs′(n, m, r̃). Using l̃cs′ to form a new l̃cs, we obtain

\begin{matrix}lcs(n, m) ≜ 1st(_{lcs}(n, m)) \\ _{lcs}(n, m) ≜ if n ≤ 0 or m ≤ 0 then < 0 > else let _{r} = _{lcs}(n − 1, m) in _{lcs}′(n − 1, m, _{r}) \end{matrix} (30)

Function l̃cs takes O(n ∗ m) time since it calls l̃cs′ O(n) times.

The auxiliary information used in this problem and the path sequence problem belongs to a special class that can be discovered using a simpler method than the general method described in this paper. Basically, the auxiliary information needed corresponds to “intermediate results” that are computed in some branches but not some other branches of the original program. Computing such results also in these other branches forms the auxiliary information. Thus, the simpler method only needs to perform (I) a modified caching transformation, extending Step A.1, (II) incrementalization, as in Step B.2, and (III) pruning, as in Step B.3.

For example, for function f in the original program (19) for the path sequence problem, the result of f(car(l), cdr(l)) is computed in the branch where arc(n, car(l)) is true but not in the alternative branch, and the truth value of the predicate arc is independent of the local variables n and l. In function f̃ in the resulting extended program (20), this result is computed and put in variable u regardless of the truth value of arc(n, car(l)), i.e., it is computed when arc(n, car(l)) is false as well, as the auxiliary information; this information is used and maintained in the resulting incremental program (21). The details of this simpler but restricted method and the derivation for the longest common subsequence problem are described in [51].

6.5 0-1 Knapsack

This problem computes the maximum value for a subset of n items of integer weight whose total weight does not exceed w [19]. It can be programmed as the following recursive function knap(n, w), where value(n) and weight(n) are the value and weight, respectively, of the nth item. The running time of knap(n, w) is O(2^{n}).

knap(n, w) ≜ if n ≤ 0 or w ≤ 0 then 0 else if weight(n) ≥ w then knap(n − 1, w) else max(value(n) + knap(n − 1, w − weight(n)), knap(n − 1, w)) (31)

Using the method in this paper, we can find a piece of auxiliary information, knap(n, w − weight(n + 1)), that is needed for computing knap(n + 1, w), where weight(n + 1) can be any integer from 0 to w. Since w − weight(n + 1) can also be any integer from 0 to w, taking knap(n, k) for any particular value k as auxiliary information is not sufficient for efficient computation of knap(n + 1, w). We need knap(n, k) for all k such that 0 ≤ k ≤ w. We can extend the method in this paper to cache an array for these values, and use the rest of the method in the same way; this yields an optimized program that runs in O(n ∗ w) time. The precise formulation of this extension needs to be studied, but if we assume that the construct for k := i to j do a[k] := e assigns expression e to the kth element of array a for each k and returns a, we can incrementalize knap under ⟨n′, w′⟩ = ⟨n, w⟩ ⊕ ⟨⟩ = ⟨n + 1, w⟩ and obtain

{ \atop knap}(n, w) ≜ \begin{matrix}for k := 0 to w do \\ _{r}[k] := knap(n, k) \end{matrix} (32)

and

{ \atop knap}′(n, w, _{r}) ≜ for k := 0 to w do r′[k] := if n + 1 ≤ 0 or w ≤ 0 then 0 else if weight(n + 1) ≥ w then _{r}[w] else max(value(n + 1) + _{r}[w − weight(n + 1)], _{r}[w]) (33)

Computing { \atop knap}(n+1, w) from scratch takes exponential time, but computing { \atop knap}′(n, w, r̃) takes only O(w) time. Note that knap(n, w) = { \atop knap}(n, w)[w]

and, if { \atop knap}(n, w) = r̃, then { \atop knap}(n + 1, w) = { \atop knap}′(n, w, r̃). Using { \atop knap}′ to

form a new { \atop knap}, we obtain

\begin{matrix}knap(n, w) ≜ { \atop knap}(n, w)[w] \\ { \atop knap}(n, w) ≜ if n ≤ 0 or w ≤ 0 then for k := 0 to w do _{r}[k] := 0 else let _{r} = { \atop knap}(n − 1, w) in { \atop knap}′(n − 1, w, _{r}) \end{matrix} (34)

The auxiliary information used here is essentially a map and is used extensively in finite differencing of set expressions [64]. It is interesting to notice that most previous methods have used arrays as the only data structures for caching, e.g., [16,19,65], and often ended up using asymptotically more space than necessary. We use arrays only when necessary, and we found that, for most problems, arrays are not needed [51].

6.6 Local neighborhood problems

Many image processing problems compute information about local neighborhoods of objects, such as pixels, rows, and regions [82,85,87,88]. For example, an algorithm to blur a picture computes the sum of a square area for every pixel. Such programs are often easily written using loops and arrays, rather than recursive functions.

To optimize these programs using incrementalization, we identify aggregate computations on arrays performed by loops as function f, identify updates to loop variables and array subscripts as operation ⊕, and incrementalize f under ⊕ [50]; incrementalization is reduced to solving constraints on loop variables and array subscripts, for which we use Omega [69].

To increase incrementality, intermediate results and a special kind of auxiliary information for aggregate array computations are also cached and used when needed [50]. The auxiliary information there corresponds to “intermediate results” that are re-associated based on associativity of operations. For example, given an array a with n_{1} rows and n_{2} columns, the row-neighborhood-summation problem computes, for each row i (0 ≤ i ≤ n_{1}−m), the sum of the m-by-n_{2} rectangle comprising rows i through i+m−1. The following straightforward program computes this in O(n_{1}n_{2}m) time.

\begin{matrix}for i := 0 to n_{1} − m do \\ s[i] := 0; \\ for k := 0 to m − 1 do \\ for l := 0 to n_{2} − 1 do \\ s[i] := s[i] + a[i + k, l] \end{matrix} (35)

Using the method in [50], we cache a special kind of auxiliary information that corresponds to the innermost for-l-loop, which computes the sum of row i + k of array a, and incrementalize the for-k-loop with respect to i′ = i + 1 of the for-i-loop, yielding

\begin{matrix}s[0] := 0; \\ for k := 0 to m − 1 do \\ s_{1}[k] := 0; \\ for l := 0 to n_{2} − 1 do \\ s_{1}[k] := s_{1}[k] + a[k, l]; \\ s[0] := s[0] + s_{1}[k]; \\ for i := 1 to n_{1} − m do \\ s_{1}[i+m−1] := 0; \\ for l := 0 to n_{2} − 1 do \\ s_{1}[i+m−1] := s_{1}[i+m−1] + a[i+m−1, l]; \\ s[i] := s[i−1] − s_{1}[i−1] + s_{1}[i+m−1] \end{matrix} (36)

Using the ideas in this paper, we can discover the above auxiliary information by directly incrementalizing the for-k-loop with respect to i′ = i + 1 of the for -i-loop without caching them specially in the first place. Observe that computing the for-k-loop for i′ = i + 1 requires adding elements in row i + m + 1 and subtracting elements in row i; this exposes these rows as auxiliary information. This approach is more general, since, even if we exchange the for-k-loop and for-l-loop in (35) to compute the same thing, this method yields the same additions and subtractions. However, if we use the method in [50], using the special kind of auxiliary information in the innermost fork-loop yields no improvement. We believe that the basic ideas for identifying general auxiliary information in this paper also apply to discovering general auxiliary information for aggregate array computations. Detailed analyses and transformations are being developed.

7 Related work and conclusion

Work related to our analysis and transformation techniques has been discussed throughout the presentation. Here, we take a closer look at related work on discovering auxiliary information for incremental computation and on strengthening invariants for efficient computation in general.

Interactive systems and reactive systems often use incremental algorithms to achieve fast response time [5,6,10,22,33,41,71,72]. Since explicit incremental algorithms are hard to write and appropriate auxiliary information is hard to discover, the general approach in this paper provides a systematic method for developing particular incremental algorithms. For example, for the dynamic incremental attribute evaluation algorithm in [73], the characteristic graph is a kind of auxiliary information that would be discovered following the general principles underlying our approach. For static incremental attribute evaluation algorithms [42,43], where no auxiliary information is needed, the approach can cache intermediate results and maintain them automatically [54].

Strength reduction [3,17,78] is a traditional compiler optimization technique that aims at computing each iteration incrementally based on the result of the previous iteration. Basically, a fixed set of strength-reduction rules for primitive operators like times and plus are used. Our method can be viewed as a principled strength reduction technique not limited to a fixed set of rules: it can be used to reduce strength of computations where no given rules apply and, furthermore, to derive or justify such rules when necessary, as shown in the integer square root example.

Finite differencing [61–64] generalizes strength reduction to set-theoretic expressions for systematic program development. Basically, rules are manually developed for differentiating set expressions. For continuous expressions, our method can derive such rules directly using properties of primitive set operations. For discontinuous set expressions, auxiliary expressions need to be discovered and rules for maintaining them derived. Based on careful study of the finite differencing rules [64], we believe that our basic ideas for discovering auxiliary information apply to discovering auxiliary set expressions as well, and once discovered, our method can also be used to derive rules that maintain these expressions. In general, such rules apply only to very-high-level languages like SETL; our method applies also to lower-level languages like Lisp or Java.

Maintaining and strengthening loop invariants has been advocated by Dijkstra, Gries, and others [21,31,32,74] for almost two decades as a standard strategy for developing loops. In order to produce efficient programs, loop invariants need to be maintained by the derived programs in an incremental fashion. To make a loop more efficient, the strategy of strengthening a loop invariant, often by introducing fresh variables, is proposed [32]. This corresponds to discovering appropriate auxiliary information and deriving incremental programs that maintain such information. Work on loop invariants stressed mental tools for programming, rather than mechanical assistance, so no systematic procedures were proposed.

Induction and generalization [11,59] are the logical foundations for recursive calls and iterative loops in deductive program synthesis [57] and constructive logics [18]. These corpora have for the most part ignored the efficiency of the programs derived, and the resulting programs “are often wantonly wasteful of time and space” [58]. In contrast, the approach in this paper is particularly concerned with the efficiency of the derived programs. Moreover, we can see that induction, whether course-of-value induction [44], structural induction [11,13], or well-founded induction [11,59], enables derived programs to use results of previous iterations in each iteration, and generalization [11,59] enables derived programs to use appropriate auxiliary information by strengthening induction hypotheses, just like strengthening loop invariants. The approach in this paper may be used for systematically constructing induction steps [44] and strengthening induction hypotheses.

The promotion and accumulation strategies are proposed by Bird [8,9] as general methods for achieving efficient transformed programs. Promotion attempts to derive a program that defines f(cons(a, x)) in terms of f(x), and accumulation generalizes a definition by including an extra argument. Thus, promotion can be regarded as deriving incremental programs, and accumulation as identifying appropriate intermediate results or auxiliary information. Bird illustrates these strategies with two examples. However, no systematic steps were given in [8]. As demonstrated with the path sequence problem, our approach can be regarded as a systematic formulation of the promotion and accumulation strategies. It helps avoid the kind of errors reported and corrected in [9].

Other work on transformational programming for improving program efficiency, including the extension technique in [20], the transformation of recursive functional programs in the CIP project [7,12,65], and the finite differencing of functional programs in the semi-automatic program development system KIDS [77], can also be further automated with our systematic approach.

In conclusion, incremental computation has widespread applications for improving the efficiency of computations throughout computing. This paper proposes a systematic approach for strengthening invariants by discovering a general class of auxiliary information for incremental computation. It is naturally combined with methods for reusing the previous result and intermediate results to form a comprehensive approach for efficient incremental computation. The modularity of the approach lets us integrate other techniques in our framework and re-use our components for other optimizations.

Although our approach is presented in terms of a first-order functional language, the underlying principles are general and apply to other languages as well. For example, the method has been used to improve imperative programs with arrays for the local neighborhood problems in image processing [50].

A prototype system, CACHET, for incrementalization has been under development. It was first developed as a semi-automatic transformation system [46] for deriving incremental programs that exploit return values [56]. Soon after, transformations and analyses for cache and prune [55,54], fully automated, were added. Later on, a subset of the transformations for exploiting return values were made fully automatic [89], and the analysis for pruning was drastically improved [52]. A separate module for optimizing aggregate array computations [50] was added most recently. All these are implemented using the Synthesizer Generator [72], making extensive use of its scripting language, STk [30], a dialect of Scheme. CACHET has been used in teaching several graduate courses on related topics. To better integrate existing functionalities in the system, and to facilitate implementation of improvements to these functionalities, we are currently planning on a redesign and cleanup of the system. This should also allow us to easily add new extensions [51,53].

References

[1] A. V. Aho, J. E. Hopcroft, and J. D. Ullman. Data Structures and Algorithms. Addison-Wesley, Reading, Mass., 1983.

[2] A. V. Aho, R. Sethi, and J. D. Ullman. Compilers, Principles, Techniques, and Tools. Addison-Wesley, Reading, Mass., 1986.

[3] F. E. Allen, J. Cocke, and K. Kennedy. Reduction of operator strength. In S. S. Muchnick and N. D. Jones, editors, Program Flow Analysis, pages 79–101. Prentice-Hall, Englewood Cliffs, N.J., 1981.

[4] B. Alpern, R. Hoover, B. Rosen, P. Sweeney, and K. Zadeck. Incremental evaluation of computational circuits. In Proceedings of the 1st Annual ACM-SIAM Symposium on Discrete Algorithms, pages 32–42. ACM, New York, Jan. 1990.

[5] R. Bahlke and G. Snelting. The PSG system: From formal language definitions to interactive programming environments. ACM Trans. Program. Lang. Syst., 8(4):547–576, Oct. 1986.

[6] R. A. Ballance, S. L. Graham, and M. L. Van De Vanter. The Pan languagebased editing system. ACM Trans. Soft. Eng. Methodol., 1(1):95–127, Jan. 1992.

[7] F. L. Bauer, B. Möller, H. Partsch, and P. Pepper. Formal program construction by transformations—Computer-aided, intuition-guided programming. IEEE Trans. Softw. Eng., 15(2):165–180, Feb. 1989.

[8] R. S. Bird. The promotion and accumulation strategies in transformational programming. ACM Trans. Program. Lang. Syst., 6(4):487–504, Oct. 1984.

[9] R. S. Bird. Addendum: The promotion and accumulation strategies in transformational programming. ACM Trans. Program. Lang. Syst., 7(3):490– 492, July 1985.

[10] P. Borras, D. Clément, T. Despeyroux, J. Incerpi, G. Kahn, B. Lang, and V. Pascual. CENTAUR: The system. In Proceedings of the ACM SIGSOFT/SIGPLAN Software Engineering Symposium on Practical Software Development Environments, pages 14–24. ACM, New York, Nov. 1988.

[11] R. S. Boyer and J. S. Moore. A Computational Logic. ACM Monograph Series. Academic Press, New York, 1979.

[12] M. Broy. Algebraic methods for program construction: The project CIP. In P. Pepper, editor, Program Transformation and Programming Environments, pages 199–222. Springer-Verlag, Berlin, 1984.

[13] R. M. Burstall. Proving properties of programs by structural induction. The Computer Journal, 12(1):41–48, 1969.

[14] R. M. Burstall and J. Darlington. A transformation system for developing recursive programs. J. ACM, 24(1):44–67, Jan. 1977.

[15] W.-N. Chin. Towards an automated tupling strategy. In Proceedings of the ACM SIGPLAN Symposium on Partial Evaluation and Semantics-Based Program Manipulation, pages 119–132. ACM, New York, June 1993.

[16] W.-N. Chin and M. Hagiya. A bounds inference method for vector-based memoization. In Proceedings of the 1997 ACM SIGPLAN International Conference on Functional Programming, pages 176–187. ACM, New York, June 1997.

[17] J. Cocke and K. Kennedy. An algorithm for reduction of operator strength. Commun. ACM, 20(11):850–856, Nov. 1977.

[18] R. L. Constable et al. Implementing Mathematics with the Nuprl Proof Development System. Prentice-Hall, Englewood Cliffs, N.J., 1986.

[19] T. H. Cormen, C. E. Leiserson, and R. L. Rivest. Introduction to Algorithms. The MIT Press/McGraw-Hill, 1990.

[20] N. Dershowitz. The Evolution of Programs, volume 5 of Progress in Computer Science. Birkhäuser, Boston, 1983.

[21] E. W. Dijkstra. A Discipline of Programming. Prentice-Hall Series in Automatic Computation. Prentice-Hall, Englewood Cliffs, N.J., 1976.

[22] V. Donzeau-Gouge, G. Huet, G. Kahn, and B. Lang. Programming environments based on structure editor: The Mentor experience. In D. R. Barstow, H. E. Shrobe, and E. Sandewall, editors, Interactive Programming Environments, pages 128–140. McGraw-Hill, New York, 1984.

[23] J. Earley. High level iterators and a method for automatically designing data structure representation. J. Comput. Lang., 1:321–342, 1976.

[24] M. S. Feather. A system for assisting program transformation. ACM Trans. Program. Lang. Syst., 4(1):1–20, Jan. 1982.

[25] J. Field and T. Teitelbaum. Incremental reduction in the lambda calculus. In Proceedings of the 1990 ACM Conference on LISP and Functional Programming, pages 307–322. ACM, New York, June 1990.

[26] I. Flores. The Logic of Computer Arithmetic. Prentice-Hall, Englewood Cliffs, N.J., 1963.

[27] Proceedings of the 4th International Conference on Functional Programming Languages and Computer Architecture. ACM, New York, Sept. 1989.

[28] J. Glanz. Mathematical logic flushes out the bugs in chip designs. Science, 267:332–333, January 20, 1995.

[29] G. Gómez and Y. A. Liu. Automatic time-bound analysis for a higher-order language. Technical Report TR 535, Computer Science Department, Indiana University, Nov. 1999.

[30] GrammaTech, Inc. Synthesizer Generator Scripting Language Supplemental Manual, Release 5.0. GrammaTech, Inc., Ithaca, New York, 1996.

[31] D. Gries. The Science of Programming. Springer-Verlag, New York, 1981.

[32] D. Gries. A note on a standard strategy for developing loop invariants and loops. Sci. Comput. Program., 2:207–214, 1984.

[33] A. N. Habermann and D. Notkin. Gandalf: Software development environments. IEEE Trans. Softw. Eng., SE-12(12):1117–1127, Dec. 1986.

[34] R. Hoover. Alphonse: Incremental computation as a programming abstraction. In Proceedings of the ACM SIGPLAN ’92 Conference on Programming Language Design and Implementation, pages 261–272. ACM, New York, June 1992.

[35] S. Horwitz and T. Teitelbaum. Generating editing environments based on relations and attributes. ACM Trans. Program. Lang. Syst., 8(4):577–608, Oct. 1986.

[36] F. Jalili and J. H. Gallier. Building friendly parsers. In Conference Record of the 9th Annual ACM Symposium on Principles of Programming Languages, pages 196–206. ACM, New York, Jan. 1982.

[37] N. D. Jones, C. K. Gomard, and P. Sestoft. Partial Evaluation and Automatic Program Generation. Prentice-Hall, Englewood Cliffs, N.J., 1993.

[38] N. D. Jones, P. Sestoft, and H. Søndergaard. An experiment in partial evaluation: The generation of a compiler generator. In J.-P. Jouannaud, editor, Rewriting Techniques and Applications, volume 202 of Lecture Notes in Computer Science, pages 124–140. Springer-Verlag, Berlin, May 1985.

[39] S. B. Jones and D. Le Métayer. Compile-time garbage collection by sharing analysis. In FPCA 1989 [27], pages 54–74.

[40] U. Jorring and W. L. Scherlis. Compilers and staging transformations. In Conference Record of the 13th Annual ACM Symposium on Principles of Programming Languages, pages 86–96. ACM, New York, Jan. 1986.

[41] G. E. Kaiser. Incremental dynamic semantics for language-based programming environments. ACM Trans. Program. Lang. Syst., 11(2):168–193, Apr. 1989.

[42] U. Kastens. Ordered attributed grammars. Acta Informatica, 13(3):229–256, 1980.

[43] T. Katayama. Translation of attribute grammars into procedures. ACM Trans. Program. Lang. Syst., 6(3):345–369, July 1984.

[44] S. C. Kleene. Introduction to Metamathematics. Van Nostrand, New York, 1952. 10th reprint, Wolters-Noordhoff Publishing, Groningen and North-Holland Publishing Company, Amsterdam, 1991.

[45] J. Launchbury. Projections for specialisation. In B. Bjørner, A. P. Ershov, and N. D. Jones, editors, Partial Evaluation and Mixed Computation, pages 299–315. North-Holland, Amsterdam, 1988.

[46] Y. A. Liu. CACHET: An interactive, incremental-attribution-based program transformation system for deriving incremental programs. In Proceedings of the 10th Knowledge-Based Software Engineering Conference, pages 19–26. IEEE CS Press, Los Alamitos, Calif., Nov. 1995.

[47] Y. A. Liu. Principled strength reduction. In R. Bird and L. Meertens, editors, Algorithmic Languages and Calculi, pages 357–381. Chapman & Hall, London, U.K., 1997.

[48] Y. A. Liu. Dependence analysis for recursive data. In Proceedings of the IEEE 1998 International Conference on Computer Languages, pages 206–215. IEEE CS Press, Los Alamitos, Calif., May 1998.

[49] Y. A. Liu and G. Gómez. Automatic accurate time-bound analysis for high-level languages. In Proceedings of the ACM SIGPLAN 1998 Workshop on Languages, Compilers, and Tools for Embedded Systems, volume 1474 of Lecture Notes in Computer Science, pages 31–40. Springer-Verlag, June 1998.

[50] Y. A. Liu and S. D. Stoller. Loop optimization for aggregate array computations. In Proceedings of the IEEE 1998 International Conference on Computer Languages, pages 262–271. IEEE CS Press, Los Alamitos, Calif., May 1998.

[51] Y. A. Liu and S. D. Stoller. Dynamic programming via static incrementalization. In Proceedings of the 8th European Symposium on Programming, volume 1576 of Lecture Notes in Computer Science, pages 288– 305. Springer-Verlag, Berlin, Mar. 1999.

[52] Y. A. Liu and S. D. Stoller. Eliminating dead code on recursive data. In Proceedings of the 6th International Static Analysis Symposium, volume 1694 of Lecture Notes in Computer Science, pages 211–231. Springer-Verlag, Berlin, Sept. 1999.

[53] Y. A. Liu and S. D. Stoller. From recursion to iteration: what are the optimizations? In Proceedings of the ACM SIGPLAN 2000 Workshop on Partial Evaluation and Semantics-Based Program Manipulation, pages 73–82. ACM, New York, Jan. 2000.

[54] Y. A. Liu, S. D. Stoller, and T. Teitelbaum. Static caching for incremental computation. ACM Trans. Program. Lang. Syst., 20(3):546–585, May 1998.

[55] Y. A. Liu and T. Teitelbaum. Caching intermediate results for program improvement. In Proceedings of the ACM SIGPLAN Symposium on Partial Evaluation and Semantics-Based Program Manipulation, pages 190–201. ACM, New York, June 1995.

[56] Y. A. Liu and T. Teitelbaum. Systematic derivation of incremental programs. Sci. Comput. Program., 24(1):1–39, Feb. 1995.

[57] Z. Manna and R. Waldinger. A deductive approach to program synthesis. ACM Trans. Program. Lang. Syst., 2(1):90–121, Jan. 1980.

[58] Z. Manna and R. Waldinger. Fundamentals of deductive program synthesis. IEEE Trans. Softw. Eng., 18(8):674–704, Aug. 1992.

[59] Z. Manna and R. Waldinger. The Deductive Foundations of Computer Programming. Addison-Wesley, Reading, Mass., 1993.

[60] J. O’Leary, M. Leeser, J. Hickey, and M. Aagaard. Non-restoring integer square root: A case study in design by principled optimization. In R. Kumar and T. Kropf, editors, Proceedings of the 2nd International Conference on Theorem Provers in Circuit Design: Theory, Practice, and Experience, volume 901 of Lecture Notes in Computer Science, pages 52–71. Springer-Verlag, Berlin, 1995.

[61] B. Paige and J. T. Schwartz. Expression continuity and the formal differentiation of algorithms. In Conference Record of the 4th Annual ACM Symposium on Principles of Programming Languages, pages 58–71. ACM, New York, Jan. 1977.

[62] R. Paige. Transformational programming—Applications to algorithms and systems. In Conference Record of the 10th Annual ACM Symposium on Principles of Programming Languages, pages 73–87. ACM, New York, Jan. 1983.

[63] R. Paige. Programming with invariants. IEEE Software, 3(1):56–69, Jan. 1986.

[64] R. Paige and S. Koenig. Finite differencing of computable expressions. ACM Trans. Program. Lang. Syst., 4(3):402–454, July 1982.

[65] H. A. Partsch. Specification and Transformation of Programs—A Formal Approach to Software Development. Springer-Verlag, Berlin, 1990.

[66] A. Pettorossi. A powerful strategy for deriving efficient programs by transformation. In Conference Record of the 1984 ACM Symposium on LISP and Functional Programming. ACM, New York, Aug. 1984.

[67] A. Pettorossi and M. Proietti. Rules and strategies for transforming functional and logic programs. ACM Comput. Surv., 28(2):360–414, June 1996.

[68] L. L. Pollock and M. L. Soffa. Incremental global reoptimization of programs. ACM Trans. Program. Lang. Syst., 14(2):173–200, Apr. 1992.

[69] W. Pugh. The Omega Test: A fast and practical integer programming algorithm for dependence analysis. Commun. ACM, 31(8):102–114, Aug. 1992.

[70] W. Pugh and T. Teitelbaum. Incremental computation via function caching. In Conference Record of the 16th Annual ACM Symposium on Principles of Programming Languages, pages 315–328. ACM, New York, Jan. 1989.

[71] S. P. Reiss. An approach to incremental compilation. In Proceedings of the ACM SIGPLAN ’84 Symposium on Compiler Construction, pages 144–156. ACM, New York, June 1984. Published as SIGPLAN Notices, 19(6).

[72] T. Reps and T. Teitelbaum. The Synthesizer Generator: A System for Constructing Language-Based Editors. Springer-Verlag, New York, 1988.

[73] T. Reps, T. Teitelbaum, and A. Demers. Incremental context-dependent analysis for language-based editors. ACM Trans. Program. Lang. Syst., 5(3):449–477, July 1983.

[74] J. C. Reynolds. The Craft of Programming. Prentice-Hall, Englewood Cliffs, N.J., 1981.

[75] M. Rosendahl. Automatic complexity analysis. In FPCA 1989 [27], pages 144– 156.

[76] B. G. Ryder and M. C. Paull. Incremental data flow analysis algorithms. ACM Trans. Program. Lang. Syst., 10(1):1–50, Jan. 1988.

[77] D. R. Smith. KIDS: A semiautomatic program development system. IEEE Trans. Softw. Eng., 16(9):1024–1043, Sept. 1990.

[78] B. Steffen, J. Knoop, and O. Rüthing. Efficient code motion and an adaption to strength reduction. In Proceedings of the 4th International Joint Conference on TAPSOFT, volume 494 of Lecture Notes in Computer Science, pages 394–415. Springer-Verlag, Berlin, 1991.

[79] R. S. Sundaresh and P. Hudak. Incremental computation via partial evaluation. In Conference Record of the 18th Annual ACM Symposium on Principles of Programming Languages, pages 1–13. ACM, New York, Jan. 1991.

[80] L. Unnikrishnan, S. D. Stoller, and Y. A. Liu. Automatic accurate stack space and heap space analysis for high-level languages. Technical Report TR 538, Computer Science Department, Indiana University, Apr. 2000.

[81] M. Wand and W. D. Clinger. Set constraints for destructive array update optimization. In Proceedings of the IEEE 1998 International Conference on Computer Languages, pages 184–193. IEEE CS Press, Los Alamitos, Calif., May 1998.

[82] J. A. Webb. Steps towards architecture-independent image processing. IEEE Comput., 25(2):21–31, Feb. 1992.

[83] B. Wegbreit. Mechanical program analysis. Commun. ACM, 18(9):528–538, Sept. 1975.

[84] M. Weiser. Program slicing. IEEE Trans. Softw. Eng., SE-10(4):352–357, July 1984.

[85] W. M. Wells, III. Efficient synthesis of Gaussian filters by cascaded uniform filters. IEEE Trans. Patt. Anal. Mach. Intell., 8(2):234–239, Mar. 1986.

[86] D. M. Yellin and R. E. Strom. INC: A language for incremental computations. ACM Trans. Program. Lang. Syst., 13(2):211–236, Apr. 1991.

[87] R. Zabih. Individuating Unknown Objects by Combining Motion and Stereo. PhD thesis, Department of Computer Science, Stanford University, Stanford, Calif., 1994.

[88] R. Zabih and J. Woodfill. Non-parametric local transforms for computing visual correspondence. In J.-O. Eklundh, editor, Proceedings of the 3rd European Conference on Computer Vision, volume 801 of Lecture Notes in Computer Science, pages 151–158. Springer-Verlag, 1994.

[89] Y. Zhang and Y. A. Liu. Automating derivation of incremental programs. In Proceedings of the 1998 ACM SIGPLAN International Conference on Functional Programming, page 350. ACM, New York, Sept. 1998.