[InstCombine] Missed optimization in math expression: sin(x) / cos(x) => tan(x)

[Quolyk]

Summary

This patch enables folding sin(x) / cos(x) -> tan(x), cos(x) / sin(x) -> 1 / tan(x) under -ffast-math flag

[hfinkel]

lib/Transforms/InstCombine/InstCombineMulDivRem.cpp

1456

Use match here?

1463

For the name here...

1. Predicate the transform on TLI->has(LibFunc_tan) (or tanf, tanl, depending on the type).
2. Use TLI->getName(LibFunc_tan) (or tanf, tanl, depending on the type).

1464

Needs

BuilderTy::FastMathFlagGuard Guard(Builder);

above this.

[Quolyk]

@hfinkel I'd like to thank you for your help and patience reviewing my patches, as I'm new to the community. I hope it's ok when I include you as a reviewer.

[hfinkel]

In D41286#960694, @Quolyk wrote:

@hfinkel I'd like to thank you for your help and patience reviewing my patches, as I'm new to the community. I hope it's ok when I include you as a reviewer.

It's certainly fine to include me as a review. Welcome to the community!

[spatel]

lib/Transforms/InstCombine/InstCombineMulDivRem.cpp

560–563

Can you move and reuse the very similar function that's currently in SimplifyLibCalls.cpp?

