static bool valid(index_type i, index_type j, matrix_type *output, std::complex<float> const &value)
 {
   vsip_cscalar_f c = vsip_cmget_f(output, i, j);
   return (equal(c.r, value.real()) && equal(c.i, value.imag()));
 }

/*!
 * \copydoc is_zero
 */
inline bool is_zero(std::complex<double> a) {
    return a.real() == 0.0 && a.imag() == 0.0;
}

 void set_iq_balance(const std::complex<double> &cor){
     _iface->poke32(REG_TX_FE_MAG_CORRECTION, fs_to_bits(cor.real(), 18));
     _iface->poke32(REG_TX_FE_PHASE_CORRECTION, fs_to_bits(cor.imag(), 18));
 }

 static inline U apply(std::complex<T> const & arg) {
     return static_cast<U>(arg.real());
 }

 arma_inline static T get(const std::complex<T>& val)
   {
   return -val.real();
   }

//=================Strassen===============================
QDebug  operator <<(QDebug qd,std::complex<double> a)
{
    qd<<"("<<QString::number(a.real())/*+","+QString::number(a.imag())*/<<")";
    return qd;
}

template <typename T> void readDouble(QDataStream &in, std::complex<T> &el)
{ double v; in >> v; el.real(v); el.imag(0); }

template <typename T> void readComplex(QDataStream &in, std::complex<T> &el)
{ double re, im; in >> re >> im; el.real(re); el.imag(im); }

double abs_sum_value( std::complex< double > const& f ) {
  using namespace std ;
  return abs(f.real()) + abs(f.imag()) ;
}

template <typename T> void dumpElement(QDataStream &out, const std::complex<T> &el)
{ out << (double)el.real() << (double)el.imag(); }

 inline static const  T  tmp_real(const std::complex<T>& X) { return X.real(); }

 inline void mumps_assign_Scalar(ZMUMPS_COMPLEX & a, std::complex<double> b)
 {
   a.r = b.real();
   a.i = b.imag();
 }

