Bill McCloskey’s Blog

I discovered a neat technique in program analysis today, so I’ll try to explain it here. The first section gives a general overview of inferring integer properties of programs, and the second part describes what I found.

Numerical Domains and Difference Constraints

Inferring numerical properties of programs is a pretty important problem in program analysis. Here’s is an example.

if (x > 0) { y = x; } else { y = 10; }

We can infer from that after this program executes, y > 0. Why? Because in the true branch, y is assigned x, which we know to be positive. In the false branch, y is assigned 10, which is also positive.

There are a variety of techniques to infer properties like this. The most common ones are based on a technique called abstract interpretation. The most well-known technique of all is called the polyhedron domain, discovered by Cousot and Halbwachs. It is quite complex to implement, but it is guaranteed to infer all possible linear inequalities satisfied by a program. Unfortunately, it’s also very slow.

A more efficient and simpler technique goes by the name of “difference constraints” or “bounded difference matrices” or the like. In this post I will call it DC. DC can infer any invariant of the form x – y <= c or x <=c or x >= c. Most of the invariants that programmers rely on are of this form, so DC is kind of a sweet spot domain. (Note that constraints like x == y can be converted to the conjunction of x -y <= 0 and y – x <= 0.)

The domain works by putting all constraints into the form x – y <= c. Constraints of the form x <= c or x >= c become x – 0 <= c or 0 – x <= -c (i.e., zero is just treated as a variable). Then it saturates the constraint system by generating all constraints that follow from the known ones by linear combination. For example, if you know that x – y <= 5 and that y – z <= -3, then you can infer that x – z <= 2. The algorithm for saturation is pretty simple: you express the constraints as a weighted directed graph (where the variables are nodes and the constraints are edges) and then use the Floyd-Warshall all-pairs shortest path algorithm.

Once the constraint system has been saturated, you can perform all sorts of operations. First, it’s easy to check if the constraints are satisfiable by looking for a constraint of the form x – x <= c, where c < 0. If you want to see if one set of constraints C1 implies another set C2, you just take every constraint x – y <= c in C2 and make sure that C1 has a constraint of the form x – y <= c’, where c’ <= c.

Probably the most common operation is assignment. If you know that some set of constraints C is true before an assignment, then you want to know the set C’ that holds after the assignment. Assume the assignment is x := e. First we add a new constraint “tmp == e” to the system (using the translation of equality above) and then we saturate the system. Then we eliminate any constraints that talk about x from the system, and rename tmp to x.

Another important operation is the join. This operation takes the constraints generated by two branches of a conditional statement (call them C1 and C2) and tries to find the best set of constraints C that abstracts both branches. For this to be true, the constraints in C1 must imply C and the constraints in C2 must also imply C. To generate C, we take all pairs of variables x and y and look for a constraint of the form x – y <= c1 in C1 and x – y <= c2 in C2. Then we add x – y <= max(c1, c2) to C.

Sparsity

The problem with DC is that saturation is expensive: it’s an O(n^3) algorithm. My goal was to fix this problem, at least in the common case. In the common case, the initial set of constraints is sparse: when you form the directed graph of constraints, there aren’t very many edges. The problem is that saturation adds a lot more edges, and those edges make later operations much more expensive.

In my experience with difference constraints, most of the extra edges added by saturation are caused by the zero node (remember that we treat zero as a variable). As an example, consider these constraints: 0 <= x <= 10, 5 <= y <= 20, 10 <= z <= 30. If we draw the graph for this system, there will be edges between x and the zero node, between y and the zero node, and between z and the zero node. But once we saturate, there will be edges between all pairs of variables saying things like y – z <= 10 and z – y <= 25.

It turns out that for some operations like assignment and checking satisfiability, there’s no need to saturate across the zero node. That is, we allow saturation to take the linear combination of x – y <= c1 and y – z <= c2 as long as y is not the zero node. It’s pretty easy to modify the Floyd-Warshall algorithm to do this. And when you do, the edges that were added in the example in the last paragraph will no longer be added, so everything will be faster and we won’t have lost any precision.

Things are a bit more complicated for join. We can still do the same trick of restricting saturation, but now we will actually lose the ability to infer some invariants. Consider the case where we know x = 1 and y = 2 on one branch and x = 5 and y = 6 on the other. We would like to infer y – x = 1. With normal saturation we will, but not with restricted saturation.

Here’s the solution. We use restricted saturation to get a joined result C from the inputs C1 and C2. Then we iterate over all pairs of variables x and y. We find the shortest path from x to y in C1 by looking for an edge directly from x to y as well as edges from x to the zero node and zero to y. Call this distance c1. Do the same in C2 to get distance c2. Now do the same in C to get c. If max(c1, c2) < c, then we simply add an edge in C directly from x to y with distance max(c1, c2). Otherwise we do nothing.

This technique achieves the same precision as using full saturation, but it adds many fewer edges in practice. My hope is that will be much faster as well, although I haven’t tested it much yet.

Here’s one way to think of this technique. Saturation produces the transitive closure of our constraint graph. However, after we perform operations like join, we’d like to take the transitive reduction so that they have a minimal representation. Unfortunately, the transitive reduction is expensive to compute and isn’t even unique for our graphs. But since we know that zero is the crucial node in most of our graphs, we can compute a sort of transitive reduction, except that the result is only reduced at the zero node.

One topic that has always interested me is performance analysis. In particular, it always frustrates me when I can’t get consistent timings for a piece of code. No matter how much work I do, the numbers are always all over the map. My friend Dave had an interesting blog post related to this subject. He showed that the running times for typical benchmarks do not fall into a Gaussian distribution as we typically expect. He estimated the actual distribution by doing lots of runs and then used it to compute things like confidence intervals, which was cool.

