#include "ModulusRemainder.h"
#include "IR.h"
#include "IROperator.h"
#include "IRPrinter.h"
#include "Simplify.h"
namespace Halide {
namespace Internal {
namespace {
// A version of mod where a % 0 == a
int64_t mod(int64_t a, int64_t b) {
if (b == 0) {
return a;
} else {
return mod_imp(a, b);
}
}
class ComputeModulusRemainder : public IRVisitor {
public:
ModulusRemainder analyze(const Expr &e);
ModulusRemainder result;
Scope<ModulusRemainder> scope;
ComputeModulusRemainder(const Scope<ModulusRemainder> *s) {
scope.set_containing_scope(s);
}
void visit(const IntImm *) override;
void visit(const UIntImm *) override;
void visit(const FloatImm *) override;
void visit(const StringImm *) override;
void visit(const Cast *) override;
void visit(const Variable *) override;
void visit(const Add *) override;
void visit(const Sub *) override;
void visit(const Mul *) override;
void visit(const Div *) override;
void visit(const Mod *) override;
void visit(const Min *) override;
void visit(const Max *) override;
void visit(const EQ *) override;
void visit(const NE *) override;
void visit(const LT *) override;
void visit(const LE *) override;
void visit(const GT *) override;
void visit(const GE *) override;
void visit(const And *) override;
void visit(const Or *) override;
void visit(const Not *) override;
void visit(const Select *) override;
void visit(const Load *) override;
void visit(const Ramp *) override;
void visit(const Broadcast *) override;
void visit(const Call *) override;
void visit(const Let *) override;
void visit(const LetStmt *) override;
void visit(const AssertStmt *) override;
void visit(const ProducerConsumer *) override;
void visit(const For *) override;
void visit(const Acquire *) override;
void visit(const Store *) override;
void visit(const Provide *) override;
void visit(const Allocate *) override;
void visit(const Realize *) override;
void visit(const Block *) override;
void visit(const Fork *) override;
void visit(const IfThenElse *) override;
void visit(const Free *) override;
void visit(const Evaluate *) override;
void visit(const Shuffle *) override;
void visit(const VectorReduce *) override;
void visit(const Prefetch *) override;
void visit(const Atomic *) override;
};
void ComputeModulusRemainder::visit(const IntImm *op) {
// Equal to op->value modulo anything. We'll use zero as the
// modulus to mark this special case. We'd better be able to
// handle zero in the rest of the code...
result = {0, op->value};
}
void ComputeModulusRemainder::visit(const UIntImm *op) {
internal_error << "modulus_remainder of uint\n";
}
void ComputeModulusRemainder::visit(const FloatImm *) {
internal_error << "modulus_remainder of float\n";
}
void ComputeModulusRemainder::visit(const StringImm *) {
internal_error << "modulus_remainder of string\n";
}
void ComputeModulusRemainder::visit(const Cast *) {
// TODO: Could probably do something reasonable for integer
// upcasts and downcasts where the modulus is a power of two.
result = ModulusRemainder{};
}
void ComputeModulusRemainder::visit(const Variable *op) {
if (scope.contains(op->name)) {
result = scope.get(op->name);
} else {
result = ModulusRemainder{};
}
}
void ComputeModulusRemainder::visit(const Add *op) {
result = analyze(op->a) + analyze(op->b);
}
void ComputeModulusRemainder::visit(const Sub *op) {
result = analyze(op->a) - analyze(op->b);
}
void ComputeModulusRemainder::visit(const Mul *op) {
result = analyze(op->a) * analyze(op->b);
}
void ComputeModulusRemainder::visit(const Div *op) {
result = analyze(op->a) / analyze(op->b);
}
void ComputeModulusRemainder::visit(const Min *op) {
result = ModulusRemainder::unify(analyze(op->a), analyze(op->b));
}
void ComputeModulusRemainder::visit(const Max *op) {
result = ModulusRemainder::unify(analyze(op->a), analyze(op->b));
}
void ComputeModulusRemainder::visit(const EQ *) {
internal_error << "modulus_remainder of bool\n";
}
void ComputeModulusRemainder::visit(const NE *) {
internal_error << "modulus_remainder of bool\n";
}
void ComputeModulusRemainder::visit(const LT *) {
internal_error << "modulus_remainder of bool\n";
}
void ComputeModulusRemainder::visit(const LE *) {
internal_error << "modulus_remainder of bool\n";
}
void ComputeModulusRemainder::visit(const GT *) {
internal_error << "modulus_remainder of bool\n";
}
void ComputeModulusRemainder::visit(const GE *) {
internal_error << "modulus_remainder of bool\n";
}
void ComputeModulusRemainder::visit(const And *) {
internal_error << "modulus_remainder of bool\n";
}
void ComputeModulusRemainder::visit(const Or *) {
internal_error << "modulus_remainder of bool\n";
}
void ComputeModulusRemainder::visit(const Not *) {
internal_error << "modulus_remainder of bool\n";
}
void ComputeModulusRemainder::visit(const Select *op) {
result = ModulusRemainder::unify(analyze(op->true_value),
analyze(op->false_value));
}
void ComputeModulusRemainder::visit(const Load *) {
result = ModulusRemainder{};
}
void ComputeModulusRemainder::visit(const Ramp *) {
internal_error << "modulus_remainder of vector\n";
}
void ComputeModulusRemainder::visit(const Broadcast *) {
internal_error << "modulus_remainder of vector\n";
}
void ComputeModulusRemainder::visit(const Call *) {
result = ModulusRemainder{};
}
void ComputeModulusRemainder::visit(const Let *op) {
if (op->value.type().is_int()) {
ScopedBinding<ModulusRemainder> bind(scope, op->name, analyze(op->value));
result = analyze(op->body);
} else {
result = analyze(op->body);
}
}
void ComputeModulusRemainder::visit(const Shuffle *op) {
// It's possible that scalar expressions are extracting a lane of
// a vector - don't fail in this case, but stop
internal_assert(op->indices.size() == 1) << "modulus_remainder of vector\n";
result = ModulusRemainder{};
}
void ComputeModulusRemainder::visit(const VectorReduce *op) {
internal_assert(op->type.is_scalar()) << "modulus_remainder of vector\n";
result = ModulusRemainder{};
}
void ComputeModulusRemainder::visit(const LetStmt *) {
internal_error << "modulus_remainder of statement\n";
}
void ComputeModulusRemainder::visit(const AssertStmt *) {
internal_error << "modulus_remainder of statement\n";
}
void ComputeModulusRemainder::visit(const ProducerConsumer *) {
internal_error << "modulus_remainder of statement\n";
}
void ComputeModulusRemainder::visit(const For *) {
internal_error << "modulus_remainder of statement\n";
}
void ComputeModulusRemainder::visit(const Acquire *) {
internal_error << "modulus_remainder of statement\n";
}
void ComputeModulusRemainder::visit(const Store *) {
internal_error << "modulus_remainder of statement\n";
}
void ComputeModulusRemainder::visit(const Provide *) {
internal_error << "modulus_remainder of statement\n";
}
void ComputeModulusRemainder::visit(const Allocate *) {
internal_error << "modulus_remainder of statement\n";
}
void ComputeModulusRemainder::visit(const Realize *) {
internal_error << "modulus_remainder of statement\n";
}
void ComputeModulusRemainder::visit(const Block *) {
internal_error << "modulus_remainder of statement\n";
}
void ComputeModulusRemainder::visit(const Fork *) {
internal_error << "modulus_remainder of statement\n";
}
void ComputeModulusRemainder::visit(const Free *) {
internal_error << "modulus_remainder of statement\n";
}
void ComputeModulusRemainder::visit(const IfThenElse *) {
internal_error << "modulus_remainder of statement\n";
}
void ComputeModulusRemainder::visit(const Evaluate *) {
internal_error << "modulus_remainder of statement\n";
}
void ComputeModulusRemainder::visit(const Prefetch *) {
internal_error << "modulus_remainder of statement\n";
}
void ComputeModulusRemainder::visit(const Atomic *) {
internal_error << "modulus_remainder of statement\n";
}
} // namespace
ModulusRemainder modulus_remainder(const Expr &e) {
ComputeModulusRemainder mr(nullptr);
return mr.analyze(e);
}
ModulusRemainder modulus_remainder(const Expr &e, const Scope<ModulusRemainder> &scope) {
ComputeModulusRemainder mr(&scope);
return mr.analyze(e);
}
bool reduce_expr_modulo(const Expr &expr, int64_t modulus, int64_t *remainder) {
ModulusRemainder result = modulus_remainder(expr);
/* As an example: If we asked for expr mod 8, and the analysis
* said that expr = 16*k + 13, then because 16 % 8 == 0, the
* result is 13 % 8 == 5. But if the analysis says that expr =
* 6*k + 3, then expr mod 8 could be 1, 3, 5, or 7, so we just
* return false.
*/
if (mod(result.modulus, modulus) == 0) {
*remainder = mod(result.remainder, modulus);
return true;
} else {
return false;
}
}
bool reduce_expr_modulo(const Expr &expr, int64_t modulus, int64_t *remainder, const Scope<ModulusRemainder> &scope) {
ModulusRemainder result = modulus_remainder(expr, scope);
if (mod(result.modulus, modulus) == 0) {
*remainder = mod(result.remainder, modulus);
return true;
} else {
return false;
}
}
ModulusRemainder ComputeModulusRemainder::analyze(const Expr &e) {
e.accept(this);
return result;
}
namespace {
void check(const Expr &e, int64_t m, int64_t r) {
ModulusRemainder result = modulus_remainder(e);
if (result.modulus != m || result.remainder != r) {
std::cerr << "Test failed for modulus_remainder:\n";
std::cerr << "Expression: " << e << "\n";
std::cerr << "Correct modulus, remainder = " << m << ", " << r << "\n";
std::cerr << "Computed modulus, remainder = "
<< result.modulus << ", "
<< result.remainder << "\n";
exit(-1);
}
}
} // namespace
void modulus_remainder_test() {
Expr x = Variable::make(Int(32), "x");
Expr y = Variable::make(Int(32), "y");
check((30 * x + 3) + (40 * y + 2), 10, 5);
check((6 * x + 3) * (4 * y + 1), 2, 1);
check(max(30 * x - 24, 40 * y + 31), 5, 1);
check(10 * x - 33 * y, 1, 0);
check(10 * x - 35 * y, 5, 0);
check(123, 0, 123);
check(Let::make("y", x * 3 + 4, y * 3 + 4), 9, 7);
// Check overflow
check((5045320 * x + 4) * (405713 * y + 3) * (8000123 * x + 4354), 1, 0);
std::cout << "modulus_remainder test passed\n";
}
int64_t gcd(int64_t a, int64_t b) {
if (a < b) {
std::swap(a, b);
}
while (b != 0) {
int64_t tmp = b;
b = a % b;
a = tmp;
}
return a;
}
int64_t lcm(int64_t a, int64_t b) {
// Remove all of the common factors from one of the operands
b /= gcd(a, b);
// Then multiply
if (mul_would_overflow(64, a, b)) {
return 0;
} else {
return a * b;
}
}
ModulusRemainder operator+(const ModulusRemainder &a, const ModulusRemainder &b) {
if (add_would_overflow(64, a.remainder, b.remainder)) {
return {1, 0};
} else {
int64_t modulus = gcd(a.modulus, b.modulus);
int64_t remainder = mod(a.remainder + b.remainder, modulus);
return {modulus, remainder};
}
}
ModulusRemainder operator-(const ModulusRemainder &a, const ModulusRemainder &b) {
if (sub_would_overflow(64, a.remainder, b.remainder)) {
return {1, 0};
} else {
int64_t modulus = gcd(a.modulus, b.modulus);
int64_t remainder = mod(a.remainder - b.remainder, modulus);
return {modulus, remainder};
}
}
ModulusRemainder operator*(const ModulusRemainder &a, const ModulusRemainder &b) {
if (a.modulus == 0) {
// a is constant
if (!mul_would_overflow(64, a.remainder, b.modulus) && !mul_would_overflow(64, a.remainder, b.remainder)) {
return {a.remainder * b.modulus, a.remainder * b.remainder};
}
} else if (b.modulus == 0) {
// b is constant
if (!mul_would_overflow(64, b.remainder, a.modulus) && !mul_would_overflow(64, a.remainder, b.remainder)) {
return {b.remainder * a.modulus, a.remainder * b.remainder};
}
} else if (a.remainder == 0 && b.remainder == 0) {
// multiple times multiple
if (!mul_would_overflow(64, a.modulus, b.modulus)) {
return {a.modulus * b.modulus, 0};
}
} else if (a.remainder == 0) {
int64_t g = gcd(b.modulus, b.remainder);
if (!mul_would_overflow(64, a.modulus, g)) {
return {a.modulus * g, 0};
}
} else if (b.remainder == 0) {
int64_t g = gcd(a.modulus, a.remainder);
if (!mul_would_overflow(64, b.modulus, g)) {
return {b.modulus * g, 0};
}
} else {
// Convert them to the same modulus and multiply
if (!mul_would_overflow(64, a.remainder, b.remainder)) {
int64_t modulus = gcd(a.modulus, b.modulus);
int64_t remainder = mod(a.remainder * b.remainder, modulus);
return {modulus, remainder};
}
}
return ModulusRemainder{};
}
ModulusRemainder operator/(const ModulusRemainder &a, const ModulusRemainder &b) {
// What can we say about:
// floor((m1 * x + r1) / (m2 * y + r2))
// If m2 is zero and m1 is a multiple of r2, then we can pull the
// varying term out of the floor div and the expression simplifies
// to:
// (m1 / r2) * x + floor(r1 / r2)
// E.g. (8x + 3) / 2 -> (4x + 1)
if (b.modulus == 0 && b.remainder != 0) {
if (mod(a.modulus, b.remainder) == 0) {
return {a.modulus / b.remainder, div_imp(a.remainder, b.remainder)};
}
}
return ModulusRemainder{};
}
ModulusRemainder ModulusRemainder::unify(const ModulusRemainder &a, const ModulusRemainder &b) {
// We don't know if we're going to get a or b, so we'd better find
// a single modulus remainder that works for both.
// For example:
// max(30*_ + 13, 40*_ + 27) ->
// max(10*_ + 3, 10*_ + 7) ->
// max(2*_ + 1, 2*_ + 1) ->
// 2*_ + 1
if (b.remainder > a.remainder) {
return unify(b, a);
}
// Reduce them to the same modulus and the same remainder
int64_t modulus = gcd(a.modulus, b.modulus);
if (sub_would_overflow(64, a.remainder, b.remainder)) {
// The modulus is not representable as an int64.
return {0, 1};
}
int64_t diff = a.remainder - b.remainder;
modulus = gcd(diff, modulus);
int64_t ra = mod(a.remainder, modulus);
internal_assert(ra == mod(b.remainder, modulus))
<< "There's a bug inside ModulusRemainder in unify_alternatives:\n"
<< "a.modulus = " << a.modulus << "\n"
<< "a.remainder = " << a.remainder << "\n"
<< "b.modulus = " << b.modulus << "\n"
<< "b.remainder = " << b.remainder << "\n"
<< "diff = " << diff << "\n"
<< "unified modulus = " << modulus << "\n"
<< "unified remainder = " << ra << "\n";
return {modulus, ra};
}
ModulusRemainder ModulusRemainder::intersect(const ModulusRemainder &a, const ModulusRemainder &b) {
// We have x == ma * y + ra == mb * z + rb
// We want to synthesize these two facts into one modulus
// remainder relationship. We are permitted to be
// conservatively-large, so it's OK if some elements of the result
// only satisfy one of the two constraints.
// For coprime ma and mb you want to use the Chinese remainder
// theorem. In our case, the moduli will almost always be
// powers of two, so we should just return the smaller of the two
// sets (usually the one with the larger modulus).
if (a.modulus == 0) {
return a;
}
if (b.modulus == 0) {
return b;
}
if (a.modulus > b.modulus) {
return a;
}
return b;
}
void ComputeModulusRemainder::visit(const Mod *op) {
result = analyze(op->a) % analyze(op->b);
}
ModulusRemainder operator%(const ModulusRemainder &a, const ModulusRemainder &b) {
// For non-zero y, we can treat x mod y as x + z*y, where we know
// nothing about z.
// (ax + b) + z (cx + d) ->
// ax + b + zcx + dz ->
// gcd(a, c, d) * w + b
// E.g:
// (8x + 5) mod (6x + 2) ->
// (8x + 5) + z (6x + 2) ->
// (8x + 6zx + 2x) + 5 ->
// 2(4x + 3zx + x) + 5 ->
// 2w + 1
int64_t modulus = gcd(a.modulus, b.modulus);
modulus = gcd(modulus, b.remainder);
int64_t remainder = mod(a.remainder, modulus);
if (b.remainder == 0 && remainder != 0) {
// b could be zero, so the result could also just be zero.
if (modulus == 0) {
remainder = 0;
} else {
// This can no longer be expressed as ax + b
remainder = 0;
modulus = 1;
}
}
return {modulus, remainder};
}
ModulusRemainder operator+(const ModulusRemainder &a, int64_t b) {
return a + ModulusRemainder(0, b);
}
ModulusRemainder operator-(const ModulusRemainder &a, int64_t b) {
return a - ModulusRemainder(0, b);
}
ModulusRemainder operator*(const ModulusRemainder &a, int64_t b) {
return a * ModulusRemainder(0, b);
}
ModulusRemainder operator/(const ModulusRemainder &a, int64_t b) {
return a / ModulusRemainder(0, b);
}
ModulusRemainder operator%(const ModulusRemainder &a, int64_t b) {
return a % ModulusRemainder(0, b);
}
} // namespace Internal
} // namespace Halide