[complex] Teach Clang to preserve different-type operands to arithmetic operators where one type is a C complex type, and to emit both the efficient and correct implementation for complex arithmetic according to C11 Annex G using this extra...

[chandlerc]

Summary

...information.

For both multiply and divide the old code was writing a long-hand
reduced version of the math without any of the special handling of inf
and NaN recommended by the standard here. Instead of putting more
complexity here, this change does what GCC does which is to emit
a libcall for the fully general case.

However, the old code also failed to do the proper minimization of the
set of operations when there was a mixed complex and real operation. In
those cases, C provides a spec for much more minimal operations that are
valid. Clang now emits the exact suggested operations. This change isn't
*just* about performance though, without minimizing these operations, we
again lose the correct handling of infinities and NaNs. It is critical
that this happen in the frontend based on assymetric type operands to
complex math operations.

The performance implications of this change aren't trivial either. I've
run a set of benchmarks in Eigen, an open source mathematics library
that makes heavy use of complex. While a few have slowed down due to the
libcall being introduce, most sped up and some by a huge amount.

TODO: In order to make all of this work, also match the algorithm in the
constant evaluator to the one in the runtime library. Currently it is a broken
port of the simplifications from C's Annex G to the long-hand formulation of
the algorithm.

Splitting this patch up is very hard because none of this works without
the AST change to preserve non-complex operands. Sorry for the enormous
change.

[chandlerc]

A couple of other points on this patch.

• We could get a lot more out of these kinds of simplification, and

generally improve the sanity of Clang's model here, by adding support for
the C11 Imaginary type and using that to build up things. But this is more
work than I can sign up for right now, and I'm not the best person to do it.

• You might wonder why not add a complex type to the IR so that the

properties here are exposed at that level. There are some compelling
reasons not to do this IMO.

• Would need to "lower" this at some point which would be ugly. No chip

(that LLVM supports) has complex registers and arithmetic operations....

• Incredibly hard to change the IR in general, so bad cost/benefit

tradeoff.

• While this patch handles the simplifications to complex arithmetic we can

do when the typesystem itself shows one operand to be a real, it doesn't
handle large other categories such as when inlining or other optimizations
propagate away the imaginary part for example.

• All of these should be easily caught by a target library optimization
• I'll work on just such a target library pass to help polish off the

last few missed optimzations.

• This has uncovered some nasty behavior with SROA. Need to fix that.....

[scanon]

Thanks for working on this, Chandler, it’s an enormous improvement over the current state of affairs.

I am a bit concerned about the performance impact of unconditionally turning mul and div into libcalls; I expect we’ll see large regressions on a few focused benchmarks from this. In the long run, we likely want CX_LIMITED_RANGE on semantics to be the default, but that will require implementing support for the pragma or some other means to control the behavior (-fcx-limited-range=on?).

Nonetheless, this is a large performance improvement for cases that are actually used in a relatively mainstream project (and it’s strictly an improvement w.r.t. correctness), so I am perfectly content to put such issues into the bug tracker to be handled separately.

[hfinkel]

I agree with Steve, this looks good. To add one thing: I'd also like to skip the libcalls when possible in fast-math mode.

lib/CodeGen/CGExprComplex.cpp

590

Indent?

[rsmith]

I'd like to see some added tests for constant folding here.

lib/AST/ExprConstant.cpp

7694–7697

How do we reach this case? The below handling for real operands calls EvaluateFloat rather than EvaluateComplex.

7883

You should assign true to these somewhere =)

7910–7911

Maybe assert(!(LHSReal && RHSReal)); here?

7933

Is changeSign right here? Should 0 - (0+0i) be 0+0i or 0-0i?

8030–8031

Indentation looks off here. Maybe it's just Phab?

8038

Does this do the right thing with negative zeros? Eg, 1 / (1 + 0i) will give 1-0i but should presumably be 1+0i.

lib/CodeGen/CGExprComplex.cpp

556–557

Same question as before regarding negative zeros.

593

My 80-col sense is tingling.

618–619

Do we need to do anything about underflow/overflow in the constant evaluator?

[chandlerc]

Completely re-work the constant folding code to correctly follow the same logic
as the library calls do, but in terms of APFloat. Add extensive compile time
folding tests. These tests have to be written in C++ despite testing the
C complex number functionality because we need C++ static assertions in order
to fold these complex things. Yes, pun intended. =]

[chandlerc]

Run clang-format over the patch to fix some lingering formatting issues.

[chandlerc]

Address Richard's comments.

[chandlerc]

All comments addressed. Answers to questions below. Fresh patches uploaded and ready for review.

The constant folding is now in good shape I deem it. Took some hacking on APFloat to get there. Added test cases as well. Had fun testing complex division with Wolfram Alpha. =] We seem to get the right answers.

