Differential D48283
[SCEV] Properly solve quadratic equations Authored by kparzysz on Jun 18 2018, 9:12 AM.
Details Commit r334318 exposed a problem in SolveQuadraticEquation: an add-rec was created with coefficients in i1, and division by 2 became equivalent to division by zero. See http://lists.llvm.org/pipermail/llvm-commits/Week-of-Mon-20180611/561033.html. This also fixes https://bugs.llvm.org/show_bug.cgi?id=38024.
Diff Detail
Event TimelineComment Actions See also https://bugs.llvm.org//show_bug.cgi?id=37211 . I'm tempted to just delete the code instead, given it's screwed up in multiple ways: we shouldn't be dividing by two at all, even if the bitwidth isn't one. (It looks like someone took the quadratic formula for real numbers, and blindly translated it without considering that SCEV expressions aren't real numbers.) Comment Actions I've no idea what the best fix is but at least the current patch avoids the crash I saw. Thanks! Comment Actions Rewrite SolveQuadraticEquation. This fixes https://bugs.llvm.org//show_bug.cgi?id=37211. Comment Actions Needs more test coverage of the overflow cases.
Comment Actions
Comment Actions Following hits an assertion failure with "opt -analyze -scalar-evolution" with the latest patch, I think due to the inconsistent handling of signed-ness. define signext i32 @func() {
entry:
br label %loop
loop:
%ivr = phi i32 [ 0, %entry ], [ %ivr1, %loop ]
%inc = phi i32 [ 0, %entry ], [ %inc1, %loop ]
%acc = phi i32 [ -1, %entry ], [ %acc1, %loop ]
%ivr1 = add i32 %ivr, %inc
%inc1 = add i32 %inc, -1 ; M = inc1 = inc + N = X + N
%acc1 = add i32 %acc, %inc ; L = acc1 = X + Y
%and = and i32 %acc1, -1 ; iW
%cond = icmp ule i32 %and, -3
br i1 %cond, label %exit, label %loop
exit:
%rv = phi i32 [ %acc1, %loop ]
ret i32 %rv
}Comment Actions Well, for one 2/2 seems to be -1... Code: APInt TT, UU; APInt::udivrem(APInt(72, 2), APInt(72, 2), TT, UU); dbgs() << "TT=" << TT << '\n'; Output: TT=-1 This is because we handle the special case of LHS==RHS via "Quotient=1", which sets the bitwidth to 1... But there is also something else going on with this test. The code in getNumIterationsInRange assumes that it's the upper bound of the range that will be crossed, while this case crosses the lower bound. We should really solve for both and eliminate the wrong solution (the one for which the assertion fails). Comment Actions
The current getNumIterationsInRange implementation is consistent with "equals 0 or wraps around 2^BW in the unsigned modulo 2^BW arithmetic" (but that's not what SolveQuadraticEquation actually computes at the moment). Comment Actions Generic comments like this are not very helpful. Please be more specific in explaining your concerns. The addrec in your example is {x,+,-1,+,-1}, where x is some starting value. In unsigned arithmetic this expression has an unsigned overflow right at the first iteration (regardless of x), but it doesn't mean that it goes out of the range in the first step. Comment Actions I can change the comment to something like "returns the smallest positive x such that q(x-1) and q(x) belong to different intervals in the form of [k*2^BW, (k+1)*2^BW)". Comment Actions Check both ends of the range in getNumIterationsInRange, change comment for SolveQuadraticEquation. Comment Actions The general problem we're trying to solve in getNumIterationsInRange is how many iterations the SCEV {L,+,M,+,N} is in the range [R,S]. There's a simple, inefficient algorithm for this: just call evaluateAtIteration repeatedly until the value isn't in range. But that's obviously way too expensive, and I don't know how to solve the general problem more efficiently. So let's take a subset of the problem. The simplest subset is SCEVs where the last value in range is before the point of the first unsigned overflow. We can solve this by normalizing the range to [R, UINT_MAX], then finding the first iteration where {L,+,M,+,N} has unsigned overflow, then checking whether the iteration after the overflow is in range. This can be computed either using a straight binary search, or by converting it to an inequality over real numbers and using algebra. I thought this is what you were aiming for, since it's very close to what you've implemented? We can generalize that a bit by adding some additional code in the case where we fail to compute a solution: apply a binary NOT to both the SCEV and the range, and try again. Another possible subset is SCEVs where the last value in range is before the point of the first signed overflow. Basically, normalize the range to [R, INT_MAX], treat the coefficients and range as signed numbers, solve the system of inequalities over real numbers, then verify the solution applies to the original modular problem. This is a little more complicated because the vertex of the parabola could be inside the range. Comment Actions Another crash: define signext i32 @func() {
entry:
br label %loop
loop:
%ivr = phi i32 [ 0, %entry ], [ %ivr1, %loop ]
%inc = phi i32 [ -100000, %entry ], [ %inc1, %loop ]
%acc = phi i32 [ 100000, %entry ], [ %acc1, %loop ]
%ivr1 = add i32 %ivr, %inc
%inc1 = add i32 %inc, 1 ; M = inc1 = inc + N = X + N
%acc1 = add i32 %acc, %inc ; L = acc1 = X + Y
%and = and i32 %acc1, -1 ; iW
%cond = icmp sgt i32 %and, 5
br i1 %cond, label %exit, label %loop
exit:
%rv = phi i32 [ %acc1, %loop ]
ret i32 %rv
}
Comment Actions My approach here, in a single sentence, is to solve the quadratic equation q(n) = 0 using the formula for real numbers. Doing modulo arithmetic is tricky and it's easy to make a false assumption. For example: in modulo arithmetic there is no linear order, or if someone insists on an order, then we lose n+1>n for each n (so it will have different properties from what holds in R or Z). Additionally, a square root will now have multiple values: in modulo 16, 2^2 = 4, but also 14^2 = 4, because -2 is now "positive" as 14. Because of that I wanted to move away from the modulo arithmetic, and instead operate on integers. I get the "simulated" integers by extending everything to a bit width that will be sufficient for all required operations. Let's assume for a moment that we're not interested in inequalities, but only in the exact solutions to q(n) = 0. In integers there can only be 0, 1 or 2 such solutions, but with the modulo arithmetic that we started with there can be infinitely many solutions. Those extra solutions show up when we have something like 16*16 = 0 (mod 256). So instead of solving q(n) = 0 we have to find all solutions to q(n) = k * 2^BW and pick the smallest non-negative one. A lot of the code in SolveQuadraticEquation tries to figure out how to rewrite q(n) = k * 2^BW as q(n) - p = 0, so we can still solve q'(n) = 0 for some updated q'. It does that by exploiting the properties of a quadratic function: decreasing before the vertex, increasing afterwards, assuming positive coefficient at n^2. It handles the case where we start with some, let's say positive value, then the value decreases (but remains positive), then goes back up and hits 2^BW eventually. I'll get back to this case later on. There is still that momentary assumption that we only solve for exact n for which q(n) = 0, but the analysis that this is a part of is also interested in inequalities. The inequality here has the form of q(n) > 0 or q(n) < 0 (or greater-than-or-equal, etc), since the addrec gets shifted by the values at the ends of the range. The interesting solution here is the first step where the values change from negative to positive, or from positive to negative. In the original (modulo) arithmetic this can happen when the values cross the 2^BW boundary (any multiple of it, including 0), or when the values cross the signed overflow (i.e. an odd multiple of 2^(BW-1)). At the time of solving the equation we don't really know which one is the correct one, but since we already have results for k * 2^BW, we can stick with that. This is where the signed overflow is ignored. The returned value may not be the correct solution, but since it is validated by the callers, it is ok to return an "educated guess". Back to the case with decreasing values---in modulo arithmetic adding -1 would actually be adding 2^BW-1. That would always be an overflow, so strictly speaking we'd have to stop the first time it happens. The problem is that this is rarely the right thing to do. To avoid that, the coefficients are treated as signed numbers.
Comment Actions Separated solving of a quadratic addrec into two cases:
Both cases now guarantee that if a valid solution is found (i.e. nullifying the addrec or leaving the range), it will be the smallest non-negative one. Updated testcase to include checks for calculating failing solutions (i.e. non-minimal, etc.). In some cases an addrec will only become 0 after a number of signed and unsigned overflows has occurred. SolveQuadraticAddRecExact will still return the location of the first unsigned overflow (which will fail validation), however separating it from SolveQuadraticAddRecRange makes it easier to change it in the future to look for an exact solution. I have run this though 100k randomly generated testcases to detect potential crashes (none have occurred). Comment Actions I like the separate entry points for equality and inequality. I was curious about exact solutions, so I did a bit of searching; https://arxiv.org/pdf/1711.03621.pdf describes the complete solution for equality.
Comment Actions Btw, earlier on I found this paper, but I decided against it because it wouldn't solve inequalities: https://arxiv.org/abs/1507.07513 (and it would be more work...).
Comment Actions Addressed review comments. Added comments, an assertion to ensure that the initial value is in range in SolveQuadraticAddRecRange, more debug messages.
Comment Actions I didn't review the math, but I think there are some easy steps we can take to make this more watertight. Firstly, the code that operates purely on APInt should be moved to APInt.h. With that, we should be able to brute-force all possible inputs for a small bitwidth (i3 or i4) and check that the solutions are correct. What do you think? Comment Actions Moved the SolveQuadraticEquation to APIntOps and renamed it to SolveQuadraticEquationWrap. Changed the comments to no longer refer to any multiplication by 2, now all actions related to that are only in ScalarEvolution.cpp. Tested all widths between 2 and 10 (inclusive) using a brute-force test from D50095. Comment Actions Now that the code is in APInt we should have some tests in unittests/ADT/APIntTest.cpp. Can we pull in some subset of D50095 that runs within a reasonable amount of time (only bitwidth 4, for instance)? I also have some "black box" comments inline on how the function is documented.
Comment Actions Only one non-trivial question inline.
Comment Actions lgtm
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||