However, I’m a perfectionist and I’d like to know exactly how long my code really takes under ideal conditions. Unfortunately, there are many sources of non-determinism nowadays. First, the CPU executes instructions out of order, so the time your code takes may depend on what instructions were in the pipeline before. Cache effects and swapping are huge sources of non-determinism if they occur. Laptops may throttle the clock to save power. The OS may schedule another process in the middle of your timing loop. And finally, an interrupt may arrive and ruin everything.

Some of these issues are not too hard to address. Pipeline non-determinism isn’t a big deal, since it usually only causes the timinings to change by a few cycles. You can mitigate cache effects by warming the cache before doing the actual timing. Your laptop can be instructed to run in “performance mode” so the clock speed never changes. And you can run your benchmark at real-time priority so that the scheduler never steals a slice from you.

But even after doing all these steps I still see timing variations in the millisecond range. Given the speed of modern CPUs, that’s a huge amount of time. By a process of elimination, it would seem that interrupts are a major source of latency.

To verify if the problem really is with interrupts, I wrote a kernel module that disables interrupts and then runs my benchmark. It’s a very simple benchmark that runs entirely in cache and does nothing but additions and comparisons. You can find it here. (Instructions for compiling your own Linux kernel module can be found in Chapter 2 of Linux Device Drivers.) The next few paragraphs assume you’ve read the code.

For the benchmark to work properly, you have to have a single core with no hyperthreading (a laptop typically works). To disable interrupts, the code acquires a spinlock. It uses the rdtsc instruction to get a cycle-accurate timestamp. The cpuid instruction is used to serialize execution, as described by Intel.

I observed a few things. First, I was able to get a very accurate measurement for the latency of the CPUID/RDTSC instruction pair. On my machine, it always takes 166 cycles. The benchmark loop itself is a little more variable. It usually takes 2,000,014,244 cycles (about 1.07 seconds), but on rare occasions (1 in 10 times or so) it varies by 10 cycles or so. Being mostly ignorant of computer architecture, I have to assume that this small bit of non-determinism is caused by pipelining, or branch prediction, or something else that I can’t control.

Even though the benchmark doesn’t use the memory at all, warming the loop turns out to be important. Without this step, it takes 2,000,025,251 cycles in the common case, and it varies by as much as 150 cycles. I guess this is due to pipelining.

Disabling interrupts has a much greater effect. If I comment out the spinlock, then the running time is considerably longer. I observed differences of nearly 4 milliseconds.

In the future, I would like to investigate a bit more. First, what effect does using a real-time kernel have? Hopefully the interrupt latency would go down significantly. Second, I would like to track down which interrupts are actually causing the problem. Finally, I would like to learn enough about computer architecture to understand the sources of variation in the interrupt-free cases above. (As an aside, while working on this today I found a neat web page about the x86 architecture.)

Finally, I’d like to say that I realize this experiment was completely useless for doing actual benchmarks. No sane person would ever run their benchmark suite in a kernel module. My real goal was to understand the sources of non-determinism in timing experiments.

For the last few weeks I’ve been wondering about how exactly the compiler implements varargs functions. Tonight, I finally got a chance to sit down and take a look at the code it generates. But first let’s look at how varargs works in theory.

To write a varargs function, you include <stdarg.h>. Typically, you’ll have a non-varying argument that will describe the number of and types of the rest of the arguments (the format string for printf, say). Then, for each variable argument, you call va_arg(T), where T is the type of the argument (va_arg is a macro that expands to a compiler intrinsic, so it can do funny things like take a type as an argument). The compiler keeps a pointer to the argument on the stack that you’re currently dealing with. Each time you call va_arg, that pointer moves down according to the size of the type that you specify. This works because the arguments all appear in sequence on the stack. Seems pretty straightforward.

To understand the source of my confusion, consider the call “printf(”%c%c”, 97, 97)”. This will print the letter ‘a’ twice. Now, the compiler doesn’t know anything about printf–it just sees one fixed argument and two varying arguments. It doesn’t know that the integer 97 is a char and not an int or a long long. This is an issue because if 97 is passed as an int while printf interprets it as a char, then the pointer into the stack that va_arg keeps will get screwed up. Why doesn’t this happen in practice?

It turns out that gcc uses an argument passing convention that manages to avoid this problem most of the time, but not always. All integral types (except long long) are passed on the stack as 4-byte integers. And printf interprets them all as 4-byte integers. Even though you write %c, printf calls va_arg(int) for that argument. It does a similar thing for floating-point numbers: they’re all promoted to doubles. And %f is actually calling va_arg(double).

This still leaves a few corner cases. First, what happens if you call printf(”%lld %lld”, 12, 12)? Note that although gcc could be conservative and promote all integral types to long long, it doesn’t go quite that far, presumably for efficiency. As a result, the printf call won’t work–it prints garbage. Luckily, gcc prints a warning because it has special knowledge of printf-style format strings and it knows that you really should have written printf(”%lld %lld”, 12LL, 12LL). Unfortunately, it can’t handle this problem if you use your own kind of format string (or some other kind of calling convention all together).

The other corner case is what happens when you call va_arg(char) or va_arg(short) or va_arg(double). Once again, the gcc developers provide a helpful warning:

test.c: In function ‘f’:
test.c:7: warning: ‘short int’ is promoted to ‘int’ when passed through ‘...’
test.c:7: note: (so you should pass ‘int’ not ‘short int’ to ‘va_arg’)
test.c:7: note: if this code is reached, the program will abort

And indeed, this code generates an illegal instruction error if you reach it.

What does this mean in practice? If you write your own function, then either declare it to take printf-style format strings or else don’t ever try passing a long long to it. Otherwise you’re setting yourself up for weird bugs.