One area I'm not testing heavily is sign propagation in the face of infinities. I'd appreciate if someone (yea, I'm looking at you Steve) could contribute more thorough testing of the sign propagation properties of these operations, and/or tighten or relax the assertions on the real vs. imaginary components to match what numerics folks actually want to enshrine in the constant folder. Hopefully my test case is good enough to make it clear how to extend things going forward.

I also have plans to address both the libcall performance concerns in general and the fast-math concerns in follow-up patches, but I'd like to get a baseline of correctness in first. =]

lib/AST/ExprConstant.cpp

7694–7697

We don't. This was a speculative edit I shouldn't have made.

7883

Indeed. As you noticed (and I tried but failed to make clear) once I knew I would have to re-implement this from the principled version in our libcalls, I stopped testing it and this was completely broken and wrong.

7910–7911

Sure. I had some of that later, but nicer to state it up front.

7933

Quite surprisingly, I believe this is correct according to the table in G.5.2p2. It stipulates that a complex value subtracted from a real value as x - (u + iv) == (x - u) - iv. Because subtraction is represented as the sum of a negative, especially for complex numbers, this indicates that 'v' must be negated here, even if that would form 0.0 - 0.0i. The 'negative zero' here which is surprising resulting from subtraction of two values with the same sign isn't a big deal because in complex form it just feeds another subtraction.

8030–8031

No, formatting was off here.

lib/CodeGen/CGExprComplex.cpp

556–557

Same answer.

590

Fixed.

593

Fixed....

BY THE POWER OF CLANG-FORMAT!

618–619

Not necessarily. Unlike with integers, such events should (IMO) simply produce the relevant infinity or zero value.

My inclination is that the constant folder should error when a full operation produces a NaN. This would let the complex arithmetic logic shield it from certain NaN conditions, but still bail when we go Too Far.

[rsmith]

lib/AST/ExprConstant.cpp

8049–8050

It looks like B is not initialized if LHSReal is true here.

[chandlerc]

Closed by commit rL219557 (authored by @chandlerc).

lib/AST/ExprConstant.cpp

public:
   bool VisitBinaryOperator(const BinaryOperator *E);
   bool VisitUnaryOperator(const UnaryOperator *E);
   bool VisitInitListExpr(const InitListExpr *E);
 };
 } // end anonymous namespace
 
 static bool EvaluateComplex(const Expr *E, ComplexValue &Result,
                             EvalInfo &Info) {
-  assert(E->isRValue() && E->getType()->isAnyComplexType());
+  assert(E->isRValue() && "Can only evaluate R-values!");
+  assert((E->getType()->isAnyComplexType() ||
+          E->getType()->isRealFloatingType()) &&
+         "Complex expressions can only be composed of complex and real "
+         "floating point types.");
   return ComplexExprEvaluator(Info, Result).Visit(E);
 }
 
 bool ComplexExprEvaluator::ZeroInitialization(const Expr *E) {
   QualType ElemTy = E->getType()->castAs<ComplexType>()->getElementType();
   if (ElemTy->isRealFloatingType()) {
     Result.makeComplexFloat();
     APFloat Zero = APFloat::getZero(Info.Ctx.getFloatTypeSemantics(ElemTy));
bool ComplexExprEvaluator::VisitCastExpr(const CastExpr *E) {
 
   llvm_unreachable("unknown cast resulting in complex value");
 }
 
 bool ComplexExprEvaluator::VisitBinaryOperator(const BinaryOperator *E) {
   if (E->isPtrMemOp() || E->isAssignmentOp() || E->getOpcode() == BO_Comma)
     return ExprEvaluatorBaseTy::VisitBinaryOperator(E);
 
-  bool LHSOK = Visit(E->getLHS());
+  // Track whether the LHS or RHS is real at the type system level. When this is
+  // the case we can simplify our evaluation strategy.
+  bool LHSReal = false, RHSReal = false;
+
+  bool LHSOK;
+  if (E->getLHS()->getType()->isRealFloatingType()) {
+    LHSReal = true;
+    APFloat &Real = Result.FloatReal;
+    LHSOK = EvaluateFloat(E->getLHS(), Real, Info);
+    if (LHSOK) {
+      Result.makeComplexFloat();
+      Result.FloatImag = APFloat(Real.getSemantics());
+    }
+  } else {
+    LHSOK = Visit(E->getLHS());
+  }
   if (!LHSOK && !Info.keepEvaluatingAfterFailure())
     return false;
 
   ComplexValue RHS;
-  if (!EvaluateComplex(E->getRHS(), RHS, Info) || !LHSOK)
+  if (E->getRHS()->getType()->isRealFloatingType()) {
+    RHSReal = true;
+    APFloat &Real = RHS.FloatReal;
+    if (!EvaluateFloat(E->getRHS(), Real, Info) || !LHSOK)
+      return false;
+    RHS.makeComplexFloat();
+    RHS.FloatImag = APFloat(Real.getSemantics());
+  } else if (!EvaluateComplex(E->getRHS(), RHS, Info) || !LHSOK)
     return false;
 
-  assert(Result.isComplexFloat() == RHS.isComplexFloat() &&
-         "Invalid operands to binary operator.");
   switch (E->getOpcode()) {
   default: return Error(E);
   case BO_Add:
     if (Result.isComplexFloat()) {
       Result.getComplexFloatReal().add(RHS.getComplexFloatReal(),
                                        APFloat::rmNearestTiesToEven);
+      if (LHSReal)
+        Result.getComplexFloatImag() = RHS.getComplexFloatImag();
+      else if (!RHSReal)
         Result.getComplexFloatImag().add(RHS.getComplexFloatImag(),
                                          APFloat::rmNearestTiesToEven);
     } else {
       Result.getComplexIntReal() += RHS.getComplexIntReal();
       Result.getComplexIntImag() += RHS.getComplexIntImag();
     }
     break;
   case BO_Sub:
     if (Result.isComplexFloat()) {
       Result.getComplexFloatReal().subtract(RHS.getComplexFloatReal(),
                                             APFloat::rmNearestTiesToEven);
+      if (LHSReal) {
+        Result.getComplexFloatImag() = RHS.getComplexFloatImag();
+        Result.getComplexFloatImag().changeSign();
+      } else if (!RHSReal) {
         Result.getComplexFloatImag().subtract(RHS.getComplexFloatImag(),
                                               APFloat::rmNearestTiesToEven);
+      }
     } else {
       Result.getComplexIntReal() -= RHS.getComplexIntReal();
       Result.getComplexIntImag() -= RHS.getComplexIntImag();
     }
     break;
   case BO_Mul:
     if (Result.isComplexFloat()) {
+      // This is an implementation of complex multiplication according to the
+      // constraints laid out in C11 Annex G. The implemantion uses the
+      // following naming scheme:
+      //   (a + ib) * (c + id)
       ComplexValue LHS = Result;
-      APFloat &LHS_r = LHS.getComplexFloatReal();
-      APFloat &LHS_i = LHS.getComplexFloatImag();
-      APFloat &RHS_r = RHS.getComplexFloatReal();
-      APFloat &RHS_i = RHS.getComplexFloatImag();
-
-      APFloat Tmp = LHS_r;
-      Tmp.multiply(RHS_r, APFloat::rmNearestTiesToEven);
-      Result.getComplexFloatReal() = Tmp;
-      Tmp = LHS_i;
-      Tmp.multiply(RHS_i, APFloat::rmNearestTiesToEven);
-      Result.getComplexFloatReal().subtract(Tmp, APFloat::rmNearestTiesToEven);
-
-      Tmp = LHS_r;
-      Tmp.multiply(RHS_i, APFloat::rmNearestTiesToEven);
-      Result.getComplexFloatImag() = Tmp;
-      Tmp = LHS_i;
-      Tmp.multiply(RHS_r, APFloat::rmNearestTiesToEven);
-      Result.getComplexFloatImag().add(Tmp, APFloat::rmNearestTiesToEven);
+      APFloat &A = LHS.getComplexFloatReal();
+      APFloat &B = LHS.getComplexFloatImag();
+      APFloat &C = RHS.getComplexFloatReal();
+      APFloat &D = RHS.getComplexFloatImag();
+      APFloat &ResR = Result.getComplexFloatReal();
+      APFloat &ResI = Result.getComplexFloatImag();
+      if (LHSReal) {
+        assert(!RHSReal && "Cannot have two real operands for a complex op!");
+        ResR = A * C;
+        ResI = A * D;
+      } else if (RHSReal) {
+        ResR = C * A;
+        ResI = C * B;
+      } else {
+        // In the fully general case, we need to handle NaNs and infinities
+        // robustly.
+        APFloat AC = A * C;
+        APFloat BD = B * D;
+        APFloat AD = A * D;
+        APFloat BC = B * C;
+        ResR = AC - BD;
+        ResI = AD + BC;
+        if (ResR.isNaN() && ResI.isNaN()) {
+          bool Recalc = false;
+          if (A.isInfinity() || B.isInfinity()) {
+            A = APFloat::copySign(
+                APFloat(A.getSemantics(), A.isInfinity() ? 1 : 0), A);
+            B = APFloat::copySign(
+                APFloat(B.getSemantics(), B.isInfinity() ? 1 : 0), B);
+            if (C.isNaN())
+              C = APFloat::copySign(APFloat(C.getSemantics()), C);
+            if (D.isNaN())
+              D = APFloat::copySign(APFloat(D.getSemantics()), D);
+            Recalc = true;
+          }
+          if (C.isInfinity() || D.isInfinity()) {
+            C = APFloat::copySign(
+                APFloat(C.getSemantics(), C.isInfinity() ? 1 : 0), C);
+            D = APFloat::copySign(
+                APFloat(D.getSemantics(), D.isInfinity() ? 1 : 0), D);
+            if (A.isNaN())
+              A = APFloat::copySign(APFloat(A.getSemantics()), A);
+            if (B.isNaN())
+              B = APFloat::copySign(APFloat(B.getSemantics()), B);
+            Recalc = true;
+          }
+          if (!Recalc && (AC.isInfinity() || BD.isInfinity() ||
+                          AD.isInfinity() || BC.isInfinity())) {
+            if (A.isNaN())
+              A = APFloat::copySign(APFloat(A.getSemantics()), A);
+            if (B.isNaN())
+              B = APFloat::copySign(APFloat(B.getSemantics()), B);
+            if (C.isNaN())
+              C = APFloat::copySign(APFloat(C.getSemantics()), C);
+            if (D.isNaN())
+              D = APFloat::copySign(APFloat(D.getSemantics()), D);
+            Recalc = true;
+          }
+          if (Recalc) {
+            ResR = APFloat::getInf(A.getSemantics()) * (A * C - B * D);
+            ResI = APFloat::getInf(A.getSemantics()) * (A * D + B * C);
+          }
+        }
+      }
     } else {
       ComplexValue LHS = Result;
       Result.getComplexIntReal() =
         (LHS.getComplexIntReal() * RHS.getComplexIntReal() -
          LHS.getComplexIntImag() * RHS.getComplexIntImag());
       Result.getComplexIntImag() =
         (LHS.getComplexIntReal() * RHS.getComplexIntImag() +
          LHS.getComplexIntImag() * RHS.getComplexIntReal());
     }
     break;
   case BO_Div:
     if (Result.isComplexFloat()) {
+      // This is an implementation of complex division according to the
+      // constraints laid out in C11 Annex G. The implemantion uses the
+      // following naming scheme:
+      //   (a + ib) / (c + id)
       ComplexValue LHS = Result;
-      APFloat &LHS_r = LHS.getComplexFloatReal();
-      APFloat &LHS_i = LHS.getComplexFloatImag();
-      APFloat &RHS_r = RHS.getComplexFloatReal();
-      APFloat &RHS_i = RHS.getComplexFloatImag();
-      APFloat &Res_r = Result.getComplexFloatReal();
-      APFloat &Res_i = Result.getComplexFloatImag();
-
-      APFloat Den = RHS_r;
-      Den.multiply(RHS_r, APFloat::rmNearestTiesToEven);
-      APFloat Tmp = RHS_i;
-      Tmp.multiply(RHS_i, APFloat::rmNearestTiesToEven);
-      Den.add(Tmp, APFloat::rmNearestTiesToEven);
-
-      Res_r = LHS_r;
-      Res_r.multiply(RHS_r, APFloat::rmNearestTiesToEven);
-      Tmp = LHS_i;
-      Tmp.multiply(RHS_i, APFloat::rmNearestTiesToEven);
-      Res_r.add(Tmp, APFloat::rmNearestTiesToEven);
-      Res_r.divide(Den, APFloat::rmNearestTiesToEven);
-
-      Res_i = LHS_i;
-      Res_i.multiply(RHS_r, APFloat::rmNearestTiesToEven);
-      Tmp = LHS_r;
-      Tmp.multiply(RHS_i, APFloat::rmNearestTiesToEven);
-      Res_i.subtract(Tmp, APFloat::rmNearestTiesToEven);
-      Res_i.divide(Den, APFloat::rmNearestTiesToEven);
+      APFloat &A = LHS.getComplexFloatReal();
+      APFloat &B = LHS.getComplexFloatImag();
+      APFloat &C = RHS.getComplexFloatReal();
+      APFloat &D = RHS.getComplexFloatImag();
+      APFloat &ResR = Result.getComplexFloatReal();
+      APFloat &ResI = Result.getComplexFloatImag();
+      if (RHSReal) {
+        ResR = A / C;
+        ResI = B / C;
+      } else {
+        int DenomLogB = 0;
+        APFloat MaxCD = maxnum(abs(C), abs(D));
+        if (MaxCD.isFinite()) {
+          DenomLogB = ilogb(MaxCD);
+          C = scalbn(C, -DenomLogB);
+          D = scalbn(D, -DenomLogB);
+        }
+        APFloat Denom = C * C + D * D;
+        ResR = scalbn((A * C + B * D) / Denom, -DenomLogB);
+        ResI = scalbn((B * C - A * D) / Denom, -DenomLogB);
+        if (ResR.isNaN() && ResI.isNaN()) {
+          if (Denom.isPosZero() && (!A.isNaN() || !B.isNaN())) {
+            ResR = APFloat::getInf(ResR.getSemantics(), C.isNegative()) * A;
+            ResI = APFloat::getInf(ResR.getSemantics(), C.isNegative()) * B;
+          } else if ((A.isInfinity() || B.isInfinity()) && C.isFinite() &&
+                     D.isFinite()) {
+            A = APFloat::copySign(
+                APFloat(A.getSemantics(), A.isInfinity() ? 1 : 0), A);
+            B = APFloat::copySign(
+                APFloat(B.getSemantics(), B.isInfinity() ? 1 : 0), B);
+            ResR = APFloat::getInf(ResR.getSemantics()) * (A * C + B * D);
+            ResI = APFloat::getInf(ResI.getSemantics()) * (B * C - A * D);
+          } else if (MaxCD.isInfinity() && A.isFinite() && B.isFinite()) {
+            C = APFloat::copySign(
+                APFloat(C.getSemantics(), C.isInfinity() ? 1 : 0), C);
+            D = APFloat::copySign(
+                APFloat(D.getSemantics(), D.isInfinity() ? 1 : 0), D);
+            ResR = APFloat::getZero(ResR.getSemantics()) * (A * C + B * D);
+            ResI = APFloat::getZero(ResI.getSemantics()) * (B * C - A * D);
+          }
+        }
+      }
     } else {
       if (RHS.getComplexIntReal() == 0 && RHS.getComplexIntImag() == 0)
         return Error(E, diag::note_expr_divide_by_zero);
 
       ComplexValue LHS = Result;
       APSInt Den = RHS.getComplexIntReal() * RHS.getComplexIntReal() +
         RHS.getComplexIntImag() * RHS.getComplexIntImag();
       Result.getComplexIntReal() =

lib/CodeGen/CGExprComplex.cpp

ComplexPairTy EmitCompoundAssign(const CompoundAssignOperator *E,
                                    ComplexPairTy (ComplexExprEmitter::*Func)
                                    (const BinOpInfo &));
 
   ComplexPairTy EmitBinAdd(const BinOpInfo &Op);
   ComplexPairTy EmitBinSub(const BinOpInfo &Op);
   ComplexPairTy EmitBinMul(const BinOpInfo &Op);
   ComplexPairTy EmitBinDiv(const BinOpInfo &Op);
 
+  ComplexPairTy EmitComplexBinOpLibCall(StringRef LibCallName,
+                                        const BinOpInfo &Op);
+
+
   ComplexPairTy VisitBinAdd(const BinaryOperator *E) {
     return EmitBinAdd(EmitBinOps(E));
   }
   ComplexPairTy VisitBinSub(const BinaryOperator *E) {
     return EmitBinSub(EmitBinOps(E));
   }
   ComplexPairTy VisitBinMul(const BinaryOperator *E) {
     return EmitBinMul(EmitBinOps(E));
ComplexPairTy ComplexExprEmitter::VisitUnaryNot(const UnaryOperator *E) {
   return ComplexPairTy(Op.first, ResI);
 }
 
 ComplexPairTy ComplexExprEmitter::EmitBinAdd(const BinOpInfo &Op) {
   llvm::Value *ResR, *ResI;
 
   if (Op.LHS.first->getType()->isFloatingPointTy()) {
     ResR = Builder.CreateFAdd(Op.LHS.first,  Op.RHS.first,  "add.r");
+    if (Op.LHS.second && Op.RHS.second)
       ResI = Builder.CreateFAdd(Op.LHS.second, Op.RHS.second, "add.i");
+    else
+      ResI = Op.LHS.second ? Op.LHS.second : Op.RHS.second;
+    assert(ResI && "Only one operand may be real!");
   } else {
     ResR = Builder.CreateAdd(Op.LHS.first,  Op.RHS.first,  "add.r");
+    assert(Op.LHS.second && Op.RHS.second &&
+           "Both operands of integer complex operators must be complex!");
     ResI = Builder.CreateAdd(Op.LHS.second, Op.RHS.second, "add.i");
   }
   return ComplexPairTy(ResR, ResI);
 }
 
 ComplexPairTy ComplexExprEmitter::EmitBinSub(const BinOpInfo &Op) {
   llvm::Value *ResR, *ResI;
   if (Op.LHS.first->getType()->isFloatingPointTy()) {
-    ResR = Builder.CreateFSub(Op.LHS.first,  Op.RHS.first,  "sub.r");
+    ResR = Builder.CreateFSub(Op.LHS.first, Op.RHS.first, "sub.r");
+    if (Op.LHS.second && Op.RHS.second)
       ResI = Builder.CreateFSub(Op.LHS.second, Op.RHS.second, "sub.i");
+    else
+      ResI = Op.LHS.second ? Op.LHS.second
+                           : Builder.CreateFNeg(Op.RHS.second, "sub.i");
+    assert(ResI && "Only one operand may be real!");
   } else {
-    ResR = Builder.CreateSub(Op.LHS.first,  Op.RHS.first,  "sub.r");
+    ResR = Builder.CreateSub(Op.LHS.first, Op.RHS.first, "sub.r");
+    assert(Op.LHS.second && Op.RHS.second &&
+           "Both operands of integer complex operators must be complex!");
     ResI = Builder.CreateSub(Op.LHS.second, Op.RHS.second, "sub.i");
   }
   return ComplexPairTy(ResR, ResI);
 }
 
+/// \brief Emit a libcall for a binary operation on complex types.
+ComplexPairTy ComplexExprEmitter::EmitComplexBinOpLibCall(StringRef LibCallName,
+                                                          const BinOpInfo &Op) {
+  CallArgList Args;
+  Args.add(RValue::get(Op.LHS.first), Op.Ty->castAs<ComplexType>()->getElementType());
+  Args.add(RValue::get(Op.LHS.second), Op.Ty->castAs<ComplexType>()->getElementType());
+  Args.add(RValue::get(Op.RHS.first), Op.Ty->castAs<ComplexType>()->getElementType());
+  Args.add(RValue::get(Op.RHS.second), Op.Ty->castAs<ComplexType>()->getElementType());
+
+  // We *must* use the full CG function call building logic here because the
+  // complex type has special ABI handling.
+  const CGFunctionInfo &FuncInfo = CGF.CGM.getTypes().arrangeFreeFunctionCall(
+      Op.Ty, Args, FunctionType::ExtInfo(), RequiredArgs::All);
+  llvm::FunctionType *FTy = CGF.CGM.getTypes().GetFunctionType(FuncInfo);
+  llvm::Constant *Func = CGF.CGM.CreateRuntimeFunction(FTy, LibCallName);
+
+  llvm::Value *ArgVals[] = {Op.LHS.first, Op.LHS.second, Op.RHS.first,
+                            Op.RHS.second};
+  llvm::Value *Result = CGF.EmitRuntimeCall(Func, ArgVals);
+
+  llvm::Value *ResR, *ResI;
+  if (Result->getType()->isVectorTy()) {
+  ResR = CGF.Builder.CreateExtractElement(Result, CGF.Builder.getInt32(0));
+  ResI = CGF.Builder.CreateExtractElement(Result, CGF.Builder.getInt32(1));
+  } else {
+    assert(Result->getType()->isAggregateType() && "Only vector and aggregate libcall returns are supported!");
+    unsigned ResRIndices[] = {0};
+    ResR = CGF.Builder.CreateExtractValue(Result, ResRIndices);
+    unsigned ResIIndices[] = {1};
+    ResI = CGF.Builder.CreateExtractValue(Result, ResIIndices);
+  }
+  return ComplexPairTy(ResR, ResI);
+}
 
+// See C11 Annex G.5.1 for the semantics of multiplicative operators on complex
+// typed values.
 ComplexPairTy ComplexExprEmitter::EmitBinMul(const BinOpInfo &Op) {
   using llvm::Value;
   Value *ResR, *ResI;
 
   if (Op.LHS.first->getType()->isFloatingPointTy()) {
-    Value *ResRl = Builder.CreateFMul(Op.LHS.first, Op.RHS.first, "mul.rl");
-    Value *ResRr = Builder.CreateFMul(Op.LHS.second, Op.RHS.second,"mul.rr");
-    ResR  = Builder.CreateFSub(ResRl, ResRr, "mul.r");
-
-    Value *ResIl = Builder.CreateFMul(Op.LHS.second, Op.RHS.first, "mul.il");
-    Value *ResIr = Builder.CreateFMul(Op.LHS.first, Op.RHS.second, "mul.ir");
-    ResI  = Builder.CreateFAdd(ResIl, ResIr, "mul.i");
+    // The general formulation is:
+    // (a + ib) * (c + id) = (a * c - b * d) + i(a * d + b * c)
+    //
+    // But we can fold away components which would be zero due to a real
+    // operand according to C11 Annex G.5.1p2.
+    // FIXME: C11 also provides for imaginary types which would allow folding
+    // still more of this within the type system.
+
+    if (Op.LHS.second && Op.RHS.second) {
+      // If both operands are complex, delegate to a libcall which works to
+      // prevent underflow and overflow.
+      StringRef LibCallName;
+      switch (Op.LHS.first->getType()->getTypeID()) {
+      default:
+        llvm_unreachable("Unsupported floating point type!");
+      case llvm::Type::HalfTyID:
+        return EmitComplexBinOpLibCall("__mulhc3", Op);
+      case llvm::Type::FloatTyID:
+        return EmitComplexBinOpLibCall("__mulsc3", Op);
+      case llvm::Type::DoubleTyID:
+        return EmitComplexBinOpLibCall("__muldc3", Op);
+      case llvm::Type::X86_FP80TyID:
+        return EmitComplexBinOpLibCall("__mulxc3", Op);
+      }
+    }
+    assert((Op.LHS.second || Op.RHS.second) &&
+           "At least one operand must be complex!");
+
+    // If either of the operands is a real rather than a complex, the
+    // imaginary component is ignored when computing the real component of the
+    // result.
+    ResR = Builder.CreateFMul(Op.LHS.first, Op.RHS.first, "mul.rl");
+
+    ResI = Op.LHS.second
+               ? Builder.CreateFMul(Op.LHS.second, Op.RHS.first, "mul.il")
+               : Builder.CreateFMul(Op.LHS.first, Op.RHS.second, "mul.ir");
   } else {
+    assert(Op.LHS.second && Op.RHS.second &&
+           "Both operands of integer complex operators must be complex!");
     Value *ResRl = Builder.CreateMul(Op.LHS.first, Op.RHS.first, "mul.rl");
-    Value *ResRr = Builder.CreateMul(Op.LHS.second, Op.RHS.second,"mul.rr");
-    ResR  = Builder.CreateSub(ResRl, ResRr, "mul.r");
+    Value *ResRr = Builder.CreateMul(Op.LHS.second, Op.RHS.second, "mul.rr");
+    ResR = Builder.CreateSub(ResRl, ResRr, "mul.r");
 
     Value *ResIl = Builder.CreateMul(Op.LHS.second, Op.RHS.first, "mul.il");
     Value *ResIr = Builder.CreateMul(Op.LHS.first, Op.RHS.second, "mul.ir");
-    ResI  = Builder.CreateAdd(ResIl, ResIr, "mul.i");
+    ResI = Builder.CreateAdd(ResIl, ResIr, "mul.i");
   }
   return ComplexPairTy(ResR, ResI);
 }
 
+// See C11 Annex G.5.1 for the semantics of multiplicative operators on complex
+// typed values.
 ComplexPairTy ComplexExprEmitter::EmitBinDiv(const BinOpInfo &Op) {
   llvm::Value *LHSr = Op.LHS.first, *LHSi = Op.LHS.second;
   llvm::Value *RHSr = Op.RHS.first, *RHSi = Op.RHS.second;
 
 
   llvm::Value *DSTr, *DSTi;
-  if (Op.LHS.first->getType()->isFloatingPointTy()) {
-    // (a+ib) / (c+id) = ((ac+bd)/(cc+dd)) + i((bc-ad)/(cc+dd))
-    llvm::Value *Tmp1 = Builder.CreateFMul(LHSr, RHSr); // a*c
-    llvm::Value *Tmp2 = Builder.CreateFMul(LHSi, RHSi); // b*d
-    llvm::Value *Tmp3 = Builder.CreateFAdd(Tmp1, Tmp2); // ac+bd
-
-    llvm::Value *Tmp4 = Builder.CreateFMul(RHSr, RHSr); // c*c
-    llvm::Value *Tmp5 = Builder.CreateFMul(RHSi, RHSi); // d*d
-    llvm::Value *Tmp6 = Builder.CreateFAdd(Tmp4, Tmp5); // cc+dd
-
-    llvm::Value *Tmp7 = Builder.CreateFMul(LHSi, RHSr); // b*c
-    llvm::Value *Tmp8 = Builder.CreateFMul(LHSr, RHSi); // a*d
-    llvm::Value *Tmp9 = Builder.CreateFSub(Tmp7, Tmp8); // bc-ad
+  if (LHSr->getType()->isFloatingPointTy()) {
+    // If we have a complex operand on the RHS, we delegate to a libcall to
+    // handle all of the complexities and minimize underflow/overflow cases.
+    //
+    // FIXME: We would be able to avoid the libcall in many places if we
+    // supported imaginary types in addition to complex types.
+    if (RHSi) {
+      BinOpInfo LibCallOp = Op;
+      // If LHS was a real, supply a null imaginary part.
+      if (!LHSi)
+        LibCallOp.LHS.second = llvm::Constant::getNullValue(LHSr->getType());
 
-    DSTr = Builder.CreateFDiv(Tmp3, Tmp6);
-    DSTi = Builder.CreateFDiv(Tmp9, Tmp6);
+      StringRef LibCallName;
+      switch (LHSr->getType()->getTypeID()) {
+      default:
+        llvm_unreachable("Unsupported floating point type!");
+      case llvm::Type::HalfTyID:
+        return EmitComplexBinOpLibCall("__divhc3", LibCallOp);
+      case llvm::Type::FloatTyID:
+        return EmitComplexBinOpLibCall("__divsc3", LibCallOp);
+      case llvm::Type::DoubleTyID:
+        return EmitComplexBinOpLibCall("__divdc3", LibCallOp);
+      case llvm::Type::X86_FP80TyID:
+        return EmitComplexBinOpLibCall("__divxc3", LibCallOp);
+      }
+    }
+    assert(LHSi && "Can have at most one non-complex operand!");
+
+    DSTr = Builder.CreateFDiv(LHSr, RHSr);
+    DSTi = Builder.CreateFDiv(LHSi, RHSr);
   } else {
+    assert(Op.LHS.second && Op.RHS.second &&
+           "Both operands of integer complex operators must be complex!");
     // (a+ib) / (c+id) = ((ac+bd)/(cc+dd)) + i((bc-ad)/(cc+dd))
     llvm::Value *Tmp1 = Builder.CreateMul(LHSr, RHSr); // a*c
     llvm::Value *Tmp2 = Builder.CreateMul(LHSi, RHSi); // b*d
     llvm::Value *Tmp3 = Builder.CreateAdd(Tmp1, Tmp2); // ac+bd
 
     llvm::Value *Tmp4 = Builder.CreateMul(RHSr, RHSr); // c*c
     llvm::Value *Tmp5 = Builder.CreateMul(RHSi, RHSi); // d*d
     llvm::Value *Tmp6 = Builder.CreateAdd(Tmp4, Tmp5); // cc+dd
ComplexPairTy ComplexExprEmitter::EmitBinDiv(const BinOpInfo &Op) {
   return ComplexPairTy(DSTr, DSTi);
 }
 
 ComplexExprEmitter::BinOpInfo
 ComplexExprEmitter::EmitBinOps(const BinaryOperator *E) {
   TestAndClearIgnoreReal();
   TestAndClearIgnoreImag();
   BinOpInfo Ops;
+  if (E->getLHS()->getType()->isRealFloatingType())
+    Ops.LHS = ComplexPairTy(CGF.EmitScalarExpr(E->getLHS()), nullptr);
+  else
     Ops.LHS = Visit(E->getLHS());
+  if (E->getRHS()->getType()->isRealFloatingType())
+    Ops.RHS = ComplexPairTy(CGF.EmitScalarExpr(E->getRHS()), nullptr);
+  else
     Ops.RHS = Visit(E->getRHS());
+
   Ops.Ty = E->getType();
   return Ops;
 }
 
 
 LValue ComplexExprEmitter::
 EmitCompoundAssignLValue(const CompoundAssignOperator *E,
           ComplexPairTy (ComplexExprEmitter::*Func)(const BinOpInfo&),
                          RValue &Val) {
   TestAndClearIgnoreReal();
   TestAndClearIgnoreImag();
   QualType LHSTy = E->getLHS()->getType();
 
   BinOpInfo OpInfo;
 
   // Load the RHS and LHS operands.
   // __block variables need to have the rhs evaluated first, plus this should
   // improve codegen a little.
   OpInfo.Ty = E->getComputationResultType();
+  QualType ComplexElementTy = cast<ComplexType>(OpInfo.Ty)->getElementType();
 
   // The RHS should have been converted to the computation type.
-  assert(OpInfo.Ty->isAnyComplexType());
+  if (E->getRHS()->getType()->isRealFloatingType()) {
+    assert(CGF.getContext().hasSameUnqualifiedType(ComplexElementTy,
+                                                   E->getRHS()->getType()));
+    OpInfo.RHS = ComplexPairTy(CGF.EmitScalarExpr(E->getRHS()), nullptr);
+  } else {
     assert(CGF.getContext().hasSameUnqualifiedType(OpInfo.Ty,
                                                    E->getRHS()->getType()));
     OpInfo.RHS = Visit(E->getRHS());
+  }
 
   LValue LHS = CGF.EmitLValue(E->getLHS());
 
   // Load from the l-value and convert it.
   if (LHSTy->isAnyComplexType()) {
     ComplexPairTy LHSVal = EmitLoadOfLValue(LHS, E->getExprLoc());
     OpInfo.LHS = EmitComplexToComplexCast(LHSVal, LHSTy, OpInfo.Ty);
   } else {
     llvm::Value *LHSVal = CGF.EmitLoadOfScalar(LHS, E->getExprLoc());
+    // For floating point real operands we can directly pass the scalar form
+    // to the binary operator emission and potentially get more efficient code.
+    if (LHSTy->isRealFloatingType()) {
+      if (!CGF.getContext().hasSameUnqualifiedType(ComplexElementTy, LHSTy))
+        LHSVal = CGF.EmitScalarConversion(LHSVal, LHSTy, ComplexElementTy);
+      OpInfo.LHS = ComplexPairTy(LHSVal, nullptr);
+    } else {
       OpInfo.LHS = EmitScalarToComplexCast(LHSVal, LHSTy, OpInfo.Ty);
     }
+  }
 
   // Expand the binary operator.
   ComplexPairTy Result = (this->*Func)(OpInfo);
 
   // Truncate the result and store it into the LHS lvalue.
   if (LHSTy->isAnyComplexType()) {
     ComplexPairTy ResVal = EmitComplexToComplexCast(Result, OpInfo.Ty, LHSTy);
     EmitStoreOfComplex(ResVal, LHS, /*isInit*/ false);

lib/Sema/SemaExpr.cpp

static bool handleIntegerToComplexFloatConversion(Sema &S, ExprResult &IntExpr,
   } else {
     assert(IntTy->isComplexIntegerType());
     IntExpr = S.ImpCastExprToType(IntExpr.get(), ComplexTy,
                                   CK_IntegralComplexToFloatingComplex);
   }
   return false;
 }
 
-/// \brief Takes two complex float types and converts them to the same type.
-/// Helper function of UsualArithmeticConversions()
-static QualType
-handleComplexFloatToComplexFloatConverstion(Sema &S, ExprResult &LHS,
-                                            ExprResult &RHS, QualType LHSType,
-                                            QualType RHSType,
-                                            bool IsCompAssign) {
-  int order = S.Context.getFloatingTypeOrder(LHSType, RHSType);
-
-  if (order < 0) {
-    // _Complex float -> _Complex double
-    if (!IsCompAssign)
-      LHS = S.ImpCastExprToType(LHS.get(), RHSType, CK_FloatingComplexCast);
-    return RHSType;
-  }
-  if (order > 0)
-    // _Complex float -> _Complex double
-    RHS = S.ImpCastExprToType(RHS.get(), LHSType, CK_FloatingComplexCast);
-  return LHSType;
-}
-
-/// \brief Converts otherExpr to complex float and promotes complexExpr if
-/// necessary.  Helper function of UsualArithmeticConversions()
-static QualType handleOtherComplexFloatConversion(Sema &S,
-                                                  ExprResult &ComplexExpr,
-                                                  ExprResult &OtherExpr,
-                                                  QualType ComplexTy,
-                                                  QualType OtherTy,
-                                                  bool ConvertComplexExpr,
-                                                  bool ConvertOtherExpr) {
-  int order = S.Context.getFloatingTypeOrder(ComplexTy, OtherTy);
-
-  // If just the complexExpr is complex, the otherExpr needs to be converted,
-  // and the complexExpr might need to be promoted.
-  if (order > 0) { // complexExpr is wider
-    // float -> _Complex double
-    if (ConvertOtherExpr) {
-      QualType fp = cast<ComplexType>(ComplexTy)->getElementType();
-      OtherExpr = S.ImpCastExprToType(OtherExpr.get(), fp, CK_FloatingCast);
-      OtherExpr = S.ImpCastExprToType(OtherExpr.get(), ComplexTy,
-                                      CK_FloatingRealToComplex);
-    }
-    return ComplexTy;
-  }
-
-  // otherTy is at least as wide.  Find its corresponding complex type.
-  QualType result = (order == 0 ? ComplexTy :
-                                  S.Context.getComplexType(OtherTy));
-
-  // double -> _Complex double
-  if (ConvertOtherExpr)
-    OtherExpr = S.ImpCastExprToType(OtherExpr.get(), result,
-                                    CK_FloatingRealToComplex);
-
-  // _Complex float -> _Complex double
-  if (ConvertComplexExpr && order < 0)
-    ComplexExpr = S.ImpCastExprToType(ComplexExpr.get(), result,
-                                      CK_FloatingComplexCast);
-
-  return result;
-}
-
 /// \brief Handle arithmetic conversion with complex types.  Helper function of
 /// UsualArithmeticConversions()
 static QualType handleComplexFloatConversion(Sema &S, ExprResult &LHS,
                                              ExprResult &RHS, QualType LHSType,
                                              QualType RHSType,
                                              bool IsCompAssign) {
   // if we have an integer operand, the result is the complex type.
   if (!handleIntegerToComplexFloatConversion(S, RHS, LHS, RHSType, LHSType,
static QualType handleComplexFloatConversion(Sema &S, ExprResult &LHS,
   // to what is done when combining two real floating-point operands.
   // The fun begins when size promotion occur across type domains.
   // From H&S 6.3.4: When one operand is complex and the other is a real
   // floating-point type, the less precise type is converted, within it's
   // real or complex domain, to the precision of the other type. For example,
   // when combining a "long double" with a "double _Complex", the
   // "double _Complex" is promoted to "long double _Complex".
 
-  bool LHSComplexFloat = LHSType->isComplexType();
-  bool RHSComplexFloat = RHSType->isComplexType();
+  // Compute the rank of the two types, regardless of whether they are complex.
+  int Order = S.Context.getFloatingTypeOrder(LHSType, RHSType);
 
-  // If both are complex, just cast to the more precise type.
-  if (LHSComplexFloat && RHSComplexFloat)
-    return handleComplexFloatToComplexFloatConverstion(S, LHS, RHS,
-                                                       LHSType, RHSType,
-                                                       IsCompAssign);
-
-  // If only one operand is complex, promote it if necessary and convert the
-  // other operand to complex.
-  if (LHSComplexFloat)
-    return handleOtherComplexFloatConversion(
-        S, LHS, RHS, LHSType, RHSType, /*convertComplexExpr*/!IsCompAssign,
-        /*convertOtherExpr*/ true);
-
-  assert(RHSComplexFloat);
-  return handleOtherComplexFloatConversion(
-      S, RHS, LHS, RHSType, LHSType, /*convertComplexExpr*/true,
-      /*convertOtherExpr*/ !IsCompAssign);
+  auto *LHSComplexType = dyn_cast<ComplexType>(LHSType);
+  auto *RHSComplexType = dyn_cast<ComplexType>(RHSType);
+  QualType LHSElementType =
+      LHSComplexType ? LHSComplexType->getElementType() : LHSType;
+  QualType RHSElementType =
+      RHSComplexType ? RHSComplexType->getElementType() : RHSType;
+
+  QualType ResultType = S.Context.getComplexType(LHSElementType);
+  if (Order < 0) {
+    // Promote the precision of the LHS if not an assignment.
+    ResultType = S.Context.getComplexType(RHSElementType);
+    if (!IsCompAssign) {
+      if (LHSComplexType)
+        LHS = S.ImpCastExprToType(LHS.get(), ResultType,
+                                  CK_FloatingComplexCast);
+      else
+        LHS = S.ImpCastExprToType(LHS.get(), RHSElementType, CK_FloatingCast);
+    }
+  } else if (Order > 0) {
+    // Promote the precision of the RHS.
+    if (RHSComplexType)
+      RHS = S.ImpCastExprToType(RHS.get(), ResultType, CK_FloatingComplexCast);
+    else
+      RHS = S.ImpCastExprToType(RHS.get(), LHSElementType, CK_FloatingCast);
+  }
+  return ResultType;
 }
 
 /// \brief Hande arithmetic conversion from integer to float.  Helper function
 /// of UsualArithmeticConversions()
 static QualType handleIntToFloatConversion(Sema &S, ExprResult &FloatExpr,
                                            ExprResult &IntExpr,
                                            QualType FloatTy, QualType IntTy,
                                            bool ConvertFloat, bool ConvertInt) {

test/CodeGen/complex-math.c

+// RUN: %clang_cc1 %s -O1 -emit-llvm -triple x86_64-unknown-unknown -o - | FileCheck %s --check-prefix=X86
+
+float _Complex add_float_rr(float a, float b) {
+  // X86-LABEL: @add_float_rr(
+  // X86: fadd
+  // X86-NOT: fadd
+  // X86: ret
+  return a + b;
+}
+float _Complex add_float_cr(float _Complex a, float b) {
+  // X86-LABEL: @add_float_cr(
+  // X86: fadd
+  // X86-NOT: fadd
+  // X86: ret
+  return a + b;
+}
+float _Complex add_float_rc(float a, float _Complex b) {
+  // X86-LABEL: @add_float_rc(
+  // X86: fadd
+  // X86-NOT: fadd
+  // X86: ret
+  return a + b;
+}
+float _Complex add_float_cc(float _Complex a, float _Complex b) {
+  // X86-LABEL: @add_float_cc(
+  // X86: fadd
+  // X86: fadd
+  // X86-NOT: fadd
+  // X86: ret
+  return a + b;
+}
+
+float _Complex sub_float_rr(float a, float b) {
+  // X86-LABEL: @sub_float_rr(
+  // X86: fsub
+  // X86-NOT: fsub
+  // X86: ret
+  return a - b;
+}
+float _Complex sub_float_cr(float _Complex a, float b) {
+  // X86-LABEL: @sub_float_cr(
+  // X86: fsub
+  // X86-NOT: fsub
+  // X86: ret
+  return a - b;
+}
+float _Complex sub_float_rc(float a, float _Complex b) {
+  // X86-LABEL: @sub_float_rc(
+  // X86: fsub
+  // X86: fsub float -0.{{0+}}e+00,
+  // X86-NOT: fsub
+  // X86: ret
+  return a - b;
+}
+float _Complex sub_float_cc(float _Complex a, float _Complex b) {
+  // X86-LABEL: @sub_float_cc(
+  // X86: fsub
+  // X86: fsub
+  // X86-NOT: fsub
+  // X86: ret
+  return a - b;
+}
+
+float _Complex mul_float_rr(float a, float b) {
+  // X86-LABEL: @mul_float_rr(
+  // X86: fmul
+  // X86-NOT: fmul
+  // X86: ret
+  return a * b;
+}
+float _Complex mul_float_cr(float _Complex a, float b) {
+  // X86-LABEL: @mul_float_cr(
+  // X86: fmul
+  // X86: fmul
+  // X86-NOT: fmul
+  // X86: ret
+  return a * b;
+}
+float _Complex mul_float_rc(float a, float _Complex b) {
+  // X86-LABEL: @mul_float_rc(
+  // X86: fmul
+  // X86: fmul
+  // X86-NOT: fmul
+  // X86: ret
+  return a * b;
+}
+float _Complex mul_float_cc(float _Complex a, float _Complex b) {
+  // X86-LABEL: @mul_float_cc(
+  // X86-NOT: fmul
+  // X86: call <2 x float> @__mulsc3(
+  // X86: ret
+  return a * b;
+}
+
+float _Complex div_float_rr(float a, float b) {
+  // X86-LABEL: @div_float_rr(
+  // X86: fdiv
+  // X86-NOT: fdiv
+  // X86: ret
+  return a / b;
+}
+float _Complex div_float_cr(float _Complex a, float b) {
+  // X86-LABEL: @div_float_cr(
+  // X86: fdiv
+  // X86: fdiv
+  // X86-NOT: fdiv
+  // X86: ret
+  return a / b;
+}
+float _Complex div_float_rc(float a, float _Complex b) {
+  // X86-LABEL: @div_float_rc(
+  // X86-NOT: fdiv
+  // X86: call <2 x float> @__divsc3(
+  // X86: ret
+  return a / b;
+}
+float _Complex div_float_cc(float _Complex a, float _Complex b) {
+  // X86-LABEL: @div_float_cc(
+  // X86-NOT: fdiv
+  // X86: call <2 x float> @__divsc3(
+  // X86: ret
+  return a / b;
+}
+
+double _Complex add_double_rr(double a, double b) {
+  // X86-LABEL: @add_double_rr(
+  // X86: fadd
+  // X86-NOT: fadd
+  // X86: ret
+  return a + b;
+}
+double _Complex add_double_cr(double _Complex a, double b) {
+  // X86-LABEL: @add_double_cr(
+  // X86: fadd
+  // X86-NOT: fadd
+  // X86: ret
+  return a + b;
+}
+double _Complex add_double_rc(double a, double _Complex b) {
+  // X86-LABEL: @add_double_rc(
+  // X86: fadd
+  // X86-NOT: fadd
+  // X86: ret
+  return a + b;
+}
+double _Complex add_double_cc(double _Complex a, double _Complex b) {
+  // X86-LABEL: @add_double_cc(
+  // X86: fadd
+  // X86: fadd
+  // X86-NOT: fadd
+  // X86: ret
+  return a + b;
+}
+
+double _Complex sub_double_rr(double a, double b) {
+  // X86-LABEL: @sub_double_rr(
+  // X86: fsub
+  // X86-NOT: fsub
+  // X86: ret
+  return a - b;
+}
+double _Complex sub_double_cr(double _Complex a, double b) {
+  // X86-LABEL: @sub_double_cr(
+  // X86: fsub
+  // X86-NOT: fsub
+  // X86: ret
+  return a - b;
+}
+double _Complex sub_double_rc(double a, double _Complex b) {
+  // X86-LABEL: @sub_double_rc(
+  // X86: fsub
+  // X86: fsub double -0.{{0+}}e+00,
+  // X86-NOT: fsub
+  // X86: ret
+  return a - b;
+}
+double _Complex sub_double_cc(double _Complex a, double _Complex b) {
+  // X86-LABEL: @sub_double_cc(
+  // X86: fsub
+  // X86: fsub
+  // X86-NOT: fsub
+  // X86: ret
+  return a - b;
+}
+
+double _Complex mul_double_rr(double a, double b) {
+  // X86-LABEL: @mul_double_rr(
+  // X86: fmul
+  // X86-NOT: fmul
+  // X86: ret
+  return a * b;
+}
+double _Complex mul_double_cr(double _Complex a, double b) {
+  // X86-LABEL: @mul_double_cr(
+  // X86: fmul
+  // X86: fmul
+  // X86-NOT: fmul
+  // X86: ret
+  return a * b;
+}
+double _Complex mul_double_rc(double a, double _Complex b) {
+  // X86-LABEL: @mul_double_rc(
+  // X86: fmul
+  // X86: fmul
+  // X86-NOT: fmul
+  // X86: ret
+  return a * b;
+}
+double _Complex mul_double_cc(double _Complex a, double _Complex b) {
+  // X86-LABEL: @mul_double_cc(
+  // X86-NOT: fmul
+  // X86: call { double, double } @__muldc3(
+  // X86: ret
+  return a * b;
+}
+
+double _Complex div_double_rr(double a, double b) {
+  // X86-LABEL: @div_double_rr(
+  // X86: fdiv
+  // X86-NOT: fdiv
+  // X86: ret
+  return a / b;
+}
+double _Complex div_double_cr(double _Complex a, double b) {
+  // X86-LABEL: @div_double_cr(
+  // X86: fdiv
+  // X86: fdiv
+  // X86-NOT: fdiv
+  // X86: ret
+  return a / b;
+}
+double _Complex div_double_rc(double a, double _Complex b) {
+  // X86-LABEL: @div_double_rc(
+  // X86-NOT: fdiv
+  // X86: call { double, double } @__divdc3(
+  // X86: ret
+  return a / b;
+}
+double _Complex div_double_cc(double _Complex a, double _Complex b) {
+  // X86-LABEL: @div_double_cc(
+  // X86-NOT: fdiv
+  // X86: call { double, double } @__divdc3(
+  // X86: ret
+  return a / b;
+}
+
+long double _Complex add_long_double_rr(long double a, long double b) {
+  // X86-LABEL: @add_long_double_rr(
+  // X86: fadd
+  // X86-NOT: fadd
+  // X86: ret
+  return a + b;
+}
+long double _Complex add_long_double_cr(long double _Complex a, long double b) {
+  // X86-LABEL: @add_long_double_cr(
+  // X86: fadd
+  // X86-NOT: fadd
+  // X86: ret
+  return a + b;
+}
+long double _Complex add_long_double_rc(long double a, long double _Complex b) {
+  // X86-LABEL: @add_long_double_rc(
+  // X86: fadd
+  // X86-NOT: fadd
+  // X86: ret
+  return a + b;
+}
+long double _Complex add_long_double_cc(long double _Complex a, long double _Complex b) {
+  // X86-LABEL: @add_long_double_cc(
+  // X86: fadd
+  // X86: fadd
+  // X86-NOT: fadd
+  // X86: ret
+  return a + b;
+}
+
+long double _Complex sub_long_double_rr(long double a, long double b) {
+  // X86-LABEL: @sub_long_double_rr(
+  // X86: fsub
+  // X86-NOT: fsub
+  // X86: ret
+  return a - b;
+}
+long double _Complex sub_long_double_cr(long double _Complex a, long double b) {
+  // X86-LABEL: @sub_long_double_cr(
+  // X86: fsub
+  // X86-NOT: fsub
+  // X86: ret
+  return a - b;
+}
+long double _Complex sub_long_double_rc(long double a, long double _Complex b) {
+  // X86-LABEL: @sub_long_double_rc(
+  // X86: fsub
+  // X86: fsub x86_fp80 0xK8{{0+}},
+  // X86-NOT: fsub
+  // X86: ret
+  return a - b;
+}
+long double _Complex sub_long_double_cc(long double _Complex a, long double _Complex b) {
+  // X86-LABEL: @sub_long_double_cc(
+  // X86: fsub
+  // X86: fsub
+  // X86-NOT: fsub
+  // X86: ret
+  return a - b;
+}
+
+long double _Complex mul_long_double_rr(long double a, long double b) {
+  // X86-LABEL: @mul_long_double_rr(
+  // X86: fmul
+  // X86-NOT: fmul
+  // X86: ret
+  return a * b;
+}
+long double _Complex mul_long_double_cr(long double _Complex a, long double b) {
+  // X86-LABEL: @mul_long_double_cr(
+  // X86: fmul
+  // X86: fmul
+  // X86-NOT: fmul
+  // X86: ret
+  return a * b;
+}
+long double _Complex mul_long_double_rc(long double a, long double _Complex b) {
+  // X86-LABEL: @mul_long_double_rc(
+  // X86: fmul
+  // X86: fmul
+  // X86-NOT: fmul
+  // X86: ret
+  return a * b;
+}
+long double _Complex mul_long_double_cc(long double _Complex a, long double _Complex b) {
+  // X86-LABEL: @mul_long_double_cc(
+  // X86-NOT: fmul
+  // X86: call { x86_fp80, x86_fp80 } @__mulxc3(
+  // X86: ret
+  return a * b;
+}
+
+long double _Complex div_long_double_rr(long double a, long double b) {
+  // X86-LABEL: @div_long_double_rr(
+  // X86: fdiv
+  // X86-NOT: fdiv
+  // X86: ret
+  return a / b;
+}
+long double _Complex div_long_double_cr(long double _Complex a, long double b) {
+  // X86-LABEL: @div_long_double_cr(
+  // X86: fdiv
+  // X86: fdiv
+  // X86-NOT: fdiv
+  // X86: ret
+  return a / b;
+}
+long double _Complex div_long_double_rc(long double a, long double _Complex b) {
+  // X86-LABEL: @div_long_double_rc(
+  // X86-NOT: fdiv
+  // X86: call { x86_fp80, x86_fp80 } @__divxc3(
+  // X86: ret
+  return a / b;
+}
+long double _Complex div_long_double_cc(long double _Complex a, long double _Complex b) {
+  // X86-LABEL: @div_long_double_cc(
+  // X86-NOT: fdiv
+  // X86: call { x86_fp80, x86_fp80 } @__divxc3(
+  // X86: ret
+  return a / b;
+}

test/SemaCXX/complex-folding.cpp

+// RUN: %clang_cc1 %s -std=c++1z -fsyntax-only -verify
+//
+// Test the constant folding of builtin complex numbers.
+
+static_assert((0.0 + 0.0j) == (0.0 + 0.0j));
+
+// Walk around the complex plane stepping between angular differences and
+// equality.
+static_assert((1.0 + 0.0j) == (0.0 + 0.0j)); // expected-error {{static_assert}}
+static_assert((1.0 + 0.0j) == (1.0 + 0.0j));
+static_assert((1.0 + 1.0j) == (1.0 + 0.0j)); // expected-error {{static_assert}}
+static_assert((1.0 + 1.0j) == (1.0 + 1.0j));
+static_assert((0.0 + 1.0j) == (1.0 + 1.0j)); // expected-error {{static_assert}}
+static_assert((0.0 + 1.0j) == (0.0 + 1.0j));
+static_assert((-1.0 + 1.0j) == (0.0 + 1.0j)); // expected-error {{static_assert}}
+static_assert((-1.0 + 1.0j) == (-1.0 + 1.0j));
+static_assert((-1.0 + 0.0j) == (-1.0 + 1.0j)); // expected-error {{static_assert}}
+static_assert((-1.0 + 0.0j) == (-1.0 + 0.0j));
+static_assert((-1.0 - 1.0j) == (-1.0 + 0.0j)); // expected-error {{static_assert}}
+static_assert((-1.0 - 1.0j) == (-1.0 - 1.0j));
+static_assert((0.0 - 1.0j) == (-1.0 - 1.0j)); // expected-error {{static_assert}}
+static_assert((0.0 - 1.0j) == (0.0 - 1.0j));
+static_assert((1.0 - 1.0j) == (0.0 - 1.0j)); // expected-error {{static_assert}}
+static_assert((1.0 - 1.0j) == (1.0 - 1.0j));
+
+// Test basic mathematical folding of both complex and real operands.
+static_assert(((1.0 + 0.5j) + (0.25 - 0.75j)) == (1.25 - 0.25j));
+static_assert(((1.0 + 0.5j) + 0.25) == (1.25 + 0.5j));
+static_assert((1.0 + (0.25 - 0.75j)) == (1.25 - 0.75j));
+
+static_assert(((1.0 + 0.5j) - (0.25 - 0.75j)) == (0.75 + 1.25j));
+static_assert(((1.0 + 0.5j) - 0.25) == (0.75 + 0.5j));
+static_assert((1.0 - (0.25 - 0.75j)) == (0.75 + 0.75j));
+
+static_assert(((1.25 + 0.5j) * (0.25 - 0.75j)) == (0.6875 - 0.8125j));
+static_assert(((1.25 + 0.5j) * 0.25) == (0.3125 + 0.125j));
+static_assert((1.25 * (0.25 - 0.75j)) == (0.3125 - 0.9375j));
+
+static_assert(((1.25 + 0.5j) / (0.25 - 0.75j)) == (-0.1 + 1.7j));
+static_assert(((1.25 + 0.5j) / 0.25) == (5.0 + 2.0j));
+static_assert((1.25 / (0.25 - 0.75j)) == (0.5 + 1.5j));
+
+// Test that infinities are preserved, don't turn into NaNs, and do form zeros
+// when the divisor.
+static_assert(__builtin_isinf_sign(__real__ ((__builtin_inf() + 1.0j) * 1.0)) == 1);
+static_assert(__builtin_isinf_sign(__imag__ ((1.0 + __builtin_inf() * 1.0j) * 1.0)) == 1);
+static_assert(__builtin_isinf_sign(__real__ (1.0 * (__builtin_inf() + 1.0j))) == 1);
+static_assert(__builtin_isinf_sign(__imag__ (1.0 * (1.0 + __builtin_inf() * 1.0j))) == 1);
+
+static_assert(__builtin_isinf_sign(__real__ ((__builtin_inf() + 1.0j) * (1.0 + 1.0j))) == 1);
+static_assert(__builtin_isinf_sign(__real__ ((1.0 + 1.0j) * (__builtin_inf() + 1.0j))) == 1);
+static_assert(__builtin_isinf_sign(__real__ ((__builtin_inf() + 1.0j) * (__builtin_inf() + 1.0j))) == 1);
+
+static_assert(__builtin_isinf_sign(__real__ ((1.0 + __builtin_inf() * 1.0j) * (1.0 + 1.0j))) == -1);
+static_assert(__builtin_isinf_sign(__imag__ ((1.0 + __builtin_inf() * 1.0j) * (1.0 + 1.0j))) == 1);
+static_assert(__builtin_isinf_sign(__real__ ((1.0 + 1.0j) * (1.0 + __builtin_inf() * 1.0j))) == -1);
+static_assert(__builtin_isinf_sign(__imag__ ((1.0 + 1.0j) * (1.0 + __builtin_inf() * 1.0j))) == 1);
+
+static_assert(__builtin_isinf_sign(__real__ ((1.0 + __builtin_inf() * 1.0j) * (1.0 + __builtin_inf() * 1.0j))) == -1);
+static_assert(__builtin_isinf_sign(__real__ ((__builtin_inf() + __builtin_inf() * 1.0j) * (__builtin_inf() + __builtin_inf() * 1.0j))) == -1);
+
+static_assert(__builtin_isinf_sign(__real__ ((__builtin_inf() + 1.0j) / (1.0 + 1.0j))) == 1);
+static_assert(__builtin_isinf_sign(__imag__ (1.0 + (__builtin_inf() * 1.0j) / (1.0 + 1.0j))) == 1);
+static_assert(__builtin_isinf_sign(__imag__ ((__builtin_inf() + __builtin_inf() * 1.0j) / (1.0 + 1.0j))) == 1);
+static_assert(__builtin_isinf_sign(__real__ ((__builtin_inf() + 1.0j) / 1.0)) == 1);
+static_assert(__builtin_isinf_sign(__imag__ (1.0 + (__builtin_inf() * 1.0j) / 1.0)) == 1);
+static_assert(__builtin_isinf_sign(__imag__ ((__builtin_inf() + __builtin_inf() * 1.0j) / 1.0)) == 1);
+
+static_assert(((1.0 + 1.0j) / (__builtin_inf() + 1.0j)) == (0.0 + 0.0j));
+static_assert(((1.0 + 1.0j) / (1.0 + __builtin_inf() * 1.0j)) == (0.0 + 0.0j));
+static_assert(((1.0 + 1.0j) / (__builtin_inf() + __builtin_inf() * 1.0j)) == (0.0 + 0.0j));
+static_assert(((1.0 + 1.0j) / __builtin_inf()) == (0.0 + 0.0j));
+
+static_assert(__builtin_isinf_sign(__real__ ((1.0 + 1.0j) / (0.0 + 0.0j))) == 1);
+static_assert(__builtin_isinf_sign(__real__ ((1.0 + 1.0j) / 0.0)) == 1);
+
+static_assert(__builtin_isinf_sign(__real__ ((__builtin_inf() + 1.0j) / (0.0 + 0.0j))) == 1);
+static_assert(__builtin_isinf_sign(__imag__ ((1.0 + __builtin_inf() * 1.0j) / (0.0 + 0.0j))) == 1);
+static_assert(__builtin_isinf_sign(__imag__ ((__builtin_inf() + __builtin_inf() * 1.0j) / (0.0 + 0.0j))) == 1);
+static_assert(__builtin_isinf_sign(__real__ ((__builtin_inf() + 1.0j) / 0.0)) == 1);
+static_assert(__builtin_isinf_sign(__imag__ ((1.0 + __builtin_inf() * 1.0j) / 0.0)) == 1);
+static_assert(__builtin_isinf_sign(__imag__ ((__builtin_inf() + __builtin_inf() * 1.0j) / 0.0)) == 1);