ScalarEvolution.cpp

Go to the documentation of this file.

//===- ScalarEvolution.cpp - Scalar Evolution Analysis ----------*- C++ -*-===//
//
//                     The LLVM Compiler Infrastructure
//
// This file is distributed under the University of Illinois Open Source
// License. See LICENSE.TXT for details.
//
//===----------------------------------------------------------------------===//
//
// This file contains the implementation of the scalar evolution analysis
// engine, which is used primarily to analyze expressions involving induction
// variables in loops.
//
// There are several aspects to this library.  First is the representation of
// scalar expressions, which are represented as subclasses of the SCEV class.
// These classes are used to represent certain types of subexpressions that we
// can handle.  These classes are reference counted, managed by the SCEVHandle
// class.  We only create one SCEV of a particular shape, so pointer-comparisons
// for equality are legal.
//
// One important aspect of the SCEV objects is that they are never cyclic, even
// if there is a cycle in the dataflow for an expression (ie, a PHI node).  If
// the PHI node is one of the idioms that we can represent (e.g., a polynomial
// recurrence) then we represent it directly as a recurrence node, otherwise we
// represent it as a SCEVUnknown node.
//
// In addition to being able to represent expressions of various types, we also
// have folders that are used to build the *canonical* representation for a
// particular expression.  These folders are capable of using a variety of
// rewrite rules to simplify the expressions.
//
// Once the folders are defined, we can implement the more interesting
// higher-level code, such as the code that recognizes PHI nodes of various
// types, computes the execution count of a loop, etc.
//
// TODO: We should use these routines and value representations to implement
// dependence analysis!
//
//===----------------------------------------------------------------------===//
//
// There are several good references for the techniques used in this analysis.
//
//  Chains of recurrences -- a method to expedite the evaluation
//  of closed-form functions
//  Olaf Bachmann, Paul S. Wang, Eugene V. Zima
//
//  On computational properties of chains of recurrences
//  Eugene V. Zima
//
//  Symbolic Evaluation of Chains of Recurrences for Loop Optimization
//  Robert A. van Engelen
//
//  Efficient Symbolic Analysis for Optimizing Compilers
//  Robert A. van Engelen
//
//  Using the chains of recurrences algebra for data dependence testing and
//  induction variable substitution
//  MS Thesis, Johnie Birch
//
//===----------------------------------------------------------------------===//

#define DEBUG_TYPE "scalar-evolution"
#include "llvm/Analysis/ScalarEvolutionExpressions.h"
#include "llvm/Constants.h"
#include "llvm/DerivedTypes.h"
#include "llvm/GlobalVariable.h"
#include "llvm/Instructions.h"
#include "llvm/Analysis/ConstantFolding.h"
#include "llvm/Analysis/LoopInfo.h"
#include "llvm/Assembly/Writer.h"
#include "llvm/Transforms/Scalar.h"
#include "llvm/Support/CFG.h"
#include "llvm/Support/CommandLine.h"
#include "llvm/Support/Compiler.h"
#include "llvm/Support/ConstantRange.h"
#include "llvm/Support/InstIterator.h"
#include "llvm/Support/ManagedStatic.h"
#include "llvm/Support/MathExtras.h"
#include "llvm/Support/Streams.h"
#include "llvm/ADT/Statistic.h"
#include <ostream>
#include <algorithm>
#include <cmath>
using namespace llvm;

STATISTIC(NumArrayLenItCounts,
          "Number of trip counts computed with array length");
STATISTIC(NumTripCountsComputed,
          "Number of loops with predictable loop counts");
STATISTIC(NumTripCountsNotComputed,
          "Number of loops without predictable loop counts");
STATISTIC(NumBruteForceTripCountsComputed,
          "Number of loops with trip counts computed by force");

static cl::opt<unsigned>
MaxBruteForceIterations("scalar-evolution-max-iterations", cl::ReallyHidden,
                        cl::desc("Maximum number of iterations SCEV will "
                                 "symbolically execute a constant derived loop"),
                        cl::init(100));

static RegisterPass<ScalarEvolution>
R("scalar-evolution", "Scalar Evolution Analysis", false, true);
char ScalarEvolution::ID = 0;

//===----------------------------------------------------------------------===//
//                           SCEV class definitions
//===----------------------------------------------------------------------===//

//===----------------------------------------------------------------------===//
// Implementation of the SCEV class.
//
SCEV::~SCEV() {}
void SCEV::dump() const {
  print(cerr);
  cerr << '\n';
}

uint32_t SCEV::getBitWidth() const {
  if (const IntegerType* ITy = dyn_cast<IntegerType>(getType()))
    return ITy->getBitWidth();
  return 0;
}

bool SCEV::isZero() const {
  if (const SCEVConstant *SC = dyn_cast<SCEVConstant>(this))
    return SC->getValue()->isZero();
  return false;
}


SCEVCouldNotCompute::SCEVCouldNotCompute() : SCEV(scCouldNotCompute) {}

bool SCEVCouldNotCompute::isLoopInvariant(const Loop *L) const {
  assert(0 && "Attempt to use a SCEVCouldNotCompute object!");
  return false;
}

