[SCEV] Properly solve quadratic equations

[kparzysz]

Summary

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.

[efriedma]

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.)

[uabelho]

I've no idea what the best fix is but at least the current patch avoids the crash I saw. Thanks!

[kparzysz]

Rewrite SolveQuadraticEquation. This fixes https://bugs.llvm.org//show_bug.cgi?id=37211.

[kparzysz]

Ping.

[efriedma]

Needs more test coverage of the overflow cases.

lib/Analysis/ScalarEvolution.cpp

8427

I don't like hardcoded numbers without an explanation... why is 32 special? Why is 1 special?

8489

I think you can overflow here, and a few other places where you multiply numbers. Not sure exactly how many more bits you need.

10760–10761

Would it make sense to just return getCouldNotCompute() here? Given the SolveQuadraticEquation() API, if RootVal is in the range, something overflowed, I think.

[kparzysz]

1. Added solving for arbitrary overflow, i.e. q(x) = 2^n for n \in Z-{0}, (or mod 2^n).
2. Expanded comments.
3. Added more testcases (incl. from PR38024).

[kparzysz]

Ping.

[efriedma]

lib/Analysis/ScalarEvolution.cpp

8347

Please clarify what happens if the addrec has signed overflow. Please document the conditions under which the algorithm fails to find a solution.

8386

The algorithm should still work even if you treat "i1 true" as -1, I think. Actually, you end up solving the exact same equation, because of the A.isNegative() check.

10760–10761

Please make sure we have a testcase that hits this "getCouldNotCompute".

10762–10763

Please make sure we have a testcase that hits this "getCouldNotCompute".

test/Analysis/ScalarEvolution/solve-quadratic-overflow.ll

3

Please integrate this into the other testcase, with the same comments/CHECK lines.

test/Analysis/ScalarEvolution/solve-quadratic.ll

2

debug-only needs "REQUIRES: asserts".

[kparzysz]

lib/Analysis/ScalarEvolution.cpp

8386

You're right. The check for A.isNegative() was added later.

10762–10763

This shouldn't really happen. I've changed it to an assert.

[kparzysz]

Addressed the comments.

[efriedma]

lib/Analysis/ScalarEvolution.cpp

8348

"A signed overflow is ignored" contradicts the implementation; you're treating the coefficients as signed values (by sign-extending them).

[kparzysz]

lib/Analysis/ScalarEvolution.cpp

8348

I'm solving q(x) = k * 2^BW. The reason why I'm sign-extending the values is to keep the behavior around 0 the same. I am avoiding any overflows in the calculations and I simulate the unsigned one via "2kR".

[efriedma]

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
}

[kparzysz]

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).

[efriedma]

The code in getNumIterationsInRange assumes that it's the upper bound of the range that will be crossed

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).

[kparzysz]

In D48283#1167248, @efriedma wrote:

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).

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.

[kparzysz]

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)".

[kparzysz]

Check both ends of the range in getNumIterationsInRange, change comment for SolveQuadraticEquation.

[efriedma]

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.

[efriedma]

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
}

lib/Analysis/ScalarEvolution.cpp

8352

Please explicitly state this is treating the coefficients L, M, and N as signed integers.

[kparzysz]

Yeah, this one is a bug in updating of the coefficients.

[kparzysz]

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.

[kparzysz]

Fixed some bugs, expanded comments, added more tests.

[efriedma]

lib/Analysis/ScalarEvolution.cpp

10784

I'm not convinced this is correct: why is the "valid" solution guaranteed to be smaller?

[kparzysz]

lib/Analysis/ScalarEvolution.cpp

10784

Are you referring to line 10720? Because it's the number of iterations. Both solutions are positive and both are validated as leaving the range, so the smaller one will happen first.

[kparzysz]

lib/Analysis/ScalarEvolution.cpp

10784

In this line, either UpperValid or LowerValid are true, but not both. If UpperValid is false, then the lower solution is.

[efriedma]

lib/Analysis/ScalarEvolution.cpp

10784

I'm specifically concerned about the case where U.ult(L) is true, and UpperValid is false (so we return a "valid" solution, but not the smallest solution). Or is that impossible for some reason?

[kparzysz]

lib/Analysis/ScalarEvolution.cpp

10784

Yeah, there is a problem.

We need to handle signed overflow as well. The validity of the solution checked here only checks that it leaves the range, but it doesn't check that it's the first one to do so. The solving functions needs to give the guarantee of being the smallest solution.

[kparzysz]

Separated solving of a quadratic addrec into two cases:

1. 1. Solving for an exact solution (i.e. calculating unsigned overflow: q(n) = k * 2^BW, where 0 is considered a special case of it).
2. 2. Solving for exiting from range (includes signed and unsigned overflow).

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).

[efriedma]

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.

lib/Analysis/ScalarEvolution.cpp

8384

negate() can overflow; if that overflow is impossible because of the way the inputs are constructed, please add an assertion.

8474

Is it possible for the truncation to fail? For an AddRec of bitwidth N, evaluateAtIteration(2^N)==evaluateAtIteration(0), so any solution must be less than 2^N. (And all users should treat the backedge-taken count as an unsigned quantity.)

8495

You might want to explicitly state that there is no solution in some cases.

8552

Do we have to return None here? Even if we can't return a root, we should still be able to find the point where it wraps.

If you think None is the correct representation, please ensure SolveQuadraticAddRecRange doesn't succeed in that case.

8576

Do we need to check which one is smaller before we check LeavesRange?

[kparzysz]

lib/Analysis/ScalarEvolution.cpp

8474

It is possible. The test test/Analysis/ScalarEvolution/solve-quadratic-i1.ll has coefficients in i1, but the solution is 2 (the loop iterates 3 times). I'm not sure if this is specific to i1 or not. It hasn't happened in any test I've run outside of i1.

8552

Yes, this case could likely be solved, but we fail to do it.

The code above that updates the C in Ax^2+Bx+C, does that in a way to find the smallest x that causes an overflow (or zero). It uses the standard "real number" analysis for that, which only assures the existence of such a solution (as a real number), but does nothing to make sure that there will be an integer that satisfies the requirements. If a single iteration skips over both zeros of the updated equation then both, the previous and the next iteration will produce positive values and this case will go "unnoticed" by the loop. The code updating C could, in theory, do that to make sure that an integer solution exists, I just haven't thought of an elegant way to do it.

The point about SolveQuadraticAddRecRange is definitely valid and I have fixed it.

8576

We don't need to, but we can.

[kparzysz]

Addressed comments.

[kparzysz]

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...).

[efriedma]

lib/Analysis/ScalarEvolution.cpp

8423

I think BitWidth + 1 should be sufficient? At least, that should be enough to convert the AddRec to an equation without losing any information. (Like you note in the comments later, the equation is L + nM + n(n-1)/2 N = 0 mod BitWidth, or 2L + 2M n + n(n-1) N = 0 mod BitWidth+1.)

8428

Do we prefer to sign-extend just so abs(A) is small in common cases? I think the choice of sign-extension vs zero-extension is arbitrary otherwise. Probably worth explicitly noting with a comment.

8474

Oh, I'm just wrong; in general, the solution for an AddRec of bitwidth B needs up to B+1 bits. Maybe worth noting that explicitly.

[kparzysz]

lib/Analysis/ScalarEvolution.cpp

8423

I wanted to avoid any overflows, the motivation was 2M-N. I suppose it cannot overflow BitWidth+1 with M and N being BitWidth wide.

8428

Are you talking only about this? There is another extension in SolveQuadraticEquation (at line 8325). Originally there was only one sign-extension, but after reorganizing this code, it has been split into two. I kept them identical to the original one (i.e. both as sign-extension).

In general the choice of sign-extension was actually to avoid "false positives" when checking for overflow. For example, take q(x) = (x+1)(x-2) = x^2 - x - 2, in i3. The zeros are -1 and 2, q(0) = q(1) = -2 (i.e. nothing exceeding the signed range of i3), so we can expect the calculated solution to be 2. With zero-extended values, the function would be r(x) = x^2 + 7x + 6 with both solutions negative: -6 and -1. The equation r(x) = 0 would then become r(x) = 8 to solve for an unsigned overflow, i.e. x^2 + 7x - 2 = 0. The root of that is approx. 0.27, so the returned solution would be 1.

[kparzysz]

Addressed further comments.

[efriedma]

lib/Analysis/ScalarEvolution.cpp

8428

I'm talking about specifically this one: the extension here will still be necessary even if we add an exact modular quadratic SCEV solver.

I guess the answer is essentially the same, but deserves its own comment.

[efriedma]

lib/Analysis/ScalarEvolution.cpp

8607

Min.hasValue() is always true here.

8611

How do you ensure the correct solution isn't some value between Min and Max?

Similarly, in the case where SolveForBoundary returns None, how can you be sure the correct solution isn't some value between the Max for this boundary and the Min of the other boundary?

[kparzysz]

lib/Analysis/ScalarEvolution.cpp

8607

So is Max.hasValue().

8611

How do you ensure the correct solution isn't some value between Min and Max?

Assuming that Min and Max are different values, one of them is when the first signed overflow happens, the other is when the first unsigned overflow happens. Crossing the range boundary is only possible via an overflow (treating 0 as a special case of it, modeling overflow as crossing k*2^W for some k).

The interesting case here is when Min was eliminated as an invalid solution, but Max was not. The argument is that if there was another overflow between Min and Max, it would also have been eliminated if it was considered.

For a given boundary, it is possible to have two overflows of the same type (signed/unsigned) without having the other type in between: this can happen when the vertex of the parabola is between the iterations corresponding to the overflows. This is only possible when the two overflows cross k*2^W for the same k. In such case, if the second one left the range (and was the first one to do so), the first overflow would have to enter the range, which would mean that either we had left the range before or that we started outside of it. Both of these cases are contradictions.

Similarly, in the case where SolveForBoundary returns None, how can you be sure the correct solution isn't some value between the Max for this boundary and the Min of the other boundary?

Assume that we had such Max_A and Min_B corresponding to range boundaries A and B and such that Max_A < Min_B. If there was a solution between Max_A and Min_B, it would have to be caused by an overflow corresponding to either A or B. It cannot correspond to B, since Min_B is the first occurrence of such an overflow. If it corresponded to A, it would have to be either a signed or an unsigned overflow that is larger than both eliminated overflows for A. But between the eliminated overflows and this overflow, the values would cover the entire value space, thus crossing the other boundary, which is a contradiction.

[kparzysz]

Addressed review comments.

