MLIR  17.0.0git
PolynomialApproximation.cpp
Go to the documentation of this file.
1 //===- PolynomialApproximation.cpp - Approximate math operations ----------===//
2 //
3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4 // See https://llvm.org/LICENSE.txt for license information.
5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6 //
7 //===----------------------------------------------------------------------===//
8 //
9 // This file implements expansion of math operations to fast approximations
10 // that do not rely on any of the library functions.
11 //
12 //===----------------------------------------------------------------------===//
13 
14 #include <climits>
15 #include <cstddef>
16 
25 #include "mlir/IR/Builders.h"
26 #include "mlir/IR/BuiltinTypes.h"
28 #include "mlir/IR/OpDefinition.h"
29 #include "mlir/IR/PatternMatch.h"
30 #include "mlir/IR/TypeUtilities.h"
33 #include "llvm/ADT/ArrayRef.h"
34 #include "llvm/ADT/STLExtras.h"
35 
36 using namespace mlir;
37 using namespace mlir::math;
38 using namespace mlir::vector;
39 
40 // Returns vector shape if the type is a vector. Returns an empty shape if it is
41 // not a vector.
43  auto vectorType = type.dyn_cast<VectorType>();
44  return vectorType ? vectorType.getShape() : ArrayRef<int64_t>();
45 }
46 
48  return vectorShape(value.getType());
49 }
50 
51 //----------------------------------------------------------------------------//
52 // Broadcast scalar types and values into vector types and values.
53 //----------------------------------------------------------------------------//
54 
55 // Broadcasts scalar type into vector type (iff shape is non-scalar).
56 static Type broadcast(Type type, ArrayRef<int64_t> shape) {
57  assert(!type.isa<VectorType>() && "must be scalar type");
58  return !shape.empty() ? VectorType::get(shape, type) : type;
59 }
60 
61 // Broadcasts scalar value into vector (iff shape is non-scalar).
62 static Value broadcast(ImplicitLocOpBuilder &builder, Value value,
63  ArrayRef<int64_t> shape) {
64  assert(!value.getType().isa<VectorType>() && "must be scalar value");
65  auto type = broadcast(value.getType(), shape);
66  return !shape.empty() ? builder.create<BroadcastOp>(type, value) : value;
67 }
68 
69 //----------------------------------------------------------------------------//
70 // Helper function to handle n-D vectors with 1-D operations.
71 //----------------------------------------------------------------------------//
72 
73 // Expands and unrolls n-D vector operands into multiple fixed size 1-D vectors
74 // and calls the compute function with 1-D vector operands. Stitches back all
75 // results into the original n-D vector result.
76 //
77 // Examples: vectorWidth = 8
78 // - vector<4x8xf32> unrolled 4 times
79 // - vector<16xf32> expanded to vector<2x8xf32> and unrolled 2 times
80 // - vector<4x16xf32> expanded to vector<4x2x8xf32> and unrolled 4*2 times
81 //
82 // Some math approximations rely on ISA-specific operations that only accept
83 // fixed size 1-D vectors (e.g. AVX expects vectors of width 8).
84 //
85 // It is the caller's responsibility to verify that the inner dimension is
86 // divisible by the vectorWidth, and that all operands have the same vector
87 // shape.
88 static Value
90  ValueRange operands, int64_t vectorWidth,
92  assert(!operands.empty() && "operands must be not empty");
93  assert(vectorWidth > 0 && "vector width must be larger than 0");
94 
95  VectorType inputType = operands[0].getType().cast<VectorType>();
96  ArrayRef<int64_t> inputShape = inputType.getShape();
97 
98  // If input shape matches target vector width, we can just call the
99  // user-provided compute function with the operands.
100  if (inputShape == llvm::ArrayRef(vectorWidth))
101  return compute(operands);
102 
103  // Check if the inner dimension has to be expanded, or we can directly iterate
104  // over the outer dimensions of the vector.
105  int64_t innerDim = inputShape.back();
106  int64_t expansionDim = innerDim / vectorWidth;
107  assert((innerDim % vectorWidth == 0) && "invalid inner dimension size");
108 
109  // Maybe expand operands to the higher rank vector shape that we'll use to
110  // iterate over and extract one dimensional vectors.
111  SmallVector<int64_t> expandedShape(inputShape.begin(), inputShape.end());
112  SmallVector<Value> expandedOperands(operands);
113 
114  if (expansionDim > 1) {
115  // Expand shape from [..., innerDim] to [..., expansionDim, vectorWidth].
116  expandedShape.insert(expandedShape.end() - 1, expansionDim);
117  expandedShape.back() = vectorWidth;
118 
119  for (unsigned i = 0; i < operands.size(); ++i) {
120  auto operand = operands[i];
121  auto eltType = operand.getType().cast<VectorType>().getElementType();
122  auto expandedType = VectorType::get(expandedShape, eltType);
123  expandedOperands[i] =
124  builder.create<vector::ShapeCastOp>(expandedType, operand);
125  }
126  }
127 
128  // Iterate over all outer dimensions of the compute shape vector type.
129  auto iterationDims = ArrayRef<int64_t>(expandedShape).drop_back();
130  int64_t maxIndex = computeMaxLinearIndex(iterationDims);
131  auto strides = computeStrides(iterationDims);
132 
133  // Compute results for each one dimensional vector.
134  SmallVector<Value> results(maxIndex);
135 
136  for (int64_t i = 0; i < maxIndex; ++i) {
137  auto offsets = delinearize(i, strides);
138 
139  SmallVector<Value> extracted(expandedOperands.size());
140  for (const auto &tuple : llvm::enumerate(expandedOperands))
141  extracted[tuple.index()] =
142  builder.create<vector::ExtractOp>(tuple.value(), offsets);
143 
144  results[i] = compute(extracted);
145  }
146 
147  // Stitch results together into one large vector.
148  Type resultEltType = results[0].getType().cast<VectorType>().getElementType();
149  Type resultExpandedType = VectorType::get(expandedShape, resultEltType);
150  Value result = builder.create<arith::ConstantOp>(
151  resultExpandedType, builder.getZeroAttr(resultExpandedType));
152 
153  for (int64_t i = 0; i < maxIndex; ++i)
154  result = builder.create<vector::InsertOp>(results[i], result,
155  delinearize(i, strides));
156 
157  // Reshape back to the original vector shape.
158  return builder.create<vector::ShapeCastOp>(
159  VectorType::get(inputShape, resultEltType), result);
160 }
161 
162 //----------------------------------------------------------------------------//
163 // Helper functions to create constants.
164 //----------------------------------------------------------------------------//
165 
166 static Value floatCst(ImplicitLocOpBuilder &builder, float value,
167  Type elementType) {
168  assert((elementType.isF16() || elementType.isF32()) &&
169  "x must be f16 or f32 type.");
170  return builder.create<arith::ConstantOp>(
171  builder.getFloatAttr(elementType, value));
172 }
173 
174 static Value f32Cst(ImplicitLocOpBuilder &builder, float value) {
175  return builder.create<arith::ConstantOp>(builder.getF32FloatAttr(value));
176 }
177 
178 static Value i32Cst(ImplicitLocOpBuilder &builder, int32_t value) {
179  return builder.create<arith::ConstantOp>(builder.getI32IntegerAttr(value));
180 }
181 
182 static Value f32FromBits(ImplicitLocOpBuilder &builder, uint32_t bits) {
183  Value i32Value = i32Cst(builder, static_cast<int32_t>(bits));
184  return builder.create<arith::BitcastOp>(builder.getF32Type(), i32Value);
185 }
186 
187 //----------------------------------------------------------------------------//
188 // Helper functions to build math functions approximations.
189 //----------------------------------------------------------------------------//
190 
191 // Return the minimum of the two values or NaN if value is NaN
192 static Value min(ImplicitLocOpBuilder &builder, Value value, Value bound) {
193  return builder.create<arith::SelectOp>(
194  builder.create<arith::CmpFOp>(arith::CmpFPredicate::ULT, value, bound),
195  value, bound);
196 }
197 
198 // Return the maximum of the two values or NaN if value is NaN
199 static Value max(ImplicitLocOpBuilder &builder, Value value, Value bound) {
200  return builder.create<arith::SelectOp>(
201  builder.create<arith::CmpFOp>(arith::CmpFPredicate::UGT, value, bound),
202  value, bound);
203 }
204 
205 // Return the clamped value or NaN if value is NaN
206 static Value clamp(ImplicitLocOpBuilder &builder, Value value, Value lowerBound,
207  Value upperBound) {
208  return max(builder, min(builder, value, upperBound), lowerBound);
209 }
210 
211 // Decomposes given floating point value `arg` into a normalized fraction and
212 // an integral power of two (see std::frexp). Returned values have float type.
213 static std::pair<Value, Value> frexp(ImplicitLocOpBuilder &builder, Value arg,
214  bool isPositive = false) {
215  assert(getElementTypeOrSelf(arg).isF32() && "arg must be f32 type");
216  ArrayRef<int64_t> shape = vectorShape(arg);
217 
218  auto bcast = [&](Value value) -> Value {
219  return broadcast(builder, value, shape);
220  };
221 
222  auto i32 = builder.getIntegerType(32);
223  auto i32Vec = broadcast(i32, shape);
224  auto f32Vec = broadcast(builder.getF32Type(), shape);
225 
226  Value cst126f = f32Cst(builder, 126.0f);
227  Value cstHalf = f32Cst(builder, 0.5f);
228  Value cstInvMantMask = f32FromBits(builder, ~0x7f800000u);
229 
230  // Bitcast to i32 for bitwise operations.
231  Value i32Half = builder.create<arith::BitcastOp>(i32, cstHalf);
232  Value i32InvMantMask = builder.create<arith::BitcastOp>(i32, cstInvMantMask);
233  Value i32Arg = builder.create<arith::BitcastOp>(i32Vec, arg);
234 
235  // Compute normalized fraction.
236  Value tmp0 = builder.create<arith::AndIOp>(i32Arg, bcast(i32InvMantMask));
237  Value tmp1 = builder.create<arith::OrIOp>(tmp0, bcast(i32Half));
238  Value normalizedFraction = builder.create<arith::BitcastOp>(f32Vec, tmp1);
239 
240  // Compute exponent.
241  Value arg0 = isPositive ? arg : builder.create<math::AbsFOp>(arg);
242  Value biasedExponentBits = builder.create<arith::ShRUIOp>(
243  builder.create<arith::BitcastOp>(i32Vec, arg0),
244  bcast(i32Cst(builder, 23)));
245  Value biasedExponent =
246  builder.create<arith::SIToFPOp>(f32Vec, biasedExponentBits);
247  Value exponent =
248  builder.create<arith::SubFOp>(biasedExponent, bcast(cst126f));
249 
250  return {normalizedFraction, exponent};
251 }
252 
253 // Computes exp2 for an i32 argument.
254 static Value exp2I32(ImplicitLocOpBuilder &builder, Value arg) {
255  assert(getElementTypeOrSelf(arg).isInteger(32) && "arg must be i32 type");
256  ArrayRef<int64_t> shape = vectorShape(arg);
257 
258  auto bcast = [&](Value value) -> Value {
259  return broadcast(builder, value, shape);
260  };
261 
262  auto f32Vec = broadcast(builder.getF32Type(), shape);
263  // The exponent of f32 located at 23-bit.
264  auto exponetBitLocation = bcast(i32Cst(builder, 23));
265  // Set the exponent bias to zero.
266  auto bias = bcast(i32Cst(builder, 127));
267 
268  Value biasedArg = builder.create<arith::AddIOp>(arg, bias);
269  Value exp2ValueInt =
270  builder.create<arith::ShLIOp>(biasedArg, exponetBitLocation);
271  Value exp2ValueF32 = builder.create<arith::BitcastOp>(f32Vec, exp2ValueInt);
272 
273  return exp2ValueF32;
274 }
275 
276 namespace {
277 Value makePolynomialCalculation(ImplicitLocOpBuilder &builder,
278  llvm::ArrayRef<Value> coeffs, Value x) {
279  Type elementType = getElementTypeOrSelf(x);
280  assert((elementType.isF32() || elementType.isF16()) &&
281  "x must be f32 or f16 type");
282  ArrayRef<int64_t> shape = vectorShape(x);
283 
284  if (coeffs.empty())
285  return broadcast(builder, floatCst(builder, 0.0f, elementType), shape);
286 
287  if (coeffs.size() == 1)
288  return coeffs[0];
289 
290  Value res = builder.create<math::FmaOp>(x, coeffs[coeffs.size() - 1],
291  coeffs[coeffs.size() - 2]);
292  for (auto i = ptrdiff_t(coeffs.size()) - 3; i >= 0; --i) {
293  res = builder.create<math::FmaOp>(x, res, coeffs[i]);
294  }
295  return res;
296 }
297 } // namespace
298 
299 //----------------------------------------------------------------------------//
300 // Helper function/pattern to insert casts for reusing F32 bit expansion.
301 //----------------------------------------------------------------------------//
302 
303 template <typename T>
305  // Conservatively only allow where the operand and result types are exactly 1.
306  Type origType = op->getResultTypes().front();
307  for (Type t : llvm::drop_begin(op->getResultTypes()))
308  if (origType != t)
309  return rewriter.notifyMatchFailure(op, "required all types to match");
310  for (Type t : op->getOperandTypes())
311  if (origType != t)
312  return rewriter.notifyMatchFailure(op, "required all types to match");
313 
314  // Skip if already F32 or larger than 32 bits.
315  if (getElementTypeOrSelf(origType).isF32() ||
316  getElementTypeOrSelf(origType).getIntOrFloatBitWidth() > 32)
317  return failure();
318 
319  // Create F32 equivalent type.
320  Type newType;
321  if (auto shaped = origType.dyn_cast<ShapedType>()) {
322  newType = shaped.clone(rewriter.getF32Type());
323  } else if (origType.isa<FloatType>()) {
324  newType = rewriter.getF32Type();
325  } else {
326  return rewriter.notifyMatchFailure(op,
327  "unable to find F32 equivalent type");
328  }
329 
330  Location loc = op->getLoc();
331  SmallVector<Value> operands;
332  for (auto operand : op->getOperands())
333  operands.push_back(rewriter.create<arith::ExtFOp>(loc, newType, operand));
334  auto result =
335  rewriter.create<T>(loc, TypeRange{newType}, operands, op->getAttrs());
336  rewriter.replaceOpWithNewOp<arith::TruncFOp>(op, origType, result);
337  return success();
338 }
339 
340 namespace {
341 // Pattern to cast to F32 to reuse F32 expansion as fallback for single-result
342 // op.
343 // TODO: Consider revising to avoid adding multiple casts for a subgraph that is
344 // all in lower precision. Currently this is only fallback support and performs
345 // simplistic casting.
346 template <typename T>
347 struct ReuseF32Expansion : public OpRewritePattern<T> {
348 public:
350  LogicalResult matchAndRewrite(T op, PatternRewriter &rewriter) const final {
351  static_assert(
352  T::template hasTrait<mlir::OpTrait::SameOperandsAndResultType>(),
353  "requires same operands and result types");
354  return insertCasts<T>(op, rewriter);
355  }
356 };
357 } // namespace
358 
359 //----------------------------------------------------------------------------//
360 // AtanOp approximation.
361 //----------------------------------------------------------------------------//
362 
363 namespace {
364 struct AtanApproximation : public OpRewritePattern<math::AtanOp> {
365 public:
367 
368  LogicalResult matchAndRewrite(math::AtanOp op,
369  PatternRewriter &rewriter) const final;
370 };
371 } // namespace
372 
374 AtanApproximation::matchAndRewrite(math::AtanOp op,
375  PatternRewriter &rewriter) const {
376  auto operand = op.getOperand();
377  if (!getElementTypeOrSelf(operand).isF32())
378  return rewriter.notifyMatchFailure(op, "unsupported operand type");
379 
380  ArrayRef<int64_t> shape = vectorShape(op.getOperand());
381 
382  ImplicitLocOpBuilder builder(op->getLoc(), rewriter);
383  auto one = broadcast(builder, f32Cst(builder, 1.0f), shape);
384 
385  // Remap the problem over [0.0, 1.0] by looking at the absolute value and the
386  // handling symmetry.
387  Value abs = builder.create<math::AbsFOp>(operand);
388  Value reciprocal = builder.create<arith::DivFOp>(one, abs);
389  Value compare =
390  builder.create<arith::CmpFOp>(arith::CmpFPredicate::OLT, abs, reciprocal);
391  Value x = builder.create<arith::SelectOp>(compare, abs, reciprocal);
392 
393  // Perform the Taylor series approximation for atan over the range
394  // [-1.0, 1.0].
395  auto n1 = broadcast(builder, f32Cst(builder, 0.14418283f), shape);
396  auto n2 = broadcast(builder, f32Cst(builder, -0.34999234f), shape);
397  auto n3 = broadcast(builder, f32Cst(builder, -0.01067831f), shape);
398  auto n4 = broadcast(builder, f32Cst(builder, 1.00209986f), shape);
399 
400  Value p = builder.create<math::FmaOp>(x, n1, n2);
401  p = builder.create<math::FmaOp>(x, p, n3);
402  p = builder.create<math::FmaOp>(x, p, n4);
403  p = builder.create<arith::MulFOp>(x, p);
404 
405  // Remap the solution for over [0.0, 1.0] to [0.0, inf]
406  auto halfPi = broadcast(builder, f32Cst(builder, 1.57079632679f), shape);
407  Value sub = builder.create<arith::SubFOp>(halfPi, p);
408  Value select = builder.create<arith::SelectOp>(compare, p, sub);
409 
410  // Correct for signing of the input.
411  rewriter.replaceOpWithNewOp<math::CopySignOp>(op, select, operand);
412  return success();
413 }
414 
415 //----------------------------------------------------------------------------//
416 // AtanOp approximation.
417 //----------------------------------------------------------------------------//
418 
419 namespace {
420 struct Atan2Approximation : public OpRewritePattern<math::Atan2Op> {
421 public:
423 
424  LogicalResult matchAndRewrite(math::Atan2Op op,
425  PatternRewriter &rewriter) const final;
426 };
427 } // namespace
428 
430 Atan2Approximation::matchAndRewrite(math::Atan2Op op,
431  PatternRewriter &rewriter) const {
432  auto y = op.getOperand(0);
433  auto x = op.getOperand(1);
434  if (!getElementTypeOrSelf(x).isF32())
435  return rewriter.notifyMatchFailure(op, "unsupported operand type");
436 
437  ImplicitLocOpBuilder builder(op->getLoc(), rewriter);
438  ArrayRef<int64_t> shape = vectorShape(op.getResult());
439 
440  // Compute atan in the valid range.
441  auto div = builder.create<arith::DivFOp>(y, x);
442  auto atan = builder.create<math::AtanOp>(div);
443 
444  // Determine what the atan would be for a 180 degree rotation.
445  auto zero = broadcast(builder, f32Cst(builder, 0.0f), shape);
446  auto pi = broadcast(builder, f32Cst(builder, 3.14159265359f), shape);
447  auto addPi = builder.create<arith::AddFOp>(atan, pi);
448  auto subPi = builder.create<arith::SubFOp>(atan, pi);
449  auto atanGt =
450  builder.create<arith::CmpFOp>(arith::CmpFPredicate::OGT, atan, zero);
451  auto flippedAtan = builder.create<arith::SelectOp>(atanGt, subPi, addPi);
452 
453  // Determine whether to directly use atan or use the 180 degree flip
454  auto xGt = builder.create<arith::CmpFOp>(arith::CmpFPredicate::OGT, x, zero);
455  Value result = builder.create<arith::SelectOp>(xGt, atan, flippedAtan);
456 
457  // Handle x = 0, y > 0
458  Value xZero =
459  builder.create<arith::CmpFOp>(arith::CmpFPredicate::OEQ, x, zero);
460  Value yGt = builder.create<arith::CmpFOp>(arith::CmpFPredicate::OGT, y, zero);
461  Value isHalfPi = builder.create<arith::AndIOp>(xZero, yGt);
462  auto halfPi = broadcast(builder, f32Cst(builder, 1.57079632679f), shape);
463  result = builder.create<arith::SelectOp>(isHalfPi, halfPi, result);
464 
465  // Handle x = 0, y < 0
466  Value yLt = builder.create<arith::CmpFOp>(arith::CmpFPredicate::OLT, y, zero);
467  Value isNegativeHalfPiPi = builder.create<arith::AndIOp>(xZero, yLt);
468  auto negativeHalfPiPi =
469  broadcast(builder, f32Cst(builder, -1.57079632679f), shape);
470  result = builder.create<arith::SelectOp>(isNegativeHalfPiPi, negativeHalfPiPi,
471  result);
472 
473  // Handle x = 0, y = 0;
474  Value yZero =
475  builder.create<arith::CmpFOp>(arith::CmpFPredicate::OEQ, y, zero);
476  Value isNan = builder.create<arith::AndIOp>(xZero, yZero);
477  Value cstNan = broadcast(builder, f32FromBits(builder, 0x7fc00000), shape);
478  result = builder.create<arith::SelectOp>(isNan, cstNan, result);
479 
480  rewriter.replaceOp(op, result);
481  return success();
482 }
483 
484 //----------------------------------------------------------------------------//
485 // TanhOp approximation.
486 //----------------------------------------------------------------------------//
487 
488 namespace {
489 struct TanhApproximation : public OpRewritePattern<math::TanhOp> {
490 public:
492 
493  LogicalResult matchAndRewrite(math::TanhOp op,
494  PatternRewriter &rewriter) const final;
495 };
496 } // namespace
497 
499 TanhApproximation::matchAndRewrite(math::TanhOp op,
500  PatternRewriter &rewriter) const {
501  if (!getElementTypeOrSelf(op.getOperand()).isF32())
502  return rewriter.notifyMatchFailure(op, "unsupported operand type");
503 
504  ArrayRef<int64_t> shape = vectorShape(op.getOperand());
505 
506  ImplicitLocOpBuilder builder(op->getLoc(), rewriter);
507  auto bcast = [&](Value value) -> Value {
508  return broadcast(builder, value, shape);
509  };
510 
511  // Clamp operand into [plusClamp, minusClamp] range.
512  Value minusClamp = bcast(f32Cst(builder, -7.99881172180175781f));
513  Value plusClamp = bcast(f32Cst(builder, 7.99881172180175781f));
514  Value x = clamp(builder, op.getOperand(), minusClamp, plusClamp);
515 
516  // Mask for tiny values that are approximated with `operand`.
517  Value tiny = bcast(f32Cst(builder, 0.0004f));
518  Value tinyMask = builder.create<arith::CmpFOp>(
519  arith::CmpFPredicate::OLT, builder.create<math::AbsFOp>(op.getOperand()),
520  tiny);
521 
522  // The monomial coefficients of the numerator polynomial (odd).
523  Value alpha1 = bcast(f32Cst(builder, 4.89352455891786e-03f));
524  Value alpha3 = bcast(f32Cst(builder, 6.37261928875436e-04f));
525  Value alpha5 = bcast(f32Cst(builder, 1.48572235717979e-05f));
526  Value alpha7 = bcast(f32Cst(builder, 5.12229709037114e-08f));
527  Value alpha9 = bcast(f32Cst(builder, -8.60467152213735e-11f));
528  Value alpha11 = bcast(f32Cst(builder, 2.00018790482477e-13f));
529  Value alpha13 = bcast(f32Cst(builder, -2.76076847742355e-16f));
530 
531  // The monomial coefficients of the denominator polynomial (even).
532  Value beta0 = bcast(f32Cst(builder, 4.89352518554385e-03f));
533  Value beta2 = bcast(f32Cst(builder, 2.26843463243900e-03f));
534  Value beta4 = bcast(f32Cst(builder, 1.18534705686654e-04f));
535  Value beta6 = bcast(f32Cst(builder, 1.19825839466702e-06f));
536 
537  // Since the polynomials are odd/even, we need x^2.
538  Value x2 = builder.create<arith::MulFOp>(x, x);
539 
540  // Evaluate the numerator polynomial p.
541  Value p = builder.create<math::FmaOp>(x2, alpha13, alpha11);
542  p = builder.create<math::FmaOp>(x2, p, alpha9);
543  p = builder.create<math::FmaOp>(x2, p, alpha7);
544  p = builder.create<math::FmaOp>(x2, p, alpha5);
545  p = builder.create<math::FmaOp>(x2, p, alpha3);
546  p = builder.create<math::FmaOp>(x2, p, alpha1);
547  p = builder.create<arith::MulFOp>(x, p);
548 
549  // Evaluate the denominator polynomial q.
550  Value q = builder.create<math::FmaOp>(x2, beta6, beta4);
551  q = builder.create<math::FmaOp>(x2, q, beta2);
552  q = builder.create<math::FmaOp>(x2, q, beta0);
553 
554  // Divide the numerator by the denominator.
555  Value res = builder.create<arith::SelectOp>(
556  tinyMask, x, builder.create<arith::DivFOp>(p, q));
557 
558  rewriter.replaceOp(op, res);
559 
560  return success();
561 }
562 
563 #define LN2_VALUE \
564  0.693147180559945309417232121458176568075500134360255254120680009493393621L
565 #define LOG2E_VALUE \
566  1.442695040888963407359924681001892137426645954152985934135449406931109219L
567 
568 //----------------------------------------------------------------------------//
569 // LogOp and Log2Op approximation.
570 //----------------------------------------------------------------------------//
571 
572 namespace {
573 template <typename Op>
574 struct LogApproximationBase : public OpRewritePattern<Op> {
576 
577  /// Base 2 if 'base2' is set; natural logarithm (base e) otherwise.
578  LogicalResult logMatchAndRewrite(Op op, PatternRewriter &rewriter,
579  bool base2) const;
580 };
581 } // namespace
582 
583 // This approximation comes from Julien Pommier's SSE math library.
584 // Link: http://gruntthepeon.free.fr/ssemath
585 template <typename Op>
587 LogApproximationBase<Op>::logMatchAndRewrite(Op op, PatternRewriter &rewriter,
588  bool base2) const {
589  if (!getElementTypeOrSelf(op.getOperand()).isF32())
590  return rewriter.notifyMatchFailure(op, "unsupported operand type");
591 
592  ArrayRef<int64_t> shape = vectorShape(op.getOperand());
593 
594  ImplicitLocOpBuilder builder(op->getLoc(), rewriter);
595  auto bcast = [&](Value value) -> Value {
596  return broadcast(builder, value, shape);
597  };
598 
599  Value cstZero = bcast(f32Cst(builder, 0.0f));
600  Value cstOne = bcast(f32Cst(builder, 1.0f));
601  Value cstNegHalf = bcast(f32Cst(builder, -0.5f));
602 
603  // The smallest non denormalized float number.
604  Value cstMinNormPos = bcast(f32FromBits(builder, 0x00800000u));
605  Value cstMinusInf = bcast(f32FromBits(builder, 0xff800000u));
606  Value cstPosInf = bcast(f32FromBits(builder, 0x7f800000u));
607  Value cstNan = bcast(f32FromBits(builder, 0x7fc00000));
608 
609  // Polynomial coefficients.
610  Value cstCephesSQRTHF = bcast(f32Cst(builder, 0.707106781186547524f));
611  Value cstCephesLogP0 = bcast(f32Cst(builder, 7.0376836292E-2f));
612  Value cstCephesLogP1 = bcast(f32Cst(builder, -1.1514610310E-1f));
613  Value cstCephesLogP2 = bcast(f32Cst(builder, 1.1676998740E-1f));
614  Value cstCephesLogP3 = bcast(f32Cst(builder, -1.2420140846E-1f));
615  Value cstCephesLogP4 = bcast(f32Cst(builder, +1.4249322787E-1f));
616  Value cstCephesLogP5 = bcast(f32Cst(builder, -1.6668057665E-1f));
617  Value cstCephesLogP6 = bcast(f32Cst(builder, +2.0000714765E-1f));
618  Value cstCephesLogP7 = bcast(f32Cst(builder, -2.4999993993E-1f));
619  Value cstCephesLogP8 = bcast(f32Cst(builder, +3.3333331174E-1f));
620 
621  Value x = op.getOperand();
622 
623  // Truncate input values to the minimum positive normal.
624  x = max(builder, x, cstMinNormPos);
625 
626  // Extract significant in the range [0.5,1) and exponent.
627  std::pair<Value, Value> pair = frexp(builder, x, /*isPositive=*/true);
628  x = pair.first;
629  Value e = pair.second;
630 
631  // Shift the inputs from the range [0.5,1) to [sqrt(1/2), sqrt(2)) and shift
632  // by -1.0. The values are then centered around 0, which improves the
633  // stability of the polynomial evaluation:
634  //
635  // if( x < SQRTHF ) {
636  // e -= 1;
637  // x = x + x - 1.0;
638  // } else { x = x - 1.0; }
639  Value mask = builder.create<arith::CmpFOp>(arith::CmpFPredicate::OLT, x,
640  cstCephesSQRTHF);
641  Value tmp = builder.create<arith::SelectOp>(mask, x, cstZero);
642 
643  x = builder.create<arith::SubFOp>(x, cstOne);
644  e = builder.create<arith::SubFOp>(
645  e, builder.create<arith::SelectOp>(mask, cstOne, cstZero));
646  x = builder.create<arith::AddFOp>(x, tmp);
647 
648  Value x2 = builder.create<arith::MulFOp>(x, x);
649  Value x3 = builder.create<arith::MulFOp>(x2, x);
650 
651  // Evaluate the polynomial approximant of degree 8 in three parts.
652  Value y0, y1, y2;
653  y0 = builder.create<math::FmaOp>(cstCephesLogP0, x, cstCephesLogP1);
654  y1 = builder.create<math::FmaOp>(cstCephesLogP3, x, cstCephesLogP4);
655  y2 = builder.create<math::FmaOp>(cstCephesLogP6, x, cstCephesLogP7);
656  y0 = builder.create<math::FmaOp>(y0, x, cstCephesLogP2);
657  y1 = builder.create<math::FmaOp>(y1, x, cstCephesLogP5);
658  y2 = builder.create<math::FmaOp>(y2, x, cstCephesLogP8);
659  y0 = builder.create<math::FmaOp>(y0, x3, y1);
660  y0 = builder.create<math::FmaOp>(y0, x3, y2);
661  y0 = builder.create<arith::MulFOp>(y0, x3);
662 
663  y0 = builder.create<math::FmaOp>(cstNegHalf, x2, y0);
664  x = builder.create<arith::AddFOp>(x, y0);
665 
666  if (base2) {
667  Value cstLog2e = bcast(f32Cst(builder, static_cast<float>(LOG2E_VALUE)));
668  x = builder.create<math::FmaOp>(x, cstLog2e, e);
669  } else {
670  Value cstLn2 = bcast(f32Cst(builder, static_cast<float>(LN2_VALUE)));
671  x = builder.create<math::FmaOp>(e, cstLn2, x);
672  }
673 
674  Value invalidMask = builder.create<arith::CmpFOp>(arith::CmpFPredicate::ULT,
675  op.getOperand(), cstZero);
676  Value zeroMask = builder.create<arith::CmpFOp>(arith::CmpFPredicate::OEQ,
677  op.getOperand(), cstZero);
678  Value posInfMask = builder.create<arith::CmpFOp>(arith::CmpFPredicate::OEQ,
679  op.getOperand(), cstPosInf);
680 
681  // Filter out invalid values:
682  // • x == 0 -> -INF
683  // • x < 0 -> NAN
684  // • x == +INF -> +INF
685  Value aproximation = builder.create<arith::SelectOp>(
686  zeroMask, cstMinusInf,
687  builder.create<arith::SelectOp>(
688  invalidMask, cstNan,
689  builder.create<arith::SelectOp>(posInfMask, cstPosInf, x)));
690 
691  rewriter.replaceOp(op, aproximation);
692 
693  return success();
694 }
695 
696 namespace {
697 struct LogApproximation : public LogApproximationBase<math::LogOp> {
698  using LogApproximationBase::LogApproximationBase;
699 
700  LogicalResult matchAndRewrite(math::LogOp op,
701  PatternRewriter &rewriter) const final {
702  return logMatchAndRewrite(op, rewriter, /*base2=*/false);
703  }
704 };
705 } // namespace
706 
707 namespace {
708 struct Log2Approximation : public LogApproximationBase<math::Log2Op> {
709  using LogApproximationBase::LogApproximationBase;
710 
711  LogicalResult matchAndRewrite(math::Log2Op op,
712  PatternRewriter &rewriter) const final {
713  return logMatchAndRewrite(op, rewriter, /*base2=*/true);
714  }
715 };
716 } // namespace
717 
718 //----------------------------------------------------------------------------//
719 // Log1p approximation.
720 //----------------------------------------------------------------------------//
721 
722 namespace {
723 struct Log1pApproximation : public OpRewritePattern<math::Log1pOp> {
724 public:
726 
727  LogicalResult matchAndRewrite(math::Log1pOp op,
728  PatternRewriter &rewriter) const final;
729 };
730 } // namespace
731 
732 // Approximate log(1+x).
734 Log1pApproximation::matchAndRewrite(math::Log1pOp op,
735  PatternRewriter &rewriter) const {
736  if (!getElementTypeOrSelf(op.getOperand()).isF32())
737  return rewriter.notifyMatchFailure(op, "unsupported operand type");
738 
739  ArrayRef<int64_t> shape = vectorShape(op.getOperand());
740 
741  ImplicitLocOpBuilder builder(op->getLoc(), rewriter);
742  auto bcast = [&](Value value) -> Value {
743  return broadcast(builder, value, shape);
744  };
745 
746  // Approximate log(1+x) using the following, due to W. Kahan:
747  // u = x + 1.0;
748  // if (u == 1.0 || u == inf) return x;
749  // return x * log(u) / (u - 1.0);
750  // ^^^^^^^^^^^^^^^^^^^^^^
751  // "logLarge" below.
752  Value cstOne = bcast(f32Cst(builder, 1.0f));
753  Value x = op.getOperand();
754  Value u = builder.create<arith::AddFOp>(x, cstOne);
755  Value uSmall =
756  builder.create<arith::CmpFOp>(arith::CmpFPredicate::OEQ, u, cstOne);
757  Value logU = builder.create<math::LogOp>(u);
758  Value uInf =
759  builder.create<arith::CmpFOp>(arith::CmpFPredicate::OEQ, u, logU);
760  Value logLarge = builder.create<arith::MulFOp>(
761  x, builder.create<arith::DivFOp>(
762  logU, builder.create<arith::SubFOp>(u, cstOne)));
763  Value approximation = builder.create<arith::SelectOp>(
764  builder.create<arith::OrIOp>(uSmall, uInf), x, logLarge);
765  rewriter.replaceOp(op, approximation);
766  return success();
767 }
768 
769 //----------------------------------------------------------------------------//
770 // Erf approximation.
771 //----------------------------------------------------------------------------//
772 
773 // Approximates erf(x) with
774 // a - P(x)/Q(x)
775 // where P and Q are polynomials of degree 4.
776 // Different coefficients are chosen based on the value of x.
777 // The approximation error is ~2.5e-07.
778 // Boost's minimax tool that utilizes the Remez method was used to find the
779 // coefficients.
782  PatternRewriter &rewriter) const {
783  Value operand = op.getOperand();
784  Type elementType = getElementTypeOrSelf(operand);
785 
786  if (!(elementType.isF32() || elementType.isF16()))
787  return rewriter.notifyMatchFailure(op,
788  "only f32 and f16 type is supported.");
789  ArrayRef<int64_t> shape = vectorShape(operand);
790 
791  ImplicitLocOpBuilder builder(op->getLoc(), rewriter);
792  auto bcast = [&](Value value) -> Value {
793  return broadcast(builder, value, shape);
794  };
795 
796  const int intervalsCount = 3;
797  const int polyDegree = 4;
798 
799  Value zero = bcast(floatCst(builder, 0, elementType));
800  Value one = bcast(floatCst(builder, 1, elementType));
801  Value pp[intervalsCount][polyDegree + 1];
802  pp[0][0] = bcast(floatCst(builder, +0.00000000000000000e+00f, elementType));
803  pp[0][1] = bcast(floatCst(builder, +1.12837916222975858e+00f, elementType));
804  pp[0][2] = bcast(floatCst(builder, -5.23018562988006470e-01f, elementType));
805  pp[0][3] = bcast(floatCst(builder, +2.09741709609267072e-01f, elementType));
806  pp[0][4] = bcast(floatCst(builder, +2.58146801602987875e-02f, elementType));
807  pp[1][0] = bcast(floatCst(builder, +0.00000000000000000e+00f, elementType));
808  pp[1][1] = bcast(floatCst(builder, +1.12750687816789140e+00f, elementType));
809  pp[1][2] = bcast(floatCst(builder, -3.64721408487825775e-01f, elementType));
810  pp[1][3] = bcast(floatCst(builder, +1.18407396425136952e-01f, elementType));
811  pp[1][4] = bcast(floatCst(builder, +3.70645533056476558e-02f, elementType));
812  pp[2][0] = bcast(floatCst(builder, -3.30093071049483172e-03f, elementType));
813  pp[2][1] = bcast(floatCst(builder, +3.51961938357697011e-03f, elementType));
814  pp[2][2] = bcast(floatCst(builder, -1.41373622814988039e-03f, elementType));
815  pp[2][3] = bcast(floatCst(builder, +2.53447094961941348e-04f, elementType));
816  pp[2][4] = bcast(floatCst(builder, -1.71048029455037401e-05f, elementType));
817 
818  Value qq[intervalsCount][polyDegree + 1];
819  qq[0][0] = bcast(floatCst(builder, +1.000000000000000000e+00f, elementType));
820  qq[0][1] = bcast(floatCst(builder, -4.635138185962547255e-01f, elementType));
821  qq[0][2] = bcast(floatCst(builder, +5.192301327279782447e-01f, elementType));
822  qq[0][3] = bcast(floatCst(builder, -1.318089722204810087e-01f, elementType));
823  qq[0][4] = bcast(floatCst(builder, +7.397964654672315005e-02f, elementType));
824  qq[1][0] = bcast(floatCst(builder, +1.00000000000000000e+00f, elementType));
825  qq[1][1] = bcast(floatCst(builder, -3.27607011824493086e-01f, elementType));
826  qq[1][2] = bcast(floatCst(builder, +4.48369090658821977e-01f, elementType));
827  qq[1][3] = bcast(floatCst(builder, -8.83462621207857930e-02f, elementType));
828  qq[1][4] = bcast(floatCst(builder, +5.72442770283176093e-02f, elementType));
829  qq[2][0] = bcast(floatCst(builder, +1.00000000000000000e+00f, elementType));
830  qq[2][1] = bcast(floatCst(builder, -2.06069165953913769e+00f, elementType));
831  qq[2][2] = bcast(floatCst(builder, +1.62705939945477759e+00f, elementType));
832  qq[2][3] = bcast(floatCst(builder, -5.83389859211130017e-01f, elementType));
833  qq[2][4] = bcast(floatCst(builder, +8.21908939856640930e-02f, elementType));
834 
835  Value offsets[intervalsCount];
836  offsets[0] = bcast(floatCst(builder, 0.0f, elementType));
837  offsets[1] = bcast(floatCst(builder, 0.0f, elementType));
838  offsets[2] = bcast(floatCst(builder, 1.0f, elementType));
839 
840  Value bounds[intervalsCount];
841  bounds[0] = bcast(floatCst(builder, 0.8f, elementType));
842  bounds[1] = bcast(floatCst(builder, 2.0f, elementType));
843  bounds[2] = bcast(floatCst(builder, 3.75f, elementType));
844 
845  Value isNegativeArg =
846  builder.create<arith::CmpFOp>(arith::CmpFPredicate::OLT, operand, zero);
847  Value negArg = builder.create<arith::NegFOp>(operand);
848  Value x = builder.create<arith::SelectOp>(isNegativeArg, negArg, operand);
849 
850  Value offset = offsets[0];
851  Value p[polyDegree + 1];
852  Value q[polyDegree + 1];
853  for (int i = 0; i <= polyDegree; ++i) {
854  p[i] = pp[0][i];
855  q[i] = qq[0][i];
856  }
857 
858  // TODO: maybe use vector stacking to reduce the number of selects.
859  Value isLessThanBound[intervalsCount];
860  for (int j = 0; j < intervalsCount - 1; ++j) {
861  isLessThanBound[j] =
862  builder.create<arith::CmpFOp>(arith::CmpFPredicate::OLT, x, bounds[j]);
863  for (int i = 0; i <= polyDegree; ++i) {
864  p[i] = builder.create<arith::SelectOp>(isLessThanBound[j], p[i],
865  pp[j + 1][i]);
866  q[i] = builder.create<arith::SelectOp>(isLessThanBound[j], q[i],
867  qq[j + 1][i]);
868  }
869  offset = builder.create<arith::SelectOp>(isLessThanBound[j], offset,
870  offsets[j + 1]);
871  }
872  isLessThanBound[intervalsCount - 1] = builder.create<arith::CmpFOp>(
873  arith::CmpFPredicate::ULT, x, bounds[intervalsCount - 1]);
874 
875  Value pPoly = makePolynomialCalculation(builder, p, x);
876  Value qPoly = makePolynomialCalculation(builder, q, x);
877  Value rationalPoly = builder.create<arith::DivFOp>(pPoly, qPoly);
878  Value formula = builder.create<arith::AddFOp>(offset, rationalPoly);
879  formula = builder.create<arith::SelectOp>(isLessThanBound[intervalsCount - 1],
880  formula, one);
881 
882  // erf is odd function: erf(x) = -erf(-x).
883  Value negFormula = builder.create<arith::NegFOp>(formula);
884  Value res =
885  builder.create<arith::SelectOp>(isNegativeArg, negFormula, formula);
886 
887  rewriter.replaceOp(op, res);
888 
889  return success();
890 }
891 
892 //----------------------------------------------------------------------------//
893 // Exp approximation.
894 //----------------------------------------------------------------------------//
895 
896 namespace {
897 
898 struct ExpApproximation : public OpRewritePattern<math::ExpOp> {
899 public:
901 
902  LogicalResult matchAndRewrite(math::ExpOp op,
903  PatternRewriter &rewriter) const final;
904 };
905 } // namespace
906 
907 // Approximate exp(x) using its reduced range exp(y) where y is in the range
908 // [0, ln(2)], let y = x - floor(x / ln(2)) * ln(2) = x - k * ln(2), exp(x)
909 // = exp(y) * 2^k. exp(y).
911 ExpApproximation::matchAndRewrite(math::ExpOp op,
912  PatternRewriter &rewriter) const {
913  if (!getElementTypeOrSelf(op.getOperand()).isF32())
914  return rewriter.notifyMatchFailure(op, "unsupported operand type");
915 
916  ArrayRef<int64_t> shape = vectorShape(op.getOperand());
917 
918  ImplicitLocOpBuilder builder(op->getLoc(), rewriter);
919 
920  // TODO: Consider a common pattern rewriter with all methods below to
921  // write the approximations.
922  auto bcast = [&](Value value) -> Value {
923  return broadcast(builder, value, shape);
924  };
925  auto fmla = [&](Value a, Value b, Value c) {
926  return builder.create<math::FmaOp>(a, b, c);
927  };
928  auto mul = [&](Value a, Value b) -> Value {
929  return builder.create<arith::MulFOp>(a, b);
930  };
931  auto sub = [&](Value a, Value b) -> Value {
932  return builder.create<arith::SubFOp>(a, b);
933  };
934  auto floor = [&](Value a) { return builder.create<math::FloorOp>(a); };
935 
936  Value cstLn2 = bcast(f32Cst(builder, static_cast<float>(LN2_VALUE)));
937  Value cstLog2E = bcast(f32Cst(builder, static_cast<float>(LOG2E_VALUE)));
938 
939  // Polynomial coefficients.
940  Value cstCephesExpP0 = bcast(f32Cst(builder, 1.0));
941  Value cstCephesExpP1 = bcast(f32Cst(builder, 1.0));
942  Value cstCephesExpP2 = bcast(f32Cst(builder, 0.49970514590562437052f));
943  Value cstCephesExpP3 = bcast(f32Cst(builder, 0.16873890085469545053f));
944  Value cstCephesExpP4 = bcast(f32Cst(builder, 0.03668965196652099192f));
945  Value cstCephesExpP5 = bcast(f32Cst(builder, 0.01314350012789660196f));
946 
947  Value x = op.getOperand();
948 
949  Value isNan = builder.create<arith::CmpFOp>(arith::CmpFPredicate::UNO, x, x);
950 
951  // Reduced y = x - floor(x / ln(2)) * ln(2) = x - k * ln(2)
952  Value xL2Inv = mul(x, cstLog2E);
953  Value kF32 = floor(xL2Inv);
954  Value kLn2 = mul(kF32, cstLn2);
955  Value y = sub(x, kLn2);
956 
957  // Use Estrin's evaluation scheme with 3 independent parts:
958  // P(y)^y : (c0 + c1 y) + (c2 + c3 y) y^2 + (c4 + c5 y) y^4
959  Value y2 = mul(y, y);
960  Value y4 = mul(y2, y2);
961 
962  Value q0 = fmla(cstCephesExpP1, y, cstCephesExpP0);
963  Value q1 = fmla(cstCephesExpP3, y, cstCephesExpP2);
964  Value q2 = fmla(cstCephesExpP5, y, cstCephesExpP4);
965  Value expY = fmla(q1, y2, q0);
966  expY = fmla(q2, y4, expY);
967 
968  auto i32Vec = broadcast(builder.getI32Type(), shape);
969 
970  // exp2(k)
971  Value k = builder.create<arith::FPToSIOp>(i32Vec, kF32);
972  Value exp2KValue = exp2I32(builder, k);
973 
974  // exp(x) = exp(y) * exp2(k)
975  expY = mul(expY, exp2KValue);
976 
977  // Handle overflow, inf and underflow of exp(x). exp(x) range is [0, inf], its
978  // partitioned as the following:
979  // exp(x) = 0, x <= -inf
980  // exp(x) = underflow (min_float), x <= -88
981  // exp(x) = inf (min_float), x >= 88
982  // Note: |k| = 127 is the value where the 8-bits exponent saturates.
983  Value zerof32Const = bcast(f32Cst(builder, 0));
984  auto constPosInfinity =
985  bcast(f32Cst(builder, std::numeric_limits<float>::infinity()));
986  auto constNegIfinity =
987  bcast(f32Cst(builder, -std::numeric_limits<float>::infinity()));
988  auto underflow = bcast(f32Cst(builder, std::numeric_limits<float>::min()));
989 
990  Value kMaxConst = bcast(i32Cst(builder, 127));
991  Value kMaxNegConst = bcast(i32Cst(builder, -127));
992  Value rightBound =
993  builder.create<arith::CmpIOp>(arith::CmpIPredicate::sle, k, kMaxConst);
994  Value leftBound =
995  builder.create<arith::CmpIOp>(arith::CmpIPredicate::sge, k, kMaxNegConst);
996 
997  Value isNegInfinityX = builder.create<arith::CmpFOp>(
998  arith::CmpFPredicate::OEQ, x, constNegIfinity);
999  Value isPosInfinityX = builder.create<arith::CmpFOp>(
1000  arith::CmpFPredicate::OEQ, x, constPosInfinity);
1001  Value isPostiveX =
1002  builder.create<arith::CmpFOp>(arith::CmpFPredicate::OGT, x, zerof32Const);
1003  Value isComputable = builder.create<arith::AndIOp>(rightBound, leftBound);
1004 
1005  expY = builder.create<arith::SelectOp>(
1006  isNan, x,
1007  builder.create<arith::SelectOp>(
1008  isNegInfinityX, zerof32Const,
1009  builder.create<arith::SelectOp>(
1010  isPosInfinityX, constPosInfinity,
1011  builder.create<arith::SelectOp>(
1012  isComputable, expY,
1013  builder.create<arith::SelectOp>(isPostiveX, constPosInfinity,
1014  underflow)))));
1015 
1016  rewriter.replaceOp(op, expY);
1017 
1018  return success();
1019 }
1020 
1021 //----------------------------------------------------------------------------//
1022 // ExpM1 approximation.
1023 //----------------------------------------------------------------------------//
1024 
1025 namespace {
1026 
1027 struct ExpM1Approximation : public OpRewritePattern<math::ExpM1Op> {
1028 public:
1030 
1031  LogicalResult matchAndRewrite(math::ExpM1Op op,
1032  PatternRewriter &rewriter) const final;
1033 };
1034 } // namespace
1035 
1037 ExpM1Approximation::matchAndRewrite(math::ExpM1Op op,
1038  PatternRewriter &rewriter) const {
1039  if (!getElementTypeOrSelf(op.getOperand()).isF32())
1040  return rewriter.notifyMatchFailure(op, "unsupported operand type");
1041 
1042  ArrayRef<int64_t> shape = vectorShape(op.getOperand());
1043 
1044  ImplicitLocOpBuilder builder(op->getLoc(), rewriter);
1045  auto bcast = [&](Value value) -> Value {
1046  return broadcast(builder, value, shape);
1047  };
1048 
1049  // expm1(x) = exp(x) - 1 = u - 1.
1050  // We have to handle it carefully when x is near 0, i.e. u ~= 1,
1051  // and when the input is ~= -inf, i.e. u - 1 ~= -1.
1052  Value cstOne = bcast(f32Cst(builder, 1.0f));
1053  Value cstNegOne = bcast(f32Cst(builder, -1.0f));
1054  Value x = op.getOperand();
1055  Value u = builder.create<math::ExpOp>(x);
1056  Value uEqOneOrNaN =
1057  builder.create<arith::CmpFOp>(arith::CmpFPredicate::UEQ, u, cstOne);
1058  Value uMinusOne = builder.create<arith::SubFOp>(u, cstOne);
1059  Value uMinusOneEqNegOne = builder.create<arith::CmpFOp>(
1060  arith::CmpFPredicate::OEQ, uMinusOne, cstNegOne);
1061  // logU = log(u) ~= x
1062  Value logU = builder.create<math::LogOp>(u);
1063 
1064  // Detect exp(x) = +inf; written this way to avoid having to form +inf.
1065  Value isInf =
1066  builder.create<arith::CmpFOp>(arith::CmpFPredicate::OEQ, logU, u);
1067 
1068  // (u - 1) * (x / ~x)
1069  Value expm1 = builder.create<arith::MulFOp>(
1070  uMinusOne, builder.create<arith::DivFOp>(x, logU));
1071  expm1 = builder.create<arith::SelectOp>(isInf, u, expm1);
1072  Value approximation = builder.create<arith::SelectOp>(
1073  uEqOneOrNaN, x,
1074  builder.create<arith::SelectOp>(uMinusOneEqNegOne, cstNegOne, expm1));
1075  rewriter.replaceOp(op, approximation);
1076  return success();
1077 }
1078 
1079 //----------------------------------------------------------------------------//
1080 // Sin and Cos approximation.
1081 //----------------------------------------------------------------------------//
1082 
1083 namespace {
1084 
1085 template <bool isSine, typename OpTy>
1086 struct SinAndCosApproximation : public OpRewritePattern<OpTy> {
1087 public:
1089 
1090  LogicalResult matchAndRewrite(OpTy op, PatternRewriter &rewriter) const final;
1091 };
1092 } // namespace
1093 
1094 #define TWO_OVER_PI \
1095  0.6366197723675813430755350534900574481378385829618257949906693762L
1096 #define PI_OVER_2 \
1097  1.5707963267948966192313216916397514420985846996875529104874722961L
1098 
1099 // Approximates sin(x) or cos(x) by finding the best approximation polynomial in
1100 // the reduced range [0, pi/2] for both sin(x) and cos(x). Then given y in the
1101 // reduced range sin(x) will be computed as sin(y), -sin(y), cos(y) or -cos(y).
1102 template <bool isSine, typename OpTy>
1103 LogicalResult SinAndCosApproximation<isSine, OpTy>::matchAndRewrite(
1104  OpTy op, PatternRewriter &rewriter) const {
1105  static_assert(
1106  llvm::is_one_of<OpTy, math::SinOp, math::CosOp>::value,
1107  "SinAndCosApproximation pattern expects math::SinOp or math::CosOp");
1108 
1109  if (!getElementTypeOrSelf(op.getOperand()).isF32())
1110  return rewriter.notifyMatchFailure(op, "unsupported operand type");
1111 
1112  ArrayRef<int64_t> shape = vectorShape(op.getOperand());
1113 
1114  ImplicitLocOpBuilder builder(op->getLoc(), rewriter);
1115  auto bcast = [&](Value value) -> Value {
1116  return broadcast(builder, value, shape);
1117  };
1118  auto mul = [&](Value a, Value b) -> Value {
1119  return builder.create<arith::MulFOp>(a, b);
1120  };
1121  auto sub = [&](Value a, Value b) -> Value {
1122  return builder.create<arith::SubFOp>(a, b);
1123  };
1124  auto floor = [&](Value a) { return builder.create<math::FloorOp>(a); };
1125 
1126  auto i32Vec = broadcast(builder.getI32Type(), shape);
1127  auto fPToSingedInteger = [&](Value a) -> Value {
1128  return builder.create<arith::FPToSIOp>(i32Vec, a);
1129  };
1130 
1131  auto modulo4 = [&](Value a) -> Value {
1132  return builder.create<arith::AndIOp>(a, bcast(i32Cst(builder, 3)));
1133  };
1134 
1135  auto isEqualTo = [&](Value a, Value b) -> Value {
1136  return builder.create<arith::CmpIOp>(arith::CmpIPredicate::eq, a, b);
1137  };
1138 
1139  auto isGreaterThan = [&](Value a, Value b) -> Value {
1140  return builder.create<arith::CmpIOp>(arith::CmpIPredicate::sgt, a, b);
1141  };
1142 
1143  auto select = [&](Value cond, Value t, Value f) -> Value {
1144  return builder.create<arith::SelectOp>(cond, t, f);
1145  };
1146 
1147  auto fmla = [&](Value a, Value b, Value c) {
1148  return builder.create<math::FmaOp>(a, b, c);
1149  };
1150 
1151  auto bitwiseOr = [&](Value a, Value b) {
1152  return builder.create<arith::OrIOp>(a, b);
1153  };
1154 
1155  Value twoOverPi = bcast(f32Cst(builder, (float)TWO_OVER_PI));
1156  Value piOverTwo = bcast(f32Cst(builder, (float)PI_OVER_2));
1157 
1158  Value x = op.getOperand();
1159 
1160  Value k = floor(mul(x, twoOverPi));
1161 
1162  Value y = sub(x, mul(k, piOverTwo));
1163 
1164  Value cstOne = bcast(f32Cst(builder, 1.0));
1165  Value cstNegativeOne = bcast(f32Cst(builder, -1.0));
1166 
1167  Value cstSC2 = bcast(f32Cst(builder, -0.16666667163372039794921875f));
1168  Value cstSC4 = bcast(f32Cst(builder, 8.333347737789154052734375e-3f));
1169  Value cstSC6 = bcast(f32Cst(builder, -1.9842604524455964565277099609375e-4f));
1170  Value cstSC8 =
1171  bcast(f32Cst(builder, 2.760012648650445044040679931640625e-6f));
1172  Value cstSC10 =
1173  bcast(f32Cst(builder, -2.50293279435709337121807038784027099609375e-8f));
1174 
1175  Value cstCC2 = bcast(f32Cst(builder, -0.5f));
1176  Value cstCC4 = bcast(f32Cst(builder, 4.166664183139801025390625e-2f));
1177  Value cstCC6 = bcast(f32Cst(builder, -1.388833043165504932403564453125e-3f));
1178  Value cstCC8 = bcast(f32Cst(builder, 2.47562347794882953166961669921875e-5f));
1179  Value cstCC10 =
1180  bcast(f32Cst(builder, -2.59630184018533327616751194000244140625e-7f));
1181 
1182  Value kMod4 = modulo4(fPToSingedInteger(k));
1183 
1184  Value kR0 = isEqualTo(kMod4, bcast(i32Cst(builder, 0)));
1185  Value kR1 = isEqualTo(kMod4, bcast(i32Cst(builder, 1)));
1186  Value kR2 = isEqualTo(kMod4, bcast(i32Cst(builder, 2)));
1187  Value kR3 = isEqualTo(kMod4, bcast(i32Cst(builder, 3)));
1188 
1189  Value sinuseCos = isSine ? bitwiseOr(kR1, kR3) : bitwiseOr(kR0, kR2);
1190  Value negativeRange = isSine ? isGreaterThan(kMod4, bcast(i32Cst(builder, 1)))
1191  : bitwiseOr(kR1, kR2);
1192 
1193  Value y2 = mul(y, y);
1194 
1195  Value base = select(sinuseCos, cstOne, y);
1196  Value cstC2 = select(sinuseCos, cstCC2, cstSC2);
1197  Value cstC4 = select(sinuseCos, cstCC4, cstSC4);
1198  Value cstC6 = select(sinuseCos, cstCC6, cstSC6);
1199  Value cstC8 = select(sinuseCos, cstCC8, cstSC8);
1200  Value cstC10 = select(sinuseCos, cstCC10, cstSC10);
1201 
1202  Value v1 = fmla(y2, cstC10, cstC8);
1203  Value v2 = fmla(y2, v1, cstC6);
1204  Value v3 = fmla(y2, v2, cstC4);
1205  Value v4 = fmla(y2, v3, cstC2);
1206  Value v5 = fmla(y2, v4, cstOne);
1207  Value v6 = mul(base, v5);
1208 
1209  Value approximation = select(negativeRange, mul(cstNegativeOne, v6), v6);
1210 
1211  rewriter.replaceOp(op, approximation);
1212 
1213  return success();
1214 }
1215 
1216 //----------------------------------------------------------------------------//
1217 // Cbrt approximation.
1218 //----------------------------------------------------------------------------//
1219 
1220 namespace {
1221 struct CbrtApproximation : public OpRewritePattern<math::CbrtOp> {
1223 
1224  LogicalResult matchAndRewrite(math::CbrtOp op,
1225  PatternRewriter &rewriter) const final;
1226 };
1227 } // namespace
1228 
1229 // Estimation of cube-root using an algorithm defined in
1230 // Hacker's Delight 2nd Edition.
1232 CbrtApproximation::matchAndRewrite(math::CbrtOp op,
1233  PatternRewriter &rewriter) const {
1234  auto operand = op.getOperand();
1235  if (!getElementTypeOrSelf(operand).isF32())
1236  return rewriter.notifyMatchFailure(op, "unsupported operand type");
1237 
1238  ImplicitLocOpBuilder b(op->getLoc(), rewriter);
1239  ArrayRef<int64_t> shape = vectorShape(operand);
1240 
1241  Type floatTy = getElementTypeOrSelf(operand.getType());
1242  Type intTy = b.getIntegerType(floatTy.getIntOrFloatBitWidth());
1243 
1244  // Convert to vector types if necessary.
1245  floatTy = broadcast(floatTy, shape);
1246  intTy = broadcast(intTy, shape);
1247 
1248  auto bconst = [&](Attribute attr) -> Value {
1249  Value value = b.create<arith::ConstantOp>(attr);
1250  return broadcast(b, value, shape);
1251  };
1252 
1253  // Declare the initial values:
1254  Value intTwo = bconst(b.getI32IntegerAttr(2));
1255  Value intFour = bconst(b.getI32IntegerAttr(4));
1256  Value intEight = bconst(b.getI32IntegerAttr(8));
1257  Value intMagic = bconst(b.getI32IntegerAttr(0x2a5137a0));
1258  Value fpThird = bconst(b.getF32FloatAttr(0.33333333f));
1259  Value fpTwo = bconst(b.getF32FloatAttr(2.0f));
1260  Value fpZero = bconst(b.getF32FloatAttr(0.0f));
1261 
1262  // Compute an approximation of one third:
1263  // union {int ix; float x;};
1264  // x = x0;
1265  // ix = ix/4 + ix/16;
1266  Value absValue = b.create<math::AbsFOp>(operand);
1267  Value intValue = b.create<arith::BitcastOp>(intTy, absValue);
1268  Value divideBy4 = b.create<arith::ShRSIOp>(intValue, intTwo);
1269  Value divideBy16 = b.create<arith::ShRSIOp>(intValue, intFour);
1270  intValue = b.create<arith::AddIOp>(divideBy4, divideBy16);
1271 
1272  // ix = ix + ix/16;
1273  divideBy16 = b.create<arith::ShRSIOp>(intValue, intFour);
1274  intValue = b.create<arith::AddIOp>(intValue, divideBy16);
1275 
1276  // ix = ix + ix/256;
1277  Value divideBy256 = b.create<arith::ShRSIOp>(intValue, intEight);
1278  intValue = b.create<arith::AddIOp>(intValue, divideBy256);
1279 
1280  // ix = 0x2a5137a0 + ix;
1281  intValue = b.create<arith::AddIOp>(intValue, intMagic);
1282 
1283  // Perform one newtons step:
1284  // x = 0.33333333f*(2.0f*x + x0/(x*x));
1285  Value floatValue = b.create<arith::BitcastOp>(floatTy, intValue);
1286  Value squared = b.create<arith::MulFOp>(floatValue, floatValue);
1287  Value mulTwo = b.create<arith::MulFOp>(floatValue, fpTwo);
1288  Value divSquared = b.create<arith::DivFOp>(absValue, squared);
1289  floatValue = b.create<arith::AddFOp>(mulTwo, divSquared);
1290  floatValue = b.create<arith::MulFOp>(floatValue, fpThird);
1291 
1292  // x = 0.33333333f*(2.0f*x + x0/(x*x));
1293  squared = b.create<arith::MulFOp>(floatValue, floatValue);
1294  mulTwo = b.create<arith::MulFOp>(floatValue, fpTwo);
1295  divSquared = b.create<arith::DivFOp>(absValue, squared);
1296  floatValue = b.create<arith::AddFOp>(mulTwo, divSquared);
1297  floatValue = b.create<arith::MulFOp>(floatValue, fpThird);
1298 
1299  // Check for zero and restore sign.
1300  Value isZero =
1301  b.create<arith::CmpFOp>(arith::CmpFPredicate::OEQ, absValue, fpZero);
1302  floatValue = b.create<arith::SelectOp>(isZero, fpZero, floatValue);
1303  floatValue = b.create<math::CopySignOp>(floatValue, operand);
1304 
1305  rewriter.replaceOp(op, floatValue);
1306  return success();
1307 }
1308 
1309 //----------------------------------------------------------------------------//
1310 // Rsqrt approximation.
1311 //----------------------------------------------------------------------------//
1312 
1313 namespace {
1314 struct RsqrtApproximation : public OpRewritePattern<math::RsqrtOp> {
1316 
1317  LogicalResult matchAndRewrite(math::RsqrtOp op,
1318  PatternRewriter &rewriter) const final;
1319 };
1320 } // namespace
1321 
1323 RsqrtApproximation::matchAndRewrite(math::RsqrtOp op,
1324  PatternRewriter &rewriter) const {
1325  if (!getElementTypeOrSelf(op.getOperand()).isF32())
1326  return rewriter.notifyMatchFailure(op, "unsupported operand type");
1327 
1328  ArrayRef<int64_t> shape = vectorShape(op.getOperand());
1329 
1330  // Only support already-vectorized rsqrt's.
1331  if (shape.empty() || shape.back() % 8 != 0)
1332  return rewriter.notifyMatchFailure(op, "unsupported operand type");
1333 
1334  ImplicitLocOpBuilder builder(op->getLoc(), rewriter);
1335  auto bcast = [&](Value value) -> Value {
1336  return broadcast(builder, value, shape);
1337  };
1338 
1339  Value cstPosInf = bcast(f32FromBits(builder, 0x7f800000u));
1340  Value cstOnePointFive = bcast(f32Cst(builder, 1.5f));
1341  Value cstNegHalf = bcast(f32Cst(builder, -0.5f));
1342  Value cstMinNormPos = bcast(f32FromBits(builder, 0x00800000u));
1343 
1344  Value negHalf = builder.create<arith::MulFOp>(op.getOperand(), cstNegHalf);
1345 
1346  // Select only the inverse sqrt of positive normals (denormals are
1347  // flushed to zero).
1348  Value ltMinMask = builder.create<arith::CmpFOp>(
1349  arith::CmpFPredicate::OLT, op.getOperand(), cstMinNormPos);
1350  Value infMask = builder.create<arith::CmpFOp>(arith::CmpFPredicate::OEQ,
1351  op.getOperand(), cstPosInf);
1352  Value notNormalFiniteMask = builder.create<arith::OrIOp>(ltMinMask, infMask);
1353 
1354  // Compute an approximate result.
1356  builder, op->getOperands(), 8, [&builder](ValueRange operands) -> Value {
1357  return builder.create<x86vector::RsqrtOp>(operands);
1358  });
1359 
1360  // Do a single step of Newton-Raphson iteration to improve the approximation.
1361  // This uses the formula y_{n+1} = y_n * (1.5 - y_n * (0.5 * x) * y_n).
1362  // It is essential to evaluate the inner term like this because forming
1363  // y_n^2 may over- or underflow.
1364  Value inner = builder.create<arith::MulFOp>(negHalf, yApprox);
1365  Value fma = builder.create<math::FmaOp>(yApprox, inner, cstOnePointFive);
1366  Value yNewton = builder.create<arith::MulFOp>(yApprox, fma);
1367 
1368  // Select the result of the Newton-Raphson step for positive normal arguments.
1369  // For other arguments, choose the output of the intrinsic. This will
1370  // return rsqrt(+inf) = 0, rsqrt(x) = NaN if x < 0, and rsqrt(x) = +inf if
1371  // x is zero or a positive denormalized float (equivalent to flushing positive
1372  // denormalized inputs to zero).
1373  Value res =
1374  builder.create<arith::SelectOp>(notNormalFiniteMask, yApprox, yNewton);
1375  rewriter.replaceOp(op, res);
1376 
1377  return success();
1378 }
1379 
1380 //----------------------------------------------------------------------------//
1381 
1383  RewritePatternSet &patterns,
1385  // Patterns for leveraging existing f32 lowerings on other data types.
1386  patterns
1387  .add<ReuseF32Expansion<math::AtanOp>, ReuseF32Expansion<math::Atan2Op>,
1388  ReuseF32Expansion<math::TanhOp>, ReuseF32Expansion<math::LogOp>,
1389  ReuseF32Expansion<math::Log2Op>, ReuseF32Expansion<math::Log1pOp>,
1390  ReuseF32Expansion<math::ErfOp>, ReuseF32Expansion<math::ExpOp>,
1391  ReuseF32Expansion<math::ExpM1Op>, ReuseF32Expansion<math::CbrtOp>,
1392  ReuseF32Expansion<math::SinOp>, ReuseF32Expansion<math::CosOp>>(
1393  patterns.getContext());
1394 
1395  patterns.add<AtanApproximation, Atan2Approximation, TanhApproximation,
1396  LogApproximation, Log2Approximation, Log1pApproximation,
1397  ErfPolynomialApproximation, ExpApproximation, ExpM1Approximation,
1398  CbrtApproximation, SinAndCosApproximation<true, math::SinOp>,
1399  SinAndCosApproximation<false, math::CosOp>>(
1400  patterns.getContext());
1401  if (options.enableAvx2) {
1402  patterns.add<RsqrtApproximation, ReuseF32Expansion<math::RsqrtOp>>(
1403  patterns.getContext());
1404  }
1405 }
static llvm::ManagedStatic< PassManagerOptions > options
static std::pair< Value, Value > frexp(ImplicitLocOpBuilder &builder, Value arg, bool isPositive=false)
#define LN2_VALUE
static Type broadcast(Type type, ArrayRef< int64_t > shape)
static Value f32Cst(ImplicitLocOpBuilder &builder, float value)
static Value exp2I32(ImplicitLocOpBuilder &builder, Value arg)
#define PI_OVER_2
#define TWO_OVER_PI
static ArrayRef< int64_t > vectorShape(Type type)
static Value floatCst(ImplicitLocOpBuilder &builder, float value, Type elementType)
static Value handleMultidimensionalVectors(ImplicitLocOpBuilder &builder, ValueRange operands, int64_t vectorWidth, llvm::function_ref< Value(ValueRange)> compute)
LogicalResult insertCasts(Operation *op, PatternRewriter &rewriter)
static Value clamp(ImplicitLocOpBuilder &builder, Value value, Value lowerBound, Value upperBound)
static Value i32Cst(ImplicitLocOpBuilder &builder, int32_t value)
#define LOG2E_VALUE
static Value max(ImplicitLocOpBuilder &builder, Value value, Value bound)
static Value f32FromBits(ImplicitLocOpBuilder &builder, uint32_t bits)
static Value min(ImplicitLocOpBuilder &builder, Value value, Value bound)
static Type getElementType(Type type, ArrayRef< int32_t > indices, function_ref< InFlightDiagnostic(StringRef)> emitErrorFn)
Walks the given type hierarchy with the given indices, potentially down to component granularity,...
Definition: SPIRVOps.cpp:698
static bool isZero(OpFoldResult v)
Definition: Tiling.cpp:47
Attributes are known-constant values of operations.
Definition: Attributes.h:25
IntegerAttr getI32IntegerAttr(int32_t value)
Definition: Builders.cpp:202
FloatType getF32Type()
Definition: Builders.cpp:60
FloatAttr getFloatAttr(Type type, double value)
Definition: Builders.cpp:247
IntegerType getI32Type()
Definition: Builders.cpp:80
IntegerType getIntegerType(unsigned width)
Definition: Builders.cpp:84
Attribute getZeroAttr(Type type)
Definition: Builders.cpp:318
FloatAttr getF32FloatAttr(float value)
Definition: Builders.cpp:239
ImplicitLocOpBuilder maintains a 'current location', allowing use of the create<> method without spec...
OpTy create(Args &&...args)
Create an operation of specific op type at the current insertion point and location.
This class defines the main interface for locations in MLIR and acts as a non-nullable wrapper around...
Definition: Location.h:63
Operation * create(const OperationState &state)
Creates an operation given the fields represented as an OperationState.
Definition: Builders.cpp:432
Location getLoc()
The source location the operation was defined or derived from.
Definition: OpDefinition.h:109
This provides public APIs that all operations should have.
Operation is the basic unit of execution within MLIR.
Definition: Operation.h:75
Location getLoc()
The source location the operation was defined or derived from.
Definition: Operation.h:207
ArrayRef< NamedAttribute > getAttrs()
Return all of the attributes on this operation.
Definition: Operation.h:418
operand_type_range getOperandTypes()
Definition: Operation.h:376
result_type_range getResultTypes()
Definition: Operation.h:407
operand_range getOperands()
Returns an iterator on the underlying Value's.
Definition: Operation.h:357
A special type of RewriterBase that coordinates the application of a rewrite pattern on the current I...
Definition: PatternMatch.h:668
MLIRContext * getContext() const
RewritePatternSet & add(ConstructorArg &&arg, ConstructorArgs &&...args)
Add an instance of each of the pattern types 'Ts' to the pattern list with the given arguments.
std::enable_if_t<!std::is_convertible< CallbackT, Twine >::value, LogicalResult > notifyMatchFailure(Location loc, CallbackT &&reasonCallback)
Used to notify the rewriter that the IR failed to be rewritten because of a match failure,...
Definition: PatternMatch.h:597
virtual void replaceOp(Operation *op, ValueRange newValues)
This method replaces the results of the operation with the specified list of values.
OpTy replaceOpWithNewOp(Operation *op, Args &&...args)
Replaces the result op with a new op that is created without verification.
Definition: PatternMatch.h:482
This class provides an abstraction over the various different ranges of value types.
Definition: TypeRange.h:36
Instances of the Type class are uniqued, have an immutable identifier and an optional mutable compone...
Definition: Types.h:74
U dyn_cast() const
Definition: Types.h:311
bool isF32() const
Definition: Types.cpp:44
bool isF16() const
Definition: Types.cpp:43
bool isa() const
Definition: Types.h:301
unsigned getIntOrFloatBitWidth() const
Return the bit width of an integer or a float type, assert failure on other types.
Definition: Types.cpp:112
This class provides an abstraction over the different types of ranges over Values.
Definition: ValueRange.h:370
type_range getType() const
Type front()
Return first type in the range.
Definition: TypeRange.h:148
This class represents an instance of an SSA value in the MLIR system, representing a computable value...
Definition: Value.h:93
Type getType() const
Return the type of this value.
Definition: Value.h:122
constexpr void enumerate(std::tuple< Tys... > &tuple, CallbackT &&callback)
Definition: Matchers.h:223
int compare(const Fraction &x, const Fraction &y)
Three-way comparison between two fractions.
Definition: Fraction.h:59
LLVM_ATTRIBUTE_ALWAYS_INLINE MPInt abs(const MPInt &x)
Definition: MPInt.h:370
MPInt floor(const Fraction &f)
Definition: Fraction.h:68
This header declares functions that assit transformations in the MemRef dialect.
LogicalResult failure(bool isFailure=true)
Utility function to generate a LogicalResult.
Definition: LogicalResult.h:62
SmallVector< int64_t > computeStrides(ArrayRef< int64_t > sizes)
Definition: IndexingUtils.h:45
SmallVector< int64_t > delinearize(int64_t linearIndex, ArrayRef< int64_t > strides)
Given the strides together with a linear index in the dimension space, return the vector-space offset...
LogicalResult success(bool isSuccess=true)
Utility function to generate a LogicalResult.
Definition: LogicalResult.h:56
Type getElementTypeOrSelf(Type type)
Return the element type or return the type itself.
int64_t computeMaxLinearIndex(ArrayRef< int64_t > basis)
Return the number of elements of basis (i.e.
void populateMathPolynomialApproximationPatterns(RewritePatternSet &patterns, const MathPolynomialApproximationOptions &options={})
This class represents an efficient way to signal success or failure.
Definition: LogicalResult.h:26
OpRewritePattern is a wrapper around RewritePattern that allows for matching and rewriting against an...
Definition: PatternMatch.h:357
OpRewritePattern(MLIRContext *context, PatternBenefit benefit=1, ArrayRef< StringRef > generatedNames={})
Patterns must specify the root operation name they match against, and can also specify the benefit of...
Definition: PatternMatch.h:361
LogicalResult matchAndRewrite(math::ErfOp op, PatternRewriter &rewriter) const final
Eliminates variable at the specified position using Fourier-Motzkin variable elimination.