const Type *SCEVCouldNotCompute::getType() const {
  assert(0 && "Attempt to use a SCEVCouldNotCompute object!");
  return 0;
}

bool SCEVCouldNotCompute::hasComputableLoopEvolution(const Loop *L) const {
  assert(0 && "Attempt to use a SCEVCouldNotCompute object!");
  return false;
}

SCEVHandle SCEVCouldNotCompute::
replaceSymbolicValuesWithConcrete(const SCEVHandle &Sym,
                                  const SCEVHandle &Conc,
                                  ScalarEvolution &SE) const {
  return this;
}

void SCEVCouldNotCompute::print(std::ostream &OS) const {
  OS << "***COULDNOTCOMPUTE***";
}

bool SCEVCouldNotCompute::classof(const SCEV *S) {
  return S->getSCEVType() == scCouldNotCompute;
}


// SCEVConstants - Only allow the creation of one SCEVConstant for any
// particular value.  Don't use a SCEVHandle here, or else the object will
// never be deleted!
static ManagedStatic<std::map<ConstantInt*, SCEVConstant*> > SCEVConstants;


SCEVConstant::~SCEVConstant() {
  SCEVConstants->erase(V);
}

SCEVHandle ScalarEvolution::getConstant(ConstantInt *V) {
  SCEVConstant *&R = (*SCEVConstants)[V];
  if (R == 0) R = new SCEVConstant(V);
  return R;
}

SCEVHandle ScalarEvolution::getConstant(const APInt& Val) {
  return getConstant(ConstantInt::get(Val));
}

const Type *SCEVConstant::getType() const { return V->getType(); }

void SCEVConstant::print(std::ostream &OS) const {
  WriteAsOperand(OS, V, false);
}

// SCEVTruncates - Only allow the creation of one SCEVTruncateExpr for any
// particular input.  Don't use a SCEVHandle here, or else the object will
// never be deleted!
static ManagedStatic<std::map<std::pair<SCEV*, const Type*>, 
                     SCEVTruncateExpr*> > SCEVTruncates;

SCEVTruncateExpr::SCEVTruncateExpr(const SCEVHandle &op, const Type *ty)
  : SCEV(scTruncate), Op(op), Ty(ty) {
  assert(Op->getType()->isInteger() && Ty->isInteger() &&
         "Cannot truncate non-integer value!");
  assert(Op->getType()->getPrimitiveSizeInBits() > Ty->getPrimitiveSizeInBits()
         && "This is not a truncating conversion!");
}

SCEVTruncateExpr::~SCEVTruncateExpr() {
  SCEVTruncates->erase(std::make_pair(Op, Ty));
}

void SCEVTruncateExpr::print(std::ostream &OS) const {
  OS << "(truncate " << *Op << " to " << *Ty << ")";
}

// SCEVZeroExtends - Only allow the creation of one SCEVZeroExtendExpr for any
// particular input.  Don't use a SCEVHandle here, or else the object will never
// be deleted!
static ManagedStatic<std::map<std::pair<SCEV*, const Type*>,
                     SCEVZeroExtendExpr*> > SCEVZeroExtends;

SCEVZeroExtendExpr::SCEVZeroExtendExpr(const SCEVHandle &op, const Type *ty)
  : SCEV(scZeroExtend), Op(op), Ty(ty) {
  assert(Op->getType()->isInteger() && Ty->isInteger() &&
         "Cannot zero extend non-integer value!");
  assert(Op->getType()->getPrimitiveSizeInBits() < Ty->getPrimitiveSizeInBits()
         && "This is not an extending conversion!");
}

SCEVZeroExtendExpr::~SCEVZeroExtendExpr() {
  SCEVZeroExtends->erase(std::make_pair(Op, Ty));
}

void SCEVZeroExtendExpr::print(std::ostream &OS) const {
  OS << "(zeroextend " << *Op << " to " << *Ty << ")";
}

// SCEVSignExtends - Only allow the creation of one SCEVSignExtendExpr for any
// particular input.  Don't use a SCEVHandle here, or else the object will never
// be deleted!
static ManagedStatic<std::map<std::pair<SCEV*, const Type*>,
                     SCEVSignExtendExpr*> > SCEVSignExtends;

SCEVSignExtendExpr::SCEVSignExtendExpr(const SCEVHandle &op, const Type *ty)
  : SCEV(scSignExtend), Op(op), Ty(ty) {
  assert(Op->getType()->isInteger() && Ty->isInteger() &&
         "Cannot sign extend non-integer value!");
  assert(Op->getType()->getPrimitiveSizeInBits() < Ty->getPrimitiveSizeInBits()
         && "This is not an extending conversion!");
}

SCEVSignExtendExpr::~SCEVSignExtendExpr() {
  SCEVSignExtends->erase(std::make_pair(Op, Ty));
}

void SCEVSignExtendExpr::print(std::ostream &OS) const {
  OS << "(signextend " << *Op << " to " << *Ty << ")";
}

