bool verify(CryptoPP::ECPPoint Q, byte *message, unsigned message_length, CryptoPP::Integer r, CryptoPP::Integer s){
		auto ec = common::ec_parameters().GetCurve();
	  auto G = common::ec_parameters().GetSubgroupGenerator();
	  auto n = common::ec_parameters().GetGroupOrder();

    Integer z = hash_m_to_int(message, message_length, n.ByteCount());
    // verify
    if (Q == ec.Identity()){
      cerr << "Q == O" << endl;
      return false;
    }
    if (!(ec.Multiply(n, Q) == ec.Identity())){
      cerr << "n x Q != O" << endl;
      return false;
    }
    if (r <= 0 || r >= n){
      cerr << "incorrect r" << endl;
      return false;
    } 
    if (s <= 0 || s >= n){
      cerr << "incorrect s" << endl;
      return false;
    } 
    Integer w = s.InverseMod(n);
    Integer u1 = a_times_b_mod_c(z, w, n);
    Integer u2 = a_times_b_mod_c(r, w, n);
    ECPPoint P2 = ec.Add(ec.Multiply(u1, G), ec.Multiply(u2, Q));

    if (P2.x != r){
      cerr << "P2.x != r" << endl;
      return false;
    }
    return true;
  }

    //--------------------------------------------------------------------------------
    /// @brief      三角形を連続でリスト描画する(ワイヤー)
    /// @param[in]  verticies       頂点配列
    /// @param[in]  indicies        インデックス配列
    /// @param[in]  triangleNum     三角形の数
    /// @return     なし
    //--------------------------------------------------------------------------------
    void Renderer::DrawWiredTriangleList( const VECTOR4 verticies[], const unsigned int indicies[], unsigned int triangleNum )
    {
        assert( NULL != m_pGraphics );

        MATRIX4x4 mWVP = m_mWorld;

        Multiply( mWVP, m_mView );
        Multiply( mWVP, m_mProjection );

        const unsigned int polygonVerticies = 3;

        for( unsigned int triangleCnt = 0; triangleCnt < triangleNum; ++triangleCnt )
        {
            VECTOR4 vertex0 = verticies[indicies[triangleCnt * polygonVerticies]];
            VECTOR4 vertex1 = verticies[indicies[triangleCnt * polygonVerticies + 1]];
            VECTOR4 vertex2 = verticies[indicies[triangleCnt * polygonVerticies + 2]];

            Transform( vertex0, mWVP );
            Transform( vertex1, mWVP );
            Transform( vertex2, mWVP );
            m_pGraphics->DrawLine( VECTOR2( vertex0.x / vertex0.w, vertex0.y / vertex0.w ), VECTOR2( vertex1.x / vertex1.w, vertex1.y / vertex1.w ), m_color );
            m_pGraphics->DrawLine( VECTOR2( vertex1.x / vertex1.w, vertex1.y / vertex1.w ), VECTOR2( vertex2.x / vertex2.w, vertex2.y / vertex2.w ), m_color );
            m_pGraphics->DrawLine( VECTOR2( vertex2.x / vertex2.w, vertex2.y / vertex2.w ), VECTOR2( vertex0.x / vertex0.w, vertex0.y / vertex0.w ), m_color );
        }
    }

 void Multiply(const JntArrayVel& src,const doubleVel& factor,JntArrayVel& dest)
 {
     Multiply(src.q,factor.grad,dest.q);
     Multiply(src.qdot,factor.t,dest.qdot);
     Add(dest.qdot,dest.q,dest.qdot);
     Multiply(src.q,factor.t,dest.q);
 }

// Pout = Rx*Ry*Rz*Pin
void Frames::Rotate(float yaw, float pitch, float roll, float* pointIn, float* pointOut)
{
	// Rz()
	Ry[0][0] = cos(yaw);
	Ry[0][1] = -sin(yaw);
	Ry[1][0] = sin(yaw);
	Ry[1][1] = cos(yaw);
	Ry[2][2] = 1.0f;

	// Rx()
	Rr[0][0] = 1.0f;
	Rr[1][1] = cos(roll);
	Rr[1][2] = -sin(roll);
	Rr[2][1] = sin(roll);
	Rr[2][2] = cos(roll);
	
	// Ry()
	Rp[0][0] = cos(pitch);
	Rp[0][2] = sin(pitch);
	Rp[1][1] = 1.0f;
	Rp[2][0] = -sin(pitch);
	Rp[2][2] = cos(pitch);

	float t[3];
	Multiply(&Ry, pointIn, pointOut);
	Multiply(&Rp, pointOut, t);
	Multiply(&Rr, t, pointOut);

}

