SHOGUN
v3.2.0
|
00001 /* 00002 * This program is free software; you can redistribute it and/or modify 00003 * it under the terms of the GNU General Public License as published by 00004 * the Free Software Foundation; either version 3 of the License, or 00005 * (at your option) any later version. 00006 * 00007 * Written (W) 2013 Roman Votyakov 00008 * 00009 * The abscissae and weights for Gauss-Kronrod rules are taken form 00010 * QUADPACK, which is in public domain. 00011 * http://www.netlib.org/quadpack/ 00012 * 00013 * See header file for which functions are adapted from GNU Octave, 00014 * file quadgk.m: Copyright (C) 2008-2012 David Bateman under GPLv3 00015 * http://www.gnu.org/software/octave/ 00016 */ 00017 00018 #include <shogun/mathematics/Integration.h> 00019 00020 #ifdef HAVE_EIGEN3 00021 00022 #include <shogun/mathematics/eigen3.h> 00023 00024 using namespace shogun; 00025 using namespace Eigen; 00026 00027 namespace shogun 00028 { 00029 00030 #ifndef DOXYGEN_SHOULD_SKIP_THIS 00031 00042 class CITransformFunction : public CFunction 00043 { 00044 public: 00049 CITransformFunction(CFunction* f) 00050 { 00051 SG_REF(f); 00052 m_f=f; 00053 } 00054 00055 virtual ~CITransformFunction() { SG_UNREF(m_f); } 00056 00064 virtual float64_t operator() (float64_t x) 00065 { 00066 float64_t hx=1.0/(1.0-CMath::sq(x)); 00067 float64_t gx=x*hx; 00068 float64_t dgx=(1.0+CMath::sq(x))*CMath::sq(hx); 00069 00070 return (*m_f)(gx)*dgx; 00071 } 00072 00073 private: 00075 CFunction* m_f; 00076 }; 00077 00093 class CILTransformFunction : public CFunction 00094 { 00095 public: 00101 CILTransformFunction(CFunction* f, float64_t b) 00102 { 00103 SG_REF(f); 00104 m_f=f; 00105 m_b=b; 00106 } 00107 00108 virtual ~CILTransformFunction() { SG_UNREF(m_f); } 00109 00117 virtual float64_t operator() (float64_t x) 00118 { 00119 float64_t hx=1.0/(1.0+x); 00120 float64_t gx=x*hx; 00121 float64_t dgx=CMath::sq(hx); 00122 00123 return -(*m_f)(m_b-CMath::sq(gx))*2*gx*dgx; 00124 } 00125 00126 private: 00128 CFunction* m_f; 00129 00131 float64_t m_b; 00132 }; 00133 00149 class CIUTransformFunction : public CFunction 00150 { 00151 public: 00157 CIUTransformFunction(CFunction* f, float64_t a) 00158 { 00159 SG_REF(f); 00160 m_f=f; 00161 m_a=a; 00162 } 00163 00164 virtual ~CIUTransformFunction() { SG_UNREF(m_f); } 00165 00173 virtual float64_t operator() (float64_t x) 00174 { 00175 float64_t hx=1.0/(1.0-x); 00176 float64_t gx=x*hx; 00177 float64_t dgx=CMath::sq(hx); 00178 00179 return (*m_f)(m_a+CMath::sq(gx))*2*gx*dgx; 00180 } 00181 00182 private: 00184 CFunction* m_f; 00185 00187 float64_t m_a; 00188 }; 00189 00200 class CTransformFunction : public CFunction 00201 { 00202 public: 00209 CTransformFunction(CFunction* f, float64_t a, float64_t b) 00210 { 00211 SG_REF(f); 00212 m_f=f; 00213 m_a=a; 00214 m_b=b; 00215 } 00216 00217 virtual ~CTransformFunction() { SG_UNREF(m_f); } 00218 00227 virtual float64_t operator() (float64_t x) 00228 { 00229 float64_t qw=(m_b-m_a)/4.0; 00230 float64_t gx=qw*(x*(3.0-CMath::sq(x)))+(m_b+m_a)/2.0; 00231 float64_t dgx=qw*3.0*(1.0-CMath::sq(x)); 00232 00233 return (*m_f)(gx)*dgx; 00234 } 00235 00236 private: 00238 CFunction* m_f; 00239 00241 float64_t m_a; 00242 00244 float64_t m_b; 00245 }; 00246 00247 #endif /* DOXYGEN_SHOULD_SKIP_THIS */ 00248 00249 float64_t CIntegration::integrate_quadgk(CFunction* f, float64_t a, 00250 float64_t b, float64_t abs_tol, float64_t rel_tol, uint32_t max_iter, 00251 index_t sn) 00252 { 00253 // check the parameters 00254 REQUIRE(f, "Integrable function should not be NULL\n") 00255 REQUIRE(abs_tol>0.0, "Absolute tolerance must be positive, but is %f\n", 00256 abs_tol) 00257 REQUIRE(rel_tol>0.0, "Relative tolerance must be positive, but is %f\n", 00258 rel_tol) 00259 REQUIRE(max_iter>0, "Maximum number of iterations must be greater than 0, " 00260 "but is %d\n", max_iter) 00261 REQUIRE(sn>0, "Initial number of subintervals must be greater than 0, " 00262 "but is %d\n", sn) 00263 00264 // integral evaluation function 00265 typedef void TQuadGKEvaluationFunction(CFunction* f, 00266 CDynamicArray<float64_t>* subs, CDynamicArray<float64_t>* q, 00267 CDynamicArray<float64_t>* err); 00268 00269 TQuadGKEvaluationFunction* evaluate_quadgk; 00270 00271 CFunction* tf; 00272 float64_t ta; 00273 float64_t tb; 00274 float64_t q_sign; 00275 00276 // negate integral value and swap a and b, if a>b 00277 if (a>b) 00278 { 00279 ta=b; 00280 tb=a; 00281 q_sign=-1.0; 00282 } 00283 else 00284 { 00285 ta=a; 00286 tb=b; 00287 q_sign=1.0; 00288 } 00289 00290 // transform integrable function and domain of integration 00291 if (a==-CMath::INFTY && b==CMath::INFTY) 00292 { 00293 tf=new CITransformFunction(f); 00294 evaluate_quadgk=&evaluate_quadgk15; 00295 ta=-1.0; 00296 tb=1.0; 00297 } 00298 else if (a==-CMath::INFTY) 00299 { 00300 tf=new CILTransformFunction(f, b); 00301 evaluate_quadgk=&evaluate_quadgk15; 00302 ta=-1.0; 00303 tb=0.0; 00304 } 00305 else if (b==CMath::INFTY) 00306 { 00307 tf=new CIUTransformFunction(f, a); 00308 evaluate_quadgk=&evaluate_quadgk15; 00309 ta=0.0; 00310 tb=1.0; 00311 } 00312 else 00313 { 00314 tf=new CTransformFunction(f, a, b); 00315 evaluate_quadgk=&evaluate_quadgk21; 00316 ta=-1.0; 00317 tb=1.0; 00318 } 00319 00320 // compute initial subintervals, by dividing domain [a, b] into sn 00321 // parts 00322 CDynamicArray<float64_t>* subs=new CDynamicArray<float64_t>(); 00323 00324 // width of each subinterval 00325 float64_t sw=(tb-ta)/sn; 00326 00327 for (index_t i=0; i<sn; i++) 00328 { 00329 subs->push_back(ta+i*sw); 00330 subs->push_back(ta+(i+1)*sw); 00331 } 00332 00333 // evaluate integrals on initial subintervals 00334 CDynamicArray<float64_t>* q_subs=new CDynamicArray<float64_t>(); 00335 CDynamicArray<float64_t>* err_subs=new CDynamicArray<float64_t>(); 00336 00337 evaluate_quadgk(tf, subs, q_subs, err_subs); 00338 00339 // compute value of integral and error on [a, b] 00340 float64_t q=0.0; 00341 float64_t err=0.0; 00342 00343 for (index_t i=0; i<q_subs->get_num_elements(); i++) 00344 q+=(*q_subs)[i]; 00345 00346 for (index_t i=0; i<err_subs->get_num_elements(); i++) 00347 err+=(*err_subs)[i]; 00348 00349 // evaluate tolerance 00350 float64_t tol=CMath::max(abs_tol, rel_tol*CMath::abs(q)); 00351 00352 // number of iterations 00353 uint32_t iter=1; 00354 00355 CDynamicArray<float64_t>* new_subs=new CDynamicArray<float64_t>(); 00356 00357 while (err>tol && iter<max_iter) 00358 { 00359 // choose and bisect subintervals with estimated error, which 00360 // is larger or equal to tolerance 00361 for (index_t i=0; i<subs->get_num_elements()/2; i++) 00362 { 00363 if (CMath::abs((*err_subs)[i])>=tol*CMath::abs((*subs)[2*i+1]- 00364 (*subs)[2*i])/(tb-ta)) 00365 { 00366 // bisect subinterval 00367 float64_t mid=((*subs)[2*i]+(*subs)[2*i+1])/2.0; 00368 00369 new_subs->push_back((*subs)[2*i]); 00370 new_subs->push_back(mid); 00371 new_subs->push_back(mid); 00372 new_subs->push_back((*subs)[2*i+1]); 00373 00374 // subtract value of the integral and error on this 00375 // subinterval from total value and error 00376 q-=(*q_subs)[i]; 00377 err-=(*err_subs)[i]; 00378 } 00379 } 00380 00381 subs->set_array(new_subs->get_array(), new_subs->get_num_elements(), 00382 new_subs->get_num_elements()); 00383 00384 new_subs->reset_array(); 00385 00386 // break if no new subintervals 00387 if (!subs->get_num_elements()) 00388 break; 00389 00390 // evaluate integrals on selected subintervals 00391 evaluate_quadgk(tf, subs, q_subs, err_subs); 00392 00393 for (index_t i=0; i<q_subs->get_num_elements(); i++) 00394 q+=(*q_subs)[i]; 00395 00396 for (index_t i=0; i<err_subs->get_num_elements(); i++) 00397 err+=(*err_subs)[i]; 00398 00399 // evaluate tolerance 00400 tol=CMath::max(abs_tol, rel_tol*CMath::abs(q)); 00401 00402 iter++; 00403 } 00404 00405 SG_UNREF(new_subs); 00406 00407 if (err>tol) 00408 { 00409 SG_SWARNING("Error tolerance not met. Estimated error is equal to %g " 00410 "after %d iterations\n", err, iter) 00411 } 00412 00413 // clean up 00414 SG_UNREF(subs); 00415 SG_UNREF(q_subs); 00416 SG_UNREF(err_subs); 00417 SG_UNREF(tf); 00418 00419 return q_sign*q; 00420 } 00421 00422 float64_t CIntegration::integrate_quadgh(CFunction* f) 00423 { 00424 SG_REF(f); 00425 00426 // evaluate integral using Gauss-Hermite 64-point rule 00427 float64_t q=evaluate_quadgh64(f); 00428 00429 SG_UNREF(f); 00430 00431 return q; 00432 } 00433 00434 void CIntegration::evaluate_quadgk(CFunction* f, CDynamicArray<float64_t>* subs, 00435 CDynamicArray<float64_t>* q, CDynamicArray<float64_t>* err, index_t n, 00436 float64_t* xgk, float64_t* wg, float64_t* wgk) 00437 { 00438 // check the parameters 00439 REQUIRE(f, "Integrable function should not be NULL\n") 00440 REQUIRE(subs, "Array of subintervals should not be NULL\n") 00441 REQUIRE(!(subs->get_array_size()%2), "Size of the array of subintervals " 00442 "should be even\n") 00443 REQUIRE(q, "Array of values of integrals should not be NULL\n") 00444 REQUIRE(err, "Array of errors should not be NULL\n") 00445 REQUIRE(n%2, "Order of Gauss-Kronrod should be odd\n") 00446 REQUIRE(xgk, "Gauss-Kronrod nodes should not be NULL\n") 00447 REQUIRE(wgk, "Gauss-Kronrod weights should not be NULL\n") 00448 REQUIRE(wg, "Gauss weights should not be NULL\n") 00449 00450 // create eigen representation of subs, xgk, wg, wgk 00451 Map<MatrixXd> eigen_subs(subs->get_array(), 2, subs->get_num_elements()/2); 00452 Map<VectorXd> eigen_xgk(xgk, n); 00453 Map<VectorXd> eigen_wg(wg, n/2); 00454 Map<VectorXd> eigen_wgk(wgk, n); 00455 00456 // compute half width and centers of each subinterval 00457 VectorXd eigen_hw=(eigen_subs.row(1)-eigen_subs.row(0))/2.0; 00458 VectorXd eigen_center=eigen_subs.colwise().sum()/2.0; 00459 00460 // compute Gauss-Kronrod nodes x for each subinterval: x=hw*xgk+center 00461 MatrixXd x=eigen_hw*eigen_xgk.adjoint()+eigen_center* 00462 (VectorXd::Ones(n)).adjoint(); 00463 00464 // compute ygk=f(x) 00465 MatrixXd ygk(x.rows(), x.cols()); 00466 00467 for (index_t i=0; i<ygk.rows(); i++) 00468 for (index_t j=0; j<ygk.cols(); j++) 00469 ygk(i,j)=(*f)(x(i,j)); 00470 00471 // compute value of definite integral on each subinterval 00472 VectorXd eigen_q=((ygk*eigen_wgk.asDiagonal()).rowwise().sum()).cwiseProduct( 00473 eigen_hw); 00474 q->set_array(eigen_q.data(), eigen_q.size()); 00475 00476 // choose function values for Gauss nodes 00477 MatrixXd yg(ygk.rows(), ygk.cols()/2); 00478 00479 for (index_t i=1, j=0; i<ygk.cols(); i+=2, j++) 00480 yg.col(j)=ygk.col(i); 00481 00482 // compute error on each subinterval 00483 VectorXd eigen_err=(((yg*eigen_wg.asDiagonal()).rowwise().sum()).cwiseProduct( 00484 eigen_hw)-eigen_q).array().abs(); 00485 err->set_array(eigen_err.data(), eigen_err.size()); 00486 } 00487 00488 void CIntegration::evaluate_quadgk15(CFunction* f, CDynamicArray<float64_t>* subs, 00489 CDynamicArray<float64_t>* q, CDynamicArray<float64_t>* err) 00490 { 00491 static const index_t n=15; 00492 00493 // Gauss-Kronrod nodes 00494 static float64_t xgk[n]= 00495 { 00496 -0.991455371120812639206854697526329, 00497 -0.949107912342758524526189684047851, 00498 -0.864864423359769072789712788640926, 00499 -0.741531185599394439863864773280788, 00500 -0.586087235467691130294144838258730, 00501 -0.405845151377397166906606412076961, 00502 -0.207784955007898467600689403773245, 00503 0.000000000000000000000000000000000, 00504 0.207784955007898467600689403773245, 00505 0.405845151377397166906606412076961, 00506 0.586087235467691130294144838258730, 00507 0.741531185599394439863864773280788, 00508 0.864864423359769072789712788640926, 00509 0.949107912342758524526189684047851, 00510 0.991455371120812639206854697526329 00511 }; 00512 00513 // Gauss weights 00514 static float64_t wg[n/2]= 00515 { 00516 0.129484966168869693270611432679082, 00517 0.279705391489276667901467771423780, 00518 0.381830050505118944950369775488975, 00519 0.417959183673469387755102040816327, 00520 0.381830050505118944950369775488975, 00521 0.279705391489276667901467771423780, 00522 0.129484966168869693270611432679082 00523 }; 00524 00525 // Gauss-Kronrod weights 00526 static float64_t wgk[n]= 00527 { 00528 0.022935322010529224963732008058970, 00529 0.063092092629978553290700663189204, 00530 0.104790010322250183839876322541518, 00531 0.140653259715525918745189590510238, 00532 0.169004726639267902826583426598550, 00533 0.190350578064785409913256402421014, 00534 0.204432940075298892414161999234649, 00535 0.209482141084727828012999174891714, 00536 0.204432940075298892414161999234649, 00537 0.190350578064785409913256402421014, 00538 0.169004726639267902826583426598550, 00539 0.140653259715525918745189590510238, 00540 0.104790010322250183839876322541518, 00541 0.063092092629978553290700663189204, 00542 0.022935322010529224963732008058970 00543 }; 00544 00545 // evaluate definite integral on each subinterval using Gauss-Kronrod rule 00546 evaluate_quadgk(f, subs, q, err, n, xgk, wg, wgk); 00547 } 00548 00549 void CIntegration::evaluate_quadgk21(CFunction* f, CDynamicArray<float64_t>* subs, 00550 CDynamicArray<float64_t>* q, CDynamicArray<float64_t>* err) 00551 { 00552 static const index_t n=21; 00553 00554 // Gauss-Kronrod nodes 00555 static float64_t xgk[n]= 00556 { 00557 -0.995657163025808080735527280689003, 00558 -0.973906528517171720077964012084452, 00559 -0.930157491355708226001207180059508, 00560 -0.865063366688984510732096688423493, 00561 -0.780817726586416897063717578345042, 00562 -0.679409568299024406234327365114874, 00563 -0.562757134668604683339000099272694, 00564 -0.433395394129247190799265943165784, 00565 -0.294392862701460198131126603103866, 00566 -0.148874338981631210884826001129720, 00567 0.000000000000000000000000000000000, 00568 0.148874338981631210884826001129720, 00569 0.294392862701460198131126603103866, 00570 0.433395394129247190799265943165784, 00571 0.562757134668604683339000099272694, 00572 0.679409568299024406234327365114874, 00573 0.780817726586416897063717578345042, 00574 0.865063366688984510732096688423493, 00575 0.930157491355708226001207180059508, 00576 0.973906528517171720077964012084452, 00577 0.995657163025808080735527280689003 00578 }; 00579 00580 // Gauss weights 00581 static float64_t wg[n/2]= 00582 { 00583 0.066671344308688137593568809893332, 00584 0.149451349150580593145776339657697, 00585 0.219086362515982043995534934228163, 00586 0.269266719309996355091226921569469, 00587 0.295524224714752870173892994651338, 00588 0.295524224714752870173892994651338, 00589 0.269266719309996355091226921569469, 00590 0.219086362515982043995534934228163, 00591 0.149451349150580593145776339657697, 00592 0.066671344308688137593568809893332 00593 }; 00594 00595 // Gauss-Kronrod weights 00596 static float64_t wgk[n]= 00597 { 00598 0.011694638867371874278064396062192, 00599 0.032558162307964727478818972459390, 00600 0.054755896574351996031381300244580, 00601 0.075039674810919952767043140916190, 00602 0.093125454583697605535065465083366, 00603 0.109387158802297641899210590325805, 00604 0.123491976262065851077958109831074, 00605 0.134709217311473325928054001771707, 00606 0.142775938577060080797094273138717, 00607 0.147739104901338491374841515972068, 00608 0.149445554002916905664936468389821, 00609 0.147739104901338491374841515972068, 00610 0.142775938577060080797094273138717, 00611 0.134709217311473325928054001771707, 00612 0.123491976262065851077958109831074, 00613 0.109387158802297641899210590325805, 00614 0.093125454583697605535065465083366, 00615 0.075039674810919952767043140916190, 00616 0.054755896574351996031381300244580, 00617 0.032558162307964727478818972459390, 00618 0.011694638867371874278064396062192 00619 }; 00620 00621 evaluate_quadgk(f, subs, q, err, n, xgk, wg, wgk); 00622 } 00623 00624 float64_t CIntegration::evaluate_quadgh(CFunction* f, index_t n, float64_t* xgh, 00625 float64_t* wgh) 00626 { 00627 // check the parameters 00628 REQUIRE(f, "Integrable function should not be NULL\n"); 00629 REQUIRE(xgh, "Gauss-Hermite nodes should not be NULL\n"); 00630 REQUIRE(wgh, "Gauss-Hermite weights should not be NULL\n"); 00631 00632 float64_t q=0.0; 00633 00634 for (index_t i=0; i<n; i++) 00635 q+=wgh[i]*(*f)(xgh[i]); 00636 00637 return q; 00638 } 00639 00640 float64_t CIntegration::evaluate_quadgh64(CFunction* f) 00641 { 00642 static const index_t n=64; 00643 00644 // Gauss-Hermite nodes 00645 static float64_t xgh[n]= 00646 { 00647 -10.52612316796054588332682628381528, 00648 -9.895287586829539021204461477159608, 00649 -9.373159549646721162545652439723862, 00650 -8.907249099964769757295972885642943, 00651 -8.477529083379863090564166344821916, 00652 -8.073687285010225225858791140758144, 00653 -7.68954016404049682844780422986949, 00654 -7.321013032780949201189569363719477, 00655 -6.965241120551107529242642193492688, 00656 -6.620112262636027379036660108937914, 00657 -6.284011228774828235418093195070243, 00658 -5.955666326799486045344567180984366, 00659 -5.634052164349972147249920483307154, 00660 -5.318325224633270857323649515199378, 00661 -5.007779602198768196443702627184136, 00662 -4.701815647407499816097538015812822, 00663 -4.399917168228137647767932535438923, 00664 -4.101634474566656714970981238455522, 00665 -3.806571513945360461165972000460225, 00666 -3.514375935740906211539950586474333, 00667 -3.224731291992035725848171110188419, 00668 -2.93735082300462180968533902619139, 00669 -2.651972435430635011005457785998431, 00670 -2.368354588632401404111511265341516, 00671 -2.086272879881762020832563302363221, 00672 -1.805517171465544918903773574186889, 00673 -1.525889140209863662948970133151528, 00674 -1.24720015694311794069356453069359, 00675 -0.9692694230711780167435414890191023, 00676 -0.6919223058100445772682192875955947, 00677 -0.4149888241210786845769291291996859, 00678 -0.1383022449870097241150497679666744, 00679 0.1383022449870097241150497679666744, 00680 0.4149888241210786845769291291996859, 00681 0.6919223058100445772682192875955947, 00682 0.9692694230711780167435414890191023, 00683 1.24720015694311794069356453069359, 00684 1.525889140209863662948970133151528, 00685 1.805517171465544918903773574186889, 00686 2.086272879881762020832563302363221, 00687 2.368354588632401404111511265341516, 00688 2.651972435430635011005457785998431, 00689 2.93735082300462180968533902619139, 00690 3.224731291992035725848171110188419, 00691 3.514375935740906211539950586474333, 00692 3.806571513945360461165972000460225, 00693 4.101634474566656714970981238455522, 00694 4.399917168228137647767932535438923, 00695 4.701815647407499816097538015812822, 00696 5.007779602198768196443702627184136, 00697 5.318325224633270857323649515199378, 00698 5.634052164349972147249920483307154, 00699 5.955666326799486045344567180984366, 00700 6.284011228774828235418093195070243, 00701 6.620112262636027379036660108937914, 00702 6.965241120551107529242642193492688, 00703 7.321013032780949201189569363719477, 00704 7.68954016404049682844780422986949, 00705 8.073687285010225225858791140758144, 00706 8.477529083379863090564166344821916, 00707 8.907249099964769757295972885642943, 00708 9.373159549646721162545652439723862, 00709 9.895287586829539021204461477159608, 00710 10.52612316796054588332682628381528 00711 }; 00712 00713 // Gauss-Hermite weights 00714 static float64_t wgh[n]= 00715 { 00716 5.535706535856942820575463300987E-49, 00717 1.6797479901081592186662883306299E-43, 00718 3.4211380112557405043272218281457E-39, 00719 1.557390624629763802309335380265E-35, 00720 2.549660899112999256604766580441E-32, 00721 1.92910359546496685030196877906707E-29, 00722 7.8617977889259103690999914962788E-27, 00723 1.911706883300642829958456965534449E-24, 00724 2.982862784279851154478700702016E-22, 00725 3.15225456650378141612134668341E-20, 00726 2.35188471067581911695767591555844E-18, 00727 1.28009339132243804163956329526337E-16, 00728 5.218623726590847522957808513052588E-15, 00729 1.628340730709720362084307081240893E-13, 00730 3.95917776694772392723644586425458E-12, 00731 7.61521725014545135331529567531937E-11, 00732 1.1736167423215493435425064670822E-9, 00733 1.465125316476109354926622003804004E-8, 00734 1.495532936727247061102461692934817E-7, 00735 1.258340251031184576157842180019028E-6, 00736 8.7884992308503591814440474067043E-6, 00737 5.125929135786274660821911412739621E-5, 00738 2.509836985130624860823620179819094E-4, 00739 0.001036329099507577663456741746283101, 00740 0.00362258697853445876066812537162265, 00741 0.01075604050987913704946517278667313, 00742 0.0272031289536889184538348212614932, 00743 0.0587399819640994345496889462518317, 00744 0.1084983493061868406330258455060973, 00745 0.1716858423490837020007279701237768, 00746 0.2329947860626780466505660293325675, 00747 0.2713774249413039779456065084184279, 00748 0.2713774249413039779456065084184279, 00749 0.2329947860626780466505660293325675, 00750 0.1716858423490837020007279701237768, 00751 0.1084983493061868406330258455060973, 00752 0.0587399819640994345496889462518317, 00753 0.0272031289536889184538348212614932, 00754 0.01075604050987913704946517278667313, 00755 0.00362258697853445876066812537162265, 00756 0.001036329099507577663456741746283101, 00757 2.509836985130624860823620179819094E-4, 00758 5.125929135786274660821911412739621E-5, 00759 8.7884992308503591814440474067043E-6, 00760 1.258340251031184576157842180019028E-6, 00761 1.495532936727247061102461692934817E-7, 00762 1.465125316476109354926622003804004E-8, 00763 1.1736167423215493435425064670822E-9, 00764 7.61521725014545135331529567531937E-11, 00765 3.95917776694772392723644586425458E-12, 00766 1.628340730709720362084307081240893E-13, 00767 5.218623726590847522957808513052588E-15, 00768 1.28009339132243804163956329526337E-16, 00769 2.35188471067581911695767591555844E-18, 00770 3.15225456650378141612134668341E-20, 00771 2.982862784279851154478700702016E-22, 00772 1.911706883300642829958456965534449E-24, 00773 7.8617977889259103690999914962788E-27, 00774 1.92910359546496685030196877906707E-29, 00775 2.549660899112999256604766580441E-32, 00776 1.557390624629763802309335380265E-35, 00777 3.4211380112557405043272218281457E-39, 00778 1.6797479901081592186662883306299E-43, 00779 5.535706535856942820575463300987E-49 00780 }; 00781 00782 return evaluate_quadgh(f, n, xgh, wgh); 00783 } 00784 } 00785 00786 #endif /* HAVE_EIGEN3 */