// SCEVCommExprs - Only allow the creation of one SCEVCommutativeExpr for any
// particular input.  Don't use a SCEVHandle here, or else the object will never
// be deleted!
static ManagedStatic<std::map<std::pair<unsigned, std::vector<SCEV*> >,
                     SCEVCommutativeExpr*> > SCEVCommExprs;

SCEVCommutativeExpr::~SCEVCommutativeExpr() {
  SCEVCommExprs->erase(std::make_pair(getSCEVType(),
                                      std::vector<SCEV*>(Operands.begin(),
                                                         Operands.end())));
}

void SCEVCommutativeExpr::print(std::ostream &OS) const {
  assert(Operands.size() > 1 && "This plus expr shouldn't exist!");
  const char *OpStr = getOperationStr();
  OS << "(" << *Operands[0];
  for (unsigned i = 1, e = Operands.size(); i != e; ++i)
    OS << OpStr << *Operands[i];
  OS << ")";
}

SCEVHandle SCEVCommutativeExpr::
replaceSymbolicValuesWithConcrete(const SCEVHandle &Sym,
                                  const SCEVHandle &Conc,
                                  ScalarEvolution &SE) const {
  for (unsigned i = 0, e = getNumOperands(); i != e; ++i) {
    SCEVHandle H =
      getOperand(i)->replaceSymbolicValuesWithConcrete(Sym, Conc, SE);
    if (H != getOperand(i)) {
      std::vector<SCEVHandle> NewOps;
      NewOps.reserve(getNumOperands());
      for (unsigned j = 0; j != i; ++j)
        NewOps.push_back(getOperand(j));
      NewOps.push_back(H);
      for (++i; i != e; ++i)
        NewOps.push_back(getOperand(i)->
                         replaceSymbolicValuesWithConcrete(Sym, Conc, SE));

      if (isa<SCEVAddExpr>(this))
        return SE.getAddExpr(NewOps);
      else if (isa<SCEVMulExpr>(this))
        return SE.getMulExpr(NewOps);
      else if (isa<SCEVSMaxExpr>(this))
        return SE.getSMaxExpr(NewOps);
      else if (isa<SCEVUMaxExpr>(this))
        return SE.getUMaxExpr(NewOps);
      else
        assert(0 && "Unknown commutative expr!");
    }
  }
  return this;
}


// SCEVUDivs - Only allow the creation of one SCEVUDivExpr for any particular
// input.  Don't use a SCEVHandle here, or else the object will never be
// deleted!
static ManagedStatic<std::map<std::pair<SCEV*, SCEV*>, 
                     SCEVUDivExpr*> > SCEVUDivs;

SCEVUDivExpr::~SCEVUDivExpr() {
  SCEVUDivs->erase(std::make_pair(LHS, RHS));
}

void SCEVUDivExpr::print(std::ostream &OS) const {
  OS << "(" << *LHS << " /u " << *RHS << ")";
}

const Type *SCEVUDivExpr::getType() const {
  return LHS->getType();
}


// SCEVSDivs - Only allow the creation of one SCEVSDivExpr for any particular
// input.  Don't use a SCEVHandle here, or else the object will never be
// deleted!
static ManagedStatic<std::map<std::pair<SCEV*, SCEV*>, 
                     SCEVSDivExpr*> > SCEVSDivs;

SCEVSDivExpr::~SCEVSDivExpr() {
  SCEVSDivs->erase(std::make_pair(LHS, RHS));
}

void SCEVSDivExpr::print(std::ostream &OS) const {
  OS << "(" << *LHS << " /s " << *RHS << ")";
}

const Type *SCEVSDivExpr::getType() const {
  return LHS->getType();
}


// SCEVAddRecExprs - Only allow the creation of one SCEVAddRecExpr for any
// particular input.  Don't use a SCEVHandle here, or else the object will never
// be deleted!
static ManagedStatic<std::map<std::pair<const Loop *, std::vector<SCEV*> >,
                     SCEVAddRecExpr*> > SCEVAddRecExprs;

SCEVAddRecExpr::~SCEVAddRecExpr() {
  SCEVAddRecExprs->erase(std::make_pair(L,
                                        std::vector<SCEV*>(Operands.begin(),
                                                           Operands.end())));
}

SCEVHandle SCEVAddRecExpr::
replaceSymbolicValuesWithConcrete(const SCEVHandle &Sym,
                                  const SCEVHandle &Conc,
                                  ScalarEvolution &SE) const {
  for (unsigned i = 0, e = getNumOperands(); i != e; ++i) {
    SCEVHandle H =
      getOperand(i)->replaceSymbolicValuesWithConcrete(Sym, Conc, SE);
    if (H != getOperand(i)) {
      std::vector<SCEVHandle> NewOps;
      NewOps.reserve(getNumOperands());
      for (unsigned j = 0; j != i; ++j)
        NewOps.push_back(getOperand(j));
      NewOps.push_back(H);
      for (++i; i != e; ++i)
        NewOps.push_back(getOperand(i)->
                         replaceSymbolicValuesWithConcrete(Sym, Conc, SE));

      return SE.getAddRecExpr(NewOps, L);
    }
  }
  return this;
}