bool ON_Matrix::Multiply( const ON_Matrix& a, const ON_Matrix& b )
{
  int i, j, k, mult_count;
  double x;
  if (a.ColCount() != b.RowCount() )
    return false;
  if ( a.RowCount() < 1 || a.ColCount() < 1 || b.ColCount() < 1 )
    return false;
  if ( this == &a ) {
    ON_Matrix tmp(a);
    return Multiply(tmp,b);
  }
  if ( this == &b ) {
    ON_Matrix tmp(b);
    return Multiply(a,tmp);
  }
  Create( a.RowCount(), b.ColCount() );
  mult_count = a.ColCount();
  double const*const* am = a.ThisM();
  double const*const* bm = b.ThisM();
  double** this_m = ThisM();
  for ( i = 0; i < m_row_count; i++ ) for ( j = 0; j < m_col_count; j++ ) {
    x = 0.0;
    for (k = 0; k < mult_count; k++ ) {
      x += am[i][k] * bm[k][j];
    }
    this_m[i][j] = x;
  }
  return true;  
}

void Render::DrawBonds(D3DMATRIX* view, D3DMATRIX* projection)
{
    // Apply the technique contained in the effect
    BYTE previousRenderedMaterial = 255;

    if (bondRenderOptions.UseSingleMaterial)
    {
        bondRenderOptions.Material.Apply();
    }

    for(int iBond = 0; iBond < bondCount; iBond++)
    {
        Bond* currentBond = &bonds[iBond];

        Material* material = bondRenderOptions.UseSingleMaterial ?
                             &bondRenderOptions.Material :
                             elementMaterials + currentBond->Material;


        // Transform along to direction
        D3DVECTOR direction;
        direction.x = currentBond->End.x - currentBond->Begin.x;
        direction.y = currentBond->End.y - currentBond->Begin.y;
        direction.z = currentBond->End.z - currentBond->Begin.z;
        float height = GetLength(direction);
        //Normalize(direction);
        D3DMATRIX alongTo = TransformAlongTo(direction);

        // Scale
        D3DXMATRIX scale = D3DMATRIXCREATE(
                               bondRenderOptions.BondSize,0,0,0,
                               0,1,0,0,
                               0,0,bondRenderOptions.BondSize,0,
                               0,0,0,1.0f);

        // Translate the cylinder
        D3DXMATRIX translate = D3DMATRIXCREATE(
                                   1.0f,0,0,0,
                                   0,1.0f,0,0,
                                   0,0,1.0f,0,
                                   currentBond->Begin.x,currentBond->Begin.y,currentBond->Begin.z,1.0f);

        D3DXMATRIX world = Multiply(Multiply(scale, alongTo), translate);
        material->SetMatrices(&world, view, projection);


        // Setup material's parameters
        if (!bondRenderOptions.UseSingleMaterial && previousRenderedMaterial != currentBond->Material)
        {
            // Skip if the params has been already presented
            previousRenderedMaterial = currentBond->Material;
            material->Apply();
        }

        // Render the mesh with the applied technique
        highPolyCylinder->Draw();
    }
}

