C++ AB函數代碼示例

/*-------------------------------------------------------------------------*
 * PL_CURRENT_MIRROR_ALT_0                                                 *
 *                                                                         *
 *-------------------------------------------------------------------------*/
Bool
Pl_Current_Mirror_Alt_0(void)
{
  /*  int stm; */
  WamWord m_stm_word;
  StmLst *m;

  Pl_Update_Choice_Point((CodePtr) Prolog_Predicate(CURRENT_MIRROR_ALT, 0), 0);

  /*  stm = AB(B, 0); */
  m_stm_word = AB(B, 1);
  m = (StmLst *) AB(B, 2);

  if (m->next)			/* non deterministic case */
    {
#if 0				/* the following data is unchanged */
      AB(B, 0) = stm;
      AB(B, 1) = m_stm_word;
#endif
      AB(B, 2) = (WamWord) m->next;
    }
  else
    Delete_Last_Choice_Point();

  return Pl_Get_Integer(m->stm, m_stm_word);
}

/*-------------------------------------------------------------------------*
 * PL_CURRENT_STREAM_ALT_0                                                 *
 *                                                                         *
 *-------------------------------------------------------------------------*/
Bool
Pl_Current_Stream_Alt_0(void)
{
  WamWord stm_word;
  int stm;

  Pl_Update_Choice_Point((CodePtr) Prolog_Predicate(CURRENT_STREAM_ALT, 0), 0);

  stm_word = AB(B, 0);
  stm = AB(B, 1);

  for (; stm <= pl_stm_last_used; stm++)
    if (pl_stm_tbl[stm])
      break;

  if (stm >= pl_stm_last_used)
    {
      Delete_Last_Choice_Point();
      if (stm > pl_stm_last_used)
	return FALSE;
    }
  else				/* non deterministic case */
    {
#if 0 /* the following data is unchanged */
      AB(B, 0) = stm_word;
#endif
      AB(B, 1) = stm + 1;
    }

  return Pl_Get_Integer(stm, stm_word);
}

/*-------------------------------------------------------------------------*
 * PL_FOREIGN_UPDATE_CHOICE                                                *
 *                                                                         *
 *-------------------------------------------------------------------------*/
void
Pl_Foreign_Update_Choice(CodePtr codep_alt, int arity, int choice_size)
{
  pl_foreign_bkt_counter = AB(B, arity) + 1;
  AB(B, arity) = pl_foreign_bkt_counter;

  pl_foreign_bkt_buffer = (char *) (&(AB(B, arity + choice_size)));

  if (pl_foreign_bkt_counter > 0)
    {
      Pl_Update_Choice_Point(codep_alt, arity);
    }
}

inline void
SingularValuesUpper
( DistMatrix<Complex<Real> >& A,
  DistMatrix<Real,VR,STAR>& s,
  double heightRatio=1.2 )
{
#ifndef RELEASE
    PushCallStack("svd::SingularValuesUpper");
    if( heightRatio <= 1.0 )
        throw std::logic_error("Nonsensical switchpoint for SingularValues");
#endif
    typedef Complex<Real> C;
    const Grid& g = A.Grid();
    const int m = A.Height();
    const int n = A.Width();
    if( m >= heightRatio*n )
    {
        DistMatrix<C,MD,STAR> t(g);
        QR( A, t );
        DistMatrix<C> AT(g),
                      AB(g);
        PartitionDown
        ( A, AT,
             AB, n );
        MakeTrapezoidal( LEFT, UPPER, 0, AT );
        SimpleSingularValuesUpper( AT, s );
    }
    else
    {
        SimpleSingularValuesUpper( A, s );
    }
#ifndef RELEASE
    PopCallStack();
#endif
}

  // Checks whether 3D points p lies inside or outside of the triangle ABC
  bool includePointTriangle(Vec3r& P,
			    Vec3r& A,
			    Vec3r& B,
			    Vec3r& C) {
    Vec3r AB(B - A);
    Vec3r BC(C - B);
    Vec3r CA(A - C);
    Vec3r AP(P - A);
    Vec3r BP(P - B);
    Vec3r CP(P - C);

    Vec3r N(AB ^ BC); // triangle's normal

    N.normalize();

    Vec3r J(AB ^ AP), K(BC ^ BP), L(CA ^ CP);
    J.normalize();
    K.normalize();
    L.normalize();

    if(J * N < 0)
      return false; // on the right of AB

    if(K * N < 0)
      return false; // on the right of BC

    if(L * N < 0)
      return false; // on the right of CA

    return true;
  }