bool SCEVAddRecExpr::isLoopInvariant(const Loop *QueryLoop) const {
  // This recurrence is invariant w.r.t to QueryLoop iff QueryLoop doesn't
  // contain L and if the start is invariant.
  return !QueryLoop->contains(L->getHeader()) &&
         getOperand(0)->isLoopInvariant(QueryLoop);
}


void SCEVAddRecExpr::print(std::ostream &OS) const {
  OS << "{" << *Operands[0];
  for (unsigned i = 1, e = Operands.size(); i != e; ++i)
    OS << ",+," << *Operands[i];
  OS << "}<" << L->getHeader()->getName() + ">";
}

// SCEVUnknowns - Only allow the creation of one SCEVUnknown for any particular
// value.  Don't use a SCEVHandle here, or else the object will never be
// deleted!
static ManagedStatic<std::map<Value*, SCEVUnknown*> > SCEVUnknowns;

SCEVUnknown::~SCEVUnknown() { SCEVUnknowns->erase(V); }

bool SCEVUnknown::isLoopInvariant(const Loop *L) const {
  // All non-instruction values are loop invariant.  All instructions are loop
  // invariant if they are not contained in the specified loop.
  if (Instruction *I = dyn_cast<Instruction>(V))
    return !L->contains(I->getParent());
  return true;
}

const Type *SCEVUnknown::getType() const {
  return V->getType();
}

void SCEVUnknown::print(std::ostream &OS) const {
  WriteAsOperand(OS, V, false);
}

//===----------------------------------------------------------------------===//
//                               SCEV Utilities
//===----------------------------------------------------------------------===//

namespace {
  /// SCEVComplexityCompare - Return true if the complexity of the LHS is less
  /// than the complexity of the RHS.  This comparator is used to canonicalize
  /// expressions.
  struct VISIBILITY_HIDDEN SCEVComplexityCompare {
    bool operator()(const SCEV *LHS, const SCEV *RHS) const {
      return LHS->getSCEVType() < RHS->getSCEVType();
    }
  };
}

/// GroupByComplexity - Given a list of SCEV objects, order them by their
/// complexity, and group objects of the same complexity together by value.
/// When this routine is finished, we know that any duplicates in the vector are
/// consecutive and that complexity is monotonically increasing.
///
/// Note that we go take special precautions to ensure that we get determinstic
/// results from this routine.  In other words, we don't want the results of
/// this to depend on where the addresses of various SCEV objects happened to
/// land in memory.
///
static void GroupByComplexity(std::vector<SCEVHandle> &Ops) {
  if (Ops.size() < 2) return;  // Noop
  if (Ops.size() == 2) {
    // This is the common case, which also happens to be trivially simple.
    // Special case it.
    if (SCEVComplexityCompare()(Ops[1], Ops[0]))
      std::swap(Ops[0], Ops[1]);
    return;
  }

  // Do the rough sort by complexity.
  std::sort(Ops.begin(), Ops.end(), SCEVComplexityCompare());

  // Now that we are sorted by complexity, group elements of the same
  // complexity.  Note that this is, at worst, N^2, but the vector is likely to
  // be extremely short in practice.  Note that we take this approach because we
  // do not want to depend on the addresses of the objects we are grouping.
  for (unsigned i = 0, e = Ops.size(); i != e-2; ++i) {
    SCEV *S = Ops[i];
    unsigned Complexity = S->getSCEVType();

    // If there are any objects of the same complexity and same value as this
    // one, group them.
    for (unsigned j = i+1; j != e && Ops[j]->getSCEVType() == Complexity; ++j) {
      if (Ops[j] == S) { // Found a duplicate.
        // Move it to immediately after i'th element.
        std::swap(Ops[i+1], Ops[j]);
        ++i;   // no need to rescan it.
        if (i == e-2) return;  // Done!
      }
    }
  }
}



//===----------------------------------------------------------------------===//
//                      Simple SCEV method implementations
//===----------------------------------------------------------------------===//

/// getIntegerSCEV - Given an integer or FP type, create a constant for the
/// specified signed integer value and return a SCEV for the constant.
SCEVHandle ScalarEvolution::getIntegerSCEV(int Val, const Type *Ty) {
  Constant *C;
  if (Val == 0)
    C = Constant::getNullValue(Ty);
  else if (Ty->isFloatingPoint())
    C = ConstantFP::get(APFloat(Ty==Type::FloatTy ? APFloat::IEEEsingle : 
                                APFloat::IEEEdouble, Val));
  else 
    C = ConstantInt::get(Ty, Val);
  return getUnknown(C);
}

/// getNegativeSCEV - Return a SCEV corresponding to -V = -1*V
///
SCEVHandle ScalarEvolution::getNegativeSCEV(const SCEVHandle &V) {
  if (SCEVConstant *VC = dyn_cast<SCEVConstant>(V))
    return getUnknown(ConstantExpr::getNeg(VC->getValue()));

  return getMulExpr(V, getConstant(ConstantInt::getAllOnesValue(V->getType())));
}