POLY Inverse_Poly ( int n, POLY M[n][n] )
{ 
   POLY P[n][n], Q[n][n], t_M1[n][n], t_M2[n][n];
   POLY ds;
   int i, j, dk, dl, dgcd, dc_gcd, dad, dbd;
   dcmplx **gcd_coeff, tmp;

   copy(n, n, M, t_M1);  

   Smith(n, n, t_M1, P, Q);
   /* printf("the smith form is:\n");  print1(n, n, t_M1); */
   ds=assign_poly(t_M1[n-1][n-1]);

   for( i=0; i<n-1; i++ )
     { 
       gcd_coeff = ExtPolyGcd(ds.d+1, ds.p, t_M1[i][i].d+1, t_M1[i][i].p, &dgcd, &dk, &dl, &dbd, &dad);  
       free(t_M1[i][i].p);
       t_M1[i][i].d = dad;
       t_M1[i][i].p = assign(dad,gcd_coeff[4]);
       for( j=0; j<5; j++)
         free(gcd_coeff[j]);
       free(gcd_coeff);
     }
   free(t_M1[n-1][n-1].p);

   t_M1[n-1][n-1].p = (dcmplx*) calloc(1, sizeof(dcmplx));
   t_M1[n-1][n-1].d = 0;
   t_M1[n-1][n-1].p[0] = one;

   /* calculate ds*Q*(t_M1's inverse) * P = ds*M  */
   Multiply( n, n, n, Q, t_M1, t_M2 );

   free_matrix ( n, n, M );
   Multiply( n, n, n, t_M2, P, M );

   if(ds.p[ds.d].re<0)
   {
     negative(ds.d, ds.p);
     neg_polymatrix( n, n, M );
   }   

   /* make the leading coefficient of ds one */
   tmp = ds.p[ds.d];
   divide_by_number(ds.d, ds.p, tmp); 

   for(i=0; i<n; i++)
     for(j=0; j<n; j++)
       {  
        divide_by_number(M[i][j].d, M[i][j].p, tmp);
       }
   /* printf("ds="); Print_Poly(ds.d , ds.p); */
  
   free_matrix ( n, n, t_M1 );
   free_matrix ( n, n, t_M2 );
   free_matrix ( n, n, P );
   free_matrix ( n, n, Q );
   return ds;
}

 static uint64_t Power(uint64_t a, int64_t b) {
     uint64_t r = 1;
     for (; b > 0; b >>= 1) {
         if (b & 1)
             r = Multiply(r, a);
         a = Multiply(a, a);
     }
     return r;
 }

 static Float4x4 VFunction Multiply(const Float4x4& matrixA, const Float4x4& matrixB)
 {
     Float4x4 result;
     result.x = Multiply(matrixA.x, matrixB);
     result.y = Multiply(matrixA.y, matrixB);
     result.z = Multiply(matrixA.z, matrixB);
     result.w = Multiply(matrixA.w, matrixB);
     return result;
 }

void Xenon::Math::Vec2 :: Interpolate (  Vec2 & Target, const Vec2 & Source, const Vec2 & B, const float Fraction )
{
	
	Vec2 Temp;
	
	Multiply ( Target, Source, Fraction );
	Multiply ( Temp, B, 1 - Fraction );
	Add ( Target, Temp );
	
}

void Vector3 :: Interpolate ( Vector3 & A, Vector3 & B, double Fraction, Vector3 & Result )
{
	
	Vector3 Tem;
	
	Multiply ( Result, Fraction );
	Multiply ( B, 1 - Fraction, Tem );
	Add ( Result, Tem );
	
};

bool LinearSystem<Real>::SolveSymmetricCG (const GMatrix<Real>& A,
    const Real* B, Real* X)
{
    // Based on the algorithm in "Matrix Computations" by Golum and Van Loan.
    assertion(A.GetNumRows() == A.GetNumColumns(), "Matrix must be square\n");
    int size = A.GetNumRows();
    Real* R = new1<Real>(size);
    Real* P = new1<Real>(size);
    Real* W = new1<Real>(size);

    // The first iteration.
    size_t numBytes = size*sizeof(Real);
    memset(X, 0, numBytes);
    memcpy(R, B, numBytes);
    Real rho0 = Dot(size, R, R);
    memcpy(P, R, numBytes);
    Multiply(A, P, W);
    Real alpha = rho0/Dot(size, P, W);
    UpdateX(size, X, alpha, P);
    UpdateR(size, R, alpha, W);
    Real rho1 = Dot(size, R, R);

    // The remaining iterations.
    const int imax = 1024;
    int i;
    for (i = 1; i < imax; ++i)
    {
        Real root0 = Math<Real>::Sqrt(rho1);
        Real norm = Dot(size, B, B);
        Real root1 = Math<Real>::Sqrt(norm);
        if (root0 <= ZeroTolerance*root1)
        {
            break;
        }

        Real beta = rho1/rho0;
        UpdateP(size, P, beta, R);
        Multiply(A, P, W);
        alpha = rho1/Dot(size, P, W);
        UpdateX(size, X, alpha, P);
        UpdateR(size, R, alpha, W);
        rho0 = rho1;
        rho1 = Dot(size, R, R);
    }

    delete1(W);
    delete1(P);
    delete1(R);

    return i < imax;
}

