Scalar evolution fails to compute loop iteration count

[efriedma-quic]

Bugzilla Link	2088
Resolution	FIXED
Resolved on	Jul 20, 2008 11:04
Version	trunk
OS	Linux
Attachments	Testcase
CC	@nlewycky,@regehr

Extended Description

To reproduce: compare the output of "lli abortbug.bc" and "opt -indvars abortbug.bc -o - | lli" with the attached testcase.

I think it's somehow illegally rewriting the induction variable z.

Original source (from X86/README.txt):
main ()
{
int i = 0;
unsigned long int z = 0;

do {
z -= 0x00004000;
i++;
if (i > 0x00040000)
abort ();
} while (z > 0);
exit (0);
}
compiled with clang.

[efriedma-quic]

Ah, I found the bug. It's actually in ScalarEvolutionsImpl::HowFarToZero.

Whoever wrote this code apparently thought that the formula for the loop trip count for a loop of the form "int a = Start; while (a) a += Step;" is "(0 - Start/Step) iff Start % Step == 0". The correct formula is "the smallest X such that 0 == Start + (Step * X), or infinite if no such X exists", which is completely different when dealing with numbers mod 2^32. And somehow, nobody managed to spot this bug since the code was added nearly four years ago!

Is there code somewhere in LLVM for computing multiplicative inverses?

[llvmbot]

You're right. Solving modulo equation is the remedy for the bug.

Unfortunately, I don't know of any code for computing multiplicative inverses. I suppose, if there was such code, it would reside in APInt or MathExtras... Are you going to write it? I will need this functionality in my project shortly. If you are willing to wait about a week, I'll prepare a patch for LLVM.

[efriedma-quic]

Well, this is sort of fixed... now we don't compute the iteration count at all.

[nlewycky]

I'd forgotten about this bug! The change is probably because of my fix in bug 2364. We should no longer be miscompiling, but we can't solve every solvable case either, specifically if the induction variable will wrap.

Let's keep this bug open as a feature enhancement for handling such cases.

[chandlerc]

Partial patch for early review
This is just a very early version of the multiplicative inverse code adapted to be APInt safe and relatively efficient. I have ideas about helping it out, but I want to get some early feedback on the design and make sure its correct first.

-Chandler

[chandlerc]

Working patch to provide general modular inverse
Anton has pointed out that there are better algorithms than the naive euclidean for doing this (specifically the montgomery one) but I'm not familiar with it, and I'm not sure its needed at this point. I'm working on utilizing this within SCEV.

[chandlerc]

Full patch
Final patch in the flurry. This one hase SCEV hook into the new functionality to correctly compute the trip count when possible. Its probably a bit rough around the edges (its late), but the test case now succeeds. I'd like to do some more testing, and clean up any really ugly aspects.

[llvmbot]

Full patch
Hi Chandler. Thanks for tackling this.

A few comments:

+// This provides the multiplicative inverse of this APInt modulo its BitWidth.

This line is no longer true.