/// getNotSCEV - Return a SCEV corresponding to ~V = -1-V
SCEVHandle ScalarEvolution::getNotSCEV(const SCEVHandle &V) {
  if (SCEVConstant *VC = dyn_cast<SCEVConstant>(V))
    return getUnknown(ConstantExpr::getNot(VC->getValue()));

  SCEVHandle AllOnes = getConstant(ConstantInt::getAllOnesValue(V->getType()));
  return getMinusSCEV(AllOnes, V);
}

/// getMinusSCEV - Return a SCEV corresponding to LHS - RHS.
///
SCEVHandle ScalarEvolution::getMinusSCEV(const SCEVHandle &LHS,
                                         const SCEVHandle &RHS) {
  // X - Y --> X + -Y
  return getAddExpr(LHS, getNegativeSCEV(RHS));
}


/// BinomialCoefficient - Compute BC(It, K).  The result has width W.
// Assume, K > 0.
static SCEVHandle BinomialCoefficient(SCEVHandle It, unsigned K,
                                      ScalarEvolution &SE,
                                      const IntegerType* ResultTy) {
  // Handle the simplest case efficiently.
  if (K == 1)
    return SE.getTruncateOrZeroExtend(It, ResultTy);

  // We are using the following formula for BC(It, K):
  //
  //   BC(It, K) = (It * (It - 1) * ... * (It - K + 1)) / K!
  //
  // Suppose, W is the bitwidth of the return value.  We must be prepared for
  // overflow.  Hence, we must assure that the result of our computation is
  // equal to the accurate one modulo 2^W.  Unfortunately, division isn't
  // safe in modular arithmetic.
  //
  // However, this code doesn't use exactly that formula; the formula it uses
  // is something like the following, where T is the number of factors of 2 in 
  // K! (i.e. trailing zeros in the binary representation of K!), and ^ is
  // exponentiation:
  //
  //   BC(It, K) = (It * (It - 1) * ... * (It - K + 1)) / 2^T / (K! / 2^T)
  //
  // This formula is trivially equivalent to the previous formula.  However,
  // this formula can be implemented much more efficiently.  The trick is that
  // K! / 2^T is odd, and exact division by an odd number *is* safe in modular
  // arithmetic.  To do exact division in modular arithmetic, all we have
  // to do is multiply by the inverse.  Therefore, this step can be done at
  // width W.
  // 
  // The next issue is how to safely do the division by 2^T.  The way this
  // is done is by doing the multiplication step at a width of at least W + T
  // bits.  This way, the bottom W+T bits of the product are accurate. Then,
  // when we perform the division by 2^T (which is equivalent to a right shift
  // by T), the bottom W bits are accurate.  Extra bits are okay; they'll get
  // truncated out after the division by 2^T.
  //
  // In comparison to just directly using the first formula, this technique
  // is much more efficient; using the first formula requires W * K bits,
  // but this formula less than W + K bits. Also, the first formula requires
  // a division step, whereas this formula only requires multiplies and shifts.
  //
  // It doesn't matter whether the subtraction step is done in the calculation
  // width or the input iteration count's width; if the subtraction overflows,
  // the result must be zero anyway.  We prefer here to do it in the width of
  // the induction variable because it helps a lot for certain cases; CodeGen
  // isn't smart enough to ignore the overflow, which leads to much less
  // efficient code if the width of the subtraction is wider than the native
  // register width.
  //
  // (It's possible to not widen at all by pulling out factors of 2 before
  // the multiplication; for example, K=2 can be calculated as
  // It/2*(It+(It*INT_MIN/INT_MIN)+-1). However, it requires
  // extra arithmetic, so it's not an obvious win, and it gets
  // much more complicated for K > 3.)

  // Protection from insane SCEVs; this bound is conservative,
  // but it probably doesn't matter.
  if (K > 1000)
    return new SCEVCouldNotCompute();

  unsigned W = ResultTy->getBitWidth();

  // Calculate K! / 2^T and T; we divide out the factors of two before
  // multiplying for calculating K! / 2^T to avoid overflow.
  // Other overflow doesn't matter because we only care about the bottom
  // W bits of the result.
  APInt OddFactorial(W, 1);
  unsigned T = 1;
  for (unsigned i = 3; i <= K; ++i) {
    APInt Mult(W, i);
    unsigned TwoFactors = Mult.countTrailingZeros();
    T += TwoFactors;
    Mult = Mult.lshr(TwoFactors);
    OddFactorial *= Mult;
  }

  // We need at least W + T bits for the multiplication step
  // FIXME: A temporary hack; we round up the bitwidths
  // to the nearest power of 2 to be nice to the code generator.
  unsigned CalculationBits = 1U << Log2_32_Ceil(W + T);
  // FIXME: Temporary hack to avoid generating integers that are too wide.
  // Although, it's not completely clear how to determine how much
  // widening is safe; for example, on X86, we can't really widen
  // beyond 64 because we need to be able to do multiplication
  // that's CalculationBits wide, but on X86-64, we can safely widen up to
  // 128 bits.
  if (CalculationBits > 64)
    return new SCEVCouldNotCompute();

  // Calcuate 2^T, at width T+W.
  APInt DivFactor = APInt(CalculationBits, 1).shl(T);

  // Calculate the multiplicative inverse of K! / 2^T;
  // this multiplication factor will perform the exact division by
  // K! / 2^T.
  APInt Mod = APInt::getSignedMinValue(W+1);
  APInt MultiplyFactor = OddFactorial.zext(W+1);
  MultiplyFactor = MultiplyFactor.multiplicativeInverse(Mod);
  MultiplyFactor = MultiplyFactor.trunc(W);

  // Calculate the product, at width T+W
  const IntegerType *CalculationTy = IntegerType::get(CalculationBits);
  SCEVHandle Dividend = SE.getTruncateOrZeroExtend(It, CalculationTy);
  for (unsigned i = 1; i != K; ++i) {
    SCEVHandle S = SE.getMinusSCEV(It, SE.getIntegerSCEV(i, It->getType()));
    Dividend = SE.getMulExpr(Dividend,
                             SE.getTruncateOrZeroExtend(S, CalculationTy));
  }

  // Divide by 2^T
  SCEVHandle DivResult = SE.getUDivExpr(Dividend, SE.getConstant(DivFactor));

  // Truncate the result, and divide by K! / 2^T.

  return SE.getMulExpr(SE.getConstant(MultiplyFactor),
                       SE.getTruncateOrZeroExtend(DivResult, ResultTy));
}