LispObject* PowerFloat(LispObject* int1, LispObject* int2, LispEnvironment& aEnvironment,int aPrecision)
{
    if (int2->Number(aPrecision)->iNumber->iExp != 0)
        throw LispErrNotInteger();

    // Raising to the power of an integer can be done fastest by squaring
    // and bitshifting: x^(a+b) = x^a*x^b . Then, regarding each bit
    // in y (seen as a binary number) as added, the algorithm becomes:
    //
    ANumber x(*int1->Number(aPrecision)->iNumber);
    ANumber y(*int2->Number(aPrecision)->iNumber);
    bool neg = y.iNegative;
    y.iNegative=false;

    // result <- 1
    ANumber result("1",aPrecision);
    // base <- x
    ANumber base(aPrecision);
    base.CopyFrom(x);

    ANumber copy(aPrecision);

    // while (y!=0)
    while (!y.IsZero())
    {
        // if (y&1 != 0)
        if ( (y[0] & 1) != 0)
        {
            // result <- result*base
            copy.CopyFrom(result);
            Multiply(result,copy,base);
        }
        // base <- base*base
        copy.CopyFrom(base);
        Multiply(base,copy,copy);
        // y <- y>>1
        BaseShiftRight(y,1);
    }

    if (neg)
    {
        ANumber one("1",aPrecision);
        ANumber dummy(10);
        copy.CopyFrom(result);
        Divide(result,dummy,one,copy);
    }

    // result
    return FloatToString(result, aEnvironment);
}

void QR::ComputeQRMatrices() {
 CDMatrix p;
 for(int i = 0; i < r_.size() - 1; ++i) {
  Normalize(i);
  ConstructP(p);
  if(q_.size() == 0) {
   q_ = p;
  } else {
   Multiply(q_, p);
  }
  Multiply(r_, p);
 } 
 TransposeQ();
}

PointTD* ProjectionOrtoDD::Rotate3D(PointTD pPointTD, CameraTD* pCamera) {
	double X[N] = {pPointTD.fX, pPointTD.fY, pPointTD.fZ, 1};
	double Xm[N] = {0, 0, 0, 1};
	Multiply(X, Ares, Xm);

	if (pCamera->iMode == CAM_CENTRAL) {
		Xm[0] = Xm[0] / (fabs(Xm[3]) + inac);
		Xm[1] = Xm[1] / (fabs(Xm[3]) + inac);
	}
	else {
		Xm[0] = Xm[0];
		Xm[1] = Xm[1];
	}
	double Tm[N] = {0, 0, 0, 1};
	if (ObjectID == "LabaPoint") {
		Tm[0] = Xm[0];
		Tm[1] = Xm[1];
	}
	if (c <= inac && pCamera->iMode == CAM_ORTO)
		return new PointTD(0, 0, 0, pPointTD.iAction, pPointTD.iType, pPointTD.sText);
	else if (c / 2.0 <= pCamera->TCheck[2] && pCamera->iMode == CAM_CENTRAL)
		return new PointTD(0, 0, 0, pPointTD.iAction, pPointTD.iType, pPointTD.sText);
	else
		return new PointTD(Xm[0], Xm[1], Xm[2], pPointTD.iAction, pPointTD.iType, pPointTD.sText);   
}

    NekMatrix<typename NekMatrix<RhsDataType, RhsMatrixType>::NumberType, StandardMatrixTag>
    Multiply(const DataType& lhs,
                const NekMatrix<RhsDataType, RhsMatrixType>& rhs)

    {
        return Multiply(rhs, lhs);
    }