void UInputAxisKeyDelegateBinding::BindToInputComponent(UInputComponent* InputComponent) const
{
	TArray<FInputAxisKeyBinding> BindsToAdd;

	for (int32 BindIndex=0; BindIndex<InputAxisKeyDelegateBindings.Num(); ++BindIndex)
	{
		const FBlueprintInputAxisKeyDelegateBinding& Binding = InputAxisKeyDelegateBindings[BindIndex];

		FInputAxisKeyBinding AB( Binding.AxisKey );
		AB.bConsumeInput = Binding.bConsumeInput;
		AB.bExecuteWhenPaused = Binding.bExecuteWhenPaused;
		AB.AxisDelegate.BindDelegate(InputComponent->GetOwner(), Binding.FunctionNameToBind);

		if (Binding.bOverrideParentBinding)
		{
			for (int32 ExistingIndex = InputComponent->AxisKeyBindings.Num() - 1; ExistingIndex >= 0; --ExistingIndex)
			{
				const FInputAxisKeyBinding& ExistingBind = InputComponent->AxisKeyBindings[ExistingIndex];
				if (ExistingBind.AxisKey == AB.AxisKey)
				{
					InputComponent->AxisKeyBindings.RemoveAt(ExistingIndex);
				}
			}
		}

		// To avoid binds in the same layer being removed by the parent override temporarily put them in this array and add later
		BindsToAdd.Add(AB);
	}

	for (int32 Index=0; Index < BindsToAdd.Num(); ++Index)
	{
		InputComponent->AxisKeyBindings.Add(BindsToAdd[Index]);
	}
}

bool Intersect(CVector A, CVector B, CVector C, CVector D, CVector *M)
{//ABCD are assumed to lie in the xy plane
	//returns true and intersection point M if AB and CD intersect inside
	//returns false and intersection point M if AB and CD intersect outside
	M->x = 0.0;
	M->y = 0.0;
	M->z = 0.0;
	CVector AB(B.x-A.x, B.y-A.y, B.z-A.z); 
	CVector CD(D.x-C.x, D.y-C.y, D.z-C.z); 

	//Cramer's rule

	double Det  = -AB.x * CD.y + CD.x * AB.y;
	if(Det==0.0) 
	{
		//vectors are parallel, no intersection
		return false;
	}
	double Det1 = -(C.x-A.x)*CD.y + (C.y-A.y)*CD.x;
	double Det2 = -(C.x-A.x)*AB.y + (C.y-A.y)*AB.x;

	double t = Det1/Det;
	double u = Det2/Det;

	M->x = A.x + t*AB.x;
	M->y = A.y + t*AB.y;

	if (0.0<=t && t<=1.0 && 0.0<=u && u<=1.0)	return true;//M is between A and B
	else										return false;//M is outside
}

	Tetra(Point3 A,Point3 B,Point3 C,Point3 D) :
		m_A(A),m_B(B),m_C(C),m_D(D)
	{
		Vector3 AB(B-A);
		Vector3 AC{ C[0]-A[0], C[1]-A[1], C[2]-A[2] };
		Vector3 AD{ D[0]-A[0], D[1]-A[1], D[2]-A[2] };
		Vector3 BC{ C[0]-B[0], C[1]-B[1], C[2]-B[2] };
		Vector3 BD{ D[0]-B[0], D[1]-B[1], D[2]-B[2] };

		m_fABC = ImplicitPlane<float,3>(cross(AB,AC),A);
		if (m_fABC(D) < 0)
			m_fABC.flip();

		m_fABD = ImplicitPlane<float,3>(cross(AB,AD),A);
		if (m_fABD(C) < 0)
			m_fABD.flip();

		m_fACD = ImplicitPlane<float,3>(cross(AC,AD),A);
		if (m_fACD(B) < 0)
			m_fACD.flip();

		m_fBCD = ImplicitPlane<float,3>(cross(BC,BD),B);
		if (m_fBCD(A) < 0)
			m_fBCD.flip();
	}

    /// Assign this matrix to a product of three other matrices
    /// @param mult3 :: Matrix multiplication helper object.
    GSLMatrix& GSLMatrix::operator=(const GSLMatrixMult3& mult3)
    {
      // sizes of the result matrix
      size_t n1 = mult3.tr1 ? mult3.m_1.size2() : mult3.m_1.size1();
      size_t n2 = mult3.tr3 ? mult3.m_3.size1() : mult3.m_3.size2();

      this->resize(n1,n2);

      // intermediate matrix 
      GSLMatrix AB( n1, mult3.m_2.size2() );

      CBLAS_TRANSPOSE tr1 = mult3.tr1? CblasTrans : CblasNoTrans;
      CBLAS_TRANSPOSE tr2 = mult3.tr2? CblasTrans : CblasNoTrans;
      CBLAS_TRANSPOSE tr3 = mult3.tr3? CblasTrans : CblasNoTrans;

      // AB = m_1 * m_2
      gsl_blas_dgemm (tr1, tr2,
                      1.0, mult3.m_1.gsl(), mult3.m_2.gsl(),
                      0.0, AB.gsl());

      // this = AB * m_3
      gsl_blas_dgemm (CblasNoTrans, tr3,
                      1.0, AB.gsl(), mult3.m_3.gsl(),
                      0.0, gsl());

      return *this;
    }