/// evaluateAtIteration - Return the value of this chain of recurrences at
/// the specified iteration number.  We can evaluate this recurrence by
/// multiplying each element in the chain by the binomial coefficient
/// corresponding to it.  In other words, we can evaluate {A,+,B,+,C,+,D} as:
///
///   A*BC(It, 0) + B*BC(It, 1) + C*BC(It, 2) + D*BC(It, 3)
///
/// where BC(It, k) stands for binomial coefficient.
///
SCEVHandle SCEVAddRecExpr::evaluateAtIteration(SCEVHandle It,
                                               ScalarEvolution &SE) const {
  SCEVHandle Result = getStart();
  for (unsigned i = 1, e = getNumOperands(); i != e; ++i) {
    // The computation is correct in the face of overflow provided that the
    // multiplication is performed _after_ the evaluation of the binomial
    // coefficient.
    SCEVHandle Coeff = BinomialCoefficient(It, i, SE,
                                           cast<IntegerType>(getType()));
    if (isa<SCEVCouldNotCompute>(Coeff))
      return Coeff;

    Result = SE.getAddExpr(Result, SE.getMulExpr(getOperand(i), Coeff));
  }
  return Result;
}

//===----------------------------------------------------------------------===//
//                    SCEV Expression folder implementations
//===----------------------------------------------------------------------===//

SCEVHandle ScalarEvolution::getTruncateExpr(const SCEVHandle &Op, const Type *Ty) {
  if (SCEVConstant *SC = dyn_cast<SCEVConstant>(Op))
    return getUnknown(
        ConstantExpr::getTrunc(SC->getValue(), Ty));

  // If the input value is a chrec scev made out of constants, truncate
  // all of the constants.
  if (SCEVAddRecExpr *AddRec = dyn_cast<SCEVAddRecExpr>(Op)) {
    std::vector<SCEVHandle> Operands;
    for (unsigned i = 0, e = AddRec->getNumOperands(); i != e; ++i)
      // FIXME: This should allow truncation of other expression types!
      if (isa<SCEVConstant>(AddRec->getOperand(i)))
        Operands.push_back(getTruncateExpr(AddRec->getOperand(i), Ty));
      else
        break;
    if (Operands.size() == AddRec->getNumOperands())
      return getAddRecExpr(Operands, AddRec->getLoop());
  }

  SCEVTruncateExpr *&Result = (*SCEVTruncates)[std::make_pair(Op, Ty)];
  if (Result == 0) Result = new SCEVTruncateExpr(Op, Ty);
  return Result;
}

SCEVHandle ScalarEvolution::getZeroExtendExpr(const SCEVHandle &Op, const Type *Ty) {
  if (SCEVConstant *SC = dyn_cast<SCEVConstant>(Op))
    return getUnknown(
        ConstantExpr::getZExt(SC->getValue(), Ty));

  // FIXME: If the input value is a chrec scev, and we can prove that the value
  // did not overflow the old, smaller, value, we can zero extend all of the
  // operands (often constants).  This would allow analysis of something like
  // this:  for (unsigned char X = 0; X < 100; ++X) { int Y = X; }

  SCEVZeroExtendExpr *&Result = (*SCEVZeroExtends)[std::make_pair(Op, Ty)];
  if (Result == 0) Result = new SCEVZeroExtendExpr(Op, Ty);
  return Result;
}