Added comments, an assertion to ensure that the initial value is in range in SolveQuadraticAddRecRange, more debug messages.

[kparzysz]

lib/Analysis/ScalarEvolution.cpp

8600

Doh. Leftover '?'.

[efriedma]

LGTM, but this is complicated code, so hopefully someone else will look as well.

[sanjoy]

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?

[kparzysz]

I can do the brute force tests.

[kparzysz]

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.

[sanjoy]

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.

include/llvm/ADT/APInt.h

2174

Would be nice to clarify the bitwidth in which you evaluate q(n) for (b). It can't be BW since any value in [0,2^BW) will either lies in exactly one such interval or all such intervals, depending on how you slice it.

2178

I'm not sure this is correct -- sub 5, 1 does not unsign-overflow right? (since zext(5) - zext(1) == zext(5 - 1))?

2183

Is (a) relevant any more? All the coeffs are APInt.

2187

Can you add a line on how this "q(n) > T, and q(n+1) < T" ties in with the return value of this function specified at the very top "n >= 1 and q(n-1) and q(n) belong to two different intervals..."?

[kparzysz]

include/llvm/ADT/APInt.h

2178

That's true, but then sub 5, 1 is not the same as add 5, -1. I'm treating sub x, y as equivalent to add x, -y. I added a comment explaining this.

[kparzysz]

Addressed comments. Added the verification code to APIntTest.

[sanjoy]

Only one non-trivial question inline.

unittests/ADT/APIntTest.cpp

2368

Why do we need this special case?

2382

s/VX/ValueAtX

2384

I don't quite understand this. Don't you need to check that the overflow bits are different from when the equation is evaluated at X-1 (and not from when the equation is evaluated at 0, which is what it looks like you're doing here)?

2387

Probably better to call this EquationToString?

[kparzysz]

unittests/ADT/APIntTest.cpp

2368

We don't really. Removed.

2384

The overflow bits should not change between 0 and X-1 inclusive. We're looking for the first X such that they do change, or that the value is 0 (whichever happens first).

[kparzysz]

Addressed comments.

[sanjoy]

lgtm

include/llvm/ADT/APInt.h

2193

Do you meant to say that "if the coeffs are positive then use a range of BW-1 to find a solution for signed overflow"?

[kparzysz]

include/llvm/ADT/APInt.h

2193

It should work with with any coefficients. I will clarify the comments.

include/llvm/ADT/APInt.h

 namespace llvm {
 class FoldingSetNodeID;
 class StringRef;
 class hash_code;
 class raw_ostream;
 
 template <typename T> class SmallVectorImpl;
 template <typename T> class ArrayRef;
+template <typename T> class Optional;
 
 class APInt;
 
 inline APInt operator-(APInt);
 
 //===----------------------------------------------------------------------===//
 //                              APInt Class
 //===----------------------------------------------------------------------===//

 }
 
 /// Return A unsign-divided by B, rounded by the given rounding mode.
 APInt RoundingUDiv(const APInt &A, const APInt &B, APInt::Rounding RM);
 
 /// Return A sign-divided by B, rounded by the given rounding mode.
 APInt RoundingSDiv(const APInt &A, const APInt &B, APInt::Rounding RM);
 
+/// Let q(n) = An^2 + Bn + C, and BW = bit width of the value range
+/// (e.g. 32 for i32).
+/// This function finds the smallest number n, such that
+/// (a) n >= 0 and q(n) = 0, or
+/// (b) n >= 1 and q(n-1) and q(n), when evaluated in the set of all
+///     integers, belong to two different intervals [Rk, Rk+R),
+///     where R = 2^BW, and k is an integer.
+/// The idea here is to find when q(n) "overflows" 2^BW, while at the
+/// same time "allowing" subtraction (treating it as equivalent to the
+/// addition of negated numbers). In unsigned modulo arithmetic a
+/// subtraction would always count as an overflow, but here we want to
+/// allow values to decrease and increase as long as they are within
+/// the same interval.
+///
+/// This function returns None if after finding k that minimizes the
+/// positive solution to q(n) = kR, both solutions are contained between
+/// two consecutive integers.
+///
+/// There are cases where q(n) > T, and q(n+1) < T (assuming evaluation
+/// in arithmetic modulo 2^BW, and treating the values as signed) by the
+/// virtue of *signed* overflow. This function will *not* find such an n,
+/// however it may find a value of n satisfying the inequalities due to
+/// an *unsigned* overflow (if the values are treated as unsigned).
+/// To find a solution for a signed overflow, use a range with of BW-1.
+///
+/// The returned value may have a different bit width from the input
+/// coefficients.
+Optional<APInt> SolveQuadraticEquationWrap(APInt A, APInt B, APInt C,
+                                           unsigned RangeWidth);
 } // End of APIntOps namespace
 
 // See friend declaration above. This additional declaration is required in
 // order to compile LLVM with IBM xlC compiler.
 hash_code hash_value(const APInt &Arg);
 } // End of llvm namespace
 
 #endif

lib/Analysis/ScalarEvolution.cpp

static const SCEV *SolveLinEquationWithOverflow(const APInt &A, const SCEV *B,
   // 4. Compute the minimum unsigned root of the equation:
   // I * (B / D) mod (N / D)
   // To simplify the computation, we factor out the divide by D:
   // (I * B mod N) / D
   const SCEV *D = SE.getConstant(APInt::getOneBitSet(BW, Mult2));
   return SE.getUDivExactExpr(SE.getMulExpr(B, SE.getConstant(I)), D);
 }
 
