C++ Poly::eval方法代码示例

// integrate polynomial from lower to upper
double Poly::integrate(double lower, double upper) const{
    
    // perform integration
    Poly integral = this->integrate();
    
    return integral.eval(upper) - integral.eval(lower);
}

int main() {
    BIGNUM *tmp = NULL;
    BIGNUM *p = NULL, *q = NULL, *N = NULL, *e_value = NULL, *d = NULL;
    BIGNUM *r0=NULL,*r1=NULL,*r2=NULL;

    bn_ctx = BN_CTX_new();
    BN_CTX_init(bn_ctx);
    BN_CTX_start(bn_ctx);

    /* Generate RSA key using safe primes */
    p = BN_CTX_get(bn_ctx);
    BN_generate_prime_ex(p, 1024, 1, NULL, NULL, NULL);
    q = BN_CTX_get(bn_ctx);
    BN_generate_prime_ex(q, 1024, 1, NULL, NULL, NULL);
    N = BN_CTX_get(bn_ctx);
    BN_mul(N, p, q, bn_ctx);
    r0 = BN_CTX_get(bn_ctx);
    r1 = BN_CTX_get(bn_ctx);
    r2 = BN_CTX_get(bn_ctx);
    e_value = BN_new();
    d = BN_CTX_get(bn_ctx);
    BN_set_word(e_value, e);

    /* calculate d */
    BN_sub(r1,p,BN_value_one());    /* p-1 */
    BN_sub(r2,q,BN_value_one());    /* q-1 */
    BN_mul(r0,r1,r2,bn_ctx);   /* (p-1)(q-1) */
    BN_mod_inverse(d,e_value,N,bn_ctx);  /* d */

    Poly *poly = new Poly(threshold-1, N);

    /* Extract secret exponent */
    poly->set_coeff(0, d);

    /* Output modulus and public exponent */
    cout << BN_bn2hex(N) << endl;
    cout << BN_bn2hex(e_value) << endl;

    for (int i=0; i<num_nodes; i++) {
        tmp = poly->eval(i + 1);
        /* Output secret share */
        cout << BN_bn2hex(tmp) << endl;
    }

    delete poly;
    BN_CTX_end(bn_ctx);     
    BN_CTX_free(bn_ctx);
    return 0;
}

//.........这里部分代码省略.........
    // stiffness matrix
    Matrix KbendY(ns, ns);
    matrixAssembly(EIy, ns, fpp, 1.0/pow(L,3), KbendY);

    
    // distributed applied loads
    Vector FbendY(ns);
    vectorAssembly(Py, ns, f, L, FbendY);

    // put into global matrix (mass and incremental stiffness are same in x an y)
    idx = {2, 3, 8, 9};
    //idx(0) = 2; idx(1) = 3; idx(2) = 8; idx(3) = 9;
    
    for (int ii=0; ii<ns; ii++) {
      F(idx[ii]) = FbendY(ii);
      for (int jj=0; jj<ns; jj++) {
	K(idx[ii], idx[jj]) = KbendY(ii,jj);
	M(idx[ii], idx[jj]) = Mbend(ii,jj);
	Ndist(idx[ii], idx[jj]) = Nbend_dist(ii,jj);
	Nconst(idx[ii], idx[jj]) = Nbend_const(ii,jj);
      }
    }
    //project(K, idx, idx) = KbendY;
    //project(M, idx, idx) = Mbend;
    //project(Ndist, idx, idx) = Nbend_dist;
    //project(Nconst, idx, idx) = Nbend_const;
    //project(F, idx) = FbendY;
    
    
    // ----------- axial ----------------
    const int nsz = 2; // number of shape functions
    
    Poly fz[nsz] = {
		    Poly(2, -1.0, 1.0),
		    Poly(2, 1.0, 0.0)
    };
    
    // derivatives of the shape function
    Poly fzp [nsz];
    Poly fzpp [nsz];
    
    for (i = 0; i < nsz; i++) {
      fzp[i] = fz[i].differentiate();
      fzpp[i] = fzp[i].differentiate();
    }
    
    // stiffness matrix
    Matrix Kaxial(nsz, nsz);
    matrixAssembly(EA, nsz, fzp, 1.0/L, Kaxial);

    
    // inertia matrix
    Matrix Maxial(nsz, nsz);
    matrixAssembly(rhoA, nsz, fz, L, Maxial);
    
    // axial loads already given (work equivalent approach not appropriate for distributed axial loads)
    Vector Faxial(nsz);
    Faxial(0) = FzfromPz.eval(0.0);
    Faxial(1) = FzfromPz.eval(1.0);
    
    // put into global matrix
    std::vector<int> idx_z = {4, 10};
    //idx_z(0) = 4; idx_z(1) = 10;
    
    for (int ii=0; ii<nsz; ii++) {
      F(idx_z[ii]) = Faxial(ii);
      for (int jj=0; jj<nsz; jj++) {
	K(idx_z[ii], idx_z[jj]) = Kaxial(ii,jj);
	M(idx_z[ii], idx_z[jj]) = Maxial(ii,jj);
      }
    }
    //project(K, idx_z, idx_z) = Kaxial;
    //project(M, idx_z, idx_z) = Maxial;
    //project(F, idx_z) = Faxial;
    
    // --------- torsion -------------
    // same shape functions as axial
    
    // stiffness matrix
    Matrix Ktorsion(nsz, nsz);
    matrixAssembly(GJ, nsz, fzp, 1.0/L, Ktorsion);

    // inertia matrix
    Matrix Mtorsion(nsz, nsz);
    matrixAssembly(rhoJ, nsz, fz, L, Mtorsion);
    
    // put into global matrix
    idx_z = {5, 11};
    //idx_z(0) = 5; idx_z(1) = 11;
    
    for (int ii=0; ii<nsz; ii++) {
      for (int jj=0; jj<nsz; jj++) {
	K(idx_z[ii], idx_z[jj]) = Ktorsion(ii,jj);
	M(idx_z[ii], idx_z[jj]) = Mtorsion(ii,jj);
      }
    }
    //project(K, idx_z, idx_z) = Ktorsion;
    //project(M, idx_z, idx_z) = Mtorsion;
    
  }