SCEVHandle ScalarEvolution::getSignExtendExpr(const SCEVHandle &Op, const Type *Ty) {
  if (SCEVConstant *SC = dyn_cast<SCEVConstant>(Op))
    return getUnknown(
        ConstantExpr::getSExt(SC->getValue(), Ty));

  // FIXME: If the input value is a chrec scev, and we can prove that the value
  // did not overflow the old, smaller, value, we can sign extend all of the
  // operands (often constants).  This would allow analysis of something like
  // this:  for (signed char X = 0; X < 100; ++X) { int Y = X; }

  SCEVSignExtendExpr *&Result = (*SCEVSignExtends)[std::make_pair(Op, Ty)];
  if (Result == 0) Result = new SCEVSignExtendExpr(Op, Ty);
  return Result;
}

/// getTruncateOrZeroExtend - Return a SCEV corresponding to a conversion
/// of the input value to the specified type.  If the type must be
/// extended, it is zero extended.
SCEVHandle ScalarEvolution::getTruncateOrZeroExtend(const SCEVHandle &V,
                                                    const Type *Ty) {
  const Type *SrcTy = V->getType();
  assert(SrcTy->isInteger() && Ty->isInteger() &&
         "Cannot truncate or zero extend with non-integer arguments!");
  if (SrcTy->getPrimitiveSizeInBits() == Ty->getPrimitiveSizeInBits())
    return V;  // No conversion
  if (SrcTy->getPrimitiveSizeInBits() > Ty->getPrimitiveSizeInBits())
    return getTruncateExpr(V, Ty);
  return getZeroExtendExpr(V, Ty);
}