inline void
Explicit( DistMatrix<F>& A, DistMatrix<F>& R, bool colPiv=false )
{
#ifndef RELEASE
    CallStackEntry cse("qr::Explicit");
#endif
    const Grid& g = A.Grid();
    DistMatrix<F,MD,STAR> t(g);
    if( colPiv )
    {
        DistMatrix<Int,VR,STAR> p(g);
        QR( A, t, p );
    }
    else
    {
        QR( A, t );
    }
    DistMatrix<F> AT(g),
                  AB(g);
    PartitionDown
    ( A, AT,
         AB, Min(A.Height(),A.Width()) );
    R = AT;
    MakeTriangular( UPPER, R );
    ExpandPackedReflectors( LOWER, VERTICAL, UNCONJUGATED, 0, A, t );
}

  Vec3f VertexOrientation3DF0D::operator()(Interface0DIterator& iter) {
    Vec3r A,C;
    Vec3r B(iter->getX(), iter->getY(), iter->getZ());
    if(iter.isBegin())
      A = Vec3r(iter->getX(), iter->getY(), iter->getZ());
    else
      {
	Interface0DIterator previous = iter;
	--previous ;
	A = Vec3r(previous->getX(), previous->getY(), previous->getZ());
      }
    Interface0DIterator next = iter;
    ++next ;
    if(next.isEnd())
      C = Vec3r(iter->getX(), iter->getY(), iter->getZ());
    else
      C = Vec3r(next->getX(), next->getY(), next->getZ());

    Vec3r AB(B-A);
    if(AB.norm() != 0)
      AB.normalize();
    Vec3r BC(C-B);
    if(BC.norm() != 0)
      BC.normalize();
    Vec3f res (AB + BC);
    if(res.norm() != 0)
      res.normalize();
    return res;
  }