//*******************************************************************
// ilEuclidNorm
//   normalizes vectors of dim=n in mxn input matrix with the
//   Euclidean norm
//*******************************************************************
void ilEuclidNorm(pMat const& src, pMat dst)
{
    int m = src->rows;
    int n = src->cols;
    int type = src->type;
    pMat ones = CreateMat(n, 1, type);
    pMat res1 = CreateMat(m, n, type);
    pMat norm = CreateMat(m, 1, type);
    SetValue(ones, cvScalar(1.0));

    //*** compute the norm
    Multiply(src, src, res1);
    MatrixMultiply(res1, ones, norm);
    PowerMatrix(norm, norm, .5);

    //*** normalize columns of src
#if 0 //*** matrix version
    pMat normmat = CreateMat(m, n, CV_32FC1);
    pMat onesrow = CreateMat(1, n, CV_32FC1);
    SetValue(onesrow, cvScalar(1.0));
    MatrixMultiply(norm, onesrow, normmat);
    Divide(src, normmat, dst);
#endif

#if 1 // *** column version
    CvMat colmat1 = cvMat(0, 0, 0, 0);
    CvMat colmat2 = cvMat(0, 0, 0, 0);
    for (int i = 0; i<n; ++i)
    {
        cvGetCol(src.get(), &colmat1, i);
        cvGetCol(dst.get(), &colmat2, i);
        Divide(dummyGuard(&colmat1), norm, dummyGuard(&colmat2));
    }
#endif
}

void Test(matrix A, matrix B, matrix Res)
/*
 * Runs a multiplication test on an array.  Calculates and prints the
 * time it takes to multiply the matrices.
 */
{
#ifndef UPPSALAWCET
   struct timeval StartTime, StopTime;
   float TotalTime;
#endif

   /* ***UPPSALA WCET***: don't print or time */
#ifndef UPPSALAWCET
   gettimeofday(&StartTime, NULL);
#endif

   Initialize(A);
   Initialize(B);

   Multiply(A, B, Res);

   /* ***UPPSALA WCET***: don't print or time */
#ifndef UPPSALAWCET
   gettimeofday(&StopTime, NULL);
   TotalTime = (1000 * (StopTime.tv_sec - StartTime.tv_sec) + (StopTime.tv_usec - StartTime.tv_usec) / 1000) / 1000.0;
   printf("    - Size of array is %d (%ld CLOCKS_PER_SEC)\n", UPPERLIMIT, CLOCKS_PER_SEC);
   printf("    - Total multiplication time is %3.3f seconds\n\n", TotalTime);
#endif
}

inline Int RegularizedSolveAfterNoPromote
( const SparseMatrix<F>& A,
  const Matrix<Base<F>>& reg,
  const Matrix<Base<F>>& d, 
  const vector<Int>& invMap,
  const ldl::NodeInfo& info,
  const ldl::Front<F>& front, 
        Matrix<F>& B,
  Base<F> relTol,
  Int maxRefineIts, 
  bool progress,
  bool time )
{
    DEBUG_CSE

    // TODO: Use time in these lambdas
    auto applyA =
      [&]( const Matrix<F>& X, Matrix<F>& Y )
      {
        Y = X;
        DiagonalScale( LEFT, NORMAL, reg, Y ); 
        Multiply( NORMAL, F(1), A, X, F(1), Y );
      };
    auto applyAInv = 
      [&]( Matrix<F>& Y )
      {
        DiagonalSolve( LEFT, NORMAL, d, Y );
        ldl::MatrixNode<F> YNodal( invMap, info, Y );
        ldl::SolveAfter( info, front, YNodal );
        YNodal.Push( invMap, info, Y );
        DiagonalSolve( LEFT, NORMAL, d, Y );
      };

    return RefinedSolve( applyA, applyAInv, B, relTol, maxRefineIts, progress );
}

GF256::Element GF256::MultiplicativeInverse(Element a) const
{
	Element result = a;
	for (int i=1; i<7; i++)
		result = Multiply(Square(result), a);
	return Square(result);
}