// get - Get a canonical add expression, or something simpler if possible.
SCEVHandle ScalarEvolution::getAddExpr(std::vector<SCEVHandle> &Ops) {
  assert(!Ops.empty() && "Cannot get empty add!");
  if (Ops.size() == 1) return Ops[0];

  // Sort by complexity, this groups all similar expression types together.
  GroupByComplexity(Ops);

  // If there are any constants, fold them together.
  unsigned Idx = 0;
  if (SCEVConstant *LHSC = dyn_cast<SCEVConstant>(Ops[0])) {
    ++Idx;
    assert(Idx < Ops.size());
    while (SCEVConstant *RHSC = dyn_cast<SCEVConstant>(Ops[Idx])) {
      // We found two constants, fold them together!
      ConstantInt *Fold = ConstantInt::get(LHSC->getValue()->getValue() + 
                                           RHSC->getValue()->getValue());
      Ops[0] = getConstant(Fold);
      Ops.erase(Ops.begin()+1);  // Erase the folded element
      if (Ops.size() == 1) return Ops[0];
      LHSC = cast<SCEVConstant>(Ops[0]);
    }

    // If we are left with a constant zero being added, strip it off.
    if (cast<SCEVConstant>(Ops[0])->getValue()->isZero()) {
      Ops.erase(Ops.begin());
      --Idx;
    }
  }

  if (Ops.size() == 1) return Ops[0];

  // Okay, check to see if the same value occurs in the operand list twice.  If
  // so, merge them together into an multiply expression.  Since we sorted the
  // list, these values are required to be adjacent.
  const Type *Ty = Ops[0]->getType();
  for (unsigned i = 0, e = Ops.size()-1; i != e; ++i)
    if (Ops[i] == Ops[i+1]) {      //  X + Y + Y  -->  X + Y*2
      // Found a match, merge the two values into a multiply, and add any
      // remaining values to the result.
      SCEVHandle Two = getIntegerSCEV(2, Ty);
      SCEVHandle Mul = getMulExpr(Ops[i], Two);
      if (Ops.size() == 2)
        return Mul;
      Ops.erase(Ops.begin()+i, Ops.begin()+i+2);
      Ops.push_back(Mul);
      return getAddExpr(Ops);
    }

  // Now we know the first non-constant operand.  Skip past any cast SCEVs.
  while (Idx < Ops.size() && Ops[Idx]->getSCEVType() < scAddExpr)
    ++Idx;

  // If there are add operands they would be next.
  if (Idx < Ops.size()) {
    bool DeletedAdd = false;
    while (SCEVAddExpr *Add = dyn_cast<SCEVAddExpr>(Ops[Idx])) {
      // If we have an add, expand the add operands onto the end of the operands
      // list.
      Ops.insert(Ops.end(), Add->op_begin(), Add->op_end());
      Ops.erase(Ops.begin()+Idx);
      DeletedAdd = true;
    }

    // If we deleted at least one add, we added operands to the end of the list,
    // and they are not necessarily sorted.  Recurse to resort and resimplify
    // any operands we just aquired.
    if (DeletedAdd)
      return getAddExpr(Ops);
  }

  // Skip over the add expression until we get to a multiply.
  while (Idx < Ops.size() && Ops[Idx]->getSCEVType() < scMulExpr)
    ++Idx;

  // If we are adding something to a multiply expression, make sure the
  // something is not already an operand of the multiply.  If so, merge it into
  // the multiply.
  for (; Idx < Ops.size() && isa<SCEVMulExpr>(Ops[Idx]); ++Idx) {
    SCEVMulExpr *Mul = cast<SCEVMulExpr>(Ops[Idx]);
    for (unsigned MulOp = 0, e = Mul->getNumOperands(); MulOp != e; ++MulOp) {
      SCEV *MulOpSCEV = Mul->getOperand(MulOp);
      for (unsigned AddOp = 0, e = Ops.size(); AddOp != e; ++AddOp)
        if (MulOpSCEV == Ops[AddOp] && !isa<SCEVConstant>(MulOpSCEV)) {
          // Fold W + X + (X * Y * Z)  -->  W + (X * ((Y*Z)+1))
          SCEVHandle InnerMul = Mul->getOperand(MulOp == 0);
          if (Mul->getNumOperands() != 2) {
            // If the multiply has more than two operands, we must get the
            // Y*Z term.
            std::vector<SCEVHandle> MulOps(Mul->op_begin(), Mul->op_end());
            MulOps.erase(MulOps.begin()+MulOp);
            InnerMul = getMulExpr(MulOps);
          }
          SCEVHandle One = getIntegerSCEV(1, Ty);
          SCEVHandle AddOne = getAddExpr(InnerMul, One);
          SCEVHandle OuterMul = getMulExpr(AddOne, Ops[AddOp]);
          if (Ops.size() == 2) return OuterMul;
          if (AddOp < Idx) {
            Ops.erase(Ops.begin()+AddOp);
            Ops.erase(Ops.begin()+Idx-1);
          } else {
            Ops.erase(Ops.begin()+Idx);
            Ops.erase(Ops.begin()+AddOp-1);
          }
          Ops.push_back(OuterMul);
          return getAddExpr(Ops);
        }

      // Check this multiply against other multiplies being added together.
      for (unsigned OtherMulIdx = Idx+1;
           OtherMulIdx < Ops.size() && isa<SCEVMulExpr>(Ops[OtherMulIdx]);
           ++OtherMulIdx) {
        SCEVMulExpr *OtherMul = cast<SCEVMulExpr>(Ops[OtherMulIdx]);
        // If MulOp occurs in OtherMul, we can fold the two multiplies
        // together.
        for (unsigned OMulOp = 0, e = OtherMul->getNumOperands();
             OMulOp != e; ++OMulOp)
          if (OtherMul->getOperand(OMulOp) == MulOpSCEV) {
            // Fold X + (A*B*C) + (A*D*E) --> X + (A*(B*C+D*E))
            SCEVHandle InnerMul1 = Mul->getOperand(MulOp == 0);
            if (Mul->getNumOperands() != 2) {
              std::vector<SCEVHandle> MulOps(Mul->op_begin(), Mul->op_end());
              MulOps.erase(MulOps.begin()+MulOp);
              InnerMul1 = getMulExpr(MulOps);
            }
            SCEVHandle InnerMul2 = OtherMul->getOperand(OMulOp == 0);
            if (OtherMul->getNumOperands() != 2) {
              std::vector<SCEVHandle> MulOps(OtherMul->op_begin(),
                                             OtherMul->op_end());
              MulOps.erase(MulOps.begin()+OMulOp);
              InnerMul2 = getMulExpr(MulOps);
            }
            SCEVHandle InnerMulSum = getAddExpr(InnerMul1,InnerMul2);
            SCEVHandle OuterMul = getMulExpr(MulOpSCEV, InnerMulSum);
            if (Ops.size() == 2) return OuterMul;
            Ops.erase(Ops.begin()+Idx);
            Ops.erase(Ops.begin()+OtherMulIdx-1);
            Ops.push_back(OuterMul);
            return getAddExpr(Ops);
          }
      }
    }
  }

  // If there are any add recurrences in the operands list, see if any other
  // added values are loop invariant.  If so, we can fold them into the
  // recurrence.
  while (Idx < Ops.size() && Ops[Idx]->getSCEVType() < scAddRecExpr)
    ++Idx;

  // Scan over all recurrences, trying to fold loop invariants into them.
  for (; Idx < Ops.size() && isa<SCEVAddRecExpr>(Ops[Idx]); ++Idx) {
    // Scan all of the other operands to this add and add them to the vector if
    // they are loop invariant w.r.t. the recurrence.
    std::vector<SCEVHandle> LIOps;
    SCEVAddRecExpr *AddRec = cast<SCEVAddRecExpr>(Ops[Idx]);
    for (unsigned i = 0, e = Ops.size(); i != e; ++i)
      if (Ops[i]->isLoopInvariant(AddRec->getLoop())) {
        LIOps.push_back(Ops[i]);
        Ops.erase(Ops.begin()+i);
        --i; --e;
      }

    // If we found some loop invariants, fold them into the recurrence.
    if (!LIOps.empty()) {
      //  NLI + LI + {Start,+,Step}  -->  NLI + {LI+Start,+,Step}
      LIOps.push_back(AddRec->getStart());

      std::vector<SCEVHandle> AddRecOps(AddRec->op_begin(), AddRec->op_end());
      AddRecOps[0] = getAddExpr(LIOps);

      SCEVHandle NewRec = getAddRecExpr(AddRecOps, AddRec->