@@ -275,6 +275,10 @@ class ComplexExprEmitter
275275 ComplexPairTy EmitBinSub (const BinOpInfo &Op);
276276 ComplexPairTy EmitBinMul (const BinOpInfo &Op);
277277 ComplexPairTy EmitBinDiv (const BinOpInfo &Op);
278+ ComplexPairTy EmitAlgebraicDiv (llvm::Value *A, llvm::Value *B, llvm::Value *C,
279+ llvm::Value *D);
280+ ComplexPairTy EmitRangeReductionDiv (llvm::Value *A, llvm::Value *B,
281+ llvm::Value *C, llvm::Value *D);
278282
279283 ComplexPairTy EmitComplexBinOpLibCall (StringRef LibCallName,
280284 const BinOpInfo &Op);
@@ -781,6 +785,10 @@ ComplexPairTy ComplexExprEmitter::EmitBinMul(const BinOpInfo &Op) {
781785 ResR = Builder.CreateFSub (AC, BD, " mul_r"
782786 ResI = Builder.CreateFAdd (AD, BC, " mul_i"
783787
788+ if (Op.FPFeatures .getComplexRange () == LangOptions::CX_Limited ||
789+ Op.FPFeatures .getComplexRange () == LangOptions::CX_Fortran)
790+ return ComplexPairTy (ResR, ResI);
791+
784792 // Emit the test for the real part becoming NaN and create a branch to
785793 // handle it. We test for NaN by comparing the number to itself.
786794 Value *IsRNaN = Builder.CreateFCmpUNO (ResR, ResR, " isnan_cmp"
@@ -846,23 +854,139 @@ ComplexPairTy ComplexExprEmitter::EmitBinMul(const BinOpInfo &Op) {
846854 return ComplexPairTy (ResR, ResI);
847855}
848856
857+ ComplexPairTy ComplexExprEmitter::EmitAlgebraicDiv (llvm::Value *LHSr,
858+ llvm::Value *LHSi,
859+ llvm::Value *RHSr,
860+ llvm::Value *RHSi) {
861+ // (a+ib) / (c+id) = ((ac+bd)/(cc+dd)) + i((bc-ad)/(cc+dd))
862+ llvm::Value *DSTr, *DSTi;
863+
864+ llvm::Value *AC = Builder.CreateFMul (LHSr, RHSr); // a*c
865+ llvm::Value *BD = Builder.CreateFMul (LHSi, RHSi); // b*d
866+ llvm::Value *ACpBD = Builder.CreateFAdd (AC, BD); // ac+bd
867+
868+ llvm::Value *CC = Builder.CreateFMul (RHSr, RHSr); // c*c
869+ llvm::Value *DD = Builder.CreateFMul (RHSi, RHSi); // d*d
870+ llvm::Value *CCpDD = Builder.CreateFAdd (CC, DD); // cc+dd
871+
872+ llvm::Value *BC = Builder.CreateFMul (LHSi, RHSr); // b*c
873+ llvm::Value *AD = Builder.CreateFMul (LHSr, RHSi); // a*d
874+ llvm::Value *BCmAD = Builder.CreateFSub (BC, AD); // bc-ad
875+
876+ DSTr = Builder.CreateFDiv (ACpBD, CCpDD);
877+ DSTi = Builder.CreateFDiv (BCmAD, CCpDD);
878+ return ComplexPairTy (DSTr, DSTi);
879+ }
880+
881+ // EmitFAbs - Emit a call to @llvm.fabs.
882+ static llvm::Value *EmitllvmFAbs (CodeGenFunction &CGF, llvm::Value *Value) {
883+ llvm::Function *Func =
884+ CGF.CGM .getIntrinsic (llvm::Intrinsic::fabs, Value->getType ());
885+ llvm::Value *Call = CGF.Builder .CreateCall (Func, Value);
886+ return Call;
887+ }
888+
889+ // EmitRangeReductionDiv - Implements Smith's algorithm for complex division.
890+ // SMITH, R. L. Algorithm 116: Complex division. Commun. ACM 5, 8 (1962).
891+ ComplexPairTy ComplexExprEmitter::EmitRangeReductionDiv (llvm::Value *LHSr,
892+ llvm::Value *LHSi,
893+ llvm::Value *RHSr,
894+ llvm::Value *RHSi) {
895+ // (a + ib) / (c + id) = (e + if)
896+ llvm::Value *FAbsRHSr = EmitllvmFAbs (CGF, RHSr); // |c|
897+ llvm::Value *FAbsRHSi = EmitllvmFAbs (CGF, RHSi); // |d|
898+ // |c| >= |d|
899+ llvm::Value *IsR = Builder.CreateFCmpUGT (FAbsRHSr, FAbsRHSi, " abs_cmp"
900+
901+ llvm::BasicBlock *TrueBB =
902+ CGF.createBasicBlock (" abs_rhsr_greater_or_equal_abs_rhsi"
903+ llvm::BasicBlock *FalseBB =
904+ CGF.createBasicBlock (" abs_rhsr_less_than_abs_rhsi"
905+ llvm::BasicBlock *ContBB = CGF.createBasicBlock (" complex_div"
906+ Builder.CreateCondBr (IsR, TrueBB, FalseBB);
907+
908+ CGF.EmitBlock (TrueBB);
909+ // abs(c) >= abs(d)
910+ // r = d/c
911+ // tmp = c + rd
912+ // e = (a + br)/tmp
913+ // f = (b - ar)/tmp
914+ llvm::Value *DdC = Builder.CreateFDiv (RHSi, RHSr); // r=d/c
915+
916+ llvm::Value *RD = Builder.CreateFMul (DdC, RHSi); // rd
917+ llvm::Value *CpRD = Builder.CreateFAdd (RHSr, RD); // tmp=c+rd
918+
919+ llvm::Value *T3 = Builder.CreateFMul (LHSi, DdC); // br
920+ llvm::Value *T4 = Builder.CreateFAdd (LHSr, T3); // a+br
921+ llvm::Value *DSTTr = Builder.CreateFDiv (T4, CpRD); // (a+br)/tmp
922+
923+ llvm::Value *T5 = Builder.CreateFMul (LHSr, DdC); // ar
924+ llvm::Value *T6 = Builder.CreateFSub (LHSi, T5); // b-ar
925+ llvm::Value *DSTTi = Builder.CreateFDiv (T6, CpRD); // (b-ar)/tmp
926+ Builder.CreateBr (ContBB);
927+
928+ CGF.EmitBlock (FalseBB);
929+ // abs(c) < abs(d)
930+ // r = c/d
931+ // tmp = d + rc
932+ // e = (ar + b)/tmp
933+ // f = (br - a)/tmp
934+ llvm::Value *CdD = Builder.CreateFDiv (RHSr, RHSi); // r=c/d
935+
936+ llvm::Value *RC = Builder.CreateFMul (CdD, RHSr); // rc
937+ llvm::Value *DpRC = Builder.CreateFAdd (RHSi, RC); // tmp=d+rc
938+
939+ llvm::Value *T7 = Builder.CreateFMul (LHSr, RC); // ar
940+ llvm::Value *T8 = Builder.CreateFAdd (T7, LHSi); // ar+b
941+ llvm::Value *DSTFr = Builder.CreateFDiv (T8, DpRC); // (ar+b)/tmp
942+
943+ llvm::Value *T9 = Builder.CreateFMul (LHSi, CdD); // br
944+ llvm::Value *T10 = Builder.CreateFSub (T9, LHSr); // br-a
945+ llvm::Value *DSTFi = Builder.CreateFDiv (T10, DpRC); // (br-a)/tmp
946+ Builder.CreateBr (ContBB);
947+
948+ // Phi together the computation paths.
949+ CGF.EmitBlock (ContBB);
950+ llvm::PHINode *VALr = Builder.CreatePHI (DSTTr->getType (), 2 );
951+ VALr->addIncoming (DSTTr, TrueBB);
952+ VALr->addIncoming (DSTFr, FalseBB);
953+ llvm::PHINode *VALi = Builder.CreatePHI (DSTTi->getType (), 2 );
954+ VALi->addIncoming (DSTTi, TrueBB);
955+ VALi->addIncoming (DSTFi, FalseBB);
956+ return ComplexPairTy (VALr, VALi);
957+ }
958+
849959// See C11 Annex G.5.1 for the semantics of multiplicative operators on complex
850960// typed values.
851961ComplexPairTy ComplexExprEmitter::EmitBinDiv (const BinOpInfo &Op) {
852962 llvm::Value *LHSr = Op.LHS .first , *LHSi = Op.LHS .second ;
853963 llvm::Value *RHSr = Op.RHS .first , *RHSi = Op.RHS .second ;
854-
855964 llvm::Value *DSTr, *DSTi;
856965 if (LHSr->getType ()->isFloatingPointTy ()) {
857- // If we have a complex operand on the RHS and FastMath is not allowed, we
858- // delegate to a libcall to handle all of the complexities and minimize
859- // underflow/overflow cases. When FastMath is allowed we construct the
860- // divide inline using the same algorithm as for integer operands.
861- //
862- // FIXME: We would be able to avoid the libcall in many places if we
863- // supported imaginary types in addition to complex types.
864966 CodeGenFunction::CGFPOptionsRAII FPOptsRAII (CGF, Op.FPFeatures );
865- if (RHSi && !CGF.getLangOpts ().FastMath ) {
967+ if (!RHSi) {
968+ assert (LHSi && " Can have at most one non-complex operand!"
969+
970+ DSTr = Builder.CreateFDiv (LHSr, RHSr);
971+ DSTi = Builder.CreateFDiv (LHSi, RHSr);
972+ return ComplexPairTy (DSTr, DSTi);
973+ }
974+ llvm::Value *OrigLHSi = LHSi;
975+ if (!LHSi)
976+ LHSi = llvm::Constant::getNullValue (RHSi->getType ());
977+ if (Op.FPFeatures .getComplexRange () == LangOptions::CX_Fortran)
978+ return EmitRangeReductionDiv (LHSr, LHSi, RHSr, RHSi);
979+ else if (Op.FPFeatures .getComplexRange () == LangOptions::CX_Limited)
980+ return EmitAlgebraicDiv (LHSr, LHSi, RHSr, RHSi);
981+ else if (!CGF.getLangOpts ().FastMath ) {
982+ LHSi = OrigLHSi;
983+ // If we have a complex operand on the RHS and FastMath is not allowed, we
984+ // delegate to a libcall to handle all of the complexities and minimize
985+ // underflow/overflow cases. When FastMath is allowed we construct the
986+ // divide inline using the same algorithm as for integer operands.
987+ //
988+ // FIXME: We would be able to avoid the libcall in many places if we
989+ // supported imaginary types in addition to complex types.
866990 BinOpInfo LibCallOp = Op;
867991 // If LHS was a real, supply a null imaginary part.
868992 if (!LHSi)
@@ -884,30 +1008,8 @@ ComplexPairTy ComplexExprEmitter::EmitBinDiv(const BinOpInfo &Op) {
8841008 case llvm::Type::FP128TyID:
8851009 return EmitComplexBinOpLibCall (" __divtc3"
8861010 }
887- } else if (RHSi) {
888- if (!LHSi)
889- LHSi = llvm::Constant::getNullValue (RHSi->getType ());
890-
891- // (a+ib) / (c+id) = ((ac+bd)/(cc+dd)) + i((bc-ad)/(cc+dd))
892- llvm::Value *AC = Builder.CreateFMul (LHSr, RHSr); // a*c
893- llvm::Value *BD = Builder.CreateFMul (LHSi, RHSi); // b*d
894- llvm::Value *ACpBD = Builder.CreateFAdd (AC, BD); // ac+bd
895-
896- llvm::Value *CC = Builder.CreateFMul (RHSr, RHSr); // c*c
897- llvm::Value *DD = Builder.CreateFMul (RHSi, RHSi); // d*d
898- llvm::Value *CCpDD = Builder.CreateFAdd (CC, DD); // cc+dd
899-
900- llvm::Value *BC = Builder.CreateFMul (LHSi, RHSr); // b*c
901- llvm::Value *AD = Builder.CreateFMul (LHSr, RHSi); // a*d
902- llvm::Value *BCmAD = Builder.CreateFSub (BC, AD); // bc-ad
903-
904- DSTr = Builder.CreateFDiv (ACpBD, CCpDD);
905- DSTi = Builder.CreateFDiv (BCmAD, CCpDD);
9061011 } else {
907- assert (LHSi && " Can have at most one non-complex operand!"
908-
909- DSTr = Builder.CreateFDiv (LHSr, RHSr);
910- DSTi = Builder.CreateFDiv (LHSi, RHSr);
1012+ return EmitAlgebraicDiv (LHSr, LHSi, RHSr, RHSi);
9111013 }
9121014 } else {
9131015 assert (Op.LHS .second && Op.RHS .second &&
0 commit comments