double QuadQuality(const Vertex & a,const Vertex &b,const Vertex &c,const Vertex &d)
{
  // calcul de 4 angles --
  R2 A((R2)a),B((R2)b),C((R2)c),D((R2)d);
  R2 AB(B-A),BC(C-B),CD(D-C),DA(A-D);
  //  Move(A),Line(B),Line(C),Line(D),Line(A);
  const Metric & Ma  = a;
  const Metric & Mb  = b;
  const Metric & Mc  = c;
  const Metric & Md  = d;
    
  double lAB=Norme2(AB);
  double lBC=Norme2(BC);
  double lCD=Norme2(CD);
  double lDA=Norme2(DA);
  AB /= lAB;  BC /= lBC;  CD /= lCD;  DA /= lDA;
  // version aniso 
  double cosDAB= Ma(DA,AB)/(Ma(DA)*Ma(AB)),sinDAB= Det(DA,AB);
  double cosABC= Mb(AB,BC)/(Mb(AB)*Mb(BC)),sinABC= Det(AB,BC);
  double cosBCD= Mc(BC,CD)/(Mc(BC)*Mc(CD)),sinBCD= Det(BC,CD);
  double cosCDA= Md(CD,DA)/(Md(CD)*Md(DA)),sinCDA= Det(CD,DA);
  double sinmin=Min(Min(sinDAB,sinABC),Min(sinBCD,sinCDA));
  // cout << A << B << C << D ;
  // cout << " sinmin " << sinmin << " " << cosDAB << " " << cosABC << " " << cosBCD << " " << cosCDA << endl;
  // rattente(1);
  if (sinmin<=0) return sinmin;
  return 1.0-Max(Max(Abs(cosDAB),Abs(cosABC)),Max(Abs(cosBCD),Abs(cosCDA)));
}

int VertexOrientation2DF0D::operator()(Interface0DIterator& iter)
{
	Vec2f A, C;
	Vec2f B(iter->getProjectedX(), iter->getProjectedY());
	if (iter.isBegin()) {
		A = Vec2f(iter->getProjectedX(), iter->getProjectedY());
	}
	else {
		Interface0DIterator previous = iter;
		--previous ;
		A = Vec2f(previous->getProjectedX(), previous->getProjectedY());
	}
	Interface0DIterator next = iter;
	++next;
	if (next.isEnd())
		C = Vec2f(iter->getProjectedX(), iter->getProjectedY());
	else
		C = Vec2f(next->getProjectedX(), next->getProjectedY());

	Vec2f AB(B - A);
	if (AB.norm() != 0)
		AB.normalize();
	Vec2f BC(C - B);
	if (BC.norm() != 0)
		BC.normalize();
	result = AB + BC;
	if (result.norm() != 0)
		result.normalize();
	return 0;
}

///////////////////////////////////////////////////////////////////
// CLOSEST IN
///////////////////////////////////////////////////////////////////
Vec3f ClosestIn(Simplex& P){
	assert(P.GetSize() > 0);
	if(P.GetSize() == 1) return P.at(0);
	else if(P.GetSize() == 2){ //We have a line.
		Vector3f AO(-P.at(1)[0], -P.at(1)[1], -P.at(1)[2]);
		Vector3f AB(P.at(0)[0] - P.at(1)[0], P.at(0)[1] - P.at(1)[1], P.at(0)[2] - P.at(1)[2]);
		Vector3f BO(-P.at(0)[0], -P.at(0)[1], -P.at(0)[2]);
		float v = (AB.dot(AO))/(AB.dot(AO) + (-AB).dot(BO));
		return (P.at(1)*(1-v) + (P.at(0))*v);
	}else{//We have a triangle
		Vector3f A(P.at(2)[0], P.at(2)[1], P.at(2)[2]);
		Vector3f B(P.at(1)[0], P.at(1)[1], P.at(1)[2]);
		Vector3f C (P.at(0)[0], P.at(0)[1], P.at(0)[2]);
		Vector3f AB = B-A;
		Vector3f AC = C-A;
		Vector3f N = AB.cross(AC);
		Vector3f Nab = N.cross(AB);
		Vector3f Nac = N.cross(AC);

		float v = (-A).dot(Nac)/((AB).dot(Nac));
		float w = (-A).dot(Nab)/((AC).dot(Nab));
		float u = 1-v-w;

		return (P.at(2)*u + P.at(1)*v + P.at(0)*w);
	}
}

MatrixXd kroen_product(MatrixXd A, MatrixXd B){
    unsigned int Ar = A.rows(), Ac = A.cols(), Br = B.rows(), Bc = B.cols();
    MatrixXd AB(Ar*Br,Ac*Bc);
    for (unsigned int i=0; i<Ar; ++i)
        for (unsigned int j=0; j<Ac; ++j)
            AB.block(i*Br,j*Bc,Br,Bc) = A(i,j)*B;
    return AB;
}