-/// Find the roots of the quadratic equation for the given quadratic chrec
-/// {L,+,M,+,N}.  This returns either the two roots (which might be the same) or
-/// two SCEVCouldNotCompute objects.
-static Optional<std::pair<const SCEVConstant *,const SCEVConstant *>>
-SolveQuadraticEquation(const SCEVAddRecExpr *AddRec, ScalarEvolution &SE) {
+/// For a given quadratic addrec, generate coefficients of the corresponding
+/// quadratic equation, multiplied by a common value to ensure that they are
+/// integers.
+/// The returned value is a tuple { A, B, C, M, BitWidth }, where
+/// Ax^2 + Bx + C is the quadratic function, M is the value that A, B and C
+/// were multiplied by, and BitWidth is the bit width of the original addrec
+/// coefficients.
+/// This function returns None if the addrec coefficients are not compile-
+/// time constants.
+static Optional<std::tuple<APInt, APInt, APInt, APInt, unsigned>>
+GetQuadraticEquation(const SCEVAddRecExpr *AddRec) {
   assert(AddRec->getNumOperands() == 3 && "This is not a quadratic chrec!");
   const SCEVConstant *LC = dyn_cast<SCEVConstant>(AddRec->getOperand(0));
   const SCEVConstant *MC = dyn_cast<SCEVConstant>(AddRec->getOperand(1));
   const SCEVConstant *NC = dyn_cast<SCEVConstant>(AddRec->getOperand(2));
+  LLVM_DEBUG(dbgs() << __func__ << ": analyzing quadratic addrec: "
+                    << *AddRec << '\n');
 
   // We currently can only solve this if the coefficients are constants.
-  if (!LC || !MC || !NC)
+  if (!LC || !MC || !NC) {
+    LLVM_DEBUG(dbgs() << __func__ << ": coefficients are not constant\n");
+    return None;
+  }
+
+  APInt L = LC->getAPInt();
+  APInt M = MC->getAPInt();
+  APInt N = NC->getAPInt();
+  assert(!N.isNullValue() && "This is not a quadratic addrec");
+
+  unsigned BitWidth = LC->getAPInt().getBitWidth();
+  unsigned NewWidth = BitWidth + 1;
+  LLVM_DEBUG(dbgs() << __func__ << ": addrec coeff bw: "
+                    << BitWidth << '\n');
+  // The sign-extension (as opposed to a zero-extension) here matches the
+  // extension used in SolveQuadraticEquationWrap (with the same motivation).
+  N = N.sext(NewWidth);
+  M = M.sext(NewWidth);
+  L = L.sext(NewWidth);
+
+  // The increments are M, M+N, M+2N, ..., so the accumulated values are
+  //   L+M, (L+M)+(M+N), (L+M)+(M+N)+(M+2N), ..., that is,
+  //   L+M, L+2M+N, L+3M+3N, ...
+  // After n iterations the accumulated value Acc is L + nM + n(n-1)/2 N.
+  //
+  // The equation Acc = 0 is then
+  //   L + nM + n(n-1)/2 N = 0,  or  2L + 2M n + n(n-1) N = 0.
+  // In a quadratic form it becomes:
+  //   N n^2 + (2M-N) n + 2L = 0.
+
+  APInt A = N;
+  APInt B = 2 * M - A;
+  APInt C = 2 * L;
+  APInt T = APInt(NewWidth, 2);
+  LLVM_DEBUG(dbgs() << __func__ << ": equation " << A << "x^2 + " << B
+                    << "x + " << C << ", coeff bw: " << NewWidth
+                    << ", multiplied by " << T << '\n');
+  return std::make_tuple(A, B, C, T, BitWidth);
+}
+
+/// Helper function to compare optional APInts:
+/// (a) if X and Y both exist, return min(X, Y),
+/// (b) if neither X nor Y exist, return None,
+/// (c) if exactly one of X and Y exists, return that value.
+static Optional<APInt> MinOptional(Optional<APInt> X, Optional<APInt> Y) {
+  if (X.hasValue() && Y.hasValue()) {
+    unsigned W = std::max(X->getBitWidth(), Y->getBitWidth());
+    APInt XW = X->sextOrSelf(W);
+    APInt YW = Y->sextOrSelf(W);
+    return XW.slt(YW) ? *X : *Y;
+  }
+  if (!X.hasValue() && !Y.hasValue())
+    return None;
+  return X.hasValue() ? *X : *Y;
+}
+
+/// Helper function to truncate an optional APInt to a given BitWidth.
+/// When solving addrec-related equations, it is preferable to return a value
+/// that has the same bit width as the original addrec's coefficients. If the
+/// solution fits in the original bit width, truncate it (except for i1).
+/// Returning a value of a different bit width may inhibit some optimizations.
+///
+/// In general, a solution to a quadratic equation generated from an addrec
+/// may require BW+1 bits, where BW is the bit width of the addrec's
+/// coefficients. The reason is that the coefficients of the quadratic
+/// equation are BW+1 bits wide (to avoid truncation when converting from
+/// the addrec to the equation).
+static Optional<APInt> TruncIfPossible(Optional<APInt> X, unsigned BitWidth) {
+  if (!X.hasValue())
+    return None;
+  unsigned W = X->getBitWidth();
+  if (BitWidth > 1 && BitWidth < W && X->isIntN(BitWidth))
+    return X->trunc(BitWidth);
+  return X;
+}
+
+/// Let c(n) be the value of the quadratic chrec {L,+,M,+,N} after n
+/// iterations. The values L, M, N are assumed to be signed, and they
+/// should all have the same bit widths.
+/// Find the least n >= 0 such that c(n) = 0 in the arithmetic modulo 2^BW,
+/// where BW is the bit width of the addrec's coefficients.
+/// If the calculated value is a BW-bit integer (for BW > 1), it will be
+/// returned as such, otherwise the bit width of the returned value may
+/// be greater than BW.
+///
+/// This function returns None if
+/// (a) the addrec coefficients are not constant, or
+/// (b) SolveQuadraticEquationWrap was unable to find a solution. For cases
+///     like x^2 = 5, no integer solutions exist, in other cases an integer
+///     solution may exist, but SolveQuadraticEquationWrap may fail to find it.
+static Optional<APInt>
+SolveQuadraticAddRecExact(const SCEVAddRecExpr *AddRec, ScalarEvolution &SE) {
+  APInt A, B, C, M;
+  unsigned BitWidth;
+  auto T = GetQuadraticEquation(AddRec);
+  if (!T.hasValue())
     return None;
 
-  uint32_t BitWidth = LC->getAPInt().getBitWidth();
-  const APInt &L = LC->getAPInt();
-  const APInt &M = MC->getAPInt();
-  const APInt &N = NC->getAPInt();
-  APInt Two(BitWidth, 2);
-
-  // Convert from chrec coefficients to polynomial coefficients AX^2+BX+C
-
-  // The A coefficient is N/2
-  APInt A = N.sdiv(Two);
-
-  // The B coefficient is M-N/2
-  APInt B = M;
-  B -= A; // A is the same as N/2.
-
-  // The C coefficient is L.
-  const APInt& C = L;
-
-  // Compute the B^2-4ac term.
-  APInt SqrtTerm = B;
-  SqrtTerm *= B;
-  SqrtTerm -= 4 * (A * C);
+  std::tie(A, B, C, M, BitWidth) = *T;
+  LLVM_DEBUG(dbgs() << __func__ << ": solving for unsigned overflow\n");
+  Optional<APInt> X = APIntOps::SolveQuadraticEquationWrap(A, B, C, BitWidth+1);
+  if (!X.hasValue())
+    return None;
 
-  if (SqrtTerm.isNegative()) {
-    // The loop is provably infinite.
+  ConstantInt *CX = ConstantInt::get(SE.getContext(), *X);
+  ConstantInt *V = EvaluateConstantChrecAtConstant(AddRec, CX, SE);
+  if (!V->isZero())
     return None;
-  }
+
+  return TruncIfPossible(X, BitWidth);
 
-  // Compute sqrt(B^2-4ac). This is guaranteed to be the nearest
-  // integer value or else APInt::sqrt() will assert.
-  APInt SqrtVal = SqrtTerm.sqrt();
-
-  // Compute the two solutions for the quadratic formula.
-  // The divisions must be performed as signed divisions.
-  APInt NegB = -std::move(B);
-  APInt TwoA = std::move(A);
-  TwoA <<= 1;
-  if (TwoA.isNullValue())
+/// Let c(n) be the value of the quadratic chrec {0,+,M,+,N} after n
+/// iterations. The values M, N are assumed to be signed, and they
+/// should all have the same bit widths.
+/// Find the least n such that c(n) does not belong to the given range,
+/// while c(n-1) does.
+///
+/// This function returns None if
+/// (a) the addrec coefficients are not constant, or
+/// (b) SolveQuadraticEquationWrap was unable to find a solution for the
+///     bounds of the range.
+static Optional<APInt>
+SolveQuadraticAddRecRange(const SCEVAddRecExpr *AddRec,
+                          const ConstantRange &Range, ScalarEvolution &SE) {
+  assert(AddRec->getOperand(0)->isZero() &&
+         "Starting value of addrec should be 0");
+  LLVM_DEBUG(dbgs() << __func__ << ": solving boundary crossing for range "
+                    << Range << ", addrec " << *AddRec << '\n');
+  // This case is handled in getNumIterationsInRange. Here we can assume that
+  // we start in the range.
+  assert(Range.contains(APInt(SE.getTypeSizeInBits(AddRec->getType()), 0)) &&
+         "Addrec's initial value should be in range");
+
+  APInt A, B, C, M;
+  unsigned BitWidth;
+  auto T = GetQuadraticEquation(AddRec);
+  if (!T.hasValue())
     return None;
 
-  LLVMContext &Context = SE.getContext();
+  // Be careful about the return value: there can be two reasons for not
+  // returning an actual number. First, if no solutions to the equations
+  // were found, and second, if the solutions don't leave the given range.
+  // The first case means that the actual solution is "unknown", the second
+  // means that it's known, but not valid. If the solution is unknown, we
+  // cannot make any conclusions.
+  // Return a pair: the optional solution and a flag indicating if the
+  // solution was found.
+  auto SolveForBoundary = [&](APInt Bound) -> std::pair<Optional<APInt>,bool> {
+    // Solve for signed overflow and unsigned overflow, pick the lower
+    // solution.
+    LLVM_DEBUG(dbgs() << "SolveQuadraticAddRecRange: checking boundary "
+                      << Bound << " (before multiplying by " << M << ")\n");
+    Bound *= M; // The quadratic equation multiplier.
+
+    Optional<APInt> SO = None;
+    if (BitWidth > 1) {
+      LLVM_DEBUG(dbgs() << "SolveQuadraticAddRecRange: solving for "
+                           "signed overflow\n");
+      SO = APIntOps::SolveQuadraticEquationWrap(A, B, -Bound, BitWidth);
+    }
+    LLVM_DEBUG(dbgs() << "SolveQuadraticAddRecRange: solving for "
+                         "unsigned overflow\n");
+    Optional<APInt> UO = APIntOps::SolveQuadraticEquationWrap(A, B, -Bound,
+                                                              BitWidth+1);
+
+    auto LeavesRange = [&] (const APInt &X) {
+      ConstantInt *C0 = ConstantInt::get(SE.getContext(), X);
+      ConstantInt *V0 = EvaluateConstantChrecAtConstant(AddRec, C0, SE);
+      if (Range.contains(V0->getValue()))
+        return false;
+      // X should be at least 1, so X-1 is non-negative.
+      ConstantInt *C1 = ConstantInt::get(SE.getContext(), X-1);
+      ConstantInt *V1 = EvaluateConstantChrecAtConstant(AddRec, C1, SE);
+      if (Range.contains(V1->getValue()))
+        return true;
+      return false;
+    };
+
+    // If SolveQuadraticEquationWrap returns None, it means that there can
+    // be a solution, but the function failed to find it. We cannot treat it
+    // as "no solution".
+    if (!SO.hasValue() || !UO.hasValue())
+      return { None, false };
+
+    // Check the smaller value first to see if it leaves the range.
+    // At this point, both SO and UO must have values.
+    Optional<APInt> Min = MinOptional(SO, UO);
+    if (LeavesRange(*Min))
+      return { Min, true };
+    Optional<APInt> Max = Min == SO ? UO : SO;
+    if (LeavesRange(*Max))
+      return { Max, true };
+
+    // Solutions were found, but were eliminated, hence the "true".
+    return { None, true };
+  };
+
+  std::tie(A, B, C, M, BitWidth) = *T;
+  // Lower bound is inclusive, subtract 1 to represent the exiting value.
+  APInt Lower = Range.getLower().sextOrSelf(A.getBitWidth()) - 1;
+  APInt Upper = Range.getUpper().sextOrSelf(A.getBitWidth());
+  auto SL = SolveForBoundary(Lower);
+  auto SU = SolveForBoundary(Upper);
+  // If any of the solutions was unknown, no meaninigful conclusions can
+  // be made.
+  if (!SL.second || !SU.second)
+    return None;
 
-  ConstantInt *Solution1 =
-    ConstantInt::get(Context, (NegB + SqrtVal).sdiv(TwoA));
-  ConstantInt *Solution2 =
-    ConstantInt::get(Context, (NegB - SqrtVal).sdiv(TwoA));
+  // Claim: The correct solution is not some value between Min and Max.
+  //
+  // Justification: Assuming that Min and Max are different values, one of
+  // them is when the first signed overflow happens, the other is when the
+  // first unsigned overflow happens. Crossing the range boundary is only
+  // possible via an overflow (treating 0 as a special case of it, modeling
+  // an overflow as crossing k*2^W for some k).
+  //
+  // The interesting case here is when Min was eliminated as an invalid
+  // solution, but Max was not. The argument is that if there was another
+  // overflow between Min and Max, it would also have been eliminated if
+  // it was considered.
+  //
+  // For a given boundary, it is possible to have two overflows of the same
+  // type (signed/unsigned) without having the other type in between: this
+  // can happen when the vertex of the parabola is between the iterations
+  // corresponding to the overflows. This is only possible when the two
+  // overflows cross k*2^W for the same k. In such case, if the second one
+  // left the range (and was the first one to do so), the first overflow
+  // would have to enter the range, which would mean that either we had left
+  // the range before or that we started outside of it. Both of these cases
+  // are contradictions.
+  //
+  // Claim: In the case where SolveForBoundary returns None, the correct
+  // solution is not some value between the Max for this boundary and the
+  // Min of the other boundary?
+  //
+  // Justification: Assume that we had such Max_A and Min_B corresponding
+  // to range boundaries A and B and such that Max_A < Min_B. If there was
+  // a solution between Max_A and Min_B, it would have to be caused by an
+  // overflow corresponding to either A or B. It cannot correspond to B,
+  // since Min_B is the first occurrence of such an overflow. If it
+  // corresponded to A, it would have to be either a signed or an unsigned
+  // overflow that is larger than both eliminated overflows for A. But
+  // between the eliminated overflows and this overflow, the values would
+  // cover the entire value space, thus crossing the other boundary, which
+  // is a contradiction.
 
-  return std::make_pair(cast<SCEVConstant>(SE.getConstant(Solution1)),
-                        cast<SCEVConstant>(SE.getConstant(Solution2)));
+  return TruncIfPossible(MinOptional(SL.first, SU.first), BitWidth);
 }
 
 ScalarEvolution::ExitLimit
 ScalarEvolution::howFarToZero(const SCEV *V, const Loop *L, bool ControlsExit,
                               bool AllowPredicates) {
 
   // This is only used for loops with a "x != y" exit test. The exit condition
   // is now expressed as a single expression, V = x-y. So the exit test is
if (!AddRec && AllowPredicates)
     AddRec = convertSCEVToAddRecWithPredicates(V, L, Predicates);
 
   if (!AddRec || AddRec->getLoop() != L)
     return getCouldNotCompute();
 
   // If this is a quadratic (3-term) AddRec {L,+,M,+,N}, find the roots of
   // the quadratic equation to solve it.
   if (AddRec->isQuadratic() && AddRec->getType()->isIntegerTy()) {
-    if (auto Roots = SolveQuadraticEquation(AddRec, *this)) {
-      const SCEVConstant *R1 = Roots->first;
-      const SCEVConstant *R2 = Roots->second;
-      // Pick the smallest positive root value.
-      if (ConstantInt *CB = dyn_cast<ConstantInt>(ConstantExpr::getICmp(
-              CmpInst::ICMP_ULT, R1->getValue(), R2->getValue()))) {
-        if (!CB->getZExtValue())
-          std::swap(R1, R2); // R1 is the minimum root now.
-
-        // We can only use this value if the chrec ends up with an exact zero
-        // value at this index.  When solving for "X*X != 5", for example, we
-        // should not accept a root of 2.
-        const SCEV *Val = AddRec->evaluateAtIteration(R1, *this);
-        if (Val->isZero())
-          // We found a quadratic root!
-          return ExitLimit(R1, R1, false, Predicates);
-      }
+    if (auto S = SolveQuadraticAddRecExact(AddRec, *this)) {
+      const auto *R = cast<SCEVConstant>(getConstant(S.getValue()));
+      return ExitLimit(R, R, false, Predicates);
     }
     return getCouldNotCompute();
   }
 
   // Otherwise we can only handle this if it is affine.
   if (!AddRec->isAffine())
     return getCouldNotCompute();
 
if (isAffine()) {
     if (Range.contains(Val->getValue()))
       return SE.getCouldNotCompute();  // Something strange happened
 
     // Ensure that the previous value is in the range.  This is a sanity check.
     assert(Range.contains(
            EvaluateConstantChrecAtConstant(this,
            ConstantInt::get(SE.getContext(), ExitVal - 1), SE)->getValue()) &&
            "Linear scev computation is off in a bad way!");
     return SE.getConstant(ExitValue);
-  } else if (isQuadratic()) {
-    // If this is a quadratic (3-term) AddRec {L,+,M,+,N}, find the roots of the
-    // quadratic equation to solve it.  To do this, we must frame our problem in
-    // terms of figuring out when zero is crossed, instead of when
-    // Range.getUpper() is crossed.
-    SmallVector<const SCEV *, 4> NewOps(op_begin(), op_end());
-    NewOps[0] = SE.getNegativeSCEV(SE.getConstant(Range.getUpper()));
-    const SCEV *NewAddRec = SE.getAddRecExpr(NewOps, getLoop(), FlagAnyWrap);
-
-    // Next, solve the constructed addrec
-    if (auto Roots =
-            SolveQuadraticEquation(cast<SCEVAddRecExpr>(NewAddRec), SE)) {
-      const SCEVConstant *R1 = Roots->first;
-      const SCEVConstant *R2 = Roots->second;
-      // Pick the smallest positive root value.
-      if (ConstantInt *CB = dyn_cast<ConstantInt>(ConstantExpr::getICmp(
-              ICmpInst::ICMP_ULT, R1->getValue(), R2->getValue()))) {
-        if (!CB->getZExtValue())
-          std::swap(R1, R2); // R1 is the minimum root now.
-
-        // Make sure the root is not off by one.  The returned iteration should
-        // not be in the range, but the previous one should be.  When solving
-        // for "X*X < 5", for example, we should not return a root of 2.
-        ConstantInt *R1Val =
-            EvaluateConstantChrecAtConstant(this, R1->getValue(), SE);
-        if (Range.contains(R1Val->getValue())) {
-          // The next iteration must be out of the range...
-          ConstantInt *NextVal =
-              ConstantInt::get(SE.getContext(), R1->getAPInt() + 1);
-
-          R1Val = EvaluateConstantChrecAtConstant(this, NextVal, SE);
-          if (!Range.contains(R1Val->getValue()))
-            return SE.getConstant(NextVal);
-          return SE.getCouldNotCompute(); // Something strange happened
-        }
 
-        // If R1 was not in the range, then it is a good return value.  Make
-        // sure that R1-1 WAS in the range though, just in case.
-        ConstantInt *NextVal =
-            ConstantInt::get(SE.getContext(), R1->getAPInt() - 1);
-        R1Val = EvaluateConstantChrecAtConstant(this, NextVal, SE);
-        if (Range.contains(R1Val->getValue()))
-          return R1;
-        return SE.getCouldNotCompute(); // Something strange happened
-      }
-    }
+  if (isQuadratic()) {
+    if (auto S = SolveQuadraticAddRecRange(this, Range, SE))
+      return SE.getConstant(S.getValue());
   }
 
   return SE.getCouldNotCompute();
 }
 
 const SCEVAddRecExpr *
 SCEVAddRecExpr::getPostIncExpr(ScalarEvolution &SE) const {
   assert(getNumOperands() > 1 && "AddRec with zero step?");
   // There is a temptation to just call getAddExpr(this, getStepRecurrence(SE)),
   // but in this case we cannot guarantee that the value returned will be an
   // AddRec because SCEV does not have a fixed point where it stops
   // simplification: it is legal to return ({rec1} + {rec2}). For example, it
   // may happen if we reach arithmetic depth limit while simplifying. So we
   // construct the returned value explicitly.
   SmallVector<const SCEV *, 3> Ops;
   // If this is {A,+,B,+,C,...,+,N}, then its step is {B,+,C,+,...,+,N}, and
   // (this + Step) is {A+B,+,B+C,+...,+,N}.
   for (unsigned i = 0, e = getNumOperands() - 1; i < e; ++i)
     Ops.push_back(SE.getAddExpr(getOperand(i), getOperand(i + 1)));
   // We know that the last operand is not a constant zero (otherwise it would
   // have been popped out earlier). This guarantees us that if the result has
   // the same last operand, then it will also not be popped out, meaning that
   // the returned value will be an AddRec.
   const SCEV *Last = getOperand(getNumOperands() - 1);
   assert(!Last->isZero() && "Recurrency with zero step?");
   Ops.push_back(Last);
   return cast<SCEVAddRecExpr>(SE.getAddRecExpr(Ops, getLoop(),

lib/Support/APInt.cpp

 // constant values and provide a variety of arithmetic operations on them.
 //
 //===----------------------------------------------------------------------===//
 
 #include "llvm/ADT/APInt.h"
 #include "llvm/ADT/ArrayRef.h"
 #include "llvm/ADT/FoldingSet.h"
 #include "llvm/ADT/Hashing.h"
+#include "llvm/ADT/Optional.h"
 #include "llvm/ADT/SmallString.h"
 #include "llvm/ADT/StringRef.h"
 #include "llvm/Config/llvm-config.h"
 #include "llvm/Support/Debug.h"
 #include "llvm/Support/ErrorHandling.h"
 #include "llvm/Support/MathExtras.h"
 #include "llvm/Support/raw_ostream.h"
 #include <climits>
case APInt::Rounding::UP: {
     return Quo + 1;
   }
   // Currently sdiv rounds twards zero.
   case APInt::Rounding::TOWARD_ZERO:
     return A.sdiv(B);
   }
   llvm_unreachable("Unknown APInt::Rounding enum");
 }
+
+Optional<APInt>
+llvm::APIntOps::SolveQuadraticEquationWrap(APInt A, APInt B, APInt C,
+                                           unsigned RangeWidth) {
+  unsigned CoeffWidth = A.getBitWidth();
+  assert(CoeffWidth == B.getBitWidth() && CoeffWidth == C.getBitWidth());
+  assert(RangeWidth <= CoeffWidth &&
+         "Value range width should be less than coefficient width");
+  assert(RangeWidth > 1 && "Value range bit width should be > 1");
+
+  LLVM_DEBUG(dbgs() << __func__ << ": solving " << A << "x^2 + " << B
+                    << "x + " << C << ", rw:" << RangeWidth << '\n');
+
+  // Identify 0 as a (non)solution immediately.
+  if (C.sextOrTrunc(RangeWidth).isNullValue() ) {
+    LLVM_DEBUG(dbgs() << __func__ << ": zero solution\n");
+    return APInt(CoeffWidth, 0);
+  }
+
+  // The result of APInt arithmetic has the same bit width as the operands,
+  // so it can actually lose high bits. A product of two n-bit integers needs
+  // 2n-1 bits to represent the full value.
+  // The operation done below (on quadratic coefficients) that can produce
+  // the largest value is the evaluation of the equation during bisection,
+  // which needs 3 times the bitwidth of the coefficient, so the total number
+  // of required bits is 3n.
+  //
+  // The purpose of this extension is to simulate the set Z of all integers,
+  // where n+1 > n for all n in Z. In Z it makes sense to talk about positive
+  // and negative numbers (not so much in a modulo arithmetic). The method
+  // used to solve the equation is based on the standard formula for real
+  // numbers, and uses the concepts of "positive" and "negative" with their
+  // usual meanings.
+  CoeffWidth *= 3;
+  A = A.sext(CoeffWidth);
+  B = B.sext(CoeffWidth);
+  C = C.sext(CoeffWidth);
+
+  // Make A > 0 for simplicity. Negate cannot overflow at this point because
+  // the bit width has increased.
+  if (A.isNegative()) {
+    A.negate();
+    B.negate();
+    C.negate();
+  }
+
+  // Solving an equation q(x) = 0 with coefficients in modular arithmetic
+  // is really solving a set of equations q(x) = kR for k = 0, 1, 2, ...,
+  // and R = 2^BitWidth.
+  // Since we're trying not only to find exact solutions, but also values
+  // that "wrap around", such a set will always have a solution, i.e. an x
+  // that satisfies at least one of the equations, or such that |q(x)|
+  // exceeds kR, while |q(x-1)| for the same k does not.
+  //
+  // We need to find a value k, such that Ax^2 + Bx + C = kR will have a
+  // positive solution n (in the above sense), and also such that the n
+  // will be the least among all solutions corresponding to k = 0, 1, ...
+  // (more precisely, the least element in the set
+  //   { n(k) | k is such that a solution n(k) exists }).
+  //
+  // Consider the parabola (over real numbers) that corresponds to the
+  // quadratic equation. Since A > 0, the arms of the parabola will point
+  // up. Picking different values of k will shift it up and down by R.
+  //
+  // We want to shift the parabola in such a way as to reduce the problem
+  // of solving q(x) = kR to solving shifted_q(x) = 0.
+  // (The interesting solutions are the ceilings of the real number
+  // solutions.)
+  APInt R = APInt::getOneBitSet(CoeffWidth, RangeWidth);
+  APInt TwoA = 2 * A;
+  APInt SqrB = B * B;
+  bool PickLow;
+
+  auto RoundUp = [] (const APInt &V, const APInt &A) {
+    assert(A.isStrictlyPositive());
+    APInt T = V.abs().urem(A);
+    if (T.isNullValue())
+      return V;
+    return V.isNegative() ? V+T : V+(A-T);
+  };
+
+  // The vertex of the parabola is at -B/2A, but since A > 0, it's negative
+  // iff B is positive.
+  if (B.isNonNegative()) {
+    // If B >= 0, the vertex it at a negative location (or at 0), so in
+    // order to have a non-negative solution we need to pick k that makes
+    // C-kR negative. To satisfy all the requirements for the solution
+    // that we are looking for, it needs to be closest to 0 of all k.
+    C = C.srem(R);
+    if (C.isStrictlyPositive())
+      C -= R;
+    // Pick the greater solution.
+    PickLow = false;
+  } else {
+    // If B < 0, the vertex is at a positive location. For any solution
+    // to exist, the discriminant must be non-negative. This means that
+    // C-kR <= B^2/4A is a necessary condition for k, i.e. there is a
+    // lower bound on values of k: kR >= C - B^2/4A.
+    APInt LowkR = C - SqrB.udiv(2*TwoA); // udiv because all values > 0.
+    // Round LowkR up (towards +inf) to the nearest kR.
+    LowkR = RoundUp(LowkR, R);
+
+    // If there exists k meeting the condition above, and such that
+    // C-kR > 0, there will be two positive real number solutions of
+    // q(x) = kR. Out of all such values of k, pick the one that makes
+    // C-kR closest to 0, (i.e. pick maximum k such that C-kR > 0).
+    // In other words, find maximum k such that LowkR <= kR < C.
+    if (C.sgt(LowkR)) {
+      // If LowkR < C, then such a k is guaranteed to exist because
+      // LowkR itself is a multiple of R.
+      C -= -RoundUp(-C, R);      // C = C - RoundDown(C, R)
+      // Pick the smaller solution.
+      PickLow = true;
+    } else {
+      // If C-kR < 0 for all potential k's, it means that one solution
+      // will be negative, while the other will be positive. The positive
+      // solution will shift towards 0 if the parabola is moved up.
+      // Pick the kR closest to the lower bound (i.e. make C-kR closest
+      // to 0, or in other words, out of all parabolas that have solutions,
+      // pick the one that is the farthest "up").
+      // Since LowkR is itself a multiple of R, simply take C-LowkR.
+      C -= LowkR;
+      // Pick the greater solution.
+      PickLow = false;
+    }
+  }
+
+  LLVM_DEBUG(dbgs() << __func__ << ": updated coefficients " << A << "x^2 + "
+                    << B << "x + " << C << ", rw:" << RangeWidth << '\n');
+
+  APInt D = SqrB - 4*A*C;
+  assert(D.isNonNegative() && "Negative discriminant");
+  APInt SQ = D.sqrt();
+
+  APInt Q = SQ * SQ;
+  bool InexactSQ = Q != D;
+  // The calculated SQ may actually be greater than the exact (non-integer)
+  // value. If that's the case, decremement SQ to get a value that is lower.
+  if (Q.sgt(D))
+    SQ -= 1;
+
+  APInt X;
+  APInt Rem;
+
+  // SQ is rounded down (i.e SQ * SQ <= D), so the roots may be inexact.
+  // When using the quadratic formula directly, the calculated low root
+  // may be greater than the exact one, since we would be subtracting SQ.
+  // To make sure that the calculated root is not greater than the exact
+  // one, subtract SQ+1 when calculating the low root (for inexact value
+  // of SQ).
+  if (PickLow)
+    APInt::sdivrem(-B - (SQ+InexactSQ), TwoA, X, Rem);
+  else
+    APInt::sdivrem(-B + SQ, TwoA, X, Rem);
+
+  // The updated coefficients should be such that the (exact) solution is
+  // positive. Since APInt division rounds towards 0, the calculated one
+  // can be 0, but cannot be negative.
+  assert(X.isNonNegative() && "Solution should be non-negative");
+
+  if (!InexactSQ && Rem.isNullValue()) {
+    LLVM_DEBUG(dbgs() << __func__ << ": solution (root): " << X << '\n');
+    return X;
+  }
+
+  assert((SQ*SQ).sle(D) && "SQ = |_sqrt(D)_|, so SQ*SQ <= D");
+  // The exact value of the square root of D should be between SQ and SQ+1.
+  // This implies that the solution should be between that corresponding to
+  // SQ (i.e. X) and that corresponding to SQ+1.
+  //
+  // The calculated X cannot be greater than the exact (real) solution.
+  // Actually it must be strictly less than the exact solution, while
+  // X+1 will be greater than or equal to it.
+
+  APInt VX = (A*X + B)*X + C;
+  APInt VY = VX + TwoA*X + A + B;
+  bool SignChange = VX.isNegative() != VY.isNegative() ||
+                    VX.isNullValue() != VY.isNullValue();
+  // If the sign did not change between X and X+1, X is not a valid solution.
+  // This could happen when the actual (exact) roots don't have an integer
+  // between them, so they would both be contained between X and X+1.
+  if (!SignChange) {
+    LLVM_DEBUG(dbgs() << __func__ << ": no valid solution\n");
+    return None;
+  }
+
+  X += 1;
+  LLVM_DEBUG(dbgs() << __func__ << ": solution (wrap): " << X << '\n');
+  return X;
+}

test/Analysis/ScalarEvolution/solve-quadratic-i1.ll

+; RUN: opt -analyze -scalar-evolution < %s | FileCheck %s
+
+target triple = "x86_64-unknown-linux-gnu"
+
+; CHECK-LABEL: Printing analysis 'Scalar Evolution Analysis' for function 'f0':
+; CHECK-NEXT: Classifying expressions for: @f0
+; CHECK-NEXT:   %v0 = phi i16 [ 2, %b0 ], [ %v2, %b1 ]
+; CHECK-NEXT:   -->  {2,+,1}<nuw><nsw><%b1> U: [2,4) S: [2,4)         Exits: 3                LoopDispositions: { %b1: Computable }
+; CHECK-NEXT:   %v1 = phi i16 [ 1, %b0 ], [ %v3, %b1 ]
+; CHECK-NEXT:   -->  {1,+,2,+,1}<%b1> U: full-set S: full-set         Exits: 3                LoopDispositions: { %b1: Computable }
+; CHECK-NEXT:   %v2 = add nsw i16 %v0, 1
+; CHECK-NEXT:   -->  {3,+,1}<nuw><nsw><%b1> U: [3,5) S: [3,5)         Exits: 4                LoopDispositions: { %b1: Computable }
+; CHECK-NEXT:   %v3 = add nsw i16 %v1, %v0
+; CHECK-NEXT:   -->  {3,+,3,+,1}<%b1> U: full-set S: full-set         Exits: 6                LoopDispositions: { %b1: Computable }
+; CHECK-NEXT:   %v4 = and i16 %v3, 1
+; CHECK-NEXT:   -->  (zext i1 {true,+,true,+,true}<%b1> to i16) U: [0,2) S: [0,2)             Exits: 0                LoopDispositions: { %b1: Computable }
+; CHECK-NEXT: Determining loop execution counts for: @f0
+; CHECK-NEXT: Loop %b1: backedge-taken count is 1
+; CHECK-NEXT: Loop %b1: max backedge-taken count is 1
+; CHECK-NEXT: Loop %b1: Predicated backedge-taken count is 1
+; CHECK-NEXT:  Predicates:
+; CHECK-EMPTY:
+; CHECK-NEXT: Loop %b1: Trip multiple is 2
+define void @f0() {
+b0:
+  br label %b1
+
+b1:                                               ; preds = %b1, %b0
+  %v0 = phi i16 [ 2, %b0 ], [ %v2, %b1 ]
+  %v1 = phi i16 [ 1, %b0 ], [ %v3, %b1 ]
+  %v2 = add nsw i16 %v0, 1
+  %v3 = add nsw i16 %v1, %v0
+  %v4 = and i16 %v3, 1
+  %v5 = icmp ne i16 %v4, 0
+  br i1 %v5, label %b1, label %b2
+
+b2:                                               ; preds = %b1
+  ret void
+}
+
+@g0 = common dso_local global i16 0, align 2
+@g1 = common dso_local global i32 0, align 4
+@g2 = common dso_local global i32* null, align 8
+
+; CHECK-LABEL: Printing analysis 'Scalar Evolution Analysis' for function 'f1':
+; CHECK-NEXT: Classifying expressions for: @f1
+; CHECK-NEXT:   %v0 = phi i16 [ 0, %b0 ], [ %v3, %b1 ]
+; CHECK-NEXT:   -->  {0,+,3,+,1}<%b1> U: full-set S: full-set         Exits: 7                LoopDispositions: { %b1: Computable }
+; CHECK-NEXT:   %v1 = phi i32 [ 3, %b0 ], [ %v6, %b1 ]
+; CHECK-NEXT:   -->  {3,+,1}<nuw><nsw><%b1> U: [3,6) S: [3,6)         Exits: 5                LoopDispositions: { %b1: Computable }
+; CHECK-NEXT:   %v2 = trunc i32 %v1 to i16
+; CHECK-NEXT:   -->  {3,+,1}<%b1> U: [3,6) S: [3,6)           Exits: 5                LoopDispositions: { %b1: Computable }
+; CHECK-NEXT:   %v3 = add i16 %v0, %v2
+; CHECK-NEXT:   -->  {3,+,4,+,1}<%b1> U: full-set S: full-set         Exits: 12               LoopDispositions: { %b1: Computable }
+; CHECK-NEXT:   %v4 = and i16 %v3, 1
+; CHECK-NEXT:   -->  (zext i1 {true,+,false,+,true}<%b1> to i16) U: [0,2) S: [0,2)            Exits: 0                LoopDispositions: { %b1: Computable }
+; CHECK-NEXT:   %v6 = add nuw nsw i32 %v1, 1
+; CHECK-NEXT:   -->  {4,+,1}<nuw><nsw><%b1> U: [4,7) S: [4,7)         Exits: 6                LoopDispositions: { %b1: Computable }
+; CHECK-NEXT:   %v7 = phi i32 [ %v1, %b1 ]
+; CHECK-NEXT:   -->  %v7 U: [3,6) S: [3,6)
+; CHECK-NEXT:   %v8 = phi i16 [ %v3, %b1 ]
+; CHECK-NEXT:   -->  %v8 U: full-set S: full-set
+; CHECK-NEXT: Determining loop execution counts for: @f1
+; CHECK-NEXT: Loop %b3: <multiple exits> Unpredictable backedge-taken count.
+; CHECK-NEXT: Loop %b3: Unpredictable max backedge-taken count.
+; CHECK-NEXT: Loop %b3: Unpredictable predicated backedge-taken count.
+; CHECK-NEXT: Loop %b1: backedge-taken count is 2
+; CHECK-NEXT: Loop %b1: max backedge-taken count is 2
+; CHECK-NEXT: Loop %b1: Predicated backedge-taken count is 2
+; CHECK-NEXT:  Predicates:
+; CHECK-EMPTY:
+; CHECK-NEXT: Loop %b1: Trip multiple is 3
+define void @f1() #0 {
+b0:
+  store i16 0, i16* @g0, align 2
+  store i32* @g1, i32** @g2, align 8
+  br label %b1
+
+b1:                                               ; preds = %b1, %b0
+  %v0 = phi i16 [ 0, %b0 ], [ %v3, %b1 ]
+  %v1 = phi i32 [ 3, %b0 ], [ %v6, %b1 ]
+  %v2 = trunc i32 %v1 to i16
+  %v3 = add i16 %v0, %v2
+  %v4 = and i16 %v3, 1
+  %v5 = icmp eq i16 %v4, 0
+  %v6 = add nuw nsw i32 %v1, 1
+  br i1 %v5, label %b2, label %b1
+
+b2:                                               ; preds = %b1
+  %v7 = phi i32 [ %v1, %b1 ]
+  %v8 = phi i16 [ %v3, %b1 ]
+  store i32 %v7, i32* @g1, align 4
+  store i16 %v8, i16* @g0, align 2
+  br label %b3
+
+b3:                                               ; preds = %b3, %b2
+  br label %b3
+}
+
+attributes #0 = { nounwind uwtable "target-cpu"="x86-64" }

test/Analysis/ScalarEvolution/solve-quadratic-overflow.ll

+; RUN: opt -analyze -scalar-evolution -S < %s | FileCheck %s
+
+; The exit value from this loop was originally calculated as 0.
+; The actual exit condition is 256*256 == 0 (in i16).
+
+; CHECK: Printing analysis 'Scalar Evolution Analysis' for function 'f0':
+; CHECK-NEXT: Classifying expressions for: @f0
+; CHECK-NEXT:   %v1 = phi i16 [ 0, %b0 ], [ %v2, %b1 ]
+; CHECK-NEXT:   -->  {0,+,-1}<%b1> U: [-255,1) S: [-255,1)            Exits: -255            LoopDispositions: { %b1: Computable }
+; CHECK-NEXT:   %v2 = add i16 %v1, -1
+; CHECK-NEXT:   -->  {-1,+,-1}<%b1> U: [-256,0) S: [-256,0)           Exits: -256            LoopDispositions: { %b1: Computable }
+; CHECK-NEXT:   %v3 = mul i16 %v2, %v2
+; CHECK-NEXT:   -->  {1,+,3,+,2}<%b1> U: full-set S: full-set         Exits: 0               LoopDispositions: { %b1: Computable }
+; CHECK-NEXT:   %v5 = phi i16 [ %v2, %b1 ]
+; CHECK-NEXT:   -->  %v5 U: [-256,0) S: [-256,0)
+; CHECK-NEXT:   %v6 = phi i16 [ %v3, %b1 ]
+; CHECK-NEXT:   -->  %v6 U: full-set S: full-set
+; CHECK-NEXT:   %v7 = sext i16 %v5 to i32
+; CHECK-NEXT:   -->  (sext i16 %v5 to i32) U: [-256,0) S: [-256,0)
+; CHECK-NEXT: Determining loop execution counts for: @f0
+; CHECK-NEXT: Loop %b1: backedge-taken count is 255
+; CHECK-NEXT: Loop %b1: max backedge-taken count is 255
+; CHECK-NEXT: Loop %b1: Predicated backedge-taken count is 255
+; CHECK-NEXT:  Predicates:
+; CHECK-EMPTY:
+; CHECK-NEXT: Loop %b1: Trip multiple is 256
+
+
+@g0 = global i32 0, align 4
+@g1 = global i16 0, align 2
+
+define signext i32 @f0() {
+b0:
+  br label %b1
+
+b1:                                               ; preds = %b1, %b0
+  %v1 = phi i16 [ 0, %b0 ], [ %v2, %b1 ]
+  %v2 = add i16 %v1, -1
+  %v3 = mul i16 %v2, %v2
+  %v4 = icmp eq i16 %v3, 0
+  br i1 %v4, label %b2, label %b1
+
+b2:                                               ; preds = %b1
+  %v5 = phi i16 [ %v2, %b1 ]
+  %v6 = phi i16 [ %v3, %b1 ]
+  %v7 = sext i16 %v5 to i32
+  store i32 %v7, i32* @g0, align 4
+  store i16 %v6, i16* @g1, align 2
+  ret i32 0
+}
+

test/Analysis/ScalarEvolution/solve-quadratic.ll

+; RUN: opt -analyze -scalar-evolution -S -debug-only=scalar-evolution,apint < %s 2>&1 | FileCheck %s
+; REQUIRES: asserts
+
+; Use the following template to get a chrec {L,+,M,+,N}.
+;
+; define signext i32 @func() {
+; entry:
+;   br label %loop
+;
+; loop:
+;   %ivr = phi i32 [ 0, %entry ], [ %ivr1, %loop ]
+;   %inc = phi i32 [ X, %entry ], [ %inc1, %loop ]
+;   %acc = phi i32 [ Y, %entry ], [ %acc1, %loop ]
+;   %ivr1 = add i32 %ivr, %inc
+;   %inc1 = add i32 %inc, Z                 ; M = inc1 = inc + N = X + N
+;   %acc1 = add i32 %acc, %inc              ; L = acc1 = X + Y
+;   %and  = and i32 %acc1, 2^W-1            ; iW
+;   %cond = icmp eq i32 %and, 0
+;   br i1 %cond, label %exit, label %loop
+;
+; exit:
+;   %rv = phi i32 [ %acc1, %loop ]
+;   ret i32 %rv
+; }
+;
+; From
+;       X + Y = L
+;       X + Z = M
+;           Z = N
+; get
+;       X = M - N
+;       Y = N - M + L
+;       Z = N
+
+; The connection between the chrec coefficients {L,+,M,+,N} and the quadratic
+; coefficients is that the quadratic equation is N x^2 + (2M-N) x + 2L = 0,
+; where the equation was multiplied by 2 to make the coefficient at x^2 an
+; integer (the actual equation is N/2 x^2 + (M-N/2) x + L = 0).
+
+; Quadratic equation: 2x^2 + 2x + 4 in i4, solution (wrap): 4
+; {14,+,14,+,14} -> X=0, Y=14, Z=14
+;
+; CHECK-LABEL: Printing analysis 'Scalar Evolution Analysis' for function 'test01'
+; CHECK: GetQuadraticEquation: analyzing quadratic addrec: {-2,+,-2,+,-2}<%loop>
+; CHECK: GetQuadraticEquation: addrec coeff bw: 4
+; CHECK: GetQuadraticEquation: equation -2x^2 + -2x + -4, coeff bw: 5, multiplied by 2
+; CHECK: SolveQuadraticAddRecExact: solving for unsigned overflow
+; CHECK: SolveQuadraticEquationWrap: solving -2x^2 + -2x + -4, rw:5
+; CHECK: SolveQuadraticEquationWrap: updated coefficients 2x^2 + 2x + -28, rw:5
+; CHECK: SolveQuadraticEquationWrap: solution (wrap): 4
+; CHECK: Loop %loop: Unpredictable backedge-taken count
+define signext i32 @test01() {
+entry:
+  br label %loop
+
+loop:
+  %ivr = phi i32 [  0, %entry ], [ %ivr1, %loop ]
+  %inc = phi i32 [  0, %entry ], [ %inc1, %loop ]
+  %acc = phi i32 [ 14, %entry ], [ %acc1, %loop ]
+  %ivr1 = add i32 %ivr, %inc
+  %inc1 = add i32 %inc, 14
+  %acc1 = add i32 %acc, %inc
+  %and  = and i32 %acc1, 15
+  %cond = icmp eq i32 %and, 0
+  br i1 %cond, label %exit, label %loop
+
+exit:
+  %rv = phi i32 [ %acc1, %loop ]
+  ret i32 %rv
+}
+
+; Quadratic equation: 1x^2 + -73x + -146 in i32, solution (wrap): 75
+; {-72,+,-36,+,1} -> X=-37, Y=-35, Z=1
+;
+; CHECK-LABEL: Printing analysis 'Scalar Evolution Analysis' for function 'test02':
+; CHECK: GetQuadraticEquation: analyzing quadratic addrec: {0,+,-36,+,1}<%loop>
+; CHECK: GetQuadraticEquation: addrec coeff bw: 32
+; CHECK: GetQuadraticEquation: equation 1x^2 + -73x + 0, coeff bw: 33, multiplied by 2
+; CHECK: SolveQuadraticAddRecRange: solving for signed overflow
+; CHECK: SolveQuadraticEquationWrap: solving 1x^2 + -73x + 4294967154, rw:32
+; CHECK: SolveQuadraticEquationWrap: updated coefficients 1x^2 + -73x + -142, rw:32
+; CHECK: SolveQuadraticEquationWrap: solution (wrap): 75
+; CHECK: SolveQuadraticAddRecRange: solving for unsigned overflow
+; CHECK: SolveQuadraticEquationWrap: solving 1x^2 + -73x + 4294967154, rw:33
+; CHECK: SolveQuadraticEquationWrap: updated coefficients 1x^2 + -73x + -4294967438, rw:33
+; CHECK: SolveQuadraticEquationWrap: solution (wrap): 65573
+; CHECK: SolveQuadraticAddRecRange: solving for signed overflow
+; CHECK: SolveQuadraticEquationWrap: solving 1x^2 + -73x + -146, rw:32
+; CHECK: SolveQuadraticEquationWrap: updated coefficients 1x^2 + -73x + -146, rw:32
+; CHECK: SolveQuadraticEquationWrap: solution (wrap): 75
+; CHECK: SolveQuadraticAddRecRange: solving for unsigned overflow
+; CHECK: SolveQuadraticEquationWrap: solving 1x^2 + -73x + -146, rw:33
+; CHECK: SolveQuadraticEquationWrap: updated coefficients 1x^2 + -73x + -146, rw:33
+; CHECK: SolveQuadraticEquationWrap: solution (wrap): 75
+; CHECK: Loop %loop: backedge-taken count is 75
+define signext i32 @test02() {
+entry:
+  br label %loop
+
+loop:
+  %ivr = phi i32 [  0, %entry ], [ %ivr1, %loop ]
+  %inc = phi i32 [ -37, %entry ], [ %inc1, %loop ]
+  %acc = phi i32 [ -35, %entry ], [ %acc1, %loop ]
+  %ivr1 = add i32 %ivr, %inc
+  %inc1 = add i32 %inc, 1
+  %acc1 = add i32 %acc, %inc
+  %and  = and i32 %acc1, -1
+  %cond = icmp sgt i32 %and, 0
+  br i1 %cond, label %exit, label %loop
+
+exit:
+  %rv = phi i32 [ %acc1, %loop ]
+  ret i32 %rv
+}
+
+; Quadratic equation: 2x^2 - 4x + 34 in i4, solution (exact): 1.
+; {17,+,-1,+,2} -> X=-3, Y=20, Z=2
+;
+; CHECK-LABEL: Printing analysis 'Scalar Evolution Analysis' for function 'test03':
+; CHECK: GetQuadraticEquation: analyzing quadratic addrec: {1,+,-1,+,2}<%loop>
+; CHECK: GetQuadraticEquation: addrec coeff bw: 4
+; CHECK: GetQuadraticEquation: equation 2x^2 + -4x + 2, coeff bw: 5, multiplied by 2
+; CHECK: SolveQuadraticAddRecExact: solving for unsigned overflow
+; CHECK: SolveQuadraticEquationWrap: solving 2x^2 + -4x + 2, rw:5
+; CHECK: SolveQuadraticEquationWrap: updated coefficients 2x^2 + -4x + 2, rw:5
+; CHECK: SolveQuadraticEquationWrap: solution (root): 1
+; CHECK: Loop %loop: backedge-taken count is 1
+define signext i32 @test03() {
+entry:
+  br label %loop
+
+loop:
+  %ivr = phi i32 [  0, %entry ], [ %ivr1, %loop ]
+  %inc = phi i32 [ -3, %entry ], [ %inc1, %loop ]
+  %acc = phi i32 [ 20, %entry ], [ %acc1, %loop ]
+  %ivr1 = add i32 %ivr, %inc
+  %inc1 = add i32 %inc, 2
+  %acc1 = add i32 %acc, %inc
+  %and  = and i32 %acc1, 15
+  %cond = icmp eq i32 %and, 0
+  br i1 %cond, label %exit, label %loop
+
+exit:
+  %rv = phi i32 [ %acc1, %loop ]
+  ret i32 %rv
+}
+
+; Quadratic equation  4x^2 + 2x + 2 in i16, solution (wrap): 181
+; {1,+,3,+,4} -> X=-1, Y=2, Z=4 (i16)
+;
+; This is an example where the returned solution is the first time an
+; unsigned wrap occurs, whereas the actual exit condition occurs much
+; later. The number of iterations returned by SolveQuadraticEquation
+; is 181, but the loop will iterate 37174 times.
+;
+; Here is a C code that corresponds to this case that calculates the number
+; of iterations:
+;
+; int test04() {
+;   int c = 0;
+;   int ivr = 0;
+;   int inc = -1;
+;   int acc = 2;
+;
+;   while (acc & 0xffff) {
+;     c++;
+;     ivr += inc;
+;     inc += 4;
+;     acc += inc;
+;   }
+;
+;   return c;
+; }
+;
+
+; CHECK-LABEL: Printing analysis 'Scalar Evolution Analysis' for function 'test04':
+; CHECK: GetQuadraticEquation: analyzing quadratic addrec: {0,+,3,+,4}<%loop>
+; CHECK: GetQuadraticEquation: addrec coeff bw: 16
+; CHECK: GetQuadraticEquation: equation 4x^2 + 2x + 0, coeff bw: 17, multiplied by 2
+; CHECK: SolveQuadraticAddRecRange: solving for signed overflow
+; CHECK: SolveQuadraticEquationWrap: solving 4x^2 + 2x + 2, rw:16
+; CHECK: SolveQuadraticEquationWrap: updated coefficients 4x^2 + 2x + -65534, rw:16
+; CHECK: SolveQuadraticEquationWrap: solution (wrap): 128
+; CHECK: SolveQuadraticAddRecRange: solving for unsigned overflow
+; CHECK: SolveQuadraticEquationWrap: solving 4x^2 + 2x + 2, rw:17
+; CHECK: SolveQuadraticEquationWrap: updated coefficients 4x^2 + 2x + -131070, rw:17
+; CHECK: SolveQuadraticEquationWrap: solution (wrap): 181
+; CHECK: SolveQuadraticAddRecRange: solving for signed overflow
+; CHECK: SolveQuadraticEquationWrap: solving 4x^2 + 2x + 2, rw:16
+; CHECK: SolveQuadraticEquationWrap: updated coefficients 4x^2 + 2x + -65534, rw:16
+; CHECK: SolveQuadraticEquationWrap: solution (wrap): 128
+; CHECK: SolveQuadraticAddRecRange: solving for unsigned overflow
+; CHECK: SolveQuadraticEquationWrap: solving 4x^2 + 2x + 2, rw:17
+; CHECK: SolveQuadraticEquationWrap: updated coefficients 4x^2 + 2x + -131070, rw:17
+; CHECK: SolveQuadraticEquationWrap: solution (wrap): 181
+; CHECK: GetQuadraticEquation: analyzing quadratic addrec: {1,+,3,+,4}<%loop>
+; CHECK: GetQuadraticEquation: addrec coeff bw: 16
+; CHECK: GetQuadraticEquation: equation 4x^2 + 2x + 2, coeff bw: 17, multiplied by 2
+; CHECK: SolveQuadraticAddRecExact: solving for unsigned overflow
+; CHECK: SolveQuadraticEquationWrap: solving 4x^2 + 2x + 2, rw:17
+; CHECK: SolveQuadraticEquationWrap: updated coefficients 4x^2 + 2x + -131070, rw:17
+; CHECK: SolveQuadraticEquationWrap: solution (wrap): 181
+; CHECK: Loop %loop: Unpredictable backedge-taken count.
+define signext i32 @test04() {
+entry:
+  br label %loop
+
+loop:
+  %ivr = phi i32 [  0, %entry ], [ %ivr1, %loop ]
+  %inc = phi i32 [ -1, %entry ], [ %inc1, %loop ]
+  %acc = phi i32 [  2, %entry ], [ %acc1, %loop ]
+  %ivr1 = add i32 %ivr, %inc
+  %inc1 = add i32 %inc, 4
+  %acc1 = add i32 %acc, %inc
+  %and  = trunc i32 %acc1 to i16
+  %cond = icmp eq i16 %and, 0
+  br i1 %cond, label %exit, label %loop
+
+exit:
+  %rv = phi i32 [ %acc1, %loop ]
+  ret i32 %rv
+}
+
+; A case with signed arithmetic, but unsigned comparison.
+
+; CHECK-LABEL: Printing analysis 'Scalar Evolution Analysis' for function 'test05':
+; CHECK: GetQuadraticEquation: analyzing quadratic addrec: {0,+,-1,+,-1}<%loop>
+; CHECK: GetQuadraticEquation: addrec coeff bw: 32
+; CHECK: GetQuadraticEquation: equation -1x^2 + -1x + 0, coeff bw: 33, multiplied by 2
+; CHECK: SolveQuadraticAddRecRange: solving for signed overflow
+; CHECK: SolveQuadraticEquationWrap: solving -1x^2 + -1x + 4, rw:32
+; CHECK: SolveQuadraticEquationWrap: updated coefficients 1x^2 + 1x + -4, rw:32
+; CHECK: SolveQuadraticEquationWrap: solution (wrap): 2
+; CHECK: SolveQuadraticAddRecRange: solving for unsigned overflow
+; CHECK: SolveQuadraticEquationWrap: solving -1x^2 + -1x + 4, rw:33
+; CHECK: SolveQuadraticEquationWrap: updated coefficients 1x^2 + 1x + -4, rw:33
+; CHECK: SolveQuadraticEquationWrap: solution (wrap): 2
+; CHECK: SolveQuadraticAddRecRange: solving for signed overflow
+; CHECK: SolveQuadraticEquationWrap: solving -1x^2 + -1x + -2, rw:32
+; CHECK: SolveQuadraticEquationWrap: updated coefficients 1x^2 + 1x + -4294967294, rw:32
+; CHECK: SolveQuadraticEquationWrap: solution (wrap): 65536
+; CHECK: SolveQuadraticAddRecRange: solving for unsigned overflow
+; CHECK: SolveQuadraticEquationWrap: solving -1x^2 + -1x + -2, rw:33
+; CHECK: SolveQuadraticEquationWrap: updated coefficients 1x^2 + 1x + -8589934590, rw:33
+; CHECK: SolveQuadraticEquationWrap: solution (wrap): 92682
+; CHECK: Loop %loop: backedge-taken count is 2
+
+define signext i32 @test05() {
+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
+  %acc1 = add i32 %acc, %inc
+  %and  = and i32 %acc1, -1
+  %cond = icmp ule i32 %and, -3
+  br i1 %cond, label %exit, label %loop
+
+exit:
+  %rv = phi i32 [ %acc1, %loop ]
+  ret i32 %rv
+}
+
+; A test that used to crash with one of the earlier versions of the code.
+
+; CHECK-LABEL: Printing analysis 'Scalar Evolution Analysis' for function 'test06':
+; CHECK: GetQuadraticEquation: analyzing quadratic addrec: {0,+,-99999,+,1}<%loop>
+; CHECK: GetQuadraticEquation: addrec coeff bw: 32
+; CHECK: GetQuadraticEquation: equation 1x^2 + -199999x + 0, coeff bw: 33, multiplied by 2
+; CHECK: SolveQuadraticAddRecRange: solving for signed overflow
+; CHECK: SolveQuadraticEquationWrap: solving 1x^2 + -199999x + -4294967294, rw:32
+; CHECK: SolveQuadraticEquationWrap: updated coefficients 1x^2 + -199999x + 2, rw:32
+; CHECK: SolveQuadraticEquationWrap: solution (wrap): 1
+; CHECK: SolveQuadraticAddRecRange: solving for unsigned overflow
+; CHECK: SolveQuadraticEquationWrap: solving 1x^2 + -199999x + -4294967294, rw:33
+; CHECK: SolveQuadraticEquationWrap: updated coefficients 1x^2 + -199999x + 4294967298, rw:33
+; CHECK: SolveQuadraticEquationWrap: solution (wrap): 24469
+; CHECK: SolveQuadraticAddRecRange: solving for signed overflow
+; CHECK: SolveQuadraticEquationWrap: solving 1x^2 + -199999x + -12, rw:32
+; CHECK: SolveQuadraticEquationWrap: updated coefficients 1x^2 + -199999x + 4294967284, rw:32
+; CHECK: SolveQuadraticEquationWrap: solution (wrap): 24469
+; CHECK: SolveQuadraticAddRecRange: solving for unsigned overflow
+; CHECK: SolveQuadraticEquationWrap: solving 1x^2 + -199999x + -12, rw:33
+; CHECK: SolveQuadraticEquationWrap: updated coefficients 1x^2 + -199999x + 8589934580, rw:33
+; CHECK: SolveQuadraticEquationWrap: solution (wrap): 62450
+; CHECK: Loop %loop: backedge-taken count is 24469
+define signext i32 @test06() {
+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
+  %acc1 = add i32 %acc, %inc
+  %and  = and i32 %acc1, -1
+  %cond = icmp sgt i32 %and, 5
+  br i1 %cond, label %exit, label %loop
+
+exit:
+  %rv = phi i32 [ %acc1, %loop ]
+  ret i32 %rv
+}
+
+; The equation
+;   532052752x^2 + -450429774x + 71188414 = 0
+; has two exact solutions (up to two decimal digits): 0.21 and 0.64.
+; Since there is no integer between them, there is no integer n that either
+; solves the equation exactly, or changes the sign of it between n and n+1.
+
+; CHECK-LABEL: Printing analysis 'Scalar Evolution Analysis' for function 'test07':
+; CHECK: GetQuadraticEquation: analyzing quadratic addrec: {0,+,40811489,+,532052752}<%loop>
+; CHECK: GetQuadraticEquation: addrec coeff bw: 32
+; CHECK: GetQuadraticEquation: equation 532052752x^2 + -450429774x + 0, coeff bw: 33, multiplied by 2
+; CHECK: SolveQuadraticAddRecRange: solving for signed overflow
+; CHECK: SolveQuadraticEquationWrap: solving 532052752x^2 + -450429774x + 71188414, rw:32
+; CHECK: SolveQuadraticEquationWrap: updated coefficients 532052752x^2 + -450429774x + 71188414, rw:32
+; CHECK: SolveQuadraticEquationWrap: no valid solution
+; CHECK: SolveQuadraticAddRecRange: solving for unsigned overflow
+; CHECK: SolveQuadraticEquationWrap: solving 532052752x^2 + -450429774x + 71188414, rw:33
+; CHECK: SolveQuadraticEquationWrap: updated coefficients 532052752x^2 + -450429774x + 71188414, rw:33
+; CHECK: SolveQuadraticEquationWrap: no valid solution
+; CHECK: SolveQuadraticAddRecRange: solving for signed overflow
+; CHECK: SolveQuadraticEquationWrap: solving 532052752x^2 + -450429774x + 71188414, rw:32
+; CHECK: SolveQuadraticEquationWrap: updated coefficients 532052752x^2 + -450429774x + 71188414, rw:32
+; CHECK: SolveQuadraticEquationWrap: no valid solution
+; CHECK: SolveQuadraticAddRecRange: solving for unsigned overflow
+; CHECK: SolveQuadraticEquationWrap: solving 532052752x^2 + -450429774x + 71188414, rw:33
+; CHECK: SolveQuadraticEquationWrap: updated coefficients 532052752x^2 + -450429774x + 71188414, rw:33
+; CHECK: SolveQuadraticEquationWrap: no valid solution
+; CHECK: GetQuadraticEquation: analyzing quadratic addrec: {35594207,+,40811489,+,532052752}<%loop>
+; CHECK: GetQuadraticEquation: addrec coeff bw: 32
+; CHECK: GetQuadraticEquation: equation 532052752x^2 + -450429774x + 71188414, coeff bw: 33, multiplied by 2
+; CHECK: SolveQuadraticAddRecExact: solving for unsigned overflow
+; CHECK: SolveQuadraticEquationWrap: solving 532052752x^2 + -450429774x + 71188414, rw:33
+; CHECK: SolveQuadraticEquationWrap: updated coefficients 532052752x^2 + -450429774x + 71188414, rw:33
+; CHECK: SolveQuadraticEquationWrap: no valid solution
+; CHECK: Loop %loop: Unpredictable backedge-taken count.
+define signext i32 @test07() {
+entry:
+  br label %loop
+
+loop:
+  %ivr = phi i32 [ 0, %entry ], [ %ivr1, %loop ]
+  %inc = phi i32 [ -491241263, %entry ], [ %inc1, %loop ]
+  %acc = phi i32 [ 526835470, %entry ], [ %acc1, %loop ]
+  %ivr1 = add i32 %ivr, %inc
+  %inc1 = add i32 %inc, 532052752
+  %acc1 = add i32 %acc, %inc
+  %and  = and i32 %acc1, -1
+  %cond = icmp eq i32 %and, 0
+  br i1 %cond, label %exit, label %loop
+
+exit:
+  %rv = phi i32 [ %acc1, %loop ]
+  ret i32 %rv
+}
+

unittests/ADT/APIntTest.cpp

 //===- llvm/unittest/ADT/APInt.cpp - APInt unit tests ---------------------===//
 //
 //                     The LLVM Compiler Infrastructure
 //
 // This file is distributed under the University of Illinois Open Source
 // License. See LICENSE.TXT for details.
 //
 //===----------------------------------------------------------------------===//
 
 #include "llvm/ADT/APInt.h"
 #include "llvm/ADT/ArrayRef.h"
 #include "llvm/ADT/SmallString.h"
+#include "llvm/ADT/Twine.h"
 #include "gtest/gtest.h"
 #include <array>
 
 using namespace llvm;
 
 namespace {
 
 TEST(APIntTest, ValueInit) {
for (uint64_t Bi = -128; Bi <= 127; Bi++) {
       {
         APInt Quo = A.sdiv(B);
         EXPECT_EQ(Quo, APIntOps::RoundingSDiv(A, B, APInt::Rounding::TOWARD_ZERO));
       }
     }
   }
 }
 
+TEST(APIntTest, SolveQuadraticEquationWrap) {
+  // Verify that "Solution" is the first non-negative integer that solves
+  // Ax^2 + Bx + C = "0 or overflow", i.e. that it is a correct solution
+  // as calculated by SolveQuadraticEquationWrap.
+  auto Validate = [] (int A, int B, int C, unsigned Width, int Solution) {
+    int Mask = (1 << Width) - 1;
+
+    if (Solution == 0) {
+      EXPECT_EQ(C & Mask, 0);
+      return;
+    }
+    // Solution should be non-negative.
+    EXPECT_GE(Solution, 0);
+
+    auto OverflowBits = [] (int64_t V, unsigned W) {
+      return V & -(1 << W);
+    };
+
+    int64_t Over0 = OverflowBits(C, Width);
+
+    auto IsZeroOrOverflow = [&] (int X) {
+      int64_t VX = A*X*X + B*X + C;
+      int64_t OverX = OverflowBits(VX, Width);
+      return (VX & Mask) == 0 || OverX != Over0;
+    };
+
+    auto EqStr = [&] (const char *X_str) {
+      return Twine(A) + Twine(X_str) + Twine("^2 + ") + Twine(B) +
+             Twine(X_str) + Twine(" + ") + Twine(C) + Twine(", bitwidth: ") +
+             Twine(Width);
+    };
+
+    auto IsSolution = [&] (const char *X_str, int X) {
+      if (IsZeroOrOverflow(X))
+        return ::testing::AssertionSuccess()
+                  << X << " is a solution of " << EqStr(X_str);
+      return ::testing::AssertionFailure()
+                << X << " is not an expected solution of " << EqStr(X_str);
+    };
+
+    auto IsNotSolution = [&] (const char *X_str, int X) {
+      if (!IsZeroOrOverflow(X))
+        return ::testing::AssertionSuccess()
+                  << X << " is not a solution of " << EqStr(X_str);
+      return ::testing::AssertionFailure()
+                << X << " is an unexpected solution of " << EqStr(X_str);
+    };
+
+    // This is the important part: make sure that there is no solution that
+    // is less than the calculated one.
+    for (int X = 1; X < Solution-1; ++X)
+      EXPECT_PRED_FORMAT1(IsNotSolution, X);
+
+    // Verify that the calculated solution is indeed a solution.
+    EXPECT_PRED_FORMAT1(IsSolution, Solution);
+  };
+
+  // Generate all possible quadratic equations with Width-bit wide integer
+  // coefficients, get the solution from SolveQuadraticEquationWrap, and
+  // verify that the solution is correct.
+  auto Iterate = [&] (unsigned Width) {
+    assert(1 < Width && Width < 32);
+    int Low = -(1 << (Width-1));
+    int High = (1 << (Width-1));
+
+    for (int A = Low; A != High; ++A) {
+      if (A == 0)
+        continue;
+      for (int B = Low; B != High; ++B) {
+        for (int C = Low; C != High; ++C) {
+          Optional<APInt> S = APIntOps::SolveQuadraticEquationWrap(
+                                APInt(Width, A), APInt(Width, B),
+                                APInt(Width, C), Width);
+          if (S.hasValue())
+            Validate(A, B, C, Width, S->getSExtValue());
+        }
+      }
+    }
+  };
+
+  // Test all widths in [2..6].
+  for (unsigned i = 2; i <= 6; ++i)
+    Iterate(i);
+}
+
 } // end anonymous namespace