MLIR 24.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 <cmath>
16#include <cstddef>
17
24#include "mlir/IR/Builders.h"
30#include "llvm/ADT/ArrayRef.h"
31#include "llvm/ADT/STLExtras.h"
32#include "llvm/Support/MathExtras.h"
33
34using namespace mlir;
35using namespace mlir::math;
36using namespace mlir::vector;
37
38// Helper to encapsulate a vector's shape (including scalable dims).
43
44// Returns vector shape if the type is a vector, otherwise return nullopt.
45static std::optional<VectorShape> vectorShape(Type type) {
46 if (auto vectorType = dyn_cast<VectorType>(type)) {
47 return VectorShape{vectorType.getShape(), vectorType.getScalableDims()};
48 }
49 return std::nullopt;
50}
51
52static std::optional<VectorShape> vectorShape(Value value) {
53 return vectorShape(value.getType());
54}
55
56//----------------------------------------------------------------------------//
57// Broadcast scalar types and values into vector types and values.
58//----------------------------------------------------------------------------//
59
60// Broadcasts scalar type into vector type (iff shape is non-scalar).
61static Type broadcast(Type type, std::optional<VectorShape> shape) {
62 assert(!isa<VectorType>(type) && "must be scalar type");
63 return shape ? VectorType::get(shape->sizes, type, shape->scalableFlags)
64 : type;
65}
66
67// Broadcasts scalar value into vector (iff shape is non-scalar).
69 std::optional<VectorShape> shape) {
70 assert(!isa<VectorType>(value.getType()) && "must be scalar value");
71 auto type = broadcast(value.getType(), shape);
72 return shape ? BroadcastOp::create(builder, type, value) : value;
73}
74
75//----------------------------------------------------------------------------//
76// Helper functions to create constants.
77//----------------------------------------------------------------------------//
78
79static Value boolCst(ImplicitLocOpBuilder &builder, bool value) {
80 return arith::ConstantOp::create(builder, builder.getBoolAttr(value));
81}
82
83static Value floatCst(ImplicitLocOpBuilder &builder, float value,
84 Type elementType) {
85 assert((elementType.isF16() || elementType.isF32()) &&
86 "x must be f16 or f32 type.");
87 return arith::ConstantOp::create(builder,
88 builder.getFloatAttr(elementType, value));
89}
90
91static Value f32Cst(ImplicitLocOpBuilder &builder, double value) {
92 return arith::ConstantOp::create(builder, builder.getF32FloatAttr(value));
93}
94
95static Value i32Cst(ImplicitLocOpBuilder &builder, int32_t value) {
96 return arith::ConstantOp::create(builder, builder.getI32IntegerAttr(value));
97}
98
99static Value f32FromBits(ImplicitLocOpBuilder &builder, uint32_t bits) {
100 Value i32Value = i32Cst(builder, static_cast<int32_t>(bits));
101 return arith::BitcastOp::create(builder, builder.getF32Type(), i32Value);
102}
103
104//----------------------------------------------------------------------------//
105// Helper functions to build math functions approximations.
106//----------------------------------------------------------------------------//
107
108// Return the minimum of the two values or NaN if value is NaN
109static Value min(ImplicitLocOpBuilder &builder, Value value, Value bound) {
110 return arith::SelectOp::create(
111 builder,
112 arith::CmpFOp::create(builder, arith::CmpFPredicate::ULT, value, bound),
113 value, bound);
114}
115
116// Return the maximum of the two values or NaN if value is NaN
117static Value max(ImplicitLocOpBuilder &builder, Value value, Value bound) {
118 return arith::SelectOp::create(
119 builder,
120 arith::CmpFOp::create(builder, arith::CmpFPredicate::UGT, value, bound),
121 value, bound);
122}
123
124// Return the clamped value or NaN if value is NaN
125static Value clamp(ImplicitLocOpBuilder &builder, Value value, Value lowerBound,
126 Value upperBound) {
127 return max(builder, min(builder, value, upperBound), lowerBound);
128}
129
130// Decomposes given floating point value `arg` into a normalized fraction and
131// an integral power of two (see std::frexp). Returned values have float type.
132static std::pair<Value, Value> frexp(ImplicitLocOpBuilder &builder, Value arg,
133 bool isPositive = false) {
134 assert(getElementTypeOrSelf(arg).isF32() && "arg must be f32 type");
135 std::optional<VectorShape> shape = vectorShape(arg);
136
137 auto bcast = [&](Value value) -> Value {
138 return broadcast(builder, value, shape);
139 };
140
141 auto i32 = builder.getIntegerType(32);
142 auto i32Vec = broadcast(i32, shape);
143 auto f32Vec = broadcast(builder.getF32Type(), shape);
144
145 Value cst126f = f32Cst(builder, 126.0f);
146 Value cstHalf = f32Cst(builder, 0.5f);
147 Value cstInvMantMask = f32FromBits(builder, ~0x7f800000u);
148
149 // Bitcast to i32 for bitwise operations.
150 Value i32Half = arith::BitcastOp::create(builder, i32, cstHalf);
151 Value i32InvMantMask = arith::BitcastOp::create(builder, i32, cstInvMantMask);
152 Value i32Arg = arith::BitcastOp::create(builder, i32Vec, arg);
153
154 // Compute normalized fraction.
155 Value tmp0 = arith::AndIOp::create(builder, i32Arg, bcast(i32InvMantMask));
156 Value tmp1 = arith::OrIOp::create(builder, tmp0, bcast(i32Half));
157 Value normalizedFraction = arith::BitcastOp::create(builder, f32Vec, tmp1);
158
159 // Compute exponent.
160 Value arg0 = isPositive ? arg : math::AbsFOp::create(builder, arg);
161 Value biasedExponentBits = arith::ShRUIOp::create(
162 builder, arith::BitcastOp::create(builder, i32Vec, arg0),
163 bcast(i32Cst(builder, 23)));
164 Value biasedExponent =
165 arith::SIToFPOp::create(builder, f32Vec, biasedExponentBits);
166 Value exponent =
167 arith::SubFOp::create(builder, biasedExponent, bcast(cst126f));
168
169 return {normalizedFraction, exponent};
170}
171
172// Computes exp2 for an i32 argument.
174 assert(getElementTypeOrSelf(arg).isInteger(32) && "arg must be i32 type");
175 std::optional<VectorShape> shape = vectorShape(arg);
176
177 auto bcast = [&](Value value) -> Value {
178 return broadcast(builder, value, shape);
179 };
180
181 auto f32Vec = broadcast(builder.getF32Type(), shape);
182 // The exponent of f32 located at 23-bit.
183 auto exponetBitLocation = bcast(i32Cst(builder, 23));
184 // Set the exponent bias to zero.
185 auto bias = bcast(i32Cst(builder, 127));
186
187 Value biasedArg = arith::AddIOp::create(builder, arg, bias);
188 Value exp2ValueInt =
189 arith::ShLIOp::create(builder, biasedArg, exponetBitLocation);
190 Value exp2ValueF32 = arith::BitcastOp::create(builder, f32Vec, exp2ValueInt);
191
192 return exp2ValueF32;
193}
194
195namespace {
196Value makePolynomialCalculation(ImplicitLocOpBuilder &builder,
197 llvm::ArrayRef<Value> coeffs, Value x) {
198 Type elementType = getElementTypeOrSelf(x);
199 assert((elementType.isF32() || elementType.isF16()) &&
200 "x must be f32 or f16 type");
201 std::optional<VectorShape> shape = vectorShape(x);
202
203 if (coeffs.empty())
204 return broadcast(builder, floatCst(builder, 0.0f, elementType), shape);
205
206 if (coeffs.size() == 1)
207 return coeffs[0];
208
209 Value res = math::FmaOp::create(builder, x, coeffs[coeffs.size() - 1],
210 coeffs[coeffs.size() - 2]);
211 for (auto i = ptrdiff_t(coeffs.size()) - 3; i >= 0; --i) {
212 res = math::FmaOp::create(builder, x, res, coeffs[i]);
213 }
214 return res;
215}
216} // namespace
217
218//----------------------------------------------------------------------------//
219// Helper function/pattern to insert casts for reusing F32 bit expansion.
220//----------------------------------------------------------------------------//
221
222template <typename T>
223LogicalResult insertCasts(Operation *op, PatternRewriter &rewriter) {
224 // Conservatively only allow where the operand and result types are exactly 1.
225 Type origType = op->getResultTypes().front();
226 for (Type t : llvm::drop_begin(op->getResultTypes()))
227 if (origType != t)
228 return rewriter.notifyMatchFailure(op, "required all types to match");
229 for (Type t : op->getOperandTypes())
230 if (origType != t)
231 return rewriter.notifyMatchFailure(op, "required all types to match");
232
233 // Skip if already F32 or larger than 32 bits.
234 if (getElementTypeOrSelf(origType).isF32() ||
235 getElementTypeOrSelf(origType).getIntOrFloatBitWidth() > 32)
236 return failure();
237
238 // Create F32 equivalent type.
239 Type newType;
240 if (auto shaped = dyn_cast<ShapedType>(origType)) {
241 newType = shaped.clone(rewriter.getF32Type());
242 } else if (isa<FloatType>(origType)) {
243 newType = rewriter.getF32Type();
244 } else {
245 return rewriter.notifyMatchFailure(op,
246 "unable to find F32 equivalent type");
247 }
248
249 Location loc = op->getLoc();
250 SmallVector<Value> operands;
251 for (auto operand : op->getOperands())
252 operands.push_back(arith::ExtFOp::create(rewriter, loc, TypeRange{newType},
253 ValueRange{operand},
254 arith::ExtFOp::Properties{}));
255 auto result = T::create(rewriter, loc, TypeRange{newType}, operands,
256 cast<T>(op).getProperties(),
257 op->getDiscardableAttrDictionary().getValue());
258 rewriter.replaceOpWithNewOp<arith::TruncFOp>(op, origType, result);
259 return success();
260}
261
262namespace {
263// Pattern to cast to F32 to reuse F32 expansion as fallback for single-result
264// op.
265// TODO: Consider revising to avoid adding multiple casts for a subgraph that is
266// all in lower precision. Currently this is only fallback support and performs
267// simplistic casting.
268template <typename T>
269struct ReuseF32Expansion : public OpRewritePattern<T> {
270public:
271 using OpRewritePattern<T>::OpRewritePattern;
272 LogicalResult matchAndRewrite(T op, PatternRewriter &rewriter) const final {
273 static_assert(
274 T::template hasTrait<mlir::OpTrait::SameOperandsAndResultType>(),
275 "requires same operands and result types");
276 return insertCasts<T>(op, rewriter);
277 }
278};
279} // namespace
280
281//----------------------------------------------------------------------------//
282// AtanOp approximation.
283//----------------------------------------------------------------------------//
284
285namespace {
286struct AtanApproximation : public OpRewritePattern<math::AtanOp> {
287public:
288 using Base::Base;
289
290 LogicalResult matchAndRewrite(math::AtanOp op,
291 PatternRewriter &rewriter) const final;
292};
293} // namespace
294
295LogicalResult
296AtanApproximation::matchAndRewrite(math::AtanOp op,
297 PatternRewriter &rewriter) const {
298 auto operand = op.getOperand();
299 if (!getElementTypeOrSelf(operand).isF32())
300 return rewriter.notifyMatchFailure(op, "unsupported operand type");
301
302 std::optional<VectorShape> shape = vectorShape(op.getOperand());
303
304 ImplicitLocOpBuilder builder(op->getLoc(), rewriter);
305 Value abs = math::AbsFOp::create(builder, operand);
306
307 auto one = broadcast(builder, f32Cst(builder, 1.0), shape);
308
309 // When 0.66 < x <= 2.41 we do (x-1) / (x+1):
310 auto twoThirds = broadcast(builder, f32Cst(builder, 0.66), shape);
311 Value cmp2 =
312 arith::CmpFOp::create(builder, arith::CmpFPredicate::OGT, abs, twoThirds);
313 Value addone = arith::AddFOp::create(builder, abs, one);
314 Value subone = arith::SubFOp::create(builder, abs, one);
315 Value xnum = arith::SelectOp::create(builder, cmp2, subone, abs);
316 Value xden = arith::SelectOp::create(builder, cmp2, addone, one);
317
318 auto bcast = [&](Value value) -> Value {
319 return broadcast(builder, value, shape);
320 };
321
322 // Break into the <= 0.66 or > 2.41 we do x or 1/x:
323 auto tan3pio8 = bcast(f32Cst(builder, 2.41421356237309504880));
324 Value cmp1 =
325 arith::CmpFOp::create(builder, arith::CmpFPredicate::OGT, abs, tan3pio8);
326 xnum = arith::SelectOp::create(builder, cmp1, one, xnum);
327 xden = arith::SelectOp::create(builder, cmp1, abs, xden);
328
329 Value x = arith::DivFOp::create(builder, xnum, xden);
330 Value xx = arith::MulFOp::create(builder, x, x);
331
332 // Perform the Taylor series approximation for atan over the range
333 // [0.0, 0.66].
334 auto p0 = bcast(f32Cst(builder, -8.750608600031904122785e-01));
335 auto p1 = bcast(f32Cst(builder, -1.615753718733365076637e+01));
336 auto p2 = bcast(f32Cst(builder, -7.500855792314704667340e+01));
337 auto p3 = bcast(f32Cst(builder, -1.228866684490136173410e+02));
338 auto p4 = bcast(f32Cst(builder, -6.485021904942025371773e+01));
339 auto q0 = bcast(f32Cst(builder, +2.485846490142306297962e+01));
340 auto q1 = bcast(f32Cst(builder, +1.650270098316988542046e+02));
341 auto q2 = bcast(f32Cst(builder, +4.328810604912902668951e+02));
342 auto q3 = bcast(f32Cst(builder, +4.853903996359136964868e+02));
343 auto q4 = bcast(f32Cst(builder, +1.945506571482613964425e+02));
344
345 // Apply the polynomial approximation for the numerator:
346 Value n = p0;
347 n = math::FmaOp::create(builder, xx, n, p1);
348 n = math::FmaOp::create(builder, xx, n, p2);
349 n = math::FmaOp::create(builder, xx, n, p3);
350 n = math::FmaOp::create(builder, xx, n, p4);
351 n = arith::MulFOp::create(builder, n, xx);
352
353 // Apply the polynomial approximation for the denominator:
354 Value d = q0;
355 d = math::FmaOp::create(builder, xx, d, q1);
356 d = math::FmaOp::create(builder, xx, d, q2);
357 d = math::FmaOp::create(builder, xx, d, q3);
358 d = math::FmaOp::create(builder, xx, d, q4);
359
360 // Compute approximation of theta:
361 Value ans0 = arith::DivFOp::create(builder, n, d);
362 ans0 = math::FmaOp::create(builder, ans0, x, x);
363
364 // Correct for the input mapping's angles:
365 Value mpi4 = bcast(f32Cst(builder, llvm::numbers::pi / 4));
366 Value ans2 = arith::AddFOp::create(builder, mpi4, ans0);
367 Value ans = arith::SelectOp::create(builder, cmp2, ans2, ans0);
368
369 Value mpi2 = bcast(f32Cst(builder, llvm::numbers::pi / 2));
370 Value ans1 = arith::SubFOp::create(builder, mpi2, ans0);
371 ans = arith::SelectOp::create(builder, cmp1, ans1, ans);
372
373 // Correct for signing of the input.
374 rewriter.replaceOpWithNewOp<math::CopySignOp>(op, ans, operand);
375 return success();
376}
377
378//----------------------------------------------------------------------------//
379// AtanOp approximation.
380//----------------------------------------------------------------------------//
381
382namespace {
383struct Atan2Approximation : public OpRewritePattern<math::Atan2Op> {
384public:
385 using Base::Base;
386
387 LogicalResult matchAndRewrite(math::Atan2Op op,
388 PatternRewriter &rewriter) const final;
389};
390} // namespace
391
392LogicalResult
393Atan2Approximation::matchAndRewrite(math::Atan2Op op,
394 PatternRewriter &rewriter) const {
395 auto y = op.getOperand(0);
396 auto x = op.getOperand(1);
397 if (!getElementTypeOrSelf(x).isF32())
398 return rewriter.notifyMatchFailure(op, "unsupported operand type");
399
400 ImplicitLocOpBuilder builder(op->getLoc(), rewriter);
401 std::optional<VectorShape> shape = vectorShape(op.getResult());
402
403 // Compute atan in the valid range.
404 auto div = arith::DivFOp::create(builder, y, x);
405 auto atan = math::AtanOp::create(builder, div);
406
407 // Determine what the atan would be for a 180 degree rotation.
408 auto zero = broadcast(builder, f32Cst(builder, 0.0f), shape);
409 auto pi = broadcast(builder, f32Cst(builder, 3.14159265359f), shape);
410 auto addPi = arith::AddFOp::create(builder, atan, pi);
411 auto subPi = arith::SubFOp::create(builder, atan, pi);
412 auto atanGt =
413 arith::CmpFOp::create(builder, arith::CmpFPredicate::OGT, atan, zero);
414 auto flippedAtan = arith::SelectOp::create(builder, atanGt, subPi, addPi);
415
416 // Determine whether to directly use atan or use the 180 degree flip
417 auto xGt = arith::CmpFOp::create(builder, arith::CmpFPredicate::OGT, x, zero);
418 Value result = arith::SelectOp::create(builder, xGt, atan, flippedAtan);
419
420 // Handle x = 0, y > 0
421 Value xZero =
422 arith::CmpFOp::create(builder, arith::CmpFPredicate::OEQ, x, zero);
423 Value yGt =
424 arith::CmpFOp::create(builder, arith::CmpFPredicate::OGT, y, zero);
425 Value isHalfPi = arith::AndIOp::create(builder, xZero, yGt);
426 auto halfPi = broadcast(builder, f32Cst(builder, 1.57079632679f), shape);
427 result = arith::SelectOp::create(builder, isHalfPi, halfPi, result);
428
429 // Handle x = 0, y < 0
430 Value yLt =
431 arith::CmpFOp::create(builder, arith::CmpFPredicate::OLT, y, zero);
432 Value isNegativeHalfPiPi = arith::AndIOp::create(builder, xZero, yLt);
433 auto negativeHalfPiPi =
434 broadcast(builder, f32Cst(builder, -1.57079632679f), shape);
435 result = arith::SelectOp::create(builder, isNegativeHalfPiPi,
436 negativeHalfPiPi, result);
437
438 // Handle x = 0, y = 0;
439 Value yZero =
440 arith::CmpFOp::create(builder, arith::CmpFPredicate::OEQ, y, zero);
441 Value isNan = arith::AndIOp::create(builder, xZero, yZero);
442 Value cstNan = broadcast(builder, f32FromBits(builder, 0x7fc00000), shape);
443 result = arith::SelectOp::create(builder, isNan, cstNan, result);
444
445 rewriter.replaceOp(op, result);
446 return success();
447}
448
449//----------------------------------------------------------------------------//
450// TanhOp approximation.
451//----------------------------------------------------------------------------//
452
453namespace {
454struct TanhApproximation : public OpRewritePattern<math::TanhOp> {
455public:
456 using Base::Base;
457
458 LogicalResult matchAndRewrite(math::TanhOp op,
459 PatternRewriter &rewriter) const final;
460};
461} // namespace
462
463LogicalResult
464TanhApproximation::matchAndRewrite(math::TanhOp op,
465 PatternRewriter &rewriter) const {
466 if (!getElementTypeOrSelf(op.getOperand()).isF32())
467 return rewriter.notifyMatchFailure(op, "unsupported operand type");
468
469 std::optional<VectorShape> shape = vectorShape(op.getOperand());
470
471 ImplicitLocOpBuilder builder(op->getLoc(), rewriter);
472 auto bcast = [&](Value value) -> Value {
473 return broadcast(builder, value, shape);
474 };
475
476 // Clamp operand into [plusClamp, minusClamp] range.
477 Value minusClamp = bcast(f32Cst(builder, -7.99881172180175781f));
478 Value plusClamp = bcast(f32Cst(builder, 7.99881172180175781f));
479 Value x = clamp(builder, op.getOperand(), minusClamp, plusClamp);
480
481 // Mask for tiny values that are approximated with `operand`.
482 Value tiny = bcast(f32Cst(builder, 0.0004f));
483 Value tinyMask = arith::CmpFOp::create(
484 builder, arith::CmpFPredicate::OLT,
485 math::AbsFOp::create(builder, op.getOperand()), tiny);
486
487 // The monomial coefficients of the numerator polynomial (odd).
488 Value alpha1 = bcast(f32Cst(builder, 4.89352455891786e-03f));
489 Value alpha3 = bcast(f32Cst(builder, 6.37261928875436e-04f));
490 Value alpha5 = bcast(f32Cst(builder, 1.48572235717979e-05f));
491 Value alpha7 = bcast(f32Cst(builder, 5.12229709037114e-08f));
492 Value alpha9 = bcast(f32Cst(builder, -8.60467152213735e-11f));
493 Value alpha11 = bcast(f32Cst(builder, 2.00018790482477e-13f));
494 Value alpha13 = bcast(f32Cst(builder, -2.76076847742355e-16f));
495
496 // The monomial coefficients of the denominator polynomial (even).
497 Value beta0 = bcast(f32Cst(builder, 4.89352518554385e-03f));
498 Value beta2 = bcast(f32Cst(builder, 2.26843463243900e-03f));
499 Value beta4 = bcast(f32Cst(builder, 1.18534705686654e-04f));
500 Value beta6 = bcast(f32Cst(builder, 1.19825839466702e-06f));
501
502 // Since the polynomials are odd/even, we need x^2.
503 Value x2 = arith::MulFOp::create(builder, x, x);
504
505 // Evaluate the numerator polynomial p.
506 Value p = math::FmaOp::create(builder, x2, alpha13, alpha11);
507 p = math::FmaOp::create(builder, x2, p, alpha9);
508 p = math::FmaOp::create(builder, x2, p, alpha7);
509 p = math::FmaOp::create(builder, x2, p, alpha5);
510 p = math::FmaOp::create(builder, x2, p, alpha3);
511 p = math::FmaOp::create(builder, x2, p, alpha1);
512 p = arith::MulFOp::create(builder, x, p);
513
514 // Evaluate the denominator polynomial q.
515 Value q = math::FmaOp::create(builder, x2, beta6, beta4);
516 q = math::FmaOp::create(builder, x2, q, beta2);
517 q = math::FmaOp::create(builder, x2, q, beta0);
518
519 // Divide the numerator by the denominator.
520 Value res = arith::SelectOp::create(builder, tinyMask, x,
521 arith::DivFOp::create(builder, p, q));
522
523 rewriter.replaceOp(op, res);
524
525 return success();
526}
527
528#define LN2_VALUE \
529 0.693147180559945309417232121458176568075500134360255254120680009493393621L
530#define LOG2E_VALUE \
531 1.442695040888963407359924681001892137426645954152985934135449406931109219L
532
533//----------------------------------------------------------------------------//
534// LogOp and Log2Op approximation.
535//----------------------------------------------------------------------------//
536
537namespace {
538template <typename Op>
539struct LogApproximationBase : public OpRewritePattern<Op> {
541
542 /// Base 2 if 'base2' is set; natural logarithm (base e) otherwise.
543 LogicalResult logMatchAndRewrite(Op op, PatternRewriter &rewriter,
544 bool base2) const;
545};
546} // namespace
547
548// This approximation comes from Julien Pommier's SSE math library.
549// Link: http://gruntthepeon.free.fr/ssemath
550template <typename Op>
551LogicalResult
552LogApproximationBase<Op>::logMatchAndRewrite(Op op, PatternRewriter &rewriter,
553 bool base2) const {
554 if (!getElementTypeOrSelf(op.getOperand()).isF32())
555 return rewriter.notifyMatchFailure(op, "unsupported operand type");
556
557 std::optional<VectorShape> shape = vectorShape(op.getOperand());
558
559 ImplicitLocOpBuilder builder(op->getLoc(), rewriter);
560 auto bcast = [&](Value value) -> Value {
561 return broadcast(builder, value, shape);
562 };
563
564 Value cstZero = bcast(f32Cst(builder, 0.0f));
565 Value cstOne = bcast(f32Cst(builder, 1.0f));
566 Value cstNegHalf = bcast(f32Cst(builder, -0.5f));
567
568 // The smallest non denormalized float number.
569 Value cstMinNormPos = bcast(f32FromBits(builder, 0x00800000u));
570 Value cstMinusInf = bcast(f32FromBits(builder, 0xff800000u));
571 Value cstPosInf = bcast(f32FromBits(builder, 0x7f800000u));
572 Value cstNan = bcast(f32FromBits(builder, 0x7fc00000));
573
574 // Polynomial coefficients.
575 Value cstCephesSQRTHF = bcast(f32Cst(builder, 0.707106781186547524f));
576 Value cstCephesLogP0 = bcast(f32Cst(builder, 7.0376836292E-2f));
577 Value cstCephesLogP1 = bcast(f32Cst(builder, -1.1514610310E-1f));
578 Value cstCephesLogP2 = bcast(f32Cst(builder, 1.1676998740E-1f));
579 Value cstCephesLogP3 = bcast(f32Cst(builder, -1.2420140846E-1f));
580 Value cstCephesLogP4 = bcast(f32Cst(builder, +1.4249322787E-1f));
581 Value cstCephesLogP5 = bcast(f32Cst(builder, -1.6668057665E-1f));
582 Value cstCephesLogP6 = bcast(f32Cst(builder, +2.0000714765E-1f));
583 Value cstCephesLogP7 = bcast(f32Cst(builder, -2.4999993993E-1f));
584 Value cstCephesLogP8 = bcast(f32Cst(builder, +3.3333331174E-1f));
585
586 Value x = op.getOperand();
587
588 // Truncate input values to the minimum positive normal.
589 x = max(builder, x, cstMinNormPos);
590
591 // Extract significant in the range [0.5,1) and exponent.
592 std::pair<Value, Value> pair = frexp(builder, x, /*isPositive=*/true);
593 x = pair.first;
594 Value e = pair.second;
595
596 // Shift the inputs from the range [0.5,1) to [sqrt(1/2), sqrt(2)) and shift
597 // by -1.0. The values are then centered around 0, which improves the
598 // stability of the polynomial evaluation:
599 //
600 // if( x < SQRTHF ) {
601 // e -= 1;
602 // x = x + x - 1.0;
603 // } else { x = x - 1.0; }
604 Value mask = arith::CmpFOp::create(builder, arith::CmpFPredicate::OLT, x,
605 cstCephesSQRTHF);
606 Value tmp = arith::SelectOp::create(builder, mask, x, cstZero);
607
608 x = arith::SubFOp::create(builder, x, cstOne);
609 e = arith::SubFOp::create(
610 builder, e, arith::SelectOp::create(builder, mask, cstOne, cstZero));
611 x = arith::AddFOp::create(builder, x, tmp);
612
613 Value x2 = arith::MulFOp::create(builder, x, x);
614 Value x3 = arith::MulFOp::create(builder, x2, x);
615
616 // Evaluate the polynomial approximant of degree 8 in three parts.
617 Value y0, y1, y2;
618 y0 = math::FmaOp::create(builder, cstCephesLogP0, x, cstCephesLogP1);
619 y1 = math::FmaOp::create(builder, cstCephesLogP3, x, cstCephesLogP4);
620 y2 = math::FmaOp::create(builder, cstCephesLogP6, x, cstCephesLogP7);
621 y0 = math::FmaOp::create(builder, y0, x, cstCephesLogP2);
622 y1 = math::FmaOp::create(builder, y1, x, cstCephesLogP5);
623 y2 = math::FmaOp::create(builder, y2, x, cstCephesLogP8);
624 y0 = math::FmaOp::create(builder, y0, x3, y1);
625 y0 = math::FmaOp::create(builder, y0, x3, y2);
626 y0 = arith::MulFOp::create(builder, y0, x3);
627
628 y0 = math::FmaOp::create(builder, cstNegHalf, x2, y0);
629 x = arith::AddFOp::create(builder, x, y0);
630
631 if (base2) {
632 Value cstLog2e = bcast(f32Cst(builder, static_cast<float>(LOG2E_VALUE)));
633 x = math::FmaOp::create(builder, x, cstLog2e, e);
634 } else {
635 Value cstLn2 = bcast(f32Cst(builder, static_cast<float>(LN2_VALUE)));
636 x = math::FmaOp::create(builder, e, cstLn2, x);
637 }
638
639 Value invalidMask = arith::CmpFOp::create(builder, arith::CmpFPredicate::ULT,
640 op.getOperand(), cstZero);
641 Value zeroMask = arith::CmpFOp::create(builder, arith::CmpFPredicate::OEQ,
642 op.getOperand(), cstZero);
643 Value posInfMask = arith::CmpFOp::create(builder, arith::CmpFPredicate::OEQ,
644 op.getOperand(), cstPosInf);
645
646 // Filter out invalid values:
647 // • x == 0 -> -INF
648 // • x < 0 -> NAN
649 // • x == +INF -> +INF
650 Value aproximation = arith::SelectOp::create(
651 builder, zeroMask, cstMinusInf,
652 arith::SelectOp::create(
653 builder, invalidMask, cstNan,
654 arith::SelectOp::create(builder, posInfMask, cstPosInf, x)));
655
656 rewriter.replaceOp(op, aproximation);
657
658 return success();
659}
660
661namespace {
662struct LogApproximation : public LogApproximationBase<math::LogOp> {
663 using LogApproximationBase::LogApproximationBase;
664
665 LogicalResult matchAndRewrite(math::LogOp op,
666 PatternRewriter &rewriter) const final {
667 return logMatchAndRewrite(op, rewriter, /*base2=*/false);
668 }
669};
670} // namespace
671
672namespace {
673struct Log2Approximation : public LogApproximationBase<math::Log2Op> {
674 using LogApproximationBase::LogApproximationBase;
675
676 LogicalResult matchAndRewrite(math::Log2Op op,
677 PatternRewriter &rewriter) const final {
678 return logMatchAndRewrite(op, rewriter, /*base2=*/true);
679 }
680};
681} // namespace
682
683//----------------------------------------------------------------------------//
684// Log1p approximation.
685//----------------------------------------------------------------------------//
686
687namespace {
688struct Log1pApproximation : public OpRewritePattern<math::Log1pOp> {
689public:
690 using Base::Base;
691
692 LogicalResult matchAndRewrite(math::Log1pOp op,
693 PatternRewriter &rewriter) const final;
694};
695} // namespace
696
697// Approximate log(1+x).
698LogicalResult
699Log1pApproximation::matchAndRewrite(math::Log1pOp op,
700 PatternRewriter &rewriter) const {
701 if (!getElementTypeOrSelf(op.getOperand()).isF32())
702 return rewriter.notifyMatchFailure(op, "unsupported operand type");
703
704 std::optional<VectorShape> shape = vectorShape(op.getOperand());
705
706 ImplicitLocOpBuilder builder(op->getLoc(), rewriter);
707 auto bcast = [&](Value value) -> Value {
708 return broadcast(builder, value, shape);
709 };
710
711 // Approximate log(1+x) using the following, due to W. Kahan:
712 // u = x + 1.0;
713 // if (u == 1.0 || u == inf) return x;
714 // return x * log(u) / (u - 1.0);
715 // ^^^^^^^^^^^^^^^^^^^^^^
716 // "logLarge" below.
717 Value cstOne = bcast(f32Cst(builder, 1.0f));
718 Value x = op.getOperand();
719 Value u = arith::AddFOp::create(builder, x, cstOne);
720 Value uSmall =
721 arith::CmpFOp::create(builder, arith::CmpFPredicate::OEQ, u, cstOne);
722 Value logU = math::LogOp::create(builder, u);
723 Value uInf =
724 arith::CmpFOp::create(builder, arith::CmpFPredicate::OEQ, u, logU);
725 Value logLarge = arith::MulFOp::create(
726 builder, x,
727 arith::DivFOp::create(builder, logU,
728 arith::SubFOp::create(builder, u, cstOne)));
729 Value approximation = arith::SelectOp::create(
730 builder, arith::OrIOp::create(builder, uSmall, uInf), x, logLarge);
731 rewriter.replaceOp(op, approximation);
732 return success();
733}
734
735//----------------------------------------------------------------------------//
736// Asin approximation.
737//----------------------------------------------------------------------------//
738
739// Approximates asin(x).
740// This approximation is based on the following stackoverflow post:
741// https://stackoverflow.com/a/42683455
742namespace {
743struct AsinPolynomialApproximation : public OpRewritePattern<math::AsinOp> {
744public:
745 using Base::Base;
746
747 LogicalResult matchAndRewrite(math::AsinOp op,
748 PatternRewriter &rewriter) const final;
749};
750} // namespace
751LogicalResult
752AsinPolynomialApproximation::matchAndRewrite(math::AsinOp op,
753 PatternRewriter &rewriter) const {
754 Value operand = op.getOperand();
755 Type elementType = getElementTypeOrSelf(operand);
756
757 if (!(elementType.isF32() || elementType.isF16()))
758 return rewriter.notifyMatchFailure(op,
759 "only f32 and f16 type is supported.");
760 std::optional<VectorShape> shape = vectorShape(operand);
761
762 ImplicitLocOpBuilder builder(op->getLoc(), rewriter);
763 auto bcast = [&](Value value) -> Value {
764 return broadcast(builder, value, shape);
765 };
766
767 auto fma = [&](Value a, Value b, Value c) -> Value {
768 return math::FmaOp::create(builder, a, b, c);
769 };
770
771 auto mul = [&](Value a, Value b) -> Value {
772 return arith::MulFOp::create(builder, a, b);
773 };
774
775 auto sub = [&](Value a, Value b) -> Value {
776 return arith::SubFOp::create(builder, a, b);
777 };
778
779 auto abs = [&](Value a) -> Value { return math::AbsFOp::create(builder, a); };
780
781 auto sqrt = [&](Value a) -> Value {
782 return math::SqrtOp::create(builder, a);
783 };
784
785 auto scopy = [&](Value a, Value b) -> Value {
786 return math::CopySignOp::create(builder, a, b);
787 };
788
789 auto sel = [&](Value a, Value b, Value c) -> Value {
790 return arith::SelectOp::create(builder, a, b, c);
791 };
792
793 Value abso = abs(operand);
794 Value aa = mul(operand, operand);
795 Value opp = sqrt(sub(bcast(floatCst(builder, 1.0, elementType)), aa));
796
797 Value gt = arith::CmpFOp::create(builder, arith::CmpFPredicate::OGT, aa,
798 bcast(floatCst(builder, 0.5, elementType)));
799
800 Value x = sel(gt, opp, abso);
801
802 // Asin(x) approximation for x = [-9/16, 9/16]:
803 Value s = mul(x, x);
804 Value q = mul(s, s);
805 Value r = bcast(floatCst(builder, 5.5579749017470502e-2, elementType));
806 Value t = bcast(floatCst(builder, -6.2027913464120114e-2, elementType));
807
808 r = fma(r, q, bcast(floatCst(builder, 5.4224464349245036e-2, elementType)));
809 t = fma(t, q, bcast(floatCst(builder, -1.1326992890324464e-2, elementType)));
810 r = fma(r, q, bcast(floatCst(builder, 1.5268872539397656e-2, elementType)));
811 t = fma(t, q, bcast(floatCst(builder, 1.0493798473372081e-2, elementType)));
812 r = fma(r, q, bcast(floatCst(builder, 1.4106045900607047e-2, elementType)));
813 t = fma(t, q, bcast(floatCst(builder, 1.7339776384962050e-2, elementType)));
814 r = fma(r, q, bcast(floatCst(builder, 2.2372961589651054e-2, elementType)));
815 t = fma(t, q, bcast(floatCst(builder, 3.0381912707941005e-2, elementType)));
816 r = fma(r, q, bcast(floatCst(builder, 4.4642857881094775e-2, elementType)));
817 t = fma(t, q, bcast(floatCst(builder, 7.4999999991367292e-2, elementType)));
818 r = fma(r, s, t);
819 r = fma(r, s, bcast(floatCst(builder, 1.6666666666670193e-1, elementType)));
820 t = mul(x, s);
821 r = fma(r, t, x);
822
823 Value rsub = sub(bcast(floatCst(builder, 1.57079632679, elementType)), r);
824 r = sel(gt, rsub, r);
825 r = scopy(r, operand);
826
827 rewriter.replaceOp(op, r);
828 return success();
829}
830
831//----------------------------------------------------------------------------//
832// Acos approximation.
833//----------------------------------------------------------------------------//
834
835// Approximates acos(x).
836// This approximation is based on the following stackoverflow post:
837// https://stackoverflow.com/a/42683455
838namespace {
839struct AcosPolynomialApproximation : public OpRewritePattern<math::AcosOp> {
840public:
841 using Base::Base;
842
843 LogicalResult matchAndRewrite(math::AcosOp op,
844 PatternRewriter &rewriter) const final;
845};
846} // namespace
847LogicalResult
848AcosPolynomialApproximation::matchAndRewrite(math::AcosOp op,
849 PatternRewriter &rewriter) const {
850 Value operand = op.getOperand();
851 Type elementType = getElementTypeOrSelf(operand);
852
853 if (!(elementType.isF32() || elementType.isF16()))
854 return rewriter.notifyMatchFailure(op,
855 "only f32 and f16 type is supported.");
856 std::optional<VectorShape> shape = vectorShape(operand);
857
858 ImplicitLocOpBuilder builder(op->getLoc(), rewriter);
859 auto bcast = [&](Value value) -> Value {
860 return broadcast(builder, value, shape);
861 };
862
863 auto fma = [&](Value a, Value b, Value c) -> Value {
864 return math::FmaOp::create(builder, a, b, c);
865 };
866
867 auto mul = [&](Value a, Value b) -> Value {
868 return arith::MulFOp::create(builder, a, b);
869 };
870
871 Value negOperand = arith::NegFOp::create(builder, operand);
872 Value zero = bcast(floatCst(builder, 0.0, elementType));
873 Value half = bcast(floatCst(builder, 0.5, elementType));
874 Value negOne = bcast(floatCst(builder, -1.0, elementType));
875 Value selR =
876 arith::CmpFOp::create(builder, arith::CmpFPredicate::OGT, operand, zero);
877 Value r = arith::SelectOp::create(builder, selR, negOperand, operand);
878 Value chkConst = bcast(floatCst(builder, -0.5625, elementType));
879 Value firstPred =
880 arith::CmpFOp::create(builder, arith::CmpFPredicate::OGT, r, chkConst);
881
882 Value trueVal =
883 fma(bcast(floatCst(builder, 9.3282184640716537e-1, elementType)),
884 bcast(floatCst(builder, 1.6839188885261840e+0, elementType)),
885 math::AsinOp::create(builder, r));
886
887 Value falseVal = math::SqrtOp::create(builder, fma(half, r, half));
888 falseVal = math::AsinOp::create(builder, falseVal);
889 falseVal = mul(bcast(floatCst(builder, 2.0, elementType)), falseVal);
890
891 r = arith::SelectOp::create(builder, firstPred, trueVal, falseVal);
892
893 // Check whether the operand lies in between [-1.0, 0.0).
894 Value greaterThanNegOne = arith::CmpFOp::create(
895 builder, arith::CmpFPredicate::OGE, operand, negOne);
896
897 Value lessThanZero =
898 arith::CmpFOp::create(builder, arith::CmpFPredicate::OLT, operand, zero);
899
900 Value betweenNegOneZero =
901 arith::AndIOp::create(builder, greaterThanNegOne, lessThanZero);
902
903 trueVal = fma(bcast(floatCst(builder, 1.8656436928143307e+0, elementType)),
904 bcast(floatCst(builder, 1.6839188885261840e+0, elementType)),
905 arith::NegFOp::create(builder, r));
906
907 Value finalVal =
908 arith::SelectOp::create(builder, betweenNegOneZero, trueVal, r);
909
910 rewriter.replaceOp(op, finalVal);
911 return success();
912}
913
914//----------------------------------------------------------------------------//
915// Erf approximation.
916//----------------------------------------------------------------------------//
917
918// Approximates erf(x) with
919// a - P(x)/Q(x)
920// where P and Q are polynomials of degree 4.
921// Different coefficients are chosen based on the value of x.
922// The approximation error is ~2.5e-07.
923// Boost's minimax tool that utilizes the Remez method was used to find the
924// coefficients.
925LogicalResult
927 PatternRewriter &rewriter) const {
928 Value operand = op.getOperand();
929 Type elementType = getElementTypeOrSelf(operand);
930
931 if (!(elementType.isF32() || elementType.isF16()))
932 return rewriter.notifyMatchFailure(op,
933 "only f32 and f16 type is supported.");
934 std::optional<VectorShape> shape = vectorShape(operand);
935
936 ImplicitLocOpBuilder builder(op->getLoc(), rewriter);
937 auto bcast = [&](Value value) -> Value {
938 return broadcast(builder, value, shape);
939 };
940
941 const int intervalsCount = 3;
942 const int polyDegree = 4;
943
944 Value zero = bcast(floatCst(builder, 0, elementType));
945 Value one = bcast(floatCst(builder, 1, elementType));
946 Value pp[intervalsCount][polyDegree + 1];
947 pp[0][0] = bcast(floatCst(builder, +0.00000000000000000e+00f, elementType));
948 pp[0][1] = bcast(floatCst(builder, +1.12837916222975858e+00f, elementType));
949 pp[0][2] = bcast(floatCst(builder, -5.23018562988006470e-01f, elementType));
950 pp[0][3] = bcast(floatCst(builder, +2.09741709609267072e-01f, elementType));
951 pp[0][4] = bcast(floatCst(builder, +2.58146801602987875e-02f, elementType));
952 pp[1][0] = bcast(floatCst(builder, +0.00000000000000000e+00f, elementType));
953 pp[1][1] = bcast(floatCst(builder, +1.12750687816789140e+00f, elementType));
954 pp[1][2] = bcast(floatCst(builder, -3.64721408487825775e-01f, elementType));
955 pp[1][3] = bcast(floatCst(builder, +1.18407396425136952e-01f, elementType));
956 pp[1][4] = bcast(floatCst(builder, +3.70645533056476558e-02f, elementType));
957 pp[2][0] = bcast(floatCst(builder, -3.30093071049483172e-03f, elementType));
958 pp[2][1] = bcast(floatCst(builder, +3.51961938357697011e-03f, elementType));
959 pp[2][2] = bcast(floatCst(builder, -1.41373622814988039e-03f, elementType));
960 pp[2][3] = bcast(floatCst(builder, +2.53447094961941348e-04f, elementType));
961 pp[2][4] = bcast(floatCst(builder, -1.71048029455037401e-05f, elementType));
962
963 Value qq[intervalsCount][polyDegree + 1];
964 qq[0][0] = bcast(floatCst(builder, +1.000000000000000000e+00f, elementType));
965 qq[0][1] = bcast(floatCst(builder, -4.635138185962547255e-01f, elementType));
966 qq[0][2] = bcast(floatCst(builder, +5.192301327279782447e-01f, elementType));
967 qq[0][3] = bcast(floatCst(builder, -1.318089722204810087e-01f, elementType));
968 qq[0][4] = bcast(floatCst(builder, +7.397964654672315005e-02f, elementType));
969 qq[1][0] = bcast(floatCst(builder, +1.00000000000000000e+00f, elementType));
970 qq[1][1] = bcast(floatCst(builder, -3.27607011824493086e-01f, elementType));
971 qq[1][2] = bcast(floatCst(builder, +4.48369090658821977e-01f, elementType));
972 qq[1][3] = bcast(floatCst(builder, -8.83462621207857930e-02f, elementType));
973 qq[1][4] = bcast(floatCst(builder, +5.72442770283176093e-02f, elementType));
974 qq[2][0] = bcast(floatCst(builder, +1.00000000000000000e+00f, elementType));
975 qq[2][1] = bcast(floatCst(builder, -2.06069165953913769e+00f, elementType));
976 qq[2][2] = bcast(floatCst(builder, +1.62705939945477759e+00f, elementType));
977 qq[2][3] = bcast(floatCst(builder, -5.83389859211130017e-01f, elementType));
978 qq[2][4] = bcast(floatCst(builder, +8.21908939856640930e-02f, elementType));
979
980 Value offsets[intervalsCount];
981 offsets[0] = bcast(floatCst(builder, 0.0f, elementType));
982 offsets[1] = bcast(floatCst(builder, 0.0f, elementType));
983 offsets[2] = bcast(floatCst(builder, 1.0f, elementType));
984
985 Value bounds[intervalsCount];
986 bounds[0] = bcast(floatCst(builder, 0.8f, elementType));
987 bounds[1] = bcast(floatCst(builder, 2.0f, elementType));
988 bounds[2] = bcast(floatCst(builder, 3.75f, elementType));
989
990 Value isNegativeArg =
991 arith::CmpFOp::create(builder, arith::CmpFPredicate::OLT, operand, zero);
992 Value negArg = arith::NegFOp::create(builder, operand);
993 Value x = arith::SelectOp::create(builder, isNegativeArg, negArg, operand);
994
995 Value offset = offsets[0];
996 Value p[polyDegree + 1];
997 Value q[polyDegree + 1];
998 for (int i = 0; i <= polyDegree; ++i) {
999 p[i] = pp[0][i];
1000 q[i] = qq[0][i];
1001 }
1002
1003 // TODO: maybe use vector stacking to reduce the number of selects.
1004 Value isLessThanBound[intervalsCount];
1005 for (int j = 0; j < intervalsCount - 1; ++j) {
1006 isLessThanBound[j] =
1007 arith::CmpFOp::create(builder, arith::CmpFPredicate::OLT, x, bounds[j]);
1008 for (int i = 0; i <= polyDegree; ++i) {
1009 p[i] = arith::SelectOp::create(builder, isLessThanBound[j], p[i],
1010 pp[j + 1][i]);
1011 q[i] = arith::SelectOp::create(builder, isLessThanBound[j], q[i],
1012 qq[j + 1][i]);
1013 }
1014 offset = arith::SelectOp::create(builder, isLessThanBound[j], offset,
1015 offsets[j + 1]);
1016 }
1017 isLessThanBound[intervalsCount - 1] = arith::CmpFOp::create(
1018 builder, arith::CmpFPredicate::ULT, x, bounds[intervalsCount - 1]);
1019
1020 Value pPoly = makePolynomialCalculation(builder, p, x);
1021 Value qPoly = makePolynomialCalculation(builder, q, x);
1022 Value rationalPoly = arith::DivFOp::create(builder, pPoly, qPoly);
1023 Value formula = arith::AddFOp::create(builder, offset, rationalPoly);
1024 formula = arith::SelectOp::create(
1025 builder, isLessThanBound[intervalsCount - 1], formula, one);
1026
1027 // erf is odd function: erf(x) = -erf(-x).
1028 Value negFormula = arith::NegFOp::create(builder, formula);
1029 Value res =
1030 arith::SelectOp::create(builder, isNegativeArg, negFormula, formula);
1031
1032 rewriter.replaceOp(op, res);
1033
1034 return success();
1035}
1036
1037// Approximates erfc(x) with p((x - 2) / (x + 2)), where p is a 9 degree
1038// polynomial.This approximation is based on the following stackoverflow post:
1039// https://stackoverflow.com/questions/35966695/vectorizable-implementation-of-complementary-error-function-erfcf
1040// The stackoverflow post is in turn based on:
1041// M. M. Shepherd and J. G. Laframboise, "Chebyshev Approximation of
1042// (1+2x)exp(x^2)erfc x in 0 <= x < INF", Mathematics of Computation, Vol. 36,
1043// No. 153, January 1981, pp. 249-253.
1044//
1045// Maximum error: 2.65 ulps
1046LogicalResult
1048 PatternRewriter &rewriter) const {
1049 Value x = op.getOperand();
1050 Type et = getElementTypeOrSelf(x);
1051
1052 if (!et.isF32())
1053 return rewriter.notifyMatchFailure(op, "only f32 type is supported.");
1054 std::optional<VectorShape> shape = vectorShape(x);
1055
1056 ImplicitLocOpBuilder builder(op->getLoc(), rewriter);
1057 auto bcast = [&](Value value) -> Value {
1058 return broadcast(builder, value, shape);
1059 };
1060
1061 Value trueValue = bcast(boolCst(builder, true));
1062 Value zero = bcast(floatCst(builder, 0.0f, et));
1063 Value one = bcast(floatCst(builder, 1.0f, et));
1064 Value onehalf = bcast(floatCst(builder, 0.5f, et));
1065 Value neg4 = bcast(floatCst(builder, -4.0f, et));
1066 Value neg2 = bcast(floatCst(builder, -2.0f, et));
1067 Value pos2 = bcast(floatCst(builder, 2.0f, et));
1068 Value posInf = bcast(floatCst(builder, INFINITY, et));
1069 Value clampVal = bcast(floatCst(builder, 10.0546875f, et));
1070
1071 Value a = math::AbsFOp::create(builder, x);
1072 Value p = arith::AddFOp::create(builder, a, pos2);
1073 Value r = arith::DivFOp::create(builder, one, p);
1074 Value q = math::FmaOp::create(builder, neg4, r, one);
1075 Value t = math::FmaOp::create(builder, arith::AddFOp::create(builder, q, one),
1076 neg2, a);
1077 Value e =
1078 math::FmaOp::create(builder, arith::NegFOp::create(builder, a), q, t);
1079 q = math::FmaOp::create(builder, r, e, q);
1080
1081 p = bcast(floatCst(builder, -0x1.a4a000p-12f, et)); // -4.01139259e-4
1082 Value c1 = bcast(floatCst(builder, -0x1.42a260p-10f, et)); // -1.23075210e-3
1083 p = math::FmaOp::create(builder, p, q, c1);
1084 Value c2 = bcast(floatCst(builder, 0x1.585714p-10f, et)); // 1.31355342e-3
1085 p = math::FmaOp::create(builder, p, q, c2);
1086 Value c3 = bcast(floatCst(builder, 0x1.1adcc4p-07f, et)); // 8.63227434e-3
1087 p = math::FmaOp::create(builder, p, q, c3);
1088 Value c4 = bcast(floatCst(builder, -0x1.081b82p-07f, et)); // -8.05991981e-3
1089 p = math::FmaOp::create(builder, p, q, c4);
1090 Value c5 = bcast(floatCst(builder, -0x1.bc0b6ap-05f, et)); // -5.42046614e-2
1091 p = math::FmaOp::create(builder, p, q, c5);
1092 Value c6 = bcast(floatCst(builder, 0x1.4ffc46p-03f, et)); // 1.64055392e-1
1093 p = math::FmaOp::create(builder, p, q, c6);
1094 Value c7 = bcast(floatCst(builder, -0x1.540840p-03f, et)); // -1.66031361e-1
1095 p = math::FmaOp::create(builder, p, q, c7);
1096 Value c8 = bcast(floatCst(builder, -0x1.7bf616p-04f, et)); // -9.27639827e-2
1097 p = math::FmaOp::create(builder, p, q, c8);
1098 Value c9 = bcast(floatCst(builder, 0x1.1ba03ap-02f, et)); // 2.76978403e-1
1099 p = math::FmaOp::create(builder, p, q, c9);
1100
1101 Value d = math::FmaOp::create(builder, pos2, a, one);
1102 r = arith::DivFOp::create(builder, one, d);
1103 q = math::FmaOp::create(builder, p, r, r);
1104 Value negfa = arith::NegFOp::create(builder, a);
1105 Value fmaqah = math::FmaOp::create(builder, q, negfa, onehalf);
1106 Value psubq = arith::SubFOp::create(builder, p, q);
1107 e = math::FmaOp::create(builder, fmaqah, pos2, psubq);
1108 r = math::FmaOp::create(builder, e, r, q);
1109
1110 Value s = arith::MulFOp::create(builder, a, a);
1111 e = math::ExpOp::create(builder, arith::NegFOp::create(builder, s));
1112
1113 t = math::FmaOp::create(builder, arith::NegFOp::create(builder, a), a, s);
1114 r = math::FmaOp::create(
1115 builder, r, e,
1116 arith::MulFOp::create(builder, arith::MulFOp::create(builder, r, e), t));
1117
1118 Value isNotLessThanInf = arith::XOrIOp::create(
1119 builder,
1120 arith::CmpFOp::create(builder, arith::CmpFPredicate::OLT, a, posInf),
1121 trueValue);
1122 r = arith::SelectOp::create(builder, isNotLessThanInf,
1123 arith::AddFOp::create(builder, x, x), r);
1124 Value isGreaterThanClamp =
1125 arith::CmpFOp::create(builder, arith::CmpFPredicate::OGT, a, clampVal);
1126 r = arith::SelectOp::create(builder, isGreaterThanClamp, zero, r);
1127
1128 Value isNegative =
1129 arith::CmpFOp::create(builder, arith::CmpFPredicate::OLT, x, zero);
1130 r = arith::SelectOp::create(builder, isNegative,
1131 arith::SubFOp::create(builder, pos2, r), r);
1132
1133 rewriter.replaceOp(op, r);
1134 return success();
1135}
1136//----------------------------------------------------------------------------//
1137// Exp approximation.
1138//----------------------------------------------------------------------------//
1139
1140namespace {
1141
1142Value clampWithNormals(ImplicitLocOpBuilder &builder,
1143 const std::optional<VectorShape> shape, Value value,
1144 float lowerBound, float upperBound) {
1145 assert(!std::isnan(lowerBound));
1146 assert(!std::isnan(upperBound));
1147
1148 auto bcast = [&](Value value) -> Value {
1149 return broadcast(builder, value, shape);
1150 };
1151
1152 auto selectCmp = [&builder](auto pred, Value value, Value bound) {
1153 return arith::SelectOp::create(
1154 builder, arith::CmpFOp::create(builder, pred, value, bound), value,
1155 bound);
1156 };
1157
1158 // Note: prefer UGE/ULE vs. UGT/ULT, since they generate vmaxps/vminps vs.
1159 // vcmpleps+vmovaps on x86_64. The latter outcome is also obtained with
1160 // arith::{Max,Min}FOp.
1161 value = selectCmp(arith::CmpFPredicate::UGE, value,
1162 bcast(f32Cst(builder, lowerBound)));
1163 value = selectCmp(arith::CmpFPredicate::ULE, value,
1164 bcast(f32Cst(builder, upperBound)));
1165 return value;
1166}
1167
1168struct ExpApproximation : public OpRewritePattern<math::ExpOp> {
1169public:
1170 using Base::Base;
1171
1172 LogicalResult matchAndRewrite(math::ExpOp op,
1173 PatternRewriter &rewriter) const final;
1174};
1175
1176LogicalResult
1177ExpApproximation::matchAndRewrite(math::ExpOp op,
1178 PatternRewriter &rewriter) const {
1179 auto shape = vectorShape(op.getOperand().getType());
1180 auto elementTy = getElementTypeOrSelf(op.getType());
1181 if (!elementTy.isF32())
1182 return rewriter.notifyMatchFailure(op, "unsupported operand type");
1183
1184 ImplicitLocOpBuilder builder(op->getLoc(), rewriter);
1185
1186 auto add = [&](Value a, Value b) -> Value {
1187 return arith::AddFOp::create(builder, a, b);
1188 };
1189 auto bcast = [&](Value value) -> Value {
1190 return broadcast(builder, value, shape);
1191 };
1192 auto floor = [&](Value a) { return math::FloorOp::create(builder, a); };
1193 auto fmla = [&](Value a, Value b, Value c) {
1194 return math::FmaOp::create(builder, a, b, c);
1195 };
1196 auto mul = [&](Value a, Value b) -> Value {
1197 return arith::MulFOp::create(builder, a, b);
1198 };
1199
1200 // Polynomial approximation from Cephes.
1201 //
1202 // To compute e^x, we re-express it as
1203 //
1204 // e^x = e^(a + b)
1205 // = e^(a + n log(2))
1206 // = e^a * 2^n.
1207 //
1208 // We choose n = round(x / log(2)), restricting the value of `a` to
1209 // (-log(2)/2, log(2)/2). We then use a polynomial to compute e^a. The
1210 // relative error between our approximation and the true value of e^a is less
1211 // than 2^-22.5 for all values of `a` within this range.
1212
1213 // Restrict input to a small range, including some values that evaluate to
1214 // +/- inf. Note that for our lower bound, we choose log(2^-126) instead of
1215 // log(F32_EPSILON). We do so because this routine always flushes denormal
1216 // floating points to 0. Therefore, we only need to worry about exponentiating
1217 // up to the smallest representable non-denormal floating point, which is
1218 // 2^-126.
1219
1220 // Constants.
1221 Value cstHalf = bcast(f32Cst(builder, 0.5f));
1222 Value cstOne = bcast(f32Cst(builder, 1.0f));
1223
1224 // 1/log(2)
1225 Value cstLog2ef = bcast(f32Cst(builder, 1.44269504088896341f));
1226
1227 Value cstExpC1 = bcast(f32Cst(builder, -0.693359375f));
1228 Value cstExpC2 = bcast(f32Cst(builder, 2.12194440e-4f));
1229 Value cstExpP0 = bcast(f32Cst(builder, 1.9875691500E-4f));
1230 Value cstExpP1 = bcast(f32Cst(builder, 1.3981999507E-3f));
1231 Value cstExpP2 = bcast(f32Cst(builder, 8.3334519073E-3f));
1232 Value cstExpP3 = bcast(f32Cst(builder, 4.1665795894E-2f));
1233 Value cstExpP4 = bcast(f32Cst(builder, 1.6666665459E-1f));
1234 Value cstExpP5 = bcast(f32Cst(builder, 5.0000001201E-1f));
1235
1236 // Our computations below aren't particularly sensitive to the exact choices
1237 // here, so we choose values a bit larger/smaller than
1238 //
1239 // log(F32_MAX) = 88.723...
1240 // log(2^-126) = -87.337...
1241 Value x = op.getOperand();
1242 x = clampWithNormals(builder, shape, x, -87.8f, 88.8f);
1243 Value n = floor(fmla(x, cstLog2ef, cstHalf));
1244
1245 // When we eventually do the multiplication in e^a * 2^n, we need to handle
1246 // the case when n > 127, the max fp32 exponent (so 2^n == inf) but e^a < 1
1247 // (so e^a * 2^n != inf). There's a similar problem for n < -126, the
1248 // smallest fp32 exponent.
1249 //
1250 // A straightforward solution would be to detect n out of range and split it
1251 // up, doing
1252 //
1253 // e^a * 2^n = e^a * 2^(n1 + n2)
1254 // = (2^n1 * e^a) * 2^n2.
1255 //
1256 // But it turns out this approach is quite slow, probably because it
1257 // manipulates subnormal values.
1258 //
1259 // The approach we use instead is to clamp n to [-127, 127]. Let n' be the
1260 // value of n clamped to [-127, 127]. In the case where n' = 127, `a` can grow
1261 // up to as large as 88.8 - 127 * log(2) which is about 0.7703. Even though
1262 // this value of `a` is outside our previously specified range, e^a will still
1263 // only have a relative error of approximately 2^-16 at worse. In practice
1264 // this seems to work well enough; it passes our exhaustive tests, breaking
1265 // only one result, and by one ulp (we return exp(88.7228394) = max-float but
1266 // we should return inf).
1267 //
1268 // In the case where n' = -127, the original input value of x is so small that
1269 // e^x, our final answer, is less than 2^-126. Since 2^-126 is the smallest
1270 // normal floating point, and since we flush denormals, we simply return 0. We
1271 // do this in a branchless way by observing that our code for constructing 2^n
1272 // produces 0 if n = -127.
1273 //
1274 // The proof that n' = -127 implies e^x < 2^-126 is as follows:
1275 //
1276 // n' = -127 implies n <= -127
1277 // implies round(x / log(2)) <= -127
1278 // implies x/log(2) < -126.5
1279 // implies x < -126.5 * log(2)
1280 // implies e^x < e^(-126.5 * log(2))
1281 // implies e^x < 2^-126.5 < 2^-126
1282 //
1283 // This proves that n' = -127 implies e^x < 2^-126.
1284 n = clampWithNormals(builder, shape, n, -127.0f, 127.0f);
1285
1286 // Computes x = x - n' * log(2), the value for `a`
1287 x = fmla(cstExpC1, n, x);
1288 x = fmla(cstExpC2, n, x);
1289
1290 // Polynomial to compute z = e^a, accurate for a in (-0.5, 0.5).
1291 Value z = fmla(x, cstExpP0, cstExpP1);
1292 z = fmla(z, x, cstExpP2);
1293 z = fmla(z, x, cstExpP3);
1294 z = fmla(z, x, cstExpP4);
1295 z = fmla(z, x, cstExpP5);
1296 z = fmla(z, mul(x, x), x);
1297 z = add(cstOne, z);
1298
1299 // Convert n' to an i32. This is safe because we clamped it above.
1300 auto i32Vec = broadcast(builder.getI32Type(), shape);
1301 Value nI32 = arith::FPToSIOp::create(builder, i32Vec, n);
1302
1303 // Creates the value 2^n' if -126 <= n' <= 127 and 0 if n' = -127.
1304 Value pow2 = exp2I32(builder, nI32);
1305
1306 // Return z * 2^n' if -126 <= n' <= 127 and 0 if n = -127.
1307 Value ret = mul(z, pow2);
1308
1309 rewriter.replaceOp(op, ret);
1310 return mlir::success();
1311}
1312
1313} // namespace
1314
1315//----------------------------------------------------------------------------//
1316// ExpM1 approximation.
1317//----------------------------------------------------------------------------//
1318
1319namespace {
1320
1321struct ExpM1Approximation : public OpRewritePattern<math::ExpM1Op> {
1322public:
1323 using Base::Base;
1324
1325 LogicalResult matchAndRewrite(math::ExpM1Op op,
1326 PatternRewriter &rewriter) const final;
1327};
1328} // namespace
1329
1330LogicalResult
1331ExpM1Approximation::matchAndRewrite(math::ExpM1Op op,
1332 PatternRewriter &rewriter) const {
1333 if (!getElementTypeOrSelf(op.getOperand()).isF32())
1334 return rewriter.notifyMatchFailure(op, "unsupported operand type");
1335
1336 std::optional<VectorShape> shape = vectorShape(op.getOperand());
1337
1338 ImplicitLocOpBuilder builder(op->getLoc(), rewriter);
1339 auto bcast = [&](Value value) -> Value {
1340 return broadcast(builder, value, shape);
1341 };
1342
1343 // expm1(x) = exp(x) - 1 = u - 1.
1344 // We have to handle it carefully when x is near 0, i.e. u ~= 1,
1345 // and when the input is ~= -inf, i.e. u - 1 ~= -1.
1346 Value cstOne = bcast(f32Cst(builder, 1.0f));
1347 Value cstNegOne = bcast(f32Cst(builder, -1.0f));
1348 Value x = op.getOperand();
1349 Value u = math::ExpOp::create(builder, x);
1350 Value uEqOneOrNaN =
1351 arith::CmpFOp::create(builder, arith::CmpFPredicate::UEQ, u, cstOne);
1352 Value uMinusOne = arith::SubFOp::create(builder, u, cstOne);
1353 Value uMinusOneEqNegOne = arith::CmpFOp::create(
1354 builder, arith::CmpFPredicate::OEQ, uMinusOne, cstNegOne);
1355 // logU = log(u) ~= x
1356 Value logU = math::LogOp::create(builder, u);
1357
1358 // Detect exp(x) = +inf; written this way to avoid having to form +inf.
1359 Value isInf =
1360 arith::CmpFOp::create(builder, arith::CmpFPredicate::OEQ, logU, u);
1361
1362 // (u - 1) * (x / ~x)
1363 Value expm1 = arith::MulFOp::create(builder, uMinusOne,
1364 arith::DivFOp::create(builder, x, logU));
1365 expm1 = arith::SelectOp::create(builder, isInf, u, expm1);
1366 Value approximation = arith::SelectOp::create(
1367 builder, uEqOneOrNaN, x,
1368 arith::SelectOp::create(builder, uMinusOneEqNegOne, cstNegOne, expm1));
1369 rewriter.replaceOp(op, approximation);
1370 return success();
1371}
1372
1373//----------------------------------------------------------------------------//
1374// Sin and Cos approximation.
1375//----------------------------------------------------------------------------//
1376
1377namespace {
1378
1379template <bool isSine, typename OpTy>
1380struct SinAndCosApproximation : public OpRewritePattern<OpTy> {
1381public:
1382 using OpRewritePattern<OpTy>::OpRewritePattern;
1383
1384 LogicalResult matchAndRewrite(OpTy op, PatternRewriter &rewriter) const final;
1385};
1386} // namespace
1387
1388#define TWO_OVER_PI \
1389 0.6366197723675813430755350534900574481378385829618257949906693762L
1390#define PI_OVER_2 \
1391 1.5707963267948966192313216916397514420985846996875529104874722961L
1392
1393// Approximates sin(x) or cos(x) by finding the best approximation polynomial in
1394// the reduced range [0, pi/2] for both sin(x) and cos(x). Then given y in the
1395// reduced range sin(x) will be computed as sin(y), -sin(y), cos(y) or -cos(y).
1396template <bool isSine, typename OpTy>
1397LogicalResult SinAndCosApproximation<isSine, OpTy>::matchAndRewrite(
1398 OpTy op, PatternRewriter &rewriter) const {
1399 static_assert(
1400 llvm::is_one_of<OpTy, math::SinOp, math::CosOp>::value,
1401 "SinAndCosApproximation pattern expects math::SinOp or math::CosOp");
1402
1403 if (!getElementTypeOrSelf(op.getOperand()).isF32())
1404 return rewriter.notifyMatchFailure(op, "unsupported operand type");
1405
1406 std::optional<VectorShape> shape = vectorShape(op.getOperand());
1407
1408 ImplicitLocOpBuilder builder(op->getLoc(), rewriter);
1409 auto bcast = [&](Value value) -> Value {
1410 return broadcast(builder, value, shape);
1411 };
1412 auto mul = [&](Value a, Value b) -> Value {
1413 return arith::MulFOp::create(builder, a, b);
1414 };
1415 auto sub = [&](Value a, Value b) -> Value {
1416 return arith::SubFOp::create(builder, a, b);
1417 };
1418 auto floor = [&](Value a) { return math::FloorOp::create(builder, a); };
1419
1420 auto i32Vec = broadcast(builder.getI32Type(), shape);
1421 auto fPToSingedInteger = [&](Value a) -> Value {
1422 return arith::FPToSIOp::create(builder, i32Vec, a);
1423 };
1424
1425 auto modulo4 = [&](Value a) -> Value {
1426 return arith::AndIOp::create(builder, a, bcast(i32Cst(builder, 3)));
1427 };
1428
1429 auto isEqualTo = [&](Value a, Value b) -> Value {
1430 return arith::CmpIOp::create(builder, arith::CmpIPredicate::eq, a, b);
1431 };
1432
1433 auto isGreaterThan = [&](Value a, Value b) -> Value {
1434 return arith::CmpIOp::create(builder, arith::CmpIPredicate::sgt, a, b);
1435 };
1436
1437 auto select = [&](Value cond, Value t, Value f) -> Value {
1438 return arith::SelectOp::create(builder, cond, t, f);
1439 };
1440
1441 auto fmla = [&](Value a, Value b, Value c) {
1442 return math::FmaOp::create(builder, a, b, c);
1443 };
1444
1445 auto bitwiseOr = [&](Value a, Value b) {
1446 return arith::OrIOp::create(builder, a, b);
1447 };
1448
1449 Value twoOverPi = bcast(f32Cst(builder, (float)TWO_OVER_PI));
1450 Value piOverTwo = bcast(f32Cst(builder, (float)PI_OVER_2));
1451
1452 Value x = op.getOperand();
1453
1454 Value k = floor(mul(x, twoOverPi));
1455
1456 Value y = sub(x, mul(k, piOverTwo));
1457
1458 Value cstOne = bcast(f32Cst(builder, 1.0));
1459 Value cstNegativeOne = bcast(f32Cst(builder, -1.0));
1460
1461 Value cstSC2 = bcast(f32Cst(builder, -0.16666667163372039794921875f));
1462 Value cstSC4 = bcast(f32Cst(builder, 8.333347737789154052734375e-3f));
1463 Value cstSC6 = bcast(f32Cst(builder, -1.9842604524455964565277099609375e-4f));
1464 Value cstSC8 =
1465 bcast(f32Cst(builder, 2.760012648650445044040679931640625e-6f));
1466 Value cstSC10 =
1467 bcast(f32Cst(builder, -2.50293279435709337121807038784027099609375e-8f));
1468
1469 Value cstCC2 = bcast(f32Cst(builder, -0.5f));
1470 Value cstCC4 = bcast(f32Cst(builder, 4.166664183139801025390625e-2f));
1471 Value cstCC6 = bcast(f32Cst(builder, -1.388833043165504932403564453125e-3f));
1472 Value cstCC8 = bcast(f32Cst(builder, 2.47562347794882953166961669921875e-5f));
1473 Value cstCC10 =
1474 bcast(f32Cst(builder, -2.59630184018533327616751194000244140625e-7f));
1475
1476 Value kMod4 = modulo4(fPToSingedInteger(k));
1477
1478 Value kR0 = isEqualTo(kMod4, bcast(i32Cst(builder, 0)));
1479 Value kR1 = isEqualTo(kMod4, bcast(i32Cst(builder, 1)));
1480 Value kR2 = isEqualTo(kMod4, bcast(i32Cst(builder, 2)));
1481 Value kR3 = isEqualTo(kMod4, bcast(i32Cst(builder, 3)));
1482
1483 Value sinuseCos = isSine ? bitwiseOr(kR1, kR3) : bitwiseOr(kR0, kR2);
1484 Value negativeRange = isSine ? isGreaterThan(kMod4, bcast(i32Cst(builder, 1)))
1485 : bitwiseOr(kR1, kR2);
1486
1487 Value y2 = mul(y, y);
1488
1489 Value base = select(sinuseCos, cstOne, y);
1490 Value cstC2 = select(sinuseCos, cstCC2, cstSC2);
1491 Value cstC4 = select(sinuseCos, cstCC4, cstSC4);
1492 Value cstC6 = select(sinuseCos, cstCC6, cstSC6);
1493 Value cstC8 = select(sinuseCos, cstCC8, cstSC8);
1494 Value cstC10 = select(sinuseCos, cstCC10, cstSC10);
1495
1496 Value v1 = fmla(y2, cstC10, cstC8);
1497 Value v2 = fmla(y2, v1, cstC6);
1498 Value v3 = fmla(y2, v2, cstC4);
1499 Value v4 = fmla(y2, v3, cstC2);
1500 Value v5 = fmla(y2, v4, cstOne);
1501 Value v6 = mul(base, v5);
1502
1503 Value approximation = select(negativeRange, mul(cstNegativeOne, v6), v6);
1504
1505 rewriter.replaceOp(op, approximation);
1506
1507 return success();
1508}
1509
1510//----------------------------------------------------------------------------//
1511// Cbrt approximation.
1512//----------------------------------------------------------------------------//
1513
1514namespace {
1515struct CbrtApproximation : public OpRewritePattern<math::CbrtOp> {
1516 using Base::Base;
1517
1518 LogicalResult matchAndRewrite(math::CbrtOp op,
1519 PatternRewriter &rewriter) const final;
1520};
1521} // namespace
1522
1523// Estimation of cube-root using an algorithm defined in
1524// Hacker's Delight 2nd Edition.
1525LogicalResult
1526CbrtApproximation::matchAndRewrite(math::CbrtOp op,
1527 PatternRewriter &rewriter) const {
1528 auto operand = op.getOperand();
1529 if (!getElementTypeOrSelf(operand).isF32())
1530 return rewriter.notifyMatchFailure(op, "unsupported operand type");
1531
1532 ImplicitLocOpBuilder b(op->getLoc(), rewriter);
1533 std::optional<VectorShape> shape = vectorShape(operand);
1534
1535 Type floatTy = getElementTypeOrSelf(operand.getType());
1536 Type intTy = b.getIntegerType(floatTy.getIntOrFloatBitWidth());
1537
1538 // Convert to vector types if necessary.
1539 floatTy = broadcast(floatTy, shape);
1540 intTy = broadcast(intTy, shape);
1541
1542 auto bconst = [&](TypedAttr attr) -> Value {
1543 Value value = arith::ConstantOp::create(b, attr);
1544 return broadcast(b, value, shape);
1545 };
1546
1547 // Declare the initial values:
1548 Value intTwo = bconst(b.getI32IntegerAttr(2));
1549 Value intFour = bconst(b.getI32IntegerAttr(4));
1550 Value intEight = bconst(b.getI32IntegerAttr(8));
1551 Value intMagic = bconst(b.getI32IntegerAttr(0x2a5137a0));
1552 Value fpThird = bconst(b.getF32FloatAttr(0.33333333f));
1553 Value fpTwo = bconst(b.getF32FloatAttr(2.0f));
1554 Value fpZero = bconst(b.getF32FloatAttr(0.0f));
1555
1556 // Compute an approximation of one third:
1557 // union {int ix; float x;};
1558 // x = x0;
1559 // ix = ix/4 + ix/16;
1560 Value absValue = math::AbsFOp::create(b, operand);
1561 Value intValue = arith::BitcastOp::create(b, intTy, absValue);
1562 Value divideBy4 = arith::ShRSIOp::create(b, intValue, intTwo);
1563 Value divideBy16 = arith::ShRSIOp::create(b, intValue, intFour);
1564 intValue = arith::AddIOp::create(b, divideBy4, divideBy16);
1565
1566 // ix = ix + ix/16;
1567 divideBy16 = arith::ShRSIOp::create(b, intValue, intFour);
1568 intValue = arith::AddIOp::create(b, intValue, divideBy16);
1569
1570 // ix = ix + ix/256;
1571 Value divideBy256 = arith::ShRSIOp::create(b, intValue, intEight);
1572 intValue = arith::AddIOp::create(b, intValue, divideBy256);
1573
1574 // ix = 0x2a5137a0 + ix;
1575 intValue = arith::AddIOp::create(b, intValue, intMagic);
1576
1577 // Perform one newtons step:
1578 // x = 0.33333333f*(2.0f*x + x0/(x*x));
1579 Value floatValue = arith::BitcastOp::create(b, floatTy, intValue);
1580 Value squared = arith::MulFOp::create(b, floatValue, floatValue);
1581 Value mulTwo = arith::MulFOp::create(b, floatValue, fpTwo);
1582 Value divSquared = arith::DivFOp::create(b, absValue, squared);
1583 floatValue = arith::AddFOp::create(b, mulTwo, divSquared);
1584 floatValue = arith::MulFOp::create(b, floatValue, fpThird);
1585
1586 // x = 0.33333333f*(2.0f*x + x0/(x*x));
1587 squared = arith::MulFOp::create(b, floatValue, floatValue);
1588 mulTwo = arith::MulFOp::create(b, floatValue, fpTwo);
1589 divSquared = arith::DivFOp::create(b, absValue, squared);
1590 floatValue = arith::AddFOp::create(b, mulTwo, divSquared);
1591 floatValue = arith::MulFOp::create(b, floatValue, fpThird);
1592
1593 // Check for zero and restore sign.
1594 Value isZero =
1595 arith::CmpFOp::create(b, arith::CmpFPredicate::OEQ, absValue, fpZero);
1596 floatValue = arith::SelectOp::create(b, isZero, fpZero, floatValue);
1597 floatValue = math::CopySignOp::create(b, floatValue, operand);
1598
1599 rewriter.replaceOp(op, floatValue);
1600 return success();
1601}
1602
1603//----------------------------------------------------------------------------//
1604
1606 RewritePatternSet &patterns) {
1607 patterns.add<TanhApproximation>(patterns.getContext());
1608}
1609
1614
1619
1620template <typename OpType>
1621static void
1623 llvm::function_ref<bool(StringRef)> predicate,
1624 PatternBenefit benefit) {
1625 if (predicate(OpType::getOperationName())) {
1626 patterns.add<ReuseF32Expansion<OpType>>(patterns.getContext(), benefit);
1627 }
1628}
1629
1631 RewritePatternSet &patterns, llvm::function_ref<bool(StringRef)> predicate,
1632 PatternBenefit benefit) {
1633 populateMathF32ExpansionPattern<math::AcosOp>(patterns, predicate, benefit);
1634 populateMathF32ExpansionPattern<math::AcoshOp>(patterns, predicate, benefit);
1635 populateMathF32ExpansionPattern<math::AsinOp>(patterns, predicate, benefit);
1636 populateMathF32ExpansionPattern<math::AsinhOp>(patterns, predicate, benefit);
1637 populateMathF32ExpansionPattern<math::AtanOp>(patterns, predicate, benefit);
1638 populateMathF32ExpansionPattern<math::Atan2Op>(patterns, predicate, benefit);
1639 populateMathF32ExpansionPattern<math::AtanhOp>(patterns, predicate, benefit);
1640 populateMathF32ExpansionPattern<math::CbrtOp>(patterns, predicate, benefit);
1641 populateMathF32ExpansionPattern<math::CosOp>(patterns, predicate, benefit);
1642 populateMathF32ExpansionPattern<math::CoshOp>(patterns, predicate, benefit);
1643 populateMathF32ExpansionPattern<math::ErfOp>(patterns, predicate, benefit);
1644 populateMathF32ExpansionPattern<math::ErfcOp>(patterns, predicate, benefit);
1645 populateMathF32ExpansionPattern<math::ExpOp>(patterns, predicate, benefit);
1646 populateMathF32ExpansionPattern<math::Exp2Op>(patterns, predicate, benefit);
1647 populateMathF32ExpansionPattern<math::ExpM1Op>(patterns, predicate, benefit);
1648 populateMathF32ExpansionPattern<math::LogOp>(patterns, predicate, benefit);
1649 populateMathF32ExpansionPattern<math::Log10Op>(patterns, predicate, benefit);
1650 populateMathF32ExpansionPattern<math::Log1pOp>(patterns, predicate, benefit);
1651 populateMathF32ExpansionPattern<math::Log2Op>(patterns, predicate, benefit);
1652 populateMathF32ExpansionPattern<math::PowFOp>(patterns, predicate, benefit);
1653 populateMathF32ExpansionPattern<math::RsqrtOp>(patterns, predicate, benefit);
1654 populateMathF32ExpansionPattern<math::SinOp>(patterns, predicate, benefit);
1655 populateMathF32ExpansionPattern<math::SinhOp>(patterns, predicate, benefit);
1656 populateMathF32ExpansionPattern<math::SqrtOp>(patterns, predicate, benefit);
1657 populateMathF32ExpansionPattern<math::TanOp>(patterns, predicate, benefit);
1658 populateMathF32ExpansionPattern<math::TanhOp>(patterns, predicate, benefit);
1659}
1660
1661template <typename OpType, typename PatternType>
1663 RewritePatternSet &patterns, llvm::function_ref<bool(StringRef)> predicate,
1664 PatternBenefit benefit) {
1665 if (predicate(OpType::getOperationName())) {
1666 patterns.add<PatternType>(patterns.getContext(), benefit);
1667 }
1668}
1669
1671 RewritePatternSet &patterns, llvm::function_ref<bool(StringRef)> predicate,
1672 PatternBenefit benefit) {
1674 AcosPolynomialApproximation>(
1675 patterns, predicate, benefit);
1677 AsinPolynomialApproximation>(
1678 patterns, predicate, benefit);
1680 patterns, predicate, benefit);
1682 patterns, predicate, benefit);
1684 patterns, predicate, benefit);
1686 CosOp, SinAndCosApproximation<false, math::CosOp>>(patterns, predicate,
1687 benefit);
1689 patterns, predicate, benefit);
1692 patterns, predicate, benefit);
1694 patterns, predicate, benefit);
1696 patterns, predicate, benefit);
1698 patterns, predicate, benefit);
1700 patterns, predicate, benefit);
1702 patterns, predicate, benefit);
1704 SinOp, SinAndCosApproximation<true, math::SinOp>>(patterns, predicate,
1705 benefit);
1707 patterns, predicate, benefit);
1708}
1709
1711 RewritePatternSet &patterns) {
1712 mlir::populateMathF32ExpansionPatterns(patterns, [](StringRef name) -> bool {
1713 return llvm::is_contained(
1714 {math::AtanOp::getOperationName(), math::Atan2Op::getOperationName(),
1715 math::TanhOp::getOperationName(), math::LogOp::getOperationName(),
1716 math::Log2Op::getOperationName(), math::Log1pOp::getOperationName(),
1717 math::ErfOp::getOperationName(), math::ErfcOp::getOperationName(),
1718 math::ExpOp::getOperationName(), math::ExpM1Op::getOperationName(),
1719 math::CbrtOp::getOperationName(), math::SinOp::getOperationName(),
1720 math::CosOp::getOperationName()},
1721 name);
1722 });
1723
1725 patterns, [](StringRef name) -> bool {
1726 return llvm::is_contained(
1727 {math::AtanOp::getOperationName(),
1728 math::Atan2Op::getOperationName(),
1729 math::TanhOp::getOperationName(), math::LogOp::getOperationName(),
1730 math::Log2Op::getOperationName(),
1731 math::Log1pOp::getOperationName(), math::ErfOp::getOperationName(),
1732 math::ErfcOp::getOperationName(), math::AsinOp::getOperationName(),
1733 math::AcosOp::getOperationName(), math::ExpOp::getOperationName(),
1734 math::ExpM1Op::getOperationName(),
1735 math::CbrtOp::getOperationName(), math::SinOp::getOperationName(),
1736 math::CosOp::getOperationName()},
1737 name);
1738 });
1739}
return success()
b
Return true if permutation is a valid permutation of the outer_dims_perm (case OuterOrInnerPerm::Oute...
#define LN2_VALUE
static Value exp2I32(ImplicitLocOpBuilder &builder, Value arg)
#define PI_OVER_2
static void populateMathF32ExpansionPattern(RewritePatternSet &patterns, llvm::function_ref< bool(StringRef)> predicate, PatternBenefit benefit)
#define TWO_OVER_PI
static Value boolCst(ImplicitLocOpBuilder &builder, bool value)
static Value floatCst(ImplicitLocOpBuilder &builder, float value, Type elementType)
LogicalResult insertCasts(Operation *op, PatternRewriter &rewriter)
static Value clamp(ImplicitLocOpBuilder &builder, Value value, Value lowerBound, Value upperBound)
static std::pair< Value, Value > frexp(ImplicitLocOpBuilder &builder, Value arg, bool isPositive=false)
static std::optional< VectorShape > vectorShape(Type type)
static Value i32Cst(ImplicitLocOpBuilder &builder, int32_t value)
static Type broadcast(Type type, std::optional< VectorShape > shape)
#define LOG2E_VALUE
static Value max(ImplicitLocOpBuilder &builder, Value value, Value bound)
static Value f32FromBits(ImplicitLocOpBuilder &builder, uint32_t bits)
static Value f32Cst(ImplicitLocOpBuilder &builder, double value)
static void populateMathPolynomialApproximationPattern(RewritePatternSet &patterns, llvm::function_ref< bool(StringRef)> predicate, PatternBenefit benefit)
static Value min(ImplicitLocOpBuilder &builder, Value value, Value bound)
static Value broadcast(Location loc, Value toBroadcast, unsigned numElements, const TypeConverter &typeConverter, ConversionPatternRewriter &rewriter)
Broadcasts the value to vector with numElements number of elements.
#define mul(a, b)
#define add(a, b)
#define div(a, b)
IntegerAttr getI32IntegerAttr(int32_t value)
Definition Builders.cpp:208
FloatType getF32Type()
Definition Builders.cpp:51
FloatAttr getFloatAttr(Type type, double value)
Definition Builders.cpp:263
IntegerType getI32Type()
Definition Builders.cpp:71
IntegerType getIntegerType(unsigned width)
Definition Builders.cpp:75
BoolAttr getBoolAttr(bool value)
Definition Builders.cpp:108
FloatAttr getF32FloatAttr(float value)
Definition Builders.cpp:255
ImplicitLocOpBuilder maintains a 'current location', allowing use of the create<> method without spec...
Definition Builders.h:632
This class defines the main interface for locations in MLIR and acts as a non-nullable wrapper around...
Definition Location.h:76
Location getLoc()
The source location the operation was defined or derived from.
This provides public APIs that all operations should have.
Operation is the basic unit of execution within MLIR.
Definition Operation.h:87
Location getLoc()
The source location the operation was defined or derived from.
Definition Operation.h:240
DictionaryAttr getDiscardableAttrDictionary()
Return all of the discardable attributes on this operation as a DictionaryAttr.
Definition Operation.h:555
operand_type_range getOperandTypes()
Definition Operation.h:424
result_type_range getResultTypes()
Definition Operation.h:455
operand_range getOperands()
Returns an iterator on the underlying Value's.
Definition Operation.h:405
This class represents the benefit of a pattern match in a unitless scheme that ranges from 0 (very li...
A special type of RewriterBase that coordinates the application of a rewrite pattern on the current I...
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.
virtual void replaceOp(Operation *op, ValueRange newValues)
Replace the results of the given (original) operation with the specified list of values (replacements...
std::enable_if_t<!std::is_convertible< CallbackT, Twine >::value, LogicalResult > notifyMatchFailure(Location loc, CallbackT &&reasonCallback)
Used to notify the listener that the IR failed to be rewritten because of a match failure,...
OpTy replaceOpWithNewOp(Operation *op, Args &&...args)
Replace the results of the given (original) op with a new op that is created without verification (re...
This class provides an abstraction over the various different ranges of value types.
Definition TypeRange.h:40
Instances of the Type class are uniqued, have an immutable identifier and an optional mutable compone...
Definition Types.h:74
bool isF32() const
Definition Types.cpp:40
bool isF16() const
Definition Types.cpp:38
unsigned getIntOrFloatBitWidth() const
Return the bit width of an integer or a float type, assert failure on other types.
Definition Types.cpp:124
This class provides an abstraction over the different types of ranges over Values.
Definition ValueRange.h:392
Type front()
Return first type in the range.
Definition TypeRange.h:164
This class represents an instance of an SSA value in the MLIR system, representing a computable value...
Definition Value.h:96
Type getType() const
Return the type of this value.
Definition Value.h:105
DynamicAPInt floor(const Fraction &f)
Definition Fraction.h:77
Fraction abs(const Fraction &f)
Definition Fraction.h:107
Include the generated interface declarations.
void populatePolynomialApproximateErfcPattern(RewritePatternSet &patterns)
void populateMathF32ExpansionPatterns(RewritePatternSet &patterns, llvm::function_ref< bool(StringRef)> predicate, PatternBenefit=1)
void populatePolynomialApproximateErfPattern(RewritePatternSet &patterns)
void populatePolynomialApproximateTanhPattern(RewritePatternSet &patterns)
Type getElementTypeOrSelf(Type type)
Return the element type or return the type itself.
void populateMathPolynomialApproximationPatterns(RewritePatternSet &patterns, llvm::function_ref< bool(StringRef)> predicate, PatternBenefit=1)
ArrayRef< int64_t > sizes
ArrayRef< bool > scalableFlags
OpRewritePattern is a wrapper around RewritePattern that allows for matching and rewriting against an...
LogicalResult matchAndRewrite(math::ErfOp op, PatternRewriter &rewriter) const final
LogicalResult matchAndRewrite(math::ErfcOp op, PatternRewriter &rewriter) const final
Eliminates variable at the specified position using Fourier-Motzkin variable elimination.