static bool hasUnaryFloatFn(const TargetLibraryInfo *TLI, Type *Ty,
                          LibFunc DoubleFn, LibFunc FloatFn,
                          LibFunc LongDoubleFn) {

See also:

/// Emit a call to the unary function named 'Name' (e.g.  'floor'). This
/// function is known to take a single of type matching 'Op' and returns one
/// value with the same type. If 'Op' is a long double, 'l' is added as the
/// suffix of name, if 'Op' is a float, we add a 'f' suffix.
Value *emitUnaryFloatFnCall(Value *Op, StringRef Name, IRBuilder<> &B,
                            const AttributeList &Attrs);

...in BuildLibCalls.h

1457–1460

Use m_Specific to simplify this:

// sin(a) / cos(a) -> tan(a)
Value *A;
if (match(Op0, m_Intrinsic<Intrinsic::sin>(m_Value(A))) &&
    match(Op1, m_Intrinsic<Intrinsic::cos>(m_Specific(A)))) {

test/Transforms/InstCombine/fdiv.ll

89–91

Please vary the fast-ness in these tests or add test(s) that show some variation. I'd prefer to see at least one test where we show the minimum case as we showed in one of the related patches - only the fdiv has relaxed math while the trig calls are strict.

[spatel]

lib/Transforms/InstCombine/InstCombineMulDivRem.cpp

1454

Should we also handle:

cos(a) / sin(a) -> 1 / tan(a)

?

test/Transforms/InstCombine/fdiv.ll

160–162

Do we need to declare tan functions for these tests?

[davide]

@scanon should sign off this.

lib/Transforms/InstCombine/InstCombineMulDivRem.cpp

1454

Please wait for @scanon opinion before implementing every possible 10th grade trigonometrical identity.

[davide]

For reference. this is what GCC generates (although it's unclear whether it's a good idea to follow them)
https://godbolt.org/g/YUUKGE

[spatel]

lib/Transforms/InstCombine/InstCombineMulDivRem.cpp

1454

Do you have some specific numerical concern here? As you've noted, this is a well-known math transform.

We can make cos(a) / sin(a) a 'TODO' if you think we should use a different transform.
https://stackoverflow.com/questions/3738384/stable-cotangent

[davide]

lib/Transforms/InstCombine/InstCombineMulDivRem.cpp

1454

The concern is that those transformation can overflow quite dramatically, even for -ffast-math.
The other concern is that I'm not a numerical expert, so I'd love to have this signed off from somebody who knows better than me.
The last concern is, again, we shouldn't pattern match every possible thing just because, it slows down the compiler without real benefit, so, do you know how this pattern is frequent?

[spatel]

lib/Transforms/InstCombine/InstCombineMulDivRem.cpp

1454

The compile-time concern is misguided.

This pattern, like every other "10th grade trigonometrical identity", should be optimized by an optimizing compiler because that's the job of an optimizing compiler.

These patterns can occur by way of templated code, inlining, or because the programmer may not be a computer performance expert. Think: scientists who are trying to model/simulate some math problem, but don't know much about perf...because again: that's the optimizing compiler's job.

If this patch or pass is causing a compile-time problem for you, please point to or file a bug. Obstructing patches like this is doubly bad when you're undermining the efforts of new contributors.

[davide]

lib/Transforms/InstCombine/InstCombineMulDivRem.cpp

1454

I'm afraid this is not entirely correct. The job of an optimizing compiler is that of making tradeoff on what to optimize, based on cost. I'm not obstructing this patch, I'm asking for a second opinion.
If you have a numerical explanation of why this patch can go in, I'll be happy to accept, otherwise I'll defer the review to @scanon.

[spatel]

lib/Transforms/InstCombine/InstCombineMulDivRem.cpp

1457–1458

We probably don't want to do this transform if both of the existing values have more than one use; we'd be trading a division for a libcall. If only one value has >1 use, it's probably still ok. Add more tests. :)

[spatel]

If you have a numerical explanation of why this patch can go in, I'll be happy to accept,

I'm not sure what you're looking for. Brute force sin/cos is equal to tan for every possible value?

#include <math.h>
#include <stdio.h>
#include <string.h>

int main() {
  unsigned int i;
  for (i = 0x00000000; i<0x7f800002; i++) {
    float f;
    memcpy(&f, &i, 4);
    float slow = sin(f) / cos(f);
    float fast = tan(f);
    if (slow != fast) {
      printf("\nsin(%f)/cos(%f) = %f, tan(%f) = %f\n", f, f, slow, f, fast);
    }
    if (i % (1024*1024*256) == 0) printf("0x%x...\n", i); 
  }

  return 0;
}

$  clang -O2 tan_checker.c ; ./a.out 
0x0...
0x10000000...
0x20000000...
0x30000000...
0x40000000...
0x50000000...
0x60000000...
0x70000000...

sin(inf)/cos(inf) = nan, tan(inf) = nan

sin(nan)/cos(nan) = nan, tan(nan) = nan

Note that I'm testing on macOS x86 with strict math (the sin and cos calls are replaced by __sincos_stret).
Double-check to make sure I haven't screwed anything up there, but this suggests we could do this transform without -ffast-math?

[davide]

In D41286#970342, @spatel wrote:

If you have a numerical explanation of why this patch can go in, I'll be happy to accept,

I'm not sure what you're looking for. Brute force sin/cos is equal to tan for every possible value?

I don't think this is feasible (for 64-bit values). -ffast-math is a little fun and complicated, e.g.

log(pow(x, y)) -> y*log(x)

which seems perfectly fine on paper, for x = -1, y = 4

log(pow(-1, 4)) -> 0
4*log(-1) -> NaN.

(courtesy of Steven)

My point is that reasoning about seemingly innocuous algebraic simplification turns out to be harder than expected, and therefore we should make a conscious choice on whether implement these after careful numerical analysis.

[spatel]

In D41286#970343, @davide wrote:

log(pow(x, y)) -> y*log(x)

How is this relevant? Log is only defined for positive numbers, while sin/cos/tan are valid across all numerical inputs.

I don't think there's anything I more I can say here; sorry @Quolyk , I tried. If @hfinkel , @efriedma , @andrew.w.kaylor or anyone else would like to comment that would be great.

[davide]

I don't feel qualified enough to say whether this can go in, somebody with fast-math experience should comment.

[davide]

In D41286#970372, @spatel wrote:

In D41286#970343, @davide wrote:

log(pow(x, y)) -> y*log(x)

How is this relevant? Log is only defined for positive numbers, while sin/cos/tan are valid across all numerical inputs.

My point is that fast-math implications can be non-trivial.

[spatel]

I changed the brute force float checker to test 1/tan(x), and it matches cos(x)/sin(x) in all cases on macOS 10.13. This is also true on Ubuntu 17.10 x86-64.

LGTM.

lib/Transforms/InstCombine/InstCombineMulDivRem.cpp

1474–1475

Could hoist this check to the first 'if' to reduce the duplication.

[Quolyk]

@spatel @davide thanks for review and comments

[spatel]

For reference, @bkramer improved (thanks!) the attribute propagation:
rL322284
rL322285

lib/Transforms/InstCombine/InstCombineMulDivRem.cpp

static bool isFMulOrFDivWithConstant(Value *V) {
   if (C0 && C1)
     return false;
 
   return (C0 && isFiniteNonZeroFp(C0)) || (C1 && isFiniteNonZeroFp(C1));
 }
 
 /// foldFMulConst() is a helper routine of InstCombiner::visitFMul().
 /// The input \p FMulOrDiv is a FMul/FDiv with one and only one operand
 /// being a constant (i.e. isFMulOrFDivWithConstant(FMulOrDiv) == true).
 /// This function is to simplify "FMulOrDiv * C" and returns the
 /// resulting expression. Note that this function could return NULL in
 /// case the constants cannot be folded into a normal floating-point.
 Value *InstCombiner::foldFMulConst(Instruction *FMulOrDiv, Constant *C,
                                    Instruction *InsertBefore) {
   assert(isFMulOrFDivWithConstant(FMulOrDiv) && "V is invalid");
 
   Value *Opnd0 = FMulOrDiv->getOperand(0);
   Value *Opnd1 = FMulOrDiv->getOperand(1);
 
   Constant *C0 = dyn_cast<Constant>(Opnd0);
if (AllowReassociate) {
     if (NewInst) {
       if (Instruction *T = dyn_cast<Instruction>(NewInst))
         T->setDebugLoc(I.getDebugLoc());
       SimpR->setFastMathFlags(I.getFastMathFlags());
       return SimpR;
     }
   }
 
+  // sin(a) / cos(a) -> tan(a)
+  if (AllowReassociate) {
+    if (IntrinsicInst *Sin = dyn_cast<IntrinsicInst>(Op0)) {
+      if (IntrinsicInst *Cos = dyn_cast<IntrinsicInst>(Op1)) {
+        if (Sin->getIntrinsicID() == Intrinsic::sin &&
+            Cos->getIntrinsicID() == Intrinsic::cos &&
+            Sin->getArgOperand(0) == Cos->getArgOperand(0)) {
+          Module *M = I.getModule();
+          Type *Ty = I.getType();
+          Value *Tan = M->getOrInsertFunction("tan", Ty, Ty);
+          Builder.setFastMathFlags(I.getFastMathFlags());
+          Value *TanCall = Builder.CreateCall(Tan, Cos->getArgOperand(0), "tan");
+          return replaceInstUsesWith(I, TanCall);
+        }
+      }
+    }
+  }
+
   Value *LHS;
   Value *RHS;
 
   // -x / -y -> x / y
   if (match(Op0, m_FNeg(m_Value(LHS))) && match(Op1, m_FNeg(m_Value(RHS)))) {
     I.setOperand(0, LHS);
     I.setOperand(1, RHS);
     return &I;
   }
 
   return nullptr;
 }

test/Transforms/InstCombine/fdiv.ll

 ; CHECK-LABEL @fdiv_fneg_fneg_fast(
 ; CHECK: %div = fdiv fast float %x, %y
 define float @fdiv_fneg_fneg_fast(float %x, float %y) {
   %x.fneg = fsub float -0.0, %x
   %y.fneg = fsub float -0.0, %y
   %div = fdiv fast float %x.fneg, %y.fneg
   ret float %div
 }
+
+; CHECK-LABEL @fdiv_sin_cos(
+; CHECK: call fast double @tan(double %a)
+define double @fdiv_sin_cos(double %a) {
+entry:
+  %0 = call fast double @llvm.sin.f64(double %a)
+  %1 = call fast double @llvm.cos.f64(double %a)
+  %div = fdiv fast double %0, %1
+  ret double %div
+}
+
+declare double @llvm.sin.f64(double) nounwind readnone speculatable
+declare double @llvm.cos.f64(double) nounwind readnone speculatable
+declare double @tan(double)