std::complex<T> atanh(const std::complex<T>& z)
{
   //
   // References:
   //
   // Eric W. Weisstein. "Inverse Hyperbolic Tangent." 
   // From MathWorld--A Wolfram Web Resource. 
   // http://mathworld.wolfram.com/InverseHyperbolicTangent.html
   //
   // Also: The Wolfram Functions Site,
   // http://functions.wolfram.com/ElementaryFunctions/ArcTanh/
   //
   // Also "Abramowitz and Stegun. Handbook of Mathematical Functions."
   // at : http://jove.prohosting.com/~skripty/toc.htm
   //
   
   static const T half_pi = static_cast<T>(1.57079632679489661923132169163975144L);
   static const T pi = static_cast<T>(3.141592653589793238462643383279502884197L);
   static const T one = static_cast<T>(1.0L);
   static const T two = static_cast<T>(2.0L);
   static const T four = static_cast<T>(4.0L);
   static const T zero = static_cast<T>(0);
   static const T a_crossover = static_cast<T>(0.3L);

   T x = std::fabs(z.real());
   T y = std::fabs(z.imag());

   T real, imag;  // our results

   T safe_upper = detail::safe_max(two);
   T safe_lower = detail::safe_min(static_cast<T>(2));

   //
   // Begin by handling the special cases specified in C99:
   //
   if(detail::test_is_nan(x))
   {
      if(detail::test_is_nan(y))
         return std::complex<T>(x, x);
      else if(std::numeric_limits<T>::has_infinity && (y == std::numeric_limits<T>::infinity()))
         return std::complex<T>(0, ((z.imag() < 0) ? -half_pi : half_pi));
      else
         return std::complex<T>(x, x);
   }
   else if(detail::test_is_nan(y))
   {
      if(x == 0)
         return std::complex<T>(x, y);
      if(std::numeric_limits<T>::has_infinity && (x == std::numeric_limits<T>::infinity()))
         return std::complex<T>(0, y);
      else
         return std::complex<T>(y, y);
   }
   else if((x > safe_lower) && (x < safe_upper) && (y > safe_lower) && (y < safe_upper))
   {

      T xx = x*x;
      T yy = y*y;
      T x2 = x * two;

      ///
      // The real part is given by:
      // 
      // real(atanh(z)) == log((1 + x^2 + y^2 + 2x) / (1 + x^2 + y^2 - 2x))
      // 
      // However, when x is either large (x > 1/E) or very small
      // (x < E) then this effectively simplifies
      // to log(1), leading to wildly inaccurate results.  
      // By dividing the above (top and bottom) by (1 + x^2 + y^2) we get:
      //
      // real(atanh(z)) == log((1 + (2x / (1 + x^2 + y^2))) / (1 - (-2x / (1 + x^2 + y^2))))
      //
      // which is much more sensitive to the value of x, when x is not near 1
      // (remember we can compute log(1+x) for small x very accurately).
      //
      // The cross-over from one method to the other has to be determined
      // experimentally, the value used below appears correct to within a 
      // factor of 2 (and there are larger errors from other parts
      // of the input domain anyway).
      //
      T alpha = two*x / (one + xx + yy);
      if(alpha < a_crossover)
      {
         real = boost::math::log1p(alpha) - boost::math::log1p(-alpha);
      }
      else
      {
         T xm1 = x - one;
         real = boost::math::log1p(x2 + xx + yy) - std::log(xm1*xm1 + yy);
      }
      real /= four;
      if(z.real() < 0)
         real = -real;

      imag = std::atan2((y * two), (one - xx - yy));
      imag /= two;
      if(z.imag() < 0)
         imag = -imag;
   }
   else
//.........这里部分代码省略.........

Expr::Expr(const std::complex<double>& c)
	: Playa::Handle<ExprBase>(new ComplexExpr(new ConstantExpr(c.real()),
      new ConstantExpr(c.imag())))
{}

void LTransform::set(
    PQIndex pq1, PQIndex pq2, std::complex<double> Cpq1pq2, std::complex<double> Cqp1pq2)
{
    assert(pq1.pqValid() && !pq1.pastOrder(_orderOut));
    assert(pq2.pqValid() && !pq2.pastOrder(_orderIn));

    take_ownership();
    const double RoundoffTolerance=1.e-15;
    std::complex<double> Cpq1qp2;

    if (pq2.needsConjugation()) {
        pq2 = pq2.swapPQ();
        std::complex<double> tmp=conj(Cqp1pq2);
        Cqp1pq2 = conj(Cpq1pq2);
        Cpq1pq2 = tmp;
    }
    if (pq1.needsConjugation()) {
        pq1 = pq1.swapPQ();
        std::complex<double> tmp=Cqp1pq2;
        Cqp1pq2 = Cpq1pq2;
        Cpq1pq2 = tmp;
    }

    int rIndex1 = pq1.rIndex();
    int rIndex2 = pq2.rIndex();
    int iIndex1 = rIndex1+1;
    int iIndex2 = rIndex2+1;

    if (pq1.isReal()) {
        if (Cpq1pq2!=Cqp1pq2) {
            FormatAndThrow<>()
                    << "Invalid LTransform elements for p1=q1, " << Cpq1pq2
                    << " != " << Cqp1pq2;
        }
        (*_m)(rIndex1,rIndex2) = Cpq1pq2.real() * (pq2.isReal()? 1. : 2.);
        if (pq2.isReal()) {
            if (std::abs(Cpq1pq2.imag()) > RoundoffTolerance) {
                FormatAndThrow<>()
                        << "Nonzero imaginary LTransform elements for p1=q1, p2=q2: "
                        << Cpq1pq2;
            }
        } else {
            (*_m)(rIndex1,iIndex2) = -2.*Cpq1pq2.imag();
        }
        return;
    } else if (pq2.isReal()) {
        // Here we know p1!=q1:
        if (norm(Cpq1pq2-conj(Cqp1pq2))>RoundoffTolerance) {
            FormatAndThrow<>()
                    << "Inputs to LTransform.set are not conjugate for p2=q2: "
                    << Cpq1pq2 << " vs " << Cqp1pq2 ;
        }
        (*_m)(rIndex1, rIndex2) = Cpq1pq2.real();
        (*_m)(iIndex1, rIndex2) = Cpq1pq2.imag();
    } else {
        // Neither pq is real:
        std::complex<double> z=Cpq1pq2 + Cqp1pq2;
        (*_m)(rIndex1, rIndex2) = z.real();
        (*_m)(rIndex1, iIndex2) = -z.imag();
        z=Cpq1pq2 - Cqp1pq2;
        (*_m)(iIndex1, rIndex2) = z.imag();
        (*_m)(iIndex1, iIndex2) = z.real();
    }
}

void microfacetNoExpFourierSeries(Float mu_o, Float mu_i, std::complex<Float> eta_,
                                  Float alpha, size_t n, Float phiMax,
                                  std::vector<Float> &result) {

    bool reflect = -mu_i * mu_o > 0;

    Float sinMu2 = math::safe_sqrt((1 - mu_i * mu_i) * (1 - mu_o * mu_o)),
          phiCritical = 0.0f;

    bool conductor = (eta_.imag() != 0.0f);
    std::complex<Float> eta =
        (-mu_i > 0 || conductor) ? eta_ : std::complex<Float>(1) / eta_;

    if (reflect) {
        if (!conductor)
            phiCritical = math::safe_acos((2*eta.real()*eta.real()-mu_i*mu_o-1)/sinMu2);
    } else if (!reflect) {
        if (conductor)
            throw std::runtime_error("lowfreqFourierSeries(): encountered refraction case for a conductor");
        Float etaDenser = (eta.real() > 1 ? eta.real() : 1 / eta.real());
        phiCritical = math::safe_acos((1 - etaDenser * mu_i * mu_o) /
                                      (etaDenser * sinMu2));
    }

    if (!conductor && phiCritical > math::Epsilon &&
        phiCritical < math::Pi - math::Epsilon &&
        phiCritical < phiMax - math::Epsilon) {
        /* Uh oh, some high frequency content leaked in the generally low frequency part.
           Increase the number of coefficients so that we can capture it. Fortunately, this
           happens very rarely. */
        n = std::max(n, (size_t) 100);
    }

    VectorX coeffs(n);
    coeffs.setZero();
    std::function<Float(Float)> integrand = std::bind(
        &microfacetNoExp, mu_o, mu_i, eta_, alpha, std::placeholders::_1);

    const int nEvals = 200;
    if (reflect) {
        if (phiCritical > math::Epsilon && phiCritical < phiMax-math::Epsilon) {
            filonIntegrate(integrand, coeffs.data(), n, nEvals, 0, phiCritical);
            filonIntegrate(integrand, coeffs.data(), n, nEvals, phiCritical, phiMax);
        } else {
            filonIntegrate(integrand, coeffs.data(), n, nEvals, 0, phiMax);
        }
    } else {
        filonIntegrate(integrand, coeffs.data(), n, nEvals, 0,
                       std::min(phiCritical, phiMax));
    }

    if (phiMax < math::Pi - math::Epsilon) {
        /* Precompute some sines and cosines */
        VectorX cosPhi(n), sinPhi(n);
        for (int i=0; i<n; ++i) {
            sinPhi[i] = std::sin(i*phiMax);
            cosPhi[i] = std::cos(i*phiMax);
        }

        /* The fit only occurs on a subset [0, phiMax], where the Fourier
           Fourier basis functions are not orthogonal anymore! The following
           then does a change of basis to proper Fourier coefficients. */
        MatrixX A(n, n);

        for (int i=0; i<n; ++i) {
            for (int j=0; j<=i; ++j) {
                if (i != j) {
                    A(i, j) = A(j, i) = (i * cosPhi[j] * sinPhi[i] -
                                         j * cosPhi[i] * sinPhi[j]) /
                                        (i * i - j * j);
                } else if (i != 0) {
                    A(i, i) = (std::sin(2 * i * phiMax) + 2 * i * phiMax) / (4 * i);
                } else {
                    A(i, i) = phiMax;
                }
            }
        }

        auto svd = A.bdcSvd(Eigen::ComputeFullU | Eigen::ComputeFullV);
        const MatrixX &U = svd.matrixU();
        const MatrixX &V = svd.matrixV();
        const VectorX &sigma = svd.singularValues();

        if (sigma[0] == 0) {
            result.clear();
            result.push_back(0);
            return;
        }

        VectorX temp = VectorX::Zero(n);
        coeffs[0] *= math::Pi;
        coeffs.tail(n-1) *= 0.5 * math::Pi;
        for (int i=0; i<n; ++i) {
            if (sigma[i] < 1e-9f * sigma[0])
                break;
            temp += V.col(i) * U.col(i).dot(coeffs) / sigma[i];
        }
        coeffs = temp;
    }

//.........这里部分代码省略.........

inline std::complex<T> asin(const std::complex<T>& z)
{
   //
   // This implementation is a transcription of the pseudo-code in:
   //
   // "Implementing the complex Arcsine and Arccosine Functions using Exception Handling."
   // T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang.
   // ACM Transactions on Mathematical Software, Vol 23, No 3, Sept 1997.
   //

   //
   // These static constants should really be in a maths constants library:
   //
   static const T one = static_cast<T>(1);
   //static const T two = static_cast<T>(2);
   static const T half = static_cast<T>(0.5L);
   static const T a_crossover = static_cast<T>(1.5L);
   static const T b_crossover = static_cast<T>(0.6417L);
   static const T s_pi = lslboost::math::constants::pi<T>();
   static const T half_pi = s_pi / 2;
   static const T log_two = lslboost::math::constants::ln_two<T>();
   static const T quarter_pi = s_pi / 4;
#ifdef BOOST_MSVC
#pragma warning(push)
#pragma warning(disable:4127)
#endif
   //
   // Get real and imaginary parts, discard the signs as we can 
   // figure out the sign of the result later:
   //
   T x = std::fabs(z.real());
   T y = std::fabs(z.imag());
   T real, imag;  // our results

   //
   // Begin by handling the special cases for infinities and nan's
   // specified in C99, most of this is handled by the regular logic
   // below, but handling it as a special case prevents overflow/underflow
   // arithmetic which may trip up some machines:
   //
   if((lslboost::math::isnan)(x))
   {
      if((lslboost::math::isnan)(y))
         return std::complex<T>(x, x);
      if((lslboost::math::isinf)(y))
      {
         real = x;
         imag = std::numeric_limits<T>::infinity();
      }
      else
         return std::complex<T>(x, x);
   }
   else if((lslboost::math::isnan)(y))
   {
      if(x == 0)
      {
         real = 0;
         imag = y;
      }
      else if((lslboost::math::isinf)(x))
      {
         real = y;
         imag = std::numeric_limits<T>::infinity();
      }
      else
         return std::complex<T>(y, y);
   }
   else if((lslboost::math::isinf)(x))
   {
      if((lslboost::math::isinf)(y))
      {
         real = quarter_pi;
         imag = std::numeric_limits<T>::infinity();
      }
      else
      {
         real = half_pi;
         imag = std::numeric_limits<T>::infinity();
      }
   }
   else if((lslboost::math::isinf)(y))
   {
      real = 0;
      imag = std::numeric_limits<T>::infinity();
   }
   else
   {
      //
      // special case for real numbers:
      //
      if((y == 0) && (x <= one))
         return std::complex<T>(std::asin(z.real()), z.imag());
      //
      // Figure out if our input is within the "safe area" identified by Hull et al.
      // This would be more efficient with portable floating point exception handling;
      // fortunately the quantities M and u identified by Hull et al (figure 3), 
      // match with the max and min methods of numeric_limits<T>.
      //
      T safe_max = detail::safe_max(static_cast<T>(8));
      T safe_min = detail::safe_min(static_cast<T>(4));
//.........这里部分代码省略.........