+// The iterative extended euclidean algorithm is used to solve for this value,
+// however we simplify it to speed up calculating a single inverse, and take
+// advantage of div+rem calculations. We also use some tricks to avoid copying
+// (potentially large) APInts around.
+APInt APInt::multiplicativeInverse(const APInt& modulo) const {

• // We make all of our variables have double the precision of max(this,modulo)
• // to allow any multiplies we do to be stored without loss of precision.
• // TODO: This is currently a fairly safe approximation that no value of x or
• // (r_i / r_{i-1}) will ever exceed max(this,modulo). This holds due to the
• // decreasing nature of the algorithm, but we are probably being far too
• // safe. There should be a sane tighter bound on this.
• uint64_t precision = 2 * std::max(BitWidth, modulo.BitWidth);
• APInt r[2], x[2], Q(precision, 0);
• r[0] = *this;
• r[0].sext(precision);
• r[1] = modulo;
• r[1].sext(precision);
• x[0] = APInt(precision, 1);
• x[1] = APInt(precision, 0);

Treating APInts as signed numbers may produce incorrect results.
Please note that the remainder computed by the sdivrem operation (a rem b)
may not be equivalent to the the result of modulo operation (a mod b)
when a or b is negative. I think all values should be treated as unsigned, so
zext should be used above and udivrem below (rem operation is equivalent
to modulo operation for unsigned numbers).

• // Skip the first iteration, as it is a noop: r[1] is 2^BitWidth, which by
• // definition is greater than n, resulting in a quotient of 0.

This is no longer true for the generalized multiplicativeInverse() method:
r[1] isn't guaranteed to be 2^BitWidth. The loop below probably needs
some correction.

• uint32_t i = 1;
• while (true) {
• // An overview of the math without the confusing bit-flipping:
• // Q = r[i-2] / r[i-1]
• // r[i] = r[i-2] % r[i-1]
• // x[i] = x[i-2] - x[i-1] * Q
• sdivrem(r[i&1], r[(~i)&1], Q, r[i&1]);
• if (r[i&1] == 0)

•  break;
  

• x[i&1] -= x[(~i)&1] * Q;
• ++i;
• }
• // If this, and modulo are not coprime, there is no multiplicative inverse,
• // so return 0. We check this by looking at the next-to-last remainder, which
• // is the gcd(this,modulo) as calculated by the euclidean algorithm.
• if (r[(~i)&1] != 1)
• return APInt(BitWidth,0);
• // We're returning x[i-1] here because we've flipped the bit even though we
• // haven't calculated x[i].
• return x[(~i)&1].trunc(BitWidth);
  +}

• uint32_t lhsValue = LHS.VAL;
• uint32_t rhsValue = RHS.VAL;

uint32_t -> uint64_t

if (AddRec->isAffine()) {
// If this is an affine expression the execution count of this branch is

• // equal to:

• // equal to the smallest X such that:
  //

• // (0 - Start/Step) iff Start % Step == 0

• // Start + (Step * X) = 0 (mod 2^N)
  //
• // Where N is the bit-width of Start. This can be solved for X:
• //
• // X = -Start * Step^-1 (mod 2^N)

I'm afraid it won't work for the general case. Suppose:
Start = -8
Step = 4
N = 8
Trying to compute the multiplicative inverse of Step (mod 2^N) will fail,
because gcd(Step, 2^N) == 4 != 1. However, the minimum root of the original
equation is 2. It seems the general method of solving linear equation should
be used.

[llvmbot]

Chandler, would you mind if I took care of this bug? I've prepared a general routine to solve linear modulo equations some time ago, but stopped working on this PR for a few months. Now I have time again and I could try to fix it.

[chandlerc]

Chandler, would you mind if I took care of this bug? I've prepared a general
routine to solve linear modulo equations some time ago, but stopped working on
this PR for a few months. Now I have time again and I could try to fix it.

Please. I ran into unusual bugs with the various ways I was expanding on the patch to sufficiently solve this case. A general purpose solution may well be a less buggy approach at this point. Let me know if I can help test it out.

-Chandler

[llvmbot]

Chandler, I've used your version of routine computing the multiplicative inverse of a given number. It seems to be more efficient than mine as it copies APInts less often. I've modified it a little bit. Please, review if you are interested:
http://lists.cs.uiuc.edu/pipermail/llvm-commits/Week-of-Mon-20080623/064032.html

(Now, I'm preparing the linear modulo equation solver)

[llvmbot]

*** Bug llvm/llvm-bugzilla-archive#2364 has been marked as a duplicate of this bug. ***

[llvmbot]

Fixed here:
http://lists.cs.uiuc.edu/pipermail/llvm-commits/Week-of-Mon-20080714/065225.html

[llvmbot]

mentioned in issue llvm/llvm-